text:
module TiPiTests

using OptimPackNextGen
using OptimPackNextGen.Algebra
using TiPi
using Base.Test

function reldif(a::Real, b::Real)
    a == b && return 0.0
    2*abs(a - b)/(abs(a) + abs(b))
end

function reldif{T<:Real}(a::Array{T}, b::Array{T})
    dn = enorm(a - b)
    dn == 0 && return 0.0
    2*dn/(enorm(a) + enorm(b))
end

function enorm{T<:Real}(a::Array{T})
    s::Cdouble = 0.0
    for ai in a
        s += convert(Cdouble, ai)^2
    end
    return sqrt(s)
end

function xtime_print(str, val, tbg)
    t, b, g = tbg
    if 0 < g
        @printf("%s %.15E [elapsed time: %.3E seconds, %d bytes allocated, %.2f%% gc time]\n",
                str, val, t/1e9, b, 100*g/t)
    else
        @printf("%s %.15E [elapsed time: %.3E seconds, %d bytes allocated]\n",
                str, val, t/1e9, b)
    end
end

macro xtime(ex)
    quote
        local b0 = Base.gc_bytes()
        local t0 = Base.time_ns()
        local g0 = Base.gc_time_ns()
        local val = $(esc(ex))
        local g1 = Base.gc_time_ns()
        local t1 = Base.time_ns()
        local b1 = Base.gc_bytes()
        val, (t1-t0, b1-b0, g1-g0)
    end
end

function testcost{T}(alpha::Real, param::HyperbolicEdgePreserving,
                     x::Array{T})
    res = @xtime(cost(alpha, param, x))
    dims = size(x)
    rank = length(dims)
    if rank == 1
        ref = @xtime begin
            mu1 = convert(T, param.mu[1])
            eps = convert(T, param.eps)
            dim1 = dims[1]
            x1 = x[1:end-1]
            x2 = x[2:end]
            alpha*(sum(sqrt((mu1)^2*(x2 - x1).^2 + eps^2)) - length(x1)*eps)
        end
    elseif rank == 2
        ref = @xtime begin
            mu1 = convert(T, param.mu[1])
            mu2 = convert(T, param.mu[2])
            eps = convert(T, param.eps)
            dim1 = dims[1]
            dim2 = dims[2]
            x1 = x[ 1:end-1 , 1:end-1 ]
            x2 = x[ 2:end   , 1:end-1 ]
            x3 = x[ 1:end-1 , 2:end   ]
            x4 = x[ 2:end   , 2:end   ]
            alpha*(sum(sqrt((mu1^2/2)*((x2 - x1).^2 + (x4 - x3).^2) +
                            (mu2^2/2)*((x3 - x1).^2 + (x4 - x2).^2) + eps^2))
                   - length(x1)*eps)
        end
    elseif rank == 3
        ref = @xtime begin
            mu1 = convert(T, param.mu[1])
            mu2 = convert(T, param.mu[2])
            mu3 = convert(T, param.mu[3])
            eps = convert(T, param.eps)
            dim1 = dims[1]
            dim2 = dims[2]
            x1 = x[ 1:end-1 , 1:end-1 , 1:end-1 ]
            x2 = x[ 2:end   , 1:end-1 , 1:end-1 ]
            x3 = x[ 1:end-1 , 2:end   , 1:end-1 ]
            x4 = x[ 2:end   , 2:end   , 1:end-1 ]
            x5 = x[ 1:end-1 , 1:end-1 , 2:end   ]
            x6 = x[ 2:end   , 1:end-1 , 2:end   ]
            x7 = x[ 1:end-1 , 2:end   , 2:end   ]
            x8 = x[ 2:end   , 2:end   , 2:end   ]
            alpha*(sum(sqrt((mu1^2/4)*((x2 - x1).^2 + (x4 - x3).^2 +
                                       (x6 - x5).^2 + (x8 - x7).^2) +
                            (mu2^2/4)*((x3 - x1).^2 + (x4 - x2).^2 +
                                       (x7 - x5).^2 + (x8 - x6).^2) +
                            (mu3^2/4)*((x5 - x1).^2 + (x6 - x2).^2 +
                                       (x7 - x3).^2 + (x8 - x4).^2) + eps^2))
                   - length(x1)*eps)
        end
    end
    err = reldif(res[1], ref[1])/(Base.eps(T)*sqrt(length(x)))
    xtime_print("   reference   =", ref...)
    xtime_print("   code result =", res...)
    @printf("   normalized error = %.3E\n", err)
    return err
end

function testcost!{T}(alpha::Real, param::HyperbolicEdgePreserving,
                      x::Array{T})
    dims = size(x)
    gx = Array(T, dims)
    gxref = Array(T, dims)
    res = @xtime(cost!(alpha, param, x, gx, true))
    vfill!(gxref, 0)
    rank = length(dims)
    if rank == 1
        ref = @xtime begin
            mu1 = convert(T, param.mu[1])
            eps = convert(T, param.eps)
            dim1 = dims[1]
            x1 = x[1:end-1]
            x2 = x[2:end]
            x2_x1 = x2 - x1
            r = sqrt(mu1^2*x2_x1.^2 + eps^2)
            fx = alpha*(sum(r) - length(x1)*eps)
            q = (alpha*mu1^2)./r.*x2_x1
            gxref[1:end-1] -= q
            gxref[2:end  ] += q
            fx
        end
    elseif rank == 2
        ref = @xtime begin
            mu1 = convert(T, param.mu[1])
            mu2 = convert(T, param.mu[2])
            eps = convert(T, param.eps)
            dim1 = dims[1]
            dim2 = dims[2]
            x1 = x[ 1:end-1 , 1:end-1 ]
            x2 = x[ 2:end   , 1:end-1 ]
            x3 = x[ 1:end-1 , 2:end   ]
            x4 = x[ 2:end   , 2:end   ]
            r = sqrt((mu1^2/2)*((x2 - x1).^2 + (x4 - x3).^2) +
                     (mu2^2/2)*((x3 - x1).^2 + (x4 - x2).^2) + eps^2)
            fx = alpha*(sum(r) - length(x1)*eps)
            q = alpha./r
            q1 = (mu1^2/2)*q
            q2 = (mu2^2/2)*q
            gxref[ 1:end-1 , 1:end-1 ] += q1.*(x1 - x2) + q2.*(x1 - x3)
            gxref[ 2:end   , 1:end-1 ] += q1.*(x2 - x1) + q2.*(x2 - x4)
            gxref[ 1:end-1 , 2:end   ] += q1.*(x3 - x4) + q2.*(x3 - x1)
            gxref[ 2:end   , 2:end   ] += q1.*(x4 - x3) + q2.*(x4 - x2)
            fx
        end
    elseif rank == 3
        ref = @xtime begin
            mu1 = convert(T, param.mu[1])
            mu2 = convert(T, param.mu[2])
            mu3 = convert(T, param.mu[3])
            eps = convert(T, param.eps)
            dim1 = dims[1]
            dim2 = dims[2]
            x1 = x[ 1:end-1 , 1:end-1 , 1:end-1 ]
            x2 = x[ 2:end   , 1:end-1 , 1:end-1 ]
            x3 = x[ 1:end-1 , 2:end   , 1:end-1 ]
            x4 = x[ 2:end   , 2:end   , 1:end-1 ]
            x5 = x[ 1:end-1 , 1:end-1 , 2:end   ]
            x6 = x[ 2:end   , 1:end-1 , 2:end   ]
            x7 = x[ 1:end-1 , 2:end   , 2:end   ]
            x8 = x[ 2:end   , 2:end   , 2:end   ]
            r = sqrt((mu1^2/4)*((x2 - x1).^2 + (x4 - x3).^2 +
                                (x6 - x5).^2 + (x8 - x7).^2) +
                     (mu2^2/4)*((x3 - x1).^2 + (x4 - x2).^2 +
                                (x7 - x5).^2 + (x8 - x6).^2) +
                     (mu3^2/4)*((x5 - x1).^2 + (x6 - x2).^2 +
                                (x7 - x3).^2 + (x8 - x4).^2) + eps^2)
            fx = alpha*(sum(r) - length(x1)*eps)
            q = alpha./r
            q1 = (mu1^2/4)*q
            q2 = (mu2^2/4)*q
            q3 = (mu3^2/4)*q
            gxref[ 1:end-1 , 1:end-1 , 1:end-1 ] += q1.*(x1 - x2) + q2.*(x1 - x3) + q3.*(x1 - x5)
            gxref[ 2:end   , 1:end-1 , 1:end-1 ] += q1.*(x2 - x1) + q2.*(x2 - x4) + q3.*(x2 - x6)
            gxref[ 1:end-1 , 2:end   , 1:end-1 ] += q1.*(x3 - x4) + q2.*(x3 - x1) + q3.*(x3 - x7)
            gxref[ 2:end   , 2:end   , 1:end-1 ] += q1.*(x4 - x3) + q2.*(x4 - x2) + q3.*(x4 - x8)
            gxref[ 1:end-1 , 1:end-1 , 2:end   ] += q1.*(x5 - x6) + q2.*(x5 - x7) + q3.*(x5 - x1)
            gxref[ 2:end   , 1:end-1 , 2:end   ] += q1.*(x6 - x5) + q2.*(x6 - x8) + q3.*(x6 - x2)
            gxref[ 1:end-1 , 2:end   , 2:end   ] += q1.*(x7 - x8) + q2.*(x7 - x5) + q3.*(x7 - x3)
            gxref[ 2:end   , 2:end   , 2:end   ] += q1.*(x8 - x7) + q2.*(x8 - x6) + q3.*(x8 - x4)
            fx
        end
    end
    fx_err = reldif(res[1], ref[1])/(Base.eps(T)*sqrt(length(x)))
    gx_err = reldif(gx, gxref)
    xtime_print("   reference   =", ref...)
    xtime_print("   code result =", res...)
    @printf("   normalized error for f(x) = %.3E\n", fx_err)
    @printf("   relative error for g(x) = %.3E\n", gx_err)
    return fx_err
end

# Weight and threshold for all cases:
alpha = 1.0 + rand()
eps = 0.02

# 1-D tests
println("\nRunning 1D tests...")
dims = (178,)
x = rand(dims) - 0.5
mu = (1.7,)
param = HyperbolicEdgePreserving(eps, mu)
@test testcost(alpha, param, x) <= 0.3
testcost!(alpha, param, x)

# 2-D tests
println("\nRunning 2D tests...")
dims = (174, 200)
x = rand(dims) - 0.5

println(" - Isotropic case:")
mu = (1.2, 1.2)
param = HyperbolicEdgePreserving(eps, mu)
@test testcost(alpha, param, x) <= 0.3
testcost!(alpha, param, x)

println(" - Anisotropic case:")
mu = (1.7, 1.2)
param = HyperbolicEdgePreserving(eps, mu)
@test testcost(alpha, param, x) <= 0.3
testcost!(alpha, param, x)

# 3-D tests
println("\nRunning 3D tests...")
dims = (174, 200, 108)
x = rand(dims) - 0.5

println(" - Isotropic case:")
mu = (1.2, 1.2, 1.2)
param = HyperbolicEdgePreserving(eps, mu)
@test testcost(alpha, param, x) <= 0.3
@test testcost!(alpha, param, x) <= 0.3

println(" - Anisotropic case:")
mu = (1.7, 1.2, 0.9)
param = HyperbolicEdgePreserving(eps, mu)
@test testcost(alpha, param, x) <= 0.3
@test testcost!(alpha, param, x) <= 0.3

nothing

end # module


text:
using DelimitedFiles

data = readdlm("./data/day01.txt")

# Part 1.
for d in data
    if (2020 - d) in data
        result = convert(Integer, d * (2020 - d))
        println("Two numbers: $result")
        break
    end
end

# Part 2.
found = false
for d in data
    for e in [i for i in data if ((d + i) < 2020)]
        if (2020 - d - e) in data
            result = convert(Integer, d * e * (2020 - d - e))
            println("Three numbers: $result")
            global found = true  # This sucks.
            break
        end
        if found
            break
        end
    end
end

println("Done.")

ext: jl


text:
module OpenTelemetryExporterOtlpProtoGrpc

export OtlpProtoGrpcTraceExporter

using gRPCClient
using URIs

import OpenTelemetryAPI
import OpenTelemetrySDK
import OpenTelemetryProto

import Logging

const API = OpenTelemetryAPI
const SDK = OpenTelemetrySDK
const Otlp = OpenTelemetryProto.OpentelemetryClients
const Trace = OpenTelemetryProto.OpentelemetryClients.opentelemetry.proto.trace.v1
const Metrics = OpenTelemetryProto.OpentelemetryClients.opentelemetry.proto.metrics.v1
const Logs = OpenTelemetryProto.OpentelemetryClients.opentelemetry.proto.logs.v1
const CollectorTrace =
    OpenTelemetryProto.OpentelemetryClients.opentelemetry.proto.collector.trace.v1
const CollectorMetrics =
    OpenTelemetryProto.OpentelemetryClients.opentelemetry.proto.collector.metrics.v1
const CollectorLogs =
    OpenTelemetryProto.OpentelemetryClients.opentelemetry.proto.collector.logs.v1
const Resource = OpenTelemetryProto.OpentelemetryClients.opentelemetry.proto.resource.v1
const Common = OpenTelemetryProto.OpentelemetryClients.opentelemetry.proto.common.v1

##### Traces

struct OtlpProtoGrpcTraceExporter{T} <: SDK.AbstractExporter
    client::T
end

"""
    OtlpProtoGrpcTraceExporter(;kw...)

## Keyword arguments

  - `scheme=http`
  - `host="localhost"`
  - `port=4317`
  - `is_blocking=true`, by default the `BlockingClient` is used.
  - Rest keyword arguments will be forward to the gRPC client.

`scheme`, `host` and `port` specifies the OTEL Collector to connect with.
"""
function OtlpProtoGrpcTraceExporter(;
    scheme = "http",
    host = "localhost",
    port = 4317,
    is_blocking = true,
    kw...,
)
    url = string(URI(; scheme = scheme, host = host, port = port))
    if is_blocking
        client = Otlp.TraceServiceBlockingClient(url; kw...)
    else
        client = Otlp.TraceServiceClient(url; kw...)
    end
    OtlpProtoGrpcTraceExporter(client)
end

function SDK.export!(se::OtlpProtoGrpcTraceExporter, sp)
    res, status = CollectorTrace.Export(
        se.client,
        convert(CollectorTrace.ExportTraceServiceRequest, sp),
    )
    if gRPCClient.gRPCCheck(status; throw_error = false)
        SDK.EXPORT_SUCCESS
    else
        SDK.EXPORT_FAILURE
    end
end

Base.convert(t::Type{CollectorTrace.ExportTraceServiceRequest}, s::API.AbstractSpan) =
    convert(t, [s])

Base.convert(
    ::Type{CollectorTrace.ExportTraceServiceRequest},
    s::CollectorTrace.ExportTraceServiceRequest,
) = s

function Base.convert(::Type{CollectorTrace.ExportTraceServiceRequest}, spans)
    r = CollectorTrace.ExportTraceServiceRequest(; resource_spans = [])
    s_res_pre = nothing
    s_ins_pre = nothing
    for s in spans
        s_res = API.resource(s)
        if s_res != s_res_pre
            push!(
                r.resource_spans,
                Trace.ResourceSpans(
                    resource = convert(Resource.Resource, s_res),
                    schema_url = s_res.schema_url,
                    instrumentation_library_spans = [],
                ),
            )
        end
        s_ins = API.tracer(s).instrumentation_info
        if s_ins != s_ins_pre
            push!(
                r.resource_spans[end].instrumentation_library_spans,
                Trace.InstrumentationLibrarySpans(
                    instrumentation_library = convert(Common.InstrumentationLibrary, s_ins),
                    schema_url = s_ins.schema_url,
                    spans = [],
                ),
            )
        end

        push!(
            r.resource_spans[end].instrumentation_library_spans[end].spans,
            convert(Trace.Span, s),
        )

        s_res_pre = s_res
        s_ins_pre = s_ins
    end
    r
end

function Base.convert(::Type{Resource.Resource}, r::API.Resource)
    Resource.Resource(
        attributes = convert(Vector{Common.KeyValue}, r.attributes),
        dropped_attributes_count = API.n_dropped(r.attributes),
    )
end

function Base.convert(::Type{Common.InstrumentationLibrary}, info::API.InstrumentationInfo)
    Common.InstrumentationLibrary(name = info.name, version = string(info.version))
end

# TODO: use method call instead of accessing field directly
function Base.convert(::Type{Trace.Span}, s::API.AbstractSpan)
    Trace.Span(
        trace_id = reinterpret(UInt8, [API.span_context(s).trace_id]),
        span_id = reinterpret(UInt8, [API.span_context(s).span_id]),
        trace_state = string(API.span_context(s).trace_state),
        parent_span_id = if isnothing(API.parent_span_context(s))
            UInt8[]
        else
            reinterpret(UInt8, [API.parent_span_context(s).span_id])
        end,
        name = API.span_name(s),
        kind = Int32(API.span_kind(s)),
        start_time_unix_nano = API.start_time(s),
        end_time_unix_nano = API.end_time(s),
        attributes = convert(Vector{Common.KeyValue}, API.attributes(s)),
        dropped_attributes_count = UInt32(API.n_dropped(API.attributes(s))),
        events = [convert(Trace.Span_Event, e) for e in API.span_events(s)],
        dropped_events_count = UInt32(API.n_dropped(API.span_events(s))),
        links = [convert(Trace.Span_Link, e) for e in API.span_links(s)],
        dropped_links_count = UInt32(API.n_dropped(API.span_links(s))),
        status = convert(Trace.Status, API.span_status(s)),
    )
end

function Base.convert(
    ::Type{Vector{Common.KeyValue}},
    attrs::Union{API.StaticAttrs,API.DynamicAttrs},
)
    [
        Common.KeyValue(key = string(k), value = convert(Common.AnyValue, v)) for
        (k, v) in pairs(attrs)
    ]
end

Base.convert(::Type{Common.AnyValue}, v::String) = Common.AnyValue(string_value = v)
Base.convert(::Type{Common.AnyValue}, v::Bool) = Common.AnyValue(bool_value = v)
Base.convert(::Type{Common.AnyValue}, v::Int) = Common.AnyValue(int_value = v)
Base.convert(::Type{Common.AnyValue}, v::Float64) = Common.AnyValue(double_value = v)
Base.convert(::Type{Common.AnyValue}, v::Vector{UInt8}) = Common.AnyValue(bytes_value = v)
Base.convert(::Type{Common.AnyValue}, v::Vector) =
    Common.AnyValue(array_value = convert(Common.ArrayValue, v))
Base.convert(::Type{Common.AnyValue}, v::Union{API.StaticAttrs,API.DynamicAttrs}) =
    Common.AnyValue(kvlist_value = convert(Common.KeyValueList, v))

Base.convert(::Type{Common.ArrayValue}, v::Vector) =
    Common.ArrayValue(values = [convert(Common.AnyValue, x) for x in v])
Base.convert(::Type{Common.KeyValueList}, v::Union{API.StaticAttrs,API.DynamicAttrs}) =
    Common.KeyValueList(values = convert(Vector{Common.KeyValue}, v))

function Base.convert(::Type{Trace.Span_Event}, event::API.Event)
    Trace.Span_Event(
        time_unix_nano = event.timestamp,
        name = event.name,
        attributes = convert(Vector{Common.KeyValue}, event.attributes),
        dropped_attributes_count = API.n_dropped(event.attributes),
    )
end

function Base.convert(::Type{Trace.Span_Link}, link::API.Link)
    Trace.Span_Link(
        trace_id = reinterpret(UInt8, [link.context.trace_id]),
        span_id = reinterpret(UInt8, [link.context.span_id]),
        trace_state = string(link.context.trace_state),
        attributes = convert(Vector{Common.KeyValue}, link.attributes),
        dropped_attributes_count = API.n_dropped(link.attributes),
    )
end

function Base.convert(::Type{Trace.Status}, status::API.SpanStatus)
    Trace.Status(code = Int32(status.code), message = something(status.description, ""))
end

##### Metrics

export OtlpProtoGrpcMetricsExporter

struct OtlpProtoGrpcMetricsExporter{T} <: SDK.AbstractExporter
    client::T
end

"""
    OtlpProtoGrpcMetricsExporter(;kw...)

This exporter is push-based.

## Keyword arguments

  - `scheme=http`
  - `host="localhost"`
  - `port=4317`
  - `is_blocking=true`, by default the `BlockingClient` is used.
  - Rest keyword arguments will be forward to the gRPC client.

`scheme`, `host` and `port` specifies the OTEL Collector to connect with.
"""
function OtlpProtoGrpcMetricsExporter(;
    scheme = "http",
    host = "localhost",
    port = 4317,
    is_blocking = true,
    kw...,
)
    url = string(URI(; scheme = scheme, host = host, port = port))
    if is_blocking
        client = Otlp.MetricsServiceBlockingClient(url; kw...)
    else
        client = Otlp.MetricsServiceClient(url; kw...)
    end
    OtlpProtoGrpcMetricsExporter(client)
end

function SDK.export!(e::OtlpProtoGrpcMetricsExporter, ms)
    res, status = CollectorMetrics.Export(
        e.client,
        convert(CollectorMetrics.ExportMetricsServiceRequest, ms),
    )
    if gRPCClient.gRPCCheck(status; throw_error = false)
        SDK.EXPORT_SUCCESS
    else
        SDK.EXPORT_FAILURE
    end
end

Base.convert(t::Type{CollectorMetrics.ExportMetricsServiceRequest}, s::SDK.Metric) =
    convert(t, [s])

Base.convert(
    t::Type{CollectorMetrics.ExportMetricsServiceRequest},
    ms::CollectorMetrics.ExportMetricsServiceRequest,
) = ms

function Base.convert(t::Type{CollectorMetrics.ExportMetricsServiceRequest}, ms)
    r = CollectorMetrics.ExportMetricsServiceRequest(; resource_metrics = [])
    m_res_pre = nothing
    m_ins_pre = nothing
    for m in ms
        m_res = API.resource(m)
        if m_res != m_res_pre
            push!(
                r.resource_metrics,
                Metrics.ResourceMetrics(
                    resource = convert(Resource.Resource, m_res),
                    schema_url = m_res.schema_url,
                    instrumentation_library_metrics = [],
                ),
            )
        end
        m_ins = m.instrument.meter.instrumentation_info
        if m_ins != m_ins_pre
            push!(
                r.resource_metrics[end].instrumentation_library_metrics,
                Metrics.InstrumentationLibraryMetrics(
                    instrumentation_library = convert(Common.InstrumentationLibrary, m_ins),
                    schema_url = m_ins.schema_url,
                    metrics = [],
                ),
            )
        end

        push!(
            r.resource_metrics[end].instrumentation_library_metrics[end].metrics,
            convert(Metrics.Metric, m),
        )

        m_res_pre = m_res
        m_ins_pre = m_ins
    end
    r
end

function Base.convert(::Type{Metrics.Metric}, m::SDK.Metric{<:SDK.LastValueAgg})
    Metrics.Metric(
        name = m.name,
        description = m.description,
        unit = m.instrument.unit,
        gauge = convert(Metrics.Gauge, m.aggregation),
    )
end

function Base.convert(::Type{Metrics.Gauge}, agg::SDK.LastValueAgg)
    Metrics.Gauge(data_points = convert(Vector{Metrics.NumberDataPoint}, agg.agg_store))
end

function Base.convert(
    ::Type{Vector{Metrics.NumberDataPoint}},
    agg_store::SDK.AggregationStore,
)
    data_points = Metrics.NumberDataPoint[]
    for (attr, data_point) in agg_store
        if data_point.value isa Integer
            push!(
                data_points,
                Metrics.NumberDataPoint(
                    attributes = convert(Vector{Common.KeyValue}, attr),
                    start_time_unix_nano = data_point.start_time_unix_nano,
                    time_unix_nano = data_point.time_unix_nano,
                    as_int = data_point.value,
                    exemplars = [], # TODO
                    flags = Metrics.DataPointFlags.FLAG_NONE,
                ),
            )
        else
            push!(
                data_points,
                Metrics.NumberDataPoint(
                    attributes = convert(Vector{Common.KeyValue}, attr),
                    # start_time_unix_nano = data_point.start_time_unix_nano,
                    time_unix_nano = data_point.time_unix_nano,
                    as_double = data_point.value,
                    exemplars = [], # TODO
                    flags = Metrics.DataPointFlags.FLAG_NONE,
                ),
            )
        end
    end
    data_points
end

function Base.convert(::Type{Metrics.Metric}, m::SDK.Metric{<:SDK.SumAgg})
    Metrics.Metric(
        name = m.name,
        description = m.description,
        unit = m.instrument.unit,
        sum = convert(Metrics.Sum, m.aggregation),
    )
end

function Base.convert(::Type{Metrics.Sum}, agg::SDK.SumAgg)
    Metrics.Sum(
        data_points = convert(Vector{Metrics.NumberDataPoint}, agg.agg_store),
        aggregation_temporality = Metrics.AggregationTemporality.AGGREGATION_TEMPORALITY_CUMULATIVE,  # TODO
        is_monotonic = false, # TODO
    )
end

function Base.convert(::Type{Metrics.Metric}, m::SDK.Metric{<:SDK.HistogramAgg})
    Metrics.Metric(
        name = m.name,
        description = m.description,
        unit = m.instrument.unit,
        histogram = convert(Metrics.Histogram, m.aggregation),
    )
end

function Base.convert(::Type{Metrics.Histogram}, agg::SDK.HistogramAgg)
    Metrics.Histogram(
        data_points = convert(Vector{Metrics.HistogramDataPoint}, agg.agg_store),
        aggregation_temporality = Metrics.AggregationTemporality.AGGREGATION_TEMPORALITY_CUMULATIVE,  # TODO
    )
end

function Base.convert(
    ::Type{Vector{Metrics.HistogramDataPoint}},
    agg_store::SDK.AggregationStore,
)
    [
        Metrics.HistogramDataPoint(
            attributes = convert(Vector{Common.KeyValue}, attr),
            start_time_unix_nano = data_point.start_time_unix_nano,
            time_unix_nano = data_point.time_unix_nano,
            count = sum(data_point.value.counts),
            sum = something(data_point.value.sum, 0),
            bucket_counts = [Float64(x) for x in data_point.value.counts],
            explicit_bounds = collect(data_point.value.boundaries),
            exemplars = [], # TODO
            flags = Metrics.DataPointFlags.FLAG_NONE,
        ) for (attr, data_point) in agg_store
    ]
end

##### Logs

export OtlpProtoGrpcLogsExporter

struct OtlpProtoGrpcLogsExporter{T} <: SDK.AbstractExporter
    client::T
end

"""
    OtlpProtoGrpcLogsExporter(;kw...)

This can be used as a **sink** in [`LoggingExtras.jl`](https://github.com/JuliaLogging/LoggingExtras.jl).

!!! note

    This logger only export one [`LogRecord`](@ref) at a time. A batched version may be provided later.

## Keyword arguments

  - `scheme=http`
  - `host="localhost"`
  - `port=4317`
  - `is_blocking=true`, by default the `BlockingClient` is used.
  - Rest keyword arguments will be forward to the gRPC client.

`scheme`, `host` and `port` specifies the OTEL Collector to connect with.
"""
function OtlpProtoGrpcLogsExporter(;
    scheme = "http",
    host = "localhost",
    port = 4317,
    is_blocking = true,
    kw...,
)
    url = string(URI(; scheme = scheme, host = host, port = port))
    if is_blocking
        client = Otlp.LogsServiceBlockingClient(url; kw...)
    else
        client = Otlp.LogsServiceClient(url; kw...)
    end
    OtlpProtoGrpcLogsExporter(client)
end

function SDK.export!(e::OtlpProtoGrpcLogsExporter, logs)
    res, status = CollectorLogs.Export(
        e.client,
        convert(CollectorLogs.ExportLogsServiceRequest, logs),
    )
    if gRPCClient.gRPCCheck(status; throw_error = false)
        SDK.EXPORT_SUCCESS
    else
        SDK.EXPORT_FAILURE
    end
end

Base.convert(t::Type{CollectorLogs.ExportLogsServiceRequest}, s::SDK.LogRecord) =
    convert(t, [s])

Base.convert(
    t::Type{CollectorLogs.ExportLogsServiceRequest},
    r::CollectorLogs.ExportLogsServiceRequest,
) = r

function Base.convert(t::Type{CollectorLogs.ExportLogsServiceRequest}, log_records)
    CollectorLogs.ExportLogsServiceRequest(
        resource_logs = [
            Logs.ResourceLogs(
                resource = convert(Resource.Resource, first(log_records).resource),
                schema_url = first(log_records).resource.schema_url,
                instrumentation_library_logs = [
                    Logs.InstrumentationLibraryLogs(
                        instrumentation_library = convert(
                            Common.InstrumentationLibrary,
                            first(log_records).instrumentation_info,
                        ),
                        log_records = [
                            Logs.LogRecord(
                                time_unix_nano = r.timestamp,
                                observed_time_unix_nano = UInt(time() * 10^9),
                                severity_number = r.severity_number,
                                severity_text = r.severity_text,
                                body = convert(Common.AnyValue, r.body),
                                attributes = convert(Vector{Common.KeyValue}, r.attributes),
                                dropped_attributes_count = API.n_dropped(r.attributes),
                                trace_id = reinterpret(UInt8, [r.trace_id]),
                                span_id = reinterpret(UInt8, [r.span_id]),
                                flags = r.trace_flags.sampled ? Int32(1) : Int32(0),
                            ) for r in log_records
                        ],
                        schema_url = first(log_records).instrumentation_info.schema_url,
                    ),
                ],
            ),
        ],
    )
end

function Logging.handle_message(
    logger::OtlpProtoGrpcLogsExporter,
    log_record, # single LogRecord or an iterator of LogRecords
)
    SDK.export!(logger, log_record)
end

Logging.shouldlog(logger::OtlpProtoGrpcLogsExporter, arg...) = true
Logging.min_enabled_level(logger::OtlpProtoGrpcLogsExporter) = Logging.BelowMinLevel
Logging.catch_exceptions(logger::OtlpProtoGrpcLogsExporter) = true
end # module


text:
a = 1
println(a)


text:
mutable struct Sigmoid <: DiffableFunction
    grad_field :: GradField
end

function _sigmoid(x::T) where {T <: Real}
    one(T) / (one(T) + exp(-x))
end

@inline forward(::Sigmoid, x) = _sigmoid(x)

function backward(f::Sigmoid, gy::Variable{T}) where {T <: Real}
    y = f.grad_field.outputs[1]
    return (y * (1 - y)) * gy
end

function backward(f::Sigmoid, gy::Variable{T}) where {T <: AbstractArray}
    y = f.grad_field.outputs[1]
    @. return (y * (1 - y)) * gy
end

sigmoid(x::Variable) = Sigmoid(GradField())(x)

Base.broadcasted(::typeof(sigmoid), x::Variable)  = Broadcasting(_sigmoid, Sigmoid(GradField()))(x)


text:
include("./1-graph.jl")
include("./2-variable.jl")
include("./3-math.jl")
include("./4-ArrayFun.jl")
include("./5-activation.jl")
include("./6-pool.jl")
include("./7-broadcast.jl")

max_stars_repo_name: sonosole/Delta

max_issues_repo_name: sonosole/Delta

max_forks_repo_name: sonosole/Delta


text:
mid(x::Real) = x
mid(x::Complex) = x

"""
    interval_isapprox(a::Interval, b::Interval; kwargs)

Checks whether the intervals ``a`` and ``b`` are approximate equal, that is both their
lower and upper bound are approximately equal.

### Keywords

Same of `Base.isapprox`

### Example

```jldoctest
julia> a = 1..2
[1, 2]

julia> b = a + 1e-10
[1, 2.00001]

julia> interval_isapprox(a, b)
true

julia> interval_isapprox(a, b; atol=1e-15)
false
```
"""
interval_isapprox(a::Interval, b::Interval; kwargs...) = isapprox(a.lo, b.lo; kwargs...) && isapprox(a.hi, b.hi; kwargs...)

"""
    interval_norm(A::AbstractMatrix{T}) where {T<:Interval}

computes the infinity norm of interval matrix ``A``.

### Examples

```jldoctest
julia> A = [2..4 -1..1; -1..1 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-1, 1]
 [-1, 1]   [2, 4]

julia> interval_norm(A)
5.0
```
"""
interval_norm(A::AbstractMatrix{T}) where {T<:Interval} = opnorm(mag.(A), Inf) # TODO: proper expand norm function and add 1-norm

"""
    interval_norm(A::AbstractVector{T}) where {T<:Interval}

computes the infinity norm of interval vector ``v``.

### Examples

```jldoctest
julia> b = [-2..2, -3..2]
2-element Vector{Interval{Float64}}:
 [-2, 2]
 [-3, 2]

julia> interval_norm(b)
3.0
```
"""
interval_norm(v::AbstractVector{T}) where {T<:Interval} = maximum(mag.(v))
# ? use manual loops instead

"""
    enclose(A::AbstractMatrix{T}, b::AbstractVector{T}) where {T<:Interval}

Computes an enclosure of the solution of the interval linear system ``Ax=b`` using the
algorithm described in sec. 5.7.1 of [[HOR19]](@ref).
"""
function enclose(A::StaticMatrix{N, N, T}, b::StaticVector{N, T}) where {N, T<:Interval}
    C = inv(mid.(A))
    A1 = Diagonal(ones(N)) - C*A
    e = interval_norm(C*b)/(1 - interval_norm(A1))
    x0 = MVector{N, T}(ntuple(_ -> -e..e, Val(N)))
    return x0
end

function enclose(A::AbstractMatrix{T}, b::AbstractVector{T}) where {T<:Interval}
    n = length(b)
    C = inv(mid.(A))
    A1 = Diagonal(ones(n)) - C*A
    e = interval_norm(C*b)/(1 - interval_norm(A1))
    x0 = fill(-e..e, n)
    return x0
end

"""
    comparison_matrix(A::AbstractMatrix{T}) where {T<:Interval}

Computes the comparison matrix ``⟨A⟩`` of the given interval matrix ``A`` according to the
definition ``⟨A⟩ᵢᵢ = mig(Aᵢᵢ)`` and ``⟨A⟩ᵢⱼ = -mag(Aᵢⱼ)`` if ``i≠j``.

### Examples

```jldoctest
julia> A = [2..4 -1..1; -1..1 2..4]
2×2 Matrix{Interval{Float64}}:
  [2, 4]  [-1, 1]
 [-1, 1]   [2, 4]

julia> comparison_matrix(A)
2×2 Matrix{Float64}:
  2.0  -1.0
 -1.0   2.0
```
"""
function comparison_matrix(A::SMatrix{N, N, T, M}) where {N, M, T<:Interval}
    n = size(A, 1)
    compA = -mag.(A)
    @inbounds for (i, idx) in enumerate(diagind(A))
        compA = setindex(compA, mig(A[i, i]), idx)
    end
    return compA
end

function comparison_matrix(A::AbstractMatrix{T}) where {T<:Interval}
    n = size(A, 1)
    compA = -mag.(A)
    @inbounds for i in 1:n
        compA[i, i] = mig(A[i, i])
    end
    return compA
end

"""
    Orthants

Iterator to go through all the ``2ⁿ`` vectors of length ``n`` with elements ``±1``.
This is equivalento to going through the orthants of an ``n``-dimensional euclidean space.

### Fields

n::Int -- dimension of the vector space

### Example

```jldoctest
julia> for or in Orthants(2)
       @show or
       end
or = [1, 1]
or = [-1, 1]
or = [1, -1]
or = [-1, -1]
```
"""
struct Orthants
    n::Int
end

Base.eltype(::Type{Orthants}) = Vector{Int}
Base.length(O::Orthants) = 2^(O.n)

function Base.iterate(O::Orthants, state=1)
    state > 2 ^ O.n && return nothing

    vec = -2*digits(state-1, base=2, pad=O.n) .+ 1
    return (vec, state+1)
end

function Base.getindex(O::Orthants, i::Int)
    1 <= i <= length(O) || throw(BoundsError(O, i))
    return -2*digits(i-1, base=2, pad=O.n) .+ 1
end

Base.firstindex(O::Orthants) = 1
Base.lastindex(O::Orthants) = length(O)

_unchecked_interval(x::Real) = Interval(x)
_unchecked_interval(x::Complex) = Interval(real(x)) + Interval(imag(x)) * im


text:
module BEASTSimulation

export TreeDiffusionModel,
       ResidualVarianceModel,
       LatentFactorModel,
       TraitSimulationModel,
       simulate,
       get_newick,
       randomFactorSimulationModel

using PhyloNetworks, LinearAlgebra, LinearAlgebra.BLAS, DataFrames, Distributions
using UnPack
using BeastUtils.MatrixUtils, BEASTTreeUtils

function force_hermitian(X::AbstractMatrix{Float64})
    return Hermitian(X)
end

function force_hermitian(X::Diagonal{Float64})
    return X
end

function force_hermitian(X::Hermitian{Float64})
    return X
end

import Distributions.MvNormal
function MvNormal(μ::AbstractVector{Float64}, Σ::Hermitian{Float64})
    return MvNormal(μ, Symmetric(Σ))
end

abstract type ModelExtension end

struct TreeDiffusionModel{T <: Union{Diagonal{Float64}, Hermitian{Float64}},
                          S <: AbstractVector{Float64}}
    tree::HybridNetwork
    Σ::T
    μ::S

    function TreeDiffusionModel(tree::HybridNetwork, Σ::AbstractMatrix{Float64},
                μ::AbstractVector{Float64})
        Σ = force_hermitian(Σ)
        if size(Σ, 1) != length(μ)
            throw(DimensionMismatch("variance is of dimension $(size(Σ)), " *
                    "but mean is of dimension $(length(μ))."))
        end

        return new{typeof(Σ), typeof(μ)}(tree, Σ, μ)
    end
end

function TreeDiffusionModel(tree::HybridNetwork,
                            Σ::AbstractArray{Float64, 2})
    p = size(Σ, 1)
    μ = zeros(p)
    return TreeDiffusionModel(tree, Σ, μ)
end

function TreeDiffusionModel(newick::String,
                            Σ::AbstractArray{Float64, 2})
    tree = readTopology(newick)
    return TreeDiffusionModel(tree, Σ)
end

function TreeDiffusionModel(tree::Union{String, HybridNetwork}, p::Int)
    return TreeDiffusionModel(tree, Diagonal(ones(p)))
end

struct ResidualVarianceModel{T <: AbstractMatrix{Float64}} <: ModelExtension
    Γ::T # residual variance

    function ResidualVarianceModel(Γ::AbstractArray{Float64, 2})
        if !issquare(Γ)
            error("Γ must be square.")
        end
        Γ = force_hermitian(Γ)
        return new{typeof(Γ)}(Γ)
    end

end

struct LatentFactorModel{S <: AbstractMatrix{Float64},
                         T <: AbstractMatrix{Float64}} <: ModelExtension
    L::S # loadings matrix
    Λ::T # residual variance

    function LatentFactorModel(L::AbstractMatrix{Float64},
                                Λ::AbstractMatrix{Float64})

        k, p = size(L)
        if !issquare(Λ)
            error("Residual variance Λ must be square.")
        end
        if size(Λ) != (p, p)
            q, s = size(Λ)
            error("Marices are non-conformable." *
                    " Matrix L has $p columns while matrix Λ is $q x $s.")
        end

        Λ = force_hermitian(Λ)

        S = typeof(L)
        T = typeof(Λ)

        return new{S, T}(L, Λ)
    end
end

function treeDimension(x::LatentFactorModel)
    return size(x.L, 1)
end

function treeDimension(x::ResidualVarianceModel)
    return size(x.Γ, 1)
end

function dataDimension(x::LatentFactorModel)
    return size(x.L, 2)
end

function dataDimension(x::ResidualVarianceModel)
    return treeDimension(x)
end

function get_loadings(x::LatentFactorModel)
    return x.L
end

function get_loadings(x::ResidualVarianceModel)
    return I
end

function get_residual_variance(x::LatentFactorModel)
    return x.Λ
end

function randomLFM(k::Int, p::Int;
                    L_dist::UnivariateDistribution = Normal(0, 1),
                    Λ_dist::UnivariateDistribution = Gamma(2, 0.25))
    L = rand(L_dist, k, p)
    λ = rand(Λ_dist, p)
    return LatentFactorModel(L, Diagonal(λ))
end

function randomFactorSimulationModel(n::Int, k::Int, p::Int)
    tree = BEASTTreeUtils.rtree(n, ultrametric = true)
    standardize_height!(tree)
    lfm = randomLFM(k, p)
    tsm = TraitSimulationModel(tipLabels(tree), tree, lfm)
    return tsm
end

function merge_models(models::ModelExtension...)
    n = length(models)
    tree_dims = treeDimension.(models)
    data_dims = dataDimension.(models)
    K = sum(tree_dims)
    P = sum(data_dims)
    L = zeros(K, P)
    Λ = zeros(P, P)
    k_offset = 0
    p_offset = 0
    for i = 1:n
        k_end = k_offset + tree_dims[i]
        p_end = p_offset + data_dims[i]
        k_rng = (k_offset + 1):k_end
        p_rng = (p_offset + 1):p_end
        L[k_rng, p_rng] .= get_loadings(models[i])
        Λ[p_rng, p_rng] .= get_residual_variance(models[i])
        k_offset = k_end
        p_offset = p_end
    end

    return LatentFactorModel(L, Λ)
end

mutable struct TraitSimulationModel
    taxa::AbstractArray{T, 1} where T <: AbstractString
    treeModel::TreeDiffusionModel
    extensionModel::Union{Nothing, ModelExtension}

    function TraitSimulationModel(taxa::AbstractArray{T, 1} where T <: AbstractString,
                                  treeModel::TreeDiffusionModel,
                                  extensionModel::ModelExtension)
        return new(taxa, treeModel, extensionModel)
    end

    function TraitSimulationModel(taxa::AbstractArray{T, 1} where T <: AbstractString,
                                  treeModel::TreeDiffusionModel)
        return new(taxa, treeModel, nothing)
    end

    function TraitSimulationModel(taxa::AbstractArray{<:AbstractString, 1},
                                  tree::Union{String, HybridNetwork},
                                  extensionModel::ModelExtension)
        treeModel = TreeDiffusionModel(tree, treeDimension(extensionModel))
        return new(taxa, treeModel, extensionModel)
    end

    function TraitSimulationModel(tree::Union{String, HybridNetwork},
                                  extensionModel::ModelExtension)
        treeModel = TreeDiffusionModel(tree, treeDimension(extensionModel))
        return new(tipLabels(tree), treeModel, extensionModel)
    end
end

function get_tree(tsm::TraitSimulationModel)
    return tsm.treeModel.tree
end

function get_newick(tsm::TraitSimulationModel)
    return writeTopology(get_tree(tsm))
end

function get_dimension(tdm::TreeDiffusionModel)
    return length(tdm.μ)
end

function simulate_on_tree(model::TreeDiffusionModel,
                        taxa::AbstractArray{T, 1}) where T <: AbstractString
    params = PhyloNetworks.ParamsMultiBM(model.μ, Matrix(model.Σ))
    trait_sim = PhyloNetworks.simulate(model.tree, params)
    sim_taxa = trait_sim.M.tipNames
    perm = indexin(taxa, sim_taxa)
    @assert sim_taxa[perm] == taxa #TODO remove

    n = length(taxa)
    p = get_dimension(model)

    sim_Y = trait_sim[:Tips]
    Y = sim_Y'[perm, :]

    return Y
end

function add_extension(data::Matrix{Float64}, rvm::ResidualVarianceModel)

    L_chol = cholesky(rvm.Γ).L

    n, p = size(data)

    for i = 1:n
        # yi = @view data[i, :]
        # gemm!('N', 'N', 1.0, L_chol, randn(p), 1.0, yi)
        data[i, :] .+= L_chol * randn(p)
    end
    return data
end

function add_extension(data::Matrix{Float64}, lfm::LatentFactorModel)
    n, k = size(data)
    p = size(lfm.L, 2)
    Y = data * lfm.L

    res_dist = MvNormal(zeros(p), lfm.Λ)

    for i = 1:n
        Y[i, :] .+= rand(res_dist)
    end

    return Y
end

import PhyloNetworks: simulate
function simulate(tsm::TraitSimulationModel)

    data = simulate_on_tree(tsm.treeModel, tsm.taxa)

    if !isnothing(tsm.extensionModel)
        data = add_extension(data, tsm.extensionModel)
    end

    return data

end

include("distributions.jl")

end


text:
# Used to load sequential data where the time index is preceeded by 0's: 0004
function numberDigits( number; maxdigits=4 )
   ndigs = 0
   for d in 0:maxdigits
       ndigs += number ÷ 10^d > 0
   end
   nzeros = ( maxdigits - ndigs );
   return ( nzeros < 0 ) ? string(number) : "0"^(nzeros)*string(number)
end

# computing mean and standard deviation with a single loop
#
# mean = sum( a )/N
# std = sum( a - mean )^2/N = sum( a*a + mean*mean - 2*a*mean )/N =
#     	                      sum( a*a )/N + N/N*mean^2 - 2*mean*sum( a )/N =
#                             sum( a*a )/N + mean^2 - 2*mean^2

function meanStd( a::AbstractArray{<:Real,N} ) where {N}

	len   = length(a)
	suma  = 0
	sumaa = 0
	@inbounds @simd for idx in 1:len
		suma  += a[idx]
		sumaa += a[idx]*a[idx]
	end
	mean = suma/len;
	return mean, sqrt( sumaa/len + mean*mean - 2*mean*mean )
end

function vectorFieldSize( datasize, pp::PIVParameters )

    step = pp.interSize .- pp.overlap
    return length.( StepRange.( pp.interSize, step, datasize ) )
end

# optimized to use SIMD, performing 4 comparisons at the same time
function maxval( a::AbstractArray{<:Real,N} ) where {N}

	len = length(a)

	if len < 4
		return maximum( a )
	end

	simdstart = len%4 + 1

    m1 = a[1]
    m2 = a[2]
    m3 = a[3]
    m4 = a[4]
	@inbounds @simd for i in simdstart:4:len
        m1 = ( a[ i ] > m1 ) ? a[ i ] : m1
        m2 = ( a[i+1] > m2 ) ? a[i+1] : m2
        m3 = ( a[i+2] > m3 ) ? a[i+2] : m3
        m4 = ( a[i+3] > m4 ) ? a[i+3] : m4
	end

    return maximum( [ m1, m2, m3, m4 ] )
end

# optimized to use SIMD, performing 4 comparisons at the same time
function maxidx( a::AbstractArray{<:Real,N} ) where {N}

	len = length(a)

	if len < 4
		val, idx = findmax( a )
		return idx
	end

	simdstart = len%4 + 1

    maxidx1 = 1
    maxidx2 = 2
    maxidx3 = 3
    maxidx4 = 4
	@inbounds @simd for i in simdstart:4:len
        maxidx1 = ( a[ i ] > a[maxidx1] ) ? i   : maxidx1
        maxidx2 = ( a[i+1] > a[maxidx2] ) ? i+1 : maxidx2
        maxidx3 = ( a[i+2] > a[maxidx3] ) ? i+2 : maxidx3
        maxidx4 = ( a[i+3] > a[maxidx4] ) ? i+3 : maxidx4
	end

    # Finding the max value among the maximum indices computed in parallel
    val, idx = findmax( [  a[ maxidx1 ], a[ maxidx2 ], a[ maxidx3], a[ maxidx4 ] ] );
    maxindex = (idx==1)*maxidx1 + (idx==2)*maxidx2 + (idx==3)*maxidx3 + (idx==4)*maxidx4;

	return maxindex
end

function firstPeak( cmat::Array{<:AbstractFloat,N} ) where {N}

    maxindex = maxidx( cmat )

    # Transforming from linear indexing to cartesian
	h, w, d = size(cmat,1), size(cmat,2), size(cmat,3)

    z = ceil( Int32, maxindex/(h*w) )
    x = ceil( Int32, (maxindex - (z-1)*h*w)/h )
    y = maxindex - (x-1)*h - (z-1)*h*w;

    return (y,x,z), cmat[ maxindex ]
end

function secondPeak( cmat::A{<:AbstractFloat,N}, width::I, idx::NTuple{N,I}) where {N}

	#  First, we store the original values around the maximum peak, and then change them to -1.0.
    ranges = [ max( 1, idx[d] - width ):min( size(cmat,d), idx[d] + width ) for d in 1:N ];
    OGvals = copy( cmat[ ranges... ] );
    cmat[ ranges... ] .= eltype(cmat)(-1.0);

	# Finding the second peak and setting the values back to normal (out of politeness)
    p2, val = firstPeak( cmat );
    cmat[ ranges... ] .= OGvals

    return p2, val
end

setTo0!( A::Array{T,N}... ) where {T<:Number,N} = setTo0!.( A )

function setTo0!( A::Array{T,N} ) where {T<:Number,N}
    zer = T(0)
    @inbounds @simd for i in 1:length(A)
        A[ i ] = zer
    end
end

function SAcoords( IAcoords1, IAcoords2, disp, SM_mp, datasize )

	SAcoords1  = max.(    1    , IAcoords1 .+ disp .- SM_mp )
	SAcoords2  = min.( datasize, IAcoords2 .+ disp .+ SM_mp )
	pad_offset = max.( 0, SAcoords1 .- ( IAcoords1 .+ disp .- SM_mp ) )
	y1, x1, z1 = SAcoords1
	y2, x2, z2 = SAcoords2

    return (y1, y2, x1, x2, z1, z2), pad_offset
end

function putWithinSearch!( search::Array{T,N}, data::Array{T,N},
                           so::NTuple{3,I}, coord::NTuple{6,I}
                         ) where {T<:AbstractFloat,N}

    start  = [ x + 1 for x in so ];
    r1, r2 = (coord[1], coord[2])
    c1, c2 = (coord[3], coord[4])
    d1, d2 = (coord[5], coord[6])

	if N == 2
		start[3] = 1;
		d1, d2 = 1,1;
	end

    @inbounds for z in 0:(d2-d1), c in 0:(c2-c1)
		@simd for r in 0:(r2-r1)
		    search[start[1]+r, start[2]+c, start[3]+z] = data[r1+r, c1+c, d1+z]
		end
    end
end

function putWithinPadded!( padA::Array{Complex{T},N}, A::Array{<:AbstractFloat,N}, mean,
                           offset::NTuple{3,I}, coords::NTuple{6,I}
                         )  where {T<:AbstractFloat,N}

    o1,o2,o3 = offset .+ ( 1, 1, 1 )
    r1, r2 = (coords[1], coords[2])
    c1, c2 = (coords[3], coords[4])
    d1, d2 = (coords[5], coords[6])

    Tzero = T(0)
    Tmean = convert( T, mean )

    for z in 0:(d2-d1), c in 0:(c2-c1)
        @simd for r in 0:(r2-r1)
            padA[ o1+r, o2+c, o3+z ] = complex( convert( T, A[r1+r,c1+c,d1+z] ) - Tmean, Tzero )
        end
    end
end

# construc the integral array of squared values
function integralArraySQ!( padS::AbstractArray{Complex{T},2}, intArr::Array{T,2}
                         ) where {T<:AbstractFloat}

	h, w = size(intArr) .- 1
	@inbounds for c in 2:w+1, r in 2:h+1
		intArr[r,c] = padS[r-1,c-1]^2 + intArr[r-1,c] + intArr[r,c-1] - intArr[r-1,c-1]
	end
end

# optimized sumarea, reusing previous operations
function integralArraySQ!( padS::AbstractArray{Complex{T},3}, intArr::Array{T,3}
                         ) where {T<:AbstractFloat}

    h, w, d = size(intArr) .- 1;
    @inbounds for z in 2:d+1, c in 2:w+1
        tmp = 0.0;
        for r in 2:h+1
            val2          = real( padS[r-1,c-1,z-1] )^2
            intArr[r,c,z] = val2 + intArr[r,c-1,z] + intArr[r,c,z-1] - intArr[r,c-1,z-1] + tmp;
            tmp          += val2;
    end end
end

function integralArea( intArr::AbstractArray{<:Real,2}, TL, BR )
	TL   = TL .+ 1;
	BR   = BR .+ 1;
	area = intArr[BR[1],BR[2]] - intArr[BR[1],TL[2]] - intArr[TL[1],BR[2]] + intArr[TL[1],TL[2]]
end

function integralArea( intArr::AbstractArray{<:Real,3}, TL, BR )
	TL = TL .+ 1;
	BR = BR .+ 1;
	area  = intArr[BR[1],BR[2],BR[3]] - intArr[TL[1],TL[2],TL[3]]
	area -= intArr[TL[1],BR[2],BR[3]] + intArr[BR[1],TL[2],BR[3]] + intArr[BR[1],BR[2],TL[3]]
	area += intArr[BR[1],TL[2],TL[3]] + intArr[TL[1],BR[2],TL[3]] + intArr[TL[1],TL[2],BR[3]]
    return area
end

max_issues_repo_name: ErichZimmer/quickPIV

max_forks_repo_name: ErichZimmer/quickPIV


text:
using Documenter
push!(LOAD_PATH, joinpath(@__DIR__, "..", "src"))
using Purses

DocMeta.setdocmeta!(Purses, :DocTestSetup, :(using Purses); recursive=true)
makedocs(
    sitename = "Purses.jl 👛",
    modules = [Purses],
    pages = [
        "Index" => "index.md",
        "introduction.md",
        "api.md",
    ],
)

deploydocs(
    repo = "github.com/dalum/Purses.jl.git",
)


text:
using Test
using Random
using JuliaDB
using Garamond

tbl = table((y = rand(10),
             x = collect(1:10),
             z = [randstring(1) for _ in 1:10]);
            pkey=:x)

env = (dbdata=tbl, id_key=:x, config_path=nothing)

QUERIES = [(input_query = "this is a query", filter_response = Dict{Symbol,Any}(), search_response = "this is a query"),
           (input_query = "search token x:10 y:1 z:\"a\"", filter_response = Dict(:x=>10, :y=>1.0, :z=>"a"), search_response = "search token"),
           (input_query = "x:[1,2,3] y:(\"a\",) z:", filter_response = Dict(:x=>[1,2,3], :y=>("a",)), search_response = ""),
           (input_query = "x:[1.0,2,3] y:(\"a\",) z:", filter_response = Dict(:x=>[1.0,2,3], :y=>("a",)), search_response = ""),
           (input_query = "z:\"a b c\" x:1 free form", filter_response = Dict(:z=>"a b c", :x=>1), search_response = "free form"),
           (input_query = "", filter_response=Dict(), search_response=""),
]

@testset "Input: noop_input_parser" begin
	for query in QUERIES
        fake_request = (query=query.input_query,)
		sq, filt = Garamond.noop_input_parser(env, fake_request)
        @test filt == Dict{Symbol, Any}()
		@test sq == query.input_query
	end
end

@testset "Input: base_input_parser" begin
	for query in QUERIES
        fake_request = (query=query.input_query, )
		sq, filt = Garamond.base_input_parser(env, fake_request; separator=":")
		@test filt == query.filter_response
		@test sq == query.search_response
	end
end

@testset "Input: pre_parser" begin
    PREPARSER_DATA = [("noop_input_parser>", Garamond.noop_input_parser),
                      ("base_input_parser>", Garamond.base_input_parser),
                      ("pre_parser> base_input_parser>", Garamond.pre_parser),
                      ("pre_parser> pre_parser  >", Garamond.pre_parser)]

    for (parser_specification, parser_function) in PREPARSER_DATA
        for query in QUERIES
            preparser_query = parser_specification * query.input_query
            fake_request = Garamond.InternalRequest(query=preparser_query,
                                                    input_parser=:pre_parser,
                                                    searchable_filters=[])
            sq, filt = parse_input(env, fake_request)
            if occursin("base_input_parser", parser_specification)
                # final parser is the base_input_parser
                @test filt == query.filter_response
                @test sq == query.search_response
            elseif occursin("noop_input_parser", parser_specification)
                # final parser is the noop_parser
                @test filt == Dict{Symbol, Any}()
		        @test sq == query.input_query
            else
                # Using DEFAULT_INPUT_PARSER_NAME as end parser
                # (it is already tested), just test that parsing works
                @test filt isa Dict{Symbol, Any}
                @test sq isa String
            end
        end
    end
end


text:
immutable MixtureModel <: Distribution
    components::Vector # Vector should be able to contain any type of
                       # distribution with comparable support
    probs::Vector{Float64}
    drawtable::DiscreteDistributionTable
    function MixtureModel(c::Vector, p::Vector{Float64})
        if length(c) != length(p)
            error("components and probs must have the same number of elements")
        end
        sump = 0.0
        for i in 1:length(p)
            if p[i] < 0.0
                error("MixtureModel: probabilities must be non-negative")
            end
            sump += p[i]
        end
        table = DiscreteDistributionTable(p)
        new(c, p ./ sump, table)
    end
end

function mean(d::MixtureModel)
    m = 0.0
    for i in 1:length(d.components)
        m += mean(d.components[i]) * d.probs[i]
    end
    return m
end

function pdf(d::MixtureModel, x::Any)
    p = 0.0
    for i in 1:length(d.components)
        p += pdf(d.components[i], x) * d.probs[i]
    end
    return p
end

function rand(d::MixtureModel)
    i = draw(d.drawtable)
    return rand(d.components[i])
end

# TODO: Correct this definition
function var(d::MixtureModel)
    m = 0.0
    for i in 1:length(d.components)
        m += var(d.components[i]) * d.probs[i]^2
    end
    return m
end


text:
"""

right-hand side associated with

∂t u = f(u,Θ), where

u = (x,l,v,hj) with

x - particle position
l - log determinant
v - accumulated transport costs
hj - accumulates least-squares error of -∂Φ(x,t) + 0.5 * |∇Φ(x,t)|^2 along trajectory
"""
function odefun(J,U::AbstractArray{R},Θ,t::R) where R <: Real
    (d,nex) = size(U)
    d = d-4;
    XT = [ U[1:d,:]; fill(t,1,nex)]
    Phi = J.Φ(XT,Θ) # run fwd prop to populate N.tmp
    gradPhi = getGradPotential(J.Φ,XT,Θ)
    trH = getTraceHess(J.Φ,XT,Θ)
    dx = -(1/J.α[1])*gradPhi[1:d,:]
    dl = -(1/J.α[1])*trH
    dv = R(0.5).*sum(dx.^2,dims=1)
    df = reshape(J.F(U,t),1,nex)
    hj = abs.(-reshape(gradPhi[end,:],1,:) -
                reshape(J.α[2].*getDeltaF(J.F,U,t),1,:) +
                J.α[1] .* dv)
    return [dx;reshape(dl,1,nex);dv;df;hj]
end

max_stars_repo_name: EmoryMLIP/MFGnet

max_issues_repo_name: EmoryMLIP/MFGnet

max_forks_repo_name: EmoryMLIP/MFGnet


text:
using Planck
using Test
using Unitful

@testset "Planck.jl" begin
    @testset "Consistent unit conversions" begin
        @test ustrip(u"W/m^3", blackbody(545u"nm", 6000u"K")) ≈ blackbody(545e-9, 6000)
        @test ustrip(u"W/m^2/Hz", blackbody(500u"THz", 6000u"K")) ≈ 3.441e-8 rtol=1e-3
    end
    @testset "struct interface" begin
        bb = Blackbody(6000)
        @test bb(545e-9) ≈ 3.079e13 rtol=1e-3
        @test bb(Float32, 545e-9) ≈ 3.079f13 rtol=1f-3

        qbb = Blackbody(6000u"K")
        @test qbb(545u"nm") ≈ 3.079e13u"W/m^3" rtol=1e-3
        @test qbb(500u"THz") ≈ 3.441e-8u"W/m^2/Hz" rtol=1e-3
    end
end


text:
module PosteriorFactorTest

using Test
using ForneyLab

using ForneyLab: nodesConnectedToExternalEdges, Cluster, condense, flatten, setTargets!

@testset "PosteriorFactor" begin
    g = FactorGraph()
    @RV m ~ GaussianMeanVariance(constant(0.0), constant(1.0))
    @RV w ~ Gamma(constant(1.0), constant(1.0))
    y = Variable[]
    for i = 1:3
        @RV y_i ~ GaussianMeanPrecision(m, w)
        placeholder(y_i, :y, index=i)
        push!(y, y_i)
    end

    pfz = PosteriorFactorization()
    q_m = PosteriorFactor(m)
    @test q_m.id == :posteriorfactor_1
    @test q_m.internal_edges == edges(m)
    @test pfz.posterior_factors[:posteriorfactor_1] === q_m

    q_w = PosteriorFactor(w)
    @test q_w.id == :posteriorfactor_2
    @test q_w.internal_edges == edges(w)
    @test pfz.posterior_factors[:posteriorfactor_2] === q_w

    # Joint factorizations
    q_m_w = PosteriorFactor([m, w])
    @test q_m_w.id == :posteriorfactor_3
    @test q_m_w.internal_edges == edges(Set([m, w]))
    @test pfz.posterior_factors[:posteriorfactor_3] === q_m_w

    q_y = PosteriorFactor(y)
    @test q_y.id == :posteriorfactor_4
    @test q_y.internal_edges == edges(Set(y))
    @test pfz.posterior_factors[:posteriorfactor_4] === q_y
end

@testset "setTargets!" begin
    fg = FactorGraph()
    @RV x ~ GaussianMeanPrecision(0.0, 1.0)
    @RV y ~ GaussianMeanPrecision(0.0, 1.0)
    @RV z = x + y
    @RV w ~ GaussianMeanPrecision(z, 1.0)
    placeholder(w, :w)

    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    setTargets!(pf, pfz, target_variables=Set{Variable}([z]))
    @test pfz.free_energy_flag == false
    @test pf.target_variables == Set{Variable}([z])
    @test pf.target_clusters == Set{Cluster}()

    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    setTargets!(pf, pfz, target_variables=Set{Variable}([z]), external_targets=true)
    @test pfz.free_energy_flag == false
    @test pf.target_variables == Set{Variable}([z])
    @test pf.target_clusters == Set{Cluster}()

    pfz = PosteriorFactorization()
    pf = PosteriorFactor(fg)
    setTargets!(pf, pfz, target_variables=Set{Variable}([z]), free_energy=true)
    @test pfz.free_energy_flag == true
    @test pf.target_variables == Set{Variable}([x, y, z])
    @test length(pf.target_clusters) == 1
    @test first(pf.target_clusters).id == :x_y
end

end # module


text:
module GraphQLClient

using GraphQLParser
using HTTP
using JSON3
using StructTypes

export query, mutate, open_subscription, Client, GQLEnum, Alias,
    full_introspection!, get_queries, get_mutations, get_subscriptions,
    introspect_object, get_introspected_type, initialise_introspected_struct,
    create_introspected_struct, list_all_introspected_objects, global_graphql_client,
    @gql_str

# Types
include("client.jl")
include("types.jl")

include("logging.jl")
include("gqlresponse.jl")
include("exceptions.jl")
include("schema_utils.jl")
include("args.jl")
include("variables.jl")
include("output_fields.jl")
include("type_construction.jl")
include("http_execution.jl")
include("queries.jl")
include("mutations.jl")
include("subscriptions.jl")
include("introspection.jl")
include("gql_string.jl")

end # module


text:
using HypothesisTests, Test
using Distributions

@testset "Wald-Wolfowitz test" begin

@testset "Many-Valued Observations" begin
    # Get set of dependent observations
    x = 1:1000
    tst = WaldWolfowitzTest(x)

    @test tst.z ≈ -31.575 atol=1e-3
    @test pvalue(tst) ≈ 8.052e-219 atol=1e-222

    # Test consistency of z-statistics
    @test pvalue(tst) == pvalue(Normal(tst.μ, tst.σ), tst.nruns)
    @test pvalue(tst, tail=:left) == pvalue(Normal(tst.μ, tst.σ), tst.nruns, tail=:left)
    @test pvalue(tst, tail=:right) == pvalue(Normal(tst.μ, tst.σ), tst.nruns, tail=:right)
    expected_output = """
    Wald-Wolfowitz Test
    -------------------
    Population details:
        parameter of interest:   Number of runs
        value under h_0:         501.0
        point estimate:          2

    Test summary:
        outcome with 95% confidence: reject h_0
        two-sided p-value:           <1e-99

    Details:
        number of runs:  2
        z-statistic:     -31.575338477995764
    """
    output = sprint(show, tst)
    @test output == expected_output
end

@testset "Two-Valued Observations" begin
    # equivalent data as above (half under median half over)
    x = [falses(500); trues(500)]
    tst = WaldWolfowitzTest(x)

    @test tst.z ≈ -31.575 atol=1e-3
    @test pvalue(tst) ≈ 8.052e-219 atol=1e-222

    # Test consistency of z-statistics
    @test pvalue(tst) == pvalue(Normal(tst.μ, tst.σ), tst.nruns)
    @test pvalue(tst, tail=:left) == pvalue(Normal(tst.μ, tst.σ), tst.nruns, tail=:left)
    @test pvalue(tst, tail=:right) == pvalue(Normal(tst.μ, tst.σ), tst.nruns, tail=:right)
    expected_output = """
    Wald-Wolfowitz Test
    -------------------
    Population details:
        parameter of interest:   Number of runs
        value under h_0:         501.0
        point estimate:          2

    Test summary:
        outcome with 95% confidence: reject h_0
        two-sided p-value:           <1e-99

    Details:
        number of runs:  2
        z-statistic:     -31.575338477995764
    """
    output = sprint(show, tst)
    @test output == expected_output
end

end


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

# Compute X <- a X + b Y.
function matfun_axpby!(X,a,b,Y)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:m
        @inbounds for j=1:n
            if (b isa ValueOne)
                X[i,j]+=Y[i,j]
            else
                X[i,j]+=b*Y[i,j]
            end
        end
    end
end

@inline function exp_native_jl_rho2_1(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_native_jl_rho2_1!(A_copy)
end

@inline function exp_native_jl_rho2_1!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=7
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: A2 A4 A6 A8 Ua U V Z X P
    # Computing A2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing A4 with operation: mult
    mul!(memslots3,memslots2,memslots2)
    # Computing A6 with operation: mult
    mul!(memslots4,memslots2,memslots3)
    # Computing A8 with operation: mult
    mul!(memslots5,memslots2,memslots4)
    # Computing Ua = x*I+x*A2+x*A4+x*A6+x*A8
    coeff1=8.8216128e9
    coeff2=3.027024e8
    coeff3=2.16216e6
    coeff4=3960.0
    coeff5=1.0
    memslots6 .= coeff2.*memslots2 .+ coeff3.*memslots3 .+ coeff4.*memslots4 .+ coeff5.*memslots5
    mul!(memslots6,true,I*coeff1,true,true)
    # Computing U with operation: mult
    mul!(memslots7,memslots6,memslots1)
    # Deallocating Ua in slot 6
    # Deallocating A in slot 1
    # Computing V = x*I+x*A2+x*A4+x*A6+x*A8
    coeff1=1.76432256e10
    coeff2=2.0756736e9
    coeff3=3.027024e7
    coeff4=110880.0
    coeff5=90.0
    # Smart lincomb recycle A2
    memslots2 .= coeff2.*memslots2 .+ coeff3.*memslots3 .+ coeff4.*memslots4 .+ coeff5.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Deallocating A4 in slot 3
    # Deallocating A6 in slot 4
    # Deallocating A8 in slot 5
    # Computing Z = x*V+x*U
    coeff1=1.0
    coeff2=-1.0
    memslots1 .= coeff1.*memslots2 .+ coeff2.*memslots7
    # Computing X = x*V+x*U
    coeff1=1.0
    coeff2=1.0
    # Smart lincomb recycle V
    memslots2 .= coeff1.*memslots2 .+ coeff2.*memslots7
    # Deallocating U in slot 7
    # Computing P with operation: ldiv
    memslots3 .=memslots1\memslots2
    # Deallocating Z in slot 1
    # Deallocating X in slot 2
    return memslots3 # Returning P
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true)
@testset "test dc tnep" begin
    @testset "6-bus case" begin
        resultDC = run_tnepopf("../test/data/tnep/case6_test.m", DCPPowerModel, cbc; setting = s)
        @test isapprox(resultDC["objective"], 26.3331; atol = 1e-1)
        @test isapprox(resultDC["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultDC["solution"]["convdc_ne"]["5"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultDC["solution"]["convdc_ne"]["6"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultDC["solution"]["branchdc_ne"]["2"]["isbuilt"], 0; atol = 1e-2)
        @test isapprox(resultDC["solution"]["branchdc_ne"]["10"]["isbuilt"], 0; atol = 1e-2)
        @test isapprox(resultDC["solution"]["branchdc_ne"]["1"]["pf"], 1.3737; atol = 1e-2)
    end

    @testset "9-bus case" begin
        resultDC = run_tnepopf("../test/data/tnep/case9_test.m", DCPPowerModel, cbc; setting = s)
        @test isapprox(resultDC["objective"], 10.7; atol = 1e-1)
        @test isapprox(resultDC["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
    end
    #
    @testset "14-bus case" begin
        resultDC =  run_tnepopf("../test/data/tnep/case14_test.m",  DCPPowerModel, cbc, setting = s)
        @test isapprox(resultDC["objective"], 15.6921; atol = 1e-1)
        @test isapprox(resultDC["solution"]["convdc_ne"]["6"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultDC["solution"]["convdc_ne"]["7"]["isbuilt"], 0; atol = 1e-2)
    end
    #
    @testset "39-bus case" begin
        resultDC =  run_tnepopf("../test/data/tnep/case39_test.m",  DCPPowerModel, cbc, setting = s)
        @test isapprox(resultDC["objective"], 25.1605; atol = 1e-1)
        @test isapprox(resultDC["solution"]["convdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultDC["solution"]["convdc_ne"]["4"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultDC["solution"]["convdc_ne"]["5"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultDC["solution"]["branchdc_ne"]["13"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultDC["solution"]["branchdc_ne"]["17"]["isbuilt"], 1; atol = 1e-2)

    end
end

if local_test == true
    @testset "test soc tnep" begin
        @testset "test SOC-BFM tnep" begin
            @testset "6-bus case" begin
                resultSOCBF = run_tnepopf_bf("../test/data/tnep/case6_test.m", SOCBFConicPowerModel, mosek; setting = s)
                @test isapprox(resultSOCBF["objective"], 31.6; atol = 1e-1)
                @test isapprox(resultSOCBF["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["convdc_ne"]["4"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["convdc_ne"]["5"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["convdc_ne"]["6"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["branchdc_ne"]["3"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["branchdc_ne"]["10"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["branchdc_ne"]["3"]["pf"], -2.29; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["busdc_ne"]["2"]["wdc_ne"], 1.18695; atol = 1e-2)
            end
            @testset "9-bus case" begin
                resultSOCBF = run_tnepopf_bf("../test/data/tnep/case9_test.m", SOCBFConicPowerModel, mosek; setting = s)
                @test isapprox(resultSOCBF["objective"], 10.7; atol = 1e-1)
                @test isapprox(resultSOCBF["solution"]["branchdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["convdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["branchdc_ne"]["1"]["pf"], 1.029695717391926; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["busdc_ne"]["2"]["wdc_ne"], 0.9607722492758961; atol = 1e-2)
            end
        end
        @testset "14-bus case" begin
            resultSOCBF =  run_tnepopf_bf("../test/data/tnep/case14_test.m",  SOCBFConicPowerModel, mosek; setting = s)
            @test isapprox(resultSOCBF["objective"], 20; atol = 1e-1)
            @test isapprox(resultSOCBF["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCBF["solution"]["convdc_ne"]["5"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCBF["solution"]["convdc_ne"]["7"]["isbuilt"], 1; atol = 1e-2)
        end
        #
        @testset "39-bus case" begin
            resultSOCBF =  run_tnepopf_bf("../test/data/tnep/case39_test.m",  SOCBFConicPowerModel, mosek; setting = s)
            @test isapprox(resultSOCBF["objective"], 30.4; atol = 1e-1)
            @test isapprox(resultSOCBF["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCBF["solution"]["convdc_ne"]["4"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCBF["solution"]["convdc_ne"]["5"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCBF["solution"]["convdc_ne"]["6"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCBF["solution"]["branchdc_ne"]["5"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCBF["solution"]["branchdc_ne"]["6"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCBF["solution"]["branchdc_ne"]["29"]["isbuilt"], 1; atol = 1e-2)
        end

    end

    @testset "test SOC-BIM tnep" begin
        @testset "6-bus case" begin
            resultSOCWR = run_tnepopf("../test/data/tnep/case6_test.m", SOCWRPowerModel, gurobi; setting = s)
            @test isapprox(resultSOCWR["objective"], 31.63; atol = 1e-1)
            @test isapprox(resultSOCWR["solution"]["convdc_ne"]["1"]["isbuilt"], 0; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["convdc_ne"]["4"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["convdc_ne"]["5"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["convdc_ne"]["6"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["branchdc_ne"]["3"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["branchdc_ne"]["10"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["branchdc_ne"]["3"]["pf"], -2.29067; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["busdc_ne"]["2"]["wdc_ne"], 1.18695; atol = 1e-2)
        end
        @testset "9-bus case" begin
            resultSOCWR = run_tnepopf("../test/data/tnep/case9_test.m", SOCWRPowerModel, gurobi; setting = s)
            @test isapprox(resultSOCWR["objective"], 10.7; atol = 1e-1)
            @test isapprox(resultSOCWR["solution"]["branchdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["convdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["branchdc_ne"]["1"]["pf"], 0.9468396226925635; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["busdc_ne"]["2"]["wdc_ne"], 0.8653240759694054; atol = 1e-2)

        end
    end

    @testset "test QC tnep" begin
        @testset "6-bus case" begin
            resultQC = run_tnepopf("../test/data/tnep/case6_test.m", QCRMPowerModel, gurobi; setting = s)
            @test isapprox(resultQC["objective"], 31.6; atol = 1e-1)
            @test isapprox(resultQC["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultQC["solution"]["convdc_ne"]["4"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultQC["solution"]["convdc_ne"]["5"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultQC["solution"]["convdc_ne"]["6"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultQC["solution"]["branchdc_ne"]["3"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultQC["solution"]["branchdc_ne"]["10"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultQC["solution"]["branchdc_ne"]["3"]["pf"], -2.29; atol = 1e-2)
            @test isapprox(resultQC["solution"]["busdc_ne"]["2"]["wdc_ne"], 1.18695; atol = 1e-2)
        end
        @testset "9-bus case" begin
            resultQC = run_tnepopf("../test/data/tnep/case9_test.m", QCRMPowerModel, gurobi; setting = s)
            @test isapprox(resultQC["objective"], 10.7; atol = 1e-1)
            @test isapprox(resultQC["solution"]["branchdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultQC["solution"]["convdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultQC["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultQC["solution"]["branchdc_ne"]["1"]["pf"], 0.9517040201768898; atol = 1e-2)
            @test isapprox(resultQC["solution"]["busdc_ne"]["2"]["wdc_ne"], 0.8510801499908695; atol = 1e-2)
        end
    end
end

@testset "test AC tnep" begin
    @testset "6-bus case" begin
        # ac solver huristics may build ac branch 2 or 10
        resultACP = run_tnepopf("../test/data/tnep/case6_test.m", ACPPowerModel, juniper; setting = s)
        @test isapprox(resultACP["objective"], 31.65; atol = 5e0)
        @test isapprox(resultACP["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultACP["solution"]["convdc_ne"]["4"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultACP["solution"]["convdc_ne"]["5"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultACP["solution"]["convdc_ne"]["6"]["isbuilt"], 1; atol = 1e-2)
        #@test isapprox(resultACP["solution"]["branchdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultACP["solution"]["branchdc_ne"]["6"]["isbuilt"], 0; atol = 1e-2)
        #@test isapprox(resultACP["solution"]["branchdc_ne"]["10"]["isbuilt"], 0; atol = 1e-2)
        #@test isapprox(resultACP["solution"]["branchdc_ne"]["6"]["pf"], -1.9834; atol = 1e-2)
        @test isapprox(resultACP["solution"]["busdc_ne"]["6"]["vm"], 1.09; atol = 1e-1)
    end
    @testset "9-bus case" begin
        resultACP = run_tnepopf("../test/data/tnep/case9_test.m", ACPPowerModel, juniper; setting = s)
        @test isapprox(resultACP["objective"], 10.7; atol = 1e-1)
        @test isapprox(resultACP["solution"]["branchdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultACP["solution"]["convdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultACP["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultACP["solution"]["branchdc_ne"]["1"]["pf"], 0.838; atol = 1e-2)
        @test isapprox(resultACP["solution"]["busdc_ne"]["2"]["vm"], 0.99; atol = 1e-2)
    end
end

@testset "test LPAC tnep" begin
    @testset "6-bus case" begin
        resultLPAC = run_tnepopf("../test/data/tnep/case6_test.m", LPACCPowerModel, juniper; setting = s)
        @test isapprox(resultLPAC["objective"], 30.5; atol = 1e-1)
        @test isapprox(resultLPAC["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["convdc_ne"]["4"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["convdc_ne"]["5"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["convdc_ne"]["6"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["branchdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["branchdc_ne"]["3"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["branchdc_ne"]["1"]["pf"], 1.63; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["busdc_ne"]["2"]["phivdcm_ne"], -0.00668065; atol = 1e-2)

    end
    @testset "9-bus case" begin
        resultLPAC = run_tnepopf("../test/data/tnep/case9_test.m", LPACCPowerModel, juniper; setting = s)
        @test isapprox(resultLPAC["objective"], 10.7; atol = 1e-1)
        @test isapprox(resultLPAC["solution"]["branchdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["convdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["branchdc_ne"]["1"]["pf"], 0.814; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["busdc_ne"]["2"]["phivdcm_ne"], -0.00345834; atol = 1e-2)
    end
end

@testset "ACDC tnep" begin
    @testset "DCP" begin
        resultDCP = run_acdctnepopf("../test/data/tnep/case4_acdc.m", DCPPowerModel, cbc; setting = s)
        @test isapprox(resultDCP["objective"], 303.9545; atol = 1e-1)
        @test isapprox(resultDCP["solution"]["branchdc_ne"]["3"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultDCP["solution"]["branchdc_ne"]["3"]["pf"], -1.009; atol = 1e-2)
        @test isapprox(resultDCP["solution"]["branchdc_ne"]["3"]["pf"], -1.009; atol = 1e-2)
        @test isapprox(resultDCP["solution"]["convdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultDCP["solution"]["convdc_ne"]["1"]["pconv"], -1; atol = 1e-2)
        @test isapprox(resultDCP["solution"]["ne_branch"]["1"]["built"], 1; atol = 1e-2)
    end
    @testset "ACP" begin
        resultACP = run_acdctnepopf("../test/data/tnep/case4_acdc.m", ACPPowerModel, juniper; setting = s)
        @test isapprox(resultACP["objective"], 321.00; atol = 1e0)
        #println(resultACP["objective"])
        #@test isapprox(resultACP["solution"]["branchdc_ne"]["3"]["isbuilt"], 1; atol = 1e-2)
        #@test isapprox(resultACP["solution"]["branchdc_ne"]["3"]["pf"], -0.1839; atol = 1e-2)
        @test isapprox(resultACP["solution"]["convdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultACP["solution"]["convdc_ne"]["1"]["pconv"], -0.860; atol = 1e-1)
        #println(resultACP["solution"]["convdc_ne"]["1"]["pconv"])
        @test isapprox(resultACP["solution"]["ne_branch"]["1"]["built"], 1; atol = 1e-2)
    end
    @testset "LPAC" begin
        resultLPAC = run_acdctnepopf("../test/data/tnep/case4_acdc.m", LPACCPowerModel, juniper; setting = s)
        @test isapprox(resultLPAC["objective"], 307.075; atol = 1e-1)
        @test isapprox(resultLPAC["solution"]["branchdc_ne"]["3"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["branchdc_ne"]["3"]["pf"], -1.2466;atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["convdc_ne"]["2"]["pconv"], -1.2355; atol = 1e-2)
        @test isapprox(resultLPAC["solution"]["ne_branch"]["1"]["built"], 1; atol = 1e-2)
    end
    if local_test == true
        @testset "SOCWR" begin
            resultSOCWR = run_acdctnepopf("../test/data/tnep/case4_acdc.m", SOCWRPowerModel, gurobi; setting = s)
            @test isapprox(resultSOCWR["objective"], 320.952; atol = 1e-1)
            @test isapprox(resultSOCWR["solution"]["branchdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["branchdc_ne"]["2"]["pf"], -0.8945; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["convdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["convdc_ne"]["1"]["pconv"], -0.8819; atol = 1e-2)
            @test isapprox(resultSOCWR["solution"]["ne_branch"]["1"]["built"], 1; atol = 1e-2)
        end
    end
end

s = Dict("output" => Dict("branch_flows" => true), "conv_losses_mp" => true, "process_data_internally" => false)
@testset "mp-tnep" begin
    @testset "dc tnep" begin
        @testset "DCP" begin
            data_dc = build_mn_data("../test/data/tnep/case4_original.m")
            resultDCP = run_mp_tnepopf(data_dc, DCPPowerModel, cbc, multinetwork=true; setting = s)
            @test isapprox(resultDCP["objective"], 8.2; atol = 1e-1)
            @test isapprox(resultDCP["solution"]["nw"]["1"]["branchdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultDCP["solution"]["nw"]["1"]["branchdc_ne"]["2"]["pf"], -2.0013; atol = 1e-2)
            @test isapprox(resultDCP["solution"]["nw"]["2"]["branchdc_ne"]["2"]["pf"], -2.00137; atol = 1e-2)
            @test isapprox(resultDCP["solution"]["nw"]["2"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultDCP["solution"]["nw"]["1"]["convdc_ne"]["2"]["pconv"], -2.0; atol = 1e-2)
        end
        @testset "LPAC" begin
            data_dc = build_mn_data("../test/data/tnep/case4_original.m")
            resultLPAC = run_mp_tnepopf(data_dc, LPACCPowerModel, juniper, multinetwork=true; setting = s)
            @test isapprox(resultLPAC["objective"], 10.2; atol = 1e-1)
            @test isapprox(resultLPAC["solution"]["nw"]["1"]["branchdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultLPAC["solution"]["nw"]["1"]["branchdc_ne"]["2"]["pf"], -1.2868; atol = 1e-2)
            @test isapprox(resultLPAC["solution"]["nw"]["2"]["branchdc_ne"]["2"]["pf"], -1.2868; atol = 1e-2)
            @test isapprox(resultLPAC["solution"]["nw"]["2"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultLPAC["solution"]["nw"]["1"]["convdc_ne"]["2"]["pconv"], -2.15833; atol = 1e-2)
        end
        if local_test == true
            @testset "ACP" begin   #Normally fine, but fails on travis.
                data_dc = build_mn_data("../test/data/tnep/case4_original.m")
                resultACP = run_mp_tnepopf(data_dc, ACPPowerModel, juniper, multinetwork=true; setting = s)
                @test isapprox(resultACP["objective"], 10.2; atol = 1e-1)
                @test isapprox(resultACP["solution"]["nw"]["1"]["branchdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultACP["solution"]["nw"]["1"]["branchdc_ne"]["2"]["pf"], -1.07017; atol = 1e-2)
                @test isapprox(resultACP["solution"]["nw"]["2"]["branchdc_ne"]["2"]["pf"], -1.07057; atol = 1e-2)
                @test isapprox(resultACP["solution"]["nw"]["2"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultACP["solution"]["nw"]["1"]["convdc_ne"]["2"]["pconv"], -2.08104; atol = 1e-2)
            end
            @testset "SOCWR" begin
                data_dc = build_mn_data("../test/data/tnep/case4_original.m")
                resultSOCWR = run_mp_tnepopf(data_dc, SOCWRPowerModel, gurobi, multinetwork=true; setting = s)
                @test isapprox(resultSOCWR["objective"], 10.2; atol = 1e-1)
                @test isapprox(resultSOCWR["solution"]["nw"]["1"]["branchdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCWR["solution"]["nw"]["1"]["branchdc_ne"]["2"]["pf"], -1.2073; atol = 1e-2)
                @test isapprox(resultSOCWR["solution"]["nw"]["2"]["branchdc_ne"]["2"]["pf"], -0.9908; atol = 1e-2)
                @test isapprox(resultSOCWR["solution"]["nw"]["2"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCWR["solution"]["nw"]["1"]["convdc_ne"]["2"]["pconv"], -2.1151; atol = 1e-2)
            end
                @testset "SOCBF" begin
                data_dc = build_mn_data("../test/data/tnep/case4_original.m")
                resultSOCBF = run_mp_tnepopf_bf(data_dc, SOCBFConicPowerModel, mosek, multinetwork=true; setting = s)
                @test isapprox(resultSOCBF["objective"], 10.2; atol = 1e-1)
                @test isapprox(resultSOCBF["solution"]["nw"]["1"]["branchdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["nw"]["1"]["branchdc_ne"]["2"]["pf"], -1.3660; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["nw"]["2"]["branchdc_ne"]["2"]["pf"], -1.3800; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["nw"]["2"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCBF["solution"]["nw"]["1"]["convdc_ne"]["2"]["pconv"], -2.1338; atol = 1e-2)
            end
        end
    end
    @testset "acdc tnep" begin
        @testset "DCP" begin
            data_acdc = build_mn_data("../test/data/tnep/case4_acdc.m")
            resultDCP = run_mp_acdctnepopf(data_acdc, DCPPowerModel, cbc, multinetwork=true; setting = s)
            @test isapprox(resultDCP["objective"], 607.909; atol = 1e-1)
            @test isapprox(resultDCP["solution"]["nw"]["1"]["branchdc_ne"]["3"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultDCP["solution"]["nw"]["1"]["branchdc_ne"]["3"]["pf"], -1.009; atol = 1e-2)
            @test isapprox(resultDCP["solution"]["nw"]["2"]["branchdc_ne"]["3"]["pf"], -1.009; atol = 1e-2)
            @test isapprox(resultDCP["solution"]["nw"]["2"]["convdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
            @test isapprox(resultDCP["solution"]["nw"]["1"]["convdc_ne"]["1"]["pconv"], -1; atol = 1e-2)
            @test isapprox(resultDCP["solution"]["nw"]["1"]["ne_branch"]["1"]["built"], 1; atol = 1e-2)
        end
        #@testset "LPAC" begin
          # data_acdc = build_mn_data("../test/data/tnep/case4_acdc.m")
          #  resultLPAC = run_mp_acdctnepopf(data_acdc, LPACCPowerModel, juniper, multinetwork=true; setting = s)
          #  @test isapprox(resultLPAC["objective"], 614.15; atol = 1e-1)
          #  @test isapprox(resultLPAC["solution"]["nw"]["1"]["branchdc_ne"]["3"]["isbuilt"], 1; atol = 1e-2)
          #  @test isapprox(resultLPAC["solution"]["nw"]["1"]["branchdc_ne"]["3"]["pf"], -1.2466; atol = 1e-2)
          #  @test isapprox(resultLPAC["solution"]["nw"]["2"]["branchdc_ne"]["3"]["pf"], -1.2466; atol = 1e-2)
          #  @test isapprox(resultLPAC["solution"]["nw"]["2"]["convdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
          #  @test isapprox(resultLPAC["solution"]["nw"]["1"]["convdc_ne"]["2"]["pconv"], -1.2355; atol = 1e-2)
          #  @test isapprox(resultLPAC["solution"]["nw"]["1"]["ne_branch"]["1"]["built"], 1; atol = 1e-2)
        #end
        if local_test == true
            @testset "ACP" begin # Normally fine but fails on travis
                data_acdc = build_mn_data("../test/data/tnep/case4_acdc.m")
                resultACP = run_mp_acdctnepopf(data_acdc, ACPPowerModel, juniper, multinetwork=true; setting = s)
                @test isapprox(resultACP["objective"], 648.318; atol = 1e-1)
                @test isapprox(resultACP["solution"]["nw"]["1"]["branchdc_ne"]["3"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultACP["solution"]["nw"]["1"]["branchdc_ne"]["3"]["pf"], -0.2327; atol = 1e-2)
                @test isapprox(resultACP["solution"]["nw"]["2"]["branchdc_ne"]["3"]["pf"], -0.2327; atol = 1e-2)
                @test isapprox(resultACP["solution"]["nw"]["2"]["convdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultACP["solution"]["nw"]["1"]["convdc_ne"]["1"]["pconv"], -0.5083; atol = 1e-2)
                @test isapprox(resultACP["solution"]["nw"]["1"]["ne_branch"]["1"]["built"], 1; atol = 1e-2)
            end
            @testset "SOCWR" begin
                data_acdc = build_mn_data("../test/data/tnep/case4_acdc.m")
                resultSOCWR = run_mp_acdctnepopf(data_acdc, SOCWRPowerModel, gurobi, multinetwork=true; setting = s)
                @test isapprox(resultSOCWR["objective"], 641.86; atol = 1e-1)
                @test isapprox(resultSOCWR["solution"]["nw"]["2"]["branchdc_ne"]["2"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCWR["solution"]["nw"]["1"]["branchdc_ne"]["2"]["pf"], -0.8352; atol = 1e-2)
                @test isapprox(resultSOCWR["solution"]["nw"]["2"]["branchdc_ne"]["2"]["pf"], -0.8563; atol = 1e-2)
                @test isapprox(resultSOCWR["solution"]["nw"]["2"]["convdc_ne"]["1"]["isbuilt"], 1; atol = 1e-2)
                @test isapprox(resultSOCWR["solution"]["nw"]["1"]["convdc_ne"]["1"]["pconv"], -0.8227; atol = 1e-2)
                @test isapprox(resultSOCWR["solution"]["nw"]["1"]["ne_branch"]["1"]["built"], 1; atol = 1e-2)
            end
        end
    end
end


text:
"""
	loop constraints
"""

struct LoopConstraints <: Constraints
	n
	ineq
	idx
	idx_t1
	idx_t2
end

function loop_constraints(model, idx, idx_t1, idx_t2)
	n = length(idx)
	ineq = (1:0)
	return LoopConstraints(n, ineq, idx, idx_t1, idx_t2)
end

function constraints!(c, Z, con::LoopConstraints, model, idx, h, T)
	n = con.n
	a = view(Z, idx.x[con.idx_t1][con.idx])
	b = view(Z, idx.x[con.idx_t2][con.idx])

	c[1:n] = b - a

    return nothing
end

function constraints_jacobian!(∇c, Z, con::LoopConstraints, model, idx, h, T)
	n = con.n
	shift = 0

	r_idx = 1:n
	c_idx = idx.x[con.idx_t2][con.idx]
	len = length(r_idx) * length(c_idx)
	∇c[shift .+ (1:len)] = vec(Diagonal(ones(n)))
	shift += len

	c_idx = idx.x[con.idx_t1][con.idx]
	len = length(r_idx) * length(c_idx)
	∇c[shift .+ (1:len)] = vec(Diagonal(-1.0 * ones(n)))
	shift += len
	return nothing
end

function constraints_sparsity(con::LoopConstraints, model, idx, T;
		shift_row = 0, shift_col = 0)
	row = []
    col = []

	n = con.n

	r_idx = shift_row .+ (1:n)

	c_idx = shift_col .+ idx.x[con.idx_t2][con.idx]
	row_col!(row, col, r_idx, c_idx)

	c_idx = shift_col .+ idx.x[con.idx_t1][con.idx]
	row_col!(row, col, r_idx, c_idx)

    return collect(zip(row, col))
end


text:
module TestObjectiveMap

using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

struct ObjectiveDummyBridge <: MOI.Bridges.Objective.AbstractBridge
    id::Int
end

function _test_empty(map)
    @test isempty(map)
    @test length(map) == 0
    @test isempty(values(map))
    @test iterate(map) === nothing
    @test MOI.Bridges.Objective.function_type(map) === nothing
end

function test_Map()
    map = MOI.Bridges.Objective.Map()
    _test_empty(map)
    @test sprint(MOI.Bridges.print_num_bridges, map) ==
          "\nwith 0 objective bridges"
    x = MOI.VariableIndex(1)
    MOI.Bridges.Objective.add_key_for_bridge(map, ObjectiveDummyBridge(1), x)
    @test MOI.Bridges.Objective.root_bridge(map) == ObjectiveDummyBridge(1)
    @test sprint(MOI.Bridges.print_num_bridges, map) ==
          "\nwith 1 objective bridge"
    func = 1.0x
    MOI.Bridges.Objective.add_key_for_bridge(map, ObjectiveDummyBridge(2), func)
    @test MOI.Bridges.Objective.root_bridge(map) == ObjectiveDummyBridge(2)
    @test sprint(MOI.Bridges.print_num_bridges, map) ==
          "\nwith 2 objective bridges"
    empty!(map)
    _test_empty(map)
    @test sprint(MOI.Bridges.print_num_bridges, map) ==
          "\nwith 0 objective bridges"
    return
end

function test_EmptyMap()
    map = MOI.Bridges.Objective.EmptyMap()
    _test_empty(map)
    empty!(map)
    _test_empty(map)
    @test sprint(MOI.Bridges.print_num_bridges, map) == ""
    return
end

end  # module

TestObjectiveMap.runtests()


text:
using LinearAlgebra

export reduce_tensor, expand_tensor

"""
    function decompose_gate!(gate_data::Array{<:Number, 4},
                             threshold::AbstractFloat=1e-15)

Function to decompose a tensor into two smaller tensors
"""
function decompose_gate(gate_data::Array{<:Number, 4},
                        threshold::AbstractFloat=1e-15)
    left_positions = [1, 3]
    right_positions = [2, 4]
    dims = size(gate_data)
    left_dims = [dims[x] for x in left_positions]
    right_dims = [dims[x] for x in right_positions]

    A = permutedims(gate_data, vcat(left_positions, right_positions))
    A = reshape(A, Tuple([prod(left_dims), prod(right_dims)]))

    # Use SVD here but QR could also be used
    F = svd(A)

    # find number of singular values above the threshold
    chi = sum(F.S .> threshold)
    s = sqrt.(F.S[1:chi])

    # assume that singular values and basis of U and V matrices are sorted
    # in descending order of singular value
    B = reshape(F.U[:, 1:chi] * Diagonal(s), Tuple(vcat(left_dims, [chi,])))
    C = reshape(Diagonal(s) * F.Vt[1:chi, :], Tuple(vcat([chi,], right_dims)))

    return B, C
end

"""
    function find_hyper_edges(A::AbstractArray{Elt, N}) where {Elt, N}

Function to identify hyper edges of tensors. Returns an array of tuples of indices
of the original tensor which can be identified
"""
function find_hyper_edges(A::AbstractArray{Elt, N}) where {Elt, N}
    if N == 2
        if isdiagonal(A)
            return [[1, 2]]
        else
            return Array{Array{Int64, 1}, 1}()
        end
    elseif N > 2
        groups = Array{Array{Int64, 1}, 1}()
        for i in 1:N-1
            for j in i+1:N
                if isdiagonal(A, Pair(i, j))
                    # check if there is overlap iwth any existing groups
                    if any(x -> length(intersect(x, (i, j))) > 0, groups)
                        # groups with intersections
                        igroups = [x for x in groups if length(intersect(x, (i, j))) > 0]
                        # filter these from main groups list
                        filter!(x -> length(intersect(x, (i, j))) == 0, groups)
                        merged_group = sort(reduce(union, igroups, init=[i, j]))
                        push!(groups, merged_group)
                    else
                        push!(groups, [i, j])
                    end
                end
            end
        end
        return groups
    end
    Array{Array{Int64, 1}, 1}()
end

"""
    isdiagonal(A::AbstractArray{Elt, N}, Pair{Int64, Int64}) where {Elt, N}

Function to check if the given matrix is diagonal along given axes
"""
function isdiagonal(A::AbstractArray{Elt, N}, p::Pair{Int64, Int64}) where {Elt, N}
    po = collect(1:N)
    function exchange!(a, e1, e2)  a[e1], a[e2] = a[e2], a[e1] end
    exchange!(po, 1, p[1])
    exchange!(po, 2, p[2])
    A = permutedims(A, po)
    d = size(A)
    A = reshape(A, (d[1], d[2], prod(d[3:end])))
    all(x -> isdiagonal(A[:, :, x]), 1:prod(d[3:end]))
end

"""
    isdiagonal(A::AbstractArray{Elt, 2}) where Elt

Function to check if the given matrix is diagonal
"""
function isdiagonal(A::AbstractArray{Elt, 2}) where Elt
    if real(eltype(A)) <: Integer
        all([A[i, j] == 0 for i in 1:size(A)[1], j in 1:size(A)[2] if i != j])
    else
        all([abs(A[i, j]) < eps(real(eltype(A))) for i in 1:size(A)[1], j in 1:size(A)[2] if i != j])
    end
end

function _reduce_tensor_dims(A::AbstractArray{Elt, N}, hyper_index_groups::Array{Array{Int64, 1}, 1}) where {Elt, N}
    tensor_dims = size(A)
    index_map = Dict{Int64, Int64}(x => x for x in 1:N)
    for group in hyper_index_groups
        @assert length(group) >= 2 "Not all hyper edge groups have minimum dimension 2"
        ref = group[1]
        for j in group[2:end]
            index_map[j] = ref
        end
    end
    remaining_dims = unique(sort(collect(values(index_map))))
    remaining_dim_sizes = [tensor_dims[x] for x in remaining_dims]
    return remaining_dims, remaining_dim_sizes
end

"""
    reduce_tensor(A::AbstractArray{Elt, N}, hyper_index_groups::Array{Int64, 1})

Function to reduce the rank of the given tensor assuming the given hyper edge groups.
For example a diagonal matrix will have a single hyper edge group with both indices [1, 2].
This function will reduce this to a vector containing only the diagonal elements. can
be seen as a generalisation of the diag function.

```jldoctest
julia> reduce_tensor([[1, 0] [0, 2]], [[1, 2]])
2-element Vector{Int64}:
 1
 2
```
"""
function reduce_tensor(A::AbstractArray{Elt, N}, hyper_index_groups::Array{Array{Int64, 1}, 1}) where {Elt, N}
    @assert all([g == sort(g) for g in hyper_index_groups]) "Hyperindex groups must be sorted in ascending rank order"
    remaining_dims, remaining_dim_sizes = _reduce_tensor_dims(A, hyper_index_groups)
    index_map = Dict{Int64, Int64}(x => x for x in 1:N)
    for group in hyper_index_groups
        @assert length(group) >= 2 "Not all hyper edge groups have minimum dimension 2"
        ref = group[1]
        for j in group[2:end]
            index_map[j] = ref
        end
    end

    # create array to store reduced tensor with correct dimensions
    inverse_index_map = Dict{Int64, Array{Int64, 1}}(x => Int64[] for x in remaining_dims)
    for (k, v) in pairs(index_map)
        push!(inverse_index_map[v], k)
    end
    function map_reduced_indices_to_full(reduced_index::Tuple)
        full_index_array = zeros(Int64, N)
        for (dim, val) in enumerate(reduced_index)
            full_dim = remaining_dims[dim]
            full_index_array[inverse_index_map[full_dim]] .= val
        end
        full_index_array
    end
    Ar = zeros(Elt, remaining_dim_sizes...)
    for i in CartesianIndices(Tuple(remaining_dim_sizes))
        Ar[i] = A[map_reduced_indices_to_full(Tuple(i))...]
    end
    Ar
end

struct BlockTensor{T,N,M} <: AbstractArray{T, N}
    size::NTuple{N, Int}
    index_map::NTuple{N, Int}
    data::AbstractArray{T,M}
end

BlockTensor(size::NTuple{N, Int64}) where N = BlockTensor{ComplexF64, N}(size)

"""Overload functions from base to make BlockTensor usable"""
Base.copy(a::BlockTensor{T, N}) where {T, N} = BlockTensor{T, N}(a.size)
Base.length(a::BlockTensor) = prod(size(a))
Base.size(a::BlockTensor) = a.size
Base.getindex(t::BlockTensor, i::Int64) = t[Tuple(CartesianIndices(t.size)[i])...]
# Base.getindex(t::BlockTensor, i::CartesianIndex) = t[Tuple(i)...]
Base.IndexStyle(::Type{<:BlockTensor}) = IndexLinear()

"""
    Base.getindex(t::BlockTensor{T, N, M}, i...) where {T, N, M}

Overload the getindex function. Returns zero if the indiex on ranks identified
as hyperindex groups differ and the relevant tensor index otherwise.
"""
function Base.getindex(t::BlockTensor{T, N, M}, i...) where {T, N, M}
    new_index = Vector{Int}(undef, M)
    for x in 1:N
        xm = findfirst(y -> y == t.index_map[x], t.index_map)
        if i[x] != i[xm]
            return convert(T, 0)
        end
        new_index[t.index_map[x]] = i[x]
    end
    return t.data[new_index...]
end

Base.show(io::IO, ::MIME"text/plain", a::BlockTensor) =
    print(io, "BlockTensor with dims $(a.size) and index map $(a.index_map)")

"""
    expand_tensor(A::AbstractArray{Elt, N}, hyper_index_groups::Array{Int64, 1})

Function to expand the rank of the given tensor assuming the given hyper edge groups.
Like a generalisation of Diagonal. For example if passed a vector and given hyperindex
groups [1,2], it will return a matrix where non diagonal elements are zero. Returns a
a BlockTensor object which does not store "off-diagonal" elements.

```
julia> expand_tensor([1, 2], [[1, 2]])
BlockTensor with dims (2, 2) and index map (1, 1)
```
"""
function expand_tensor(A::AbstractArray{Elt, N}, hyper_index_groups::Array{Array{Int64, 1}, 1}) where {Elt, N}
    full_dim = N + sum(map(x -> length(x) - 1, hyper_index_groups))
    index_map = OrderedDict(x => x for x in 1:full_dim)
    for g in hyper_index_groups
        g = sort(g)
        for i in g[2:end]
            index_map[i] = g[1]
        end
    end
    index_order_map = Dict(x => i for (i, x) in enumerate(sort(unique(values(index_map)))))
    index_map = Tuple(map(x -> index_order_map[x], values(index_map)))
    new_size = Tuple(map(x -> size(A)[x], values(index_map)))
    BlockTensor{Elt, full_dim, N}(new_size, index_map, A)
end


text:
# ---
# title: 673. Number of Longest Increasing Subsequence
# id: problem673
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/number-of-longest-increasing-subsequence/description/>
# hidden: true
# ---
#
# Given an integer array `nums`, return _the number of longest increasing
# subsequences._
#
# **Notice** that the sequence has to be **strictly** increasing.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,3,5,4,7]
#     Output: 2
#     Explanation: The two longest increasing subsequences are [1, 3, 4, 7] and [1, 3, 5, 7].
#
#
# **Example 2:**
#
#
#
#     Input: nums = [2,2,2,2,2]
#     Output: 5
#     Explanation: The length of longest continuous increasing subsequence is 1, and there are 5 subsequences' length is 1, so output 5.
#
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 2000`
#   * `-106 <= nums[i] <= 106`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end