text:
using LinearAlgebra
using Random
using Test

n = 5
Xorig = randn(n, n)
A = UpperTriangular(rand(n, n))
B = randn(n, n)
B = B'*B
ws = LyapdWs(n)
for row = 1:n
    X = copy(Xorig)
    LinearRationalExpectations.solve_one_row!(X, A, B, n, row, ws)
    @test X[row, 1:row] ≈ -(A[row, row]*A[1:row, 1:row] - I(row))\B[row, 1:row]
    @test X[1:row, row] == X[row, 1:row]
end

n = 5
Xorig = randn(n, n)
A = zeros(n, n) .+ UpperTriangular(rand(n, n))
B = randn(n, n)
B = B'*B
ws = LyapdWs(n)
for row = 2:n
    A[row, row - 1] = 0.5
    X = copy(Xorig)
    Borig = copy(B)
    LinearRationalExpectations.solve_two_rows!(X, A, B, n, row, ws)
    @test B == Borig
    @test ws.AA2[1:2*row, 1:2*row] == kron(A[1:row, 1:row], A[(row - 1):row, (row - 1):row]) - I
    @test vec(X[(row-1):row, 1:row]) ≈ -(kron(A[1:row, 1:row], A[(row - 1):row, (row - 1):row]) - I(2*row))\vec(B[(row-1):row, 1:row])
    @test X[1:row, row-1] == X[row-1, 1:row]
    @test X[1:row, row] == X[row, 1:row]
    A[row, row - 1] = 0
end

for i = 1:100
    s = Int64(floor(1000*rand()))
    Random.seed!(s)
    local n = 5
    X1 = randn(n, n)
    X = X1'*X1
    local A = randn(n, n)
    F = eigen(A)
    if maximum(abs.(F.values)) > 1.0
        continue
    end
    local B = X - A*X*A'
    local ws = LyapdWs(n)
    Σ = similar(X)
    extended_lyapd!(Σ, A, B, ws)
    @test Σ ≈ X
end

# with unit roots
for i = 1:100
    s = Int64(floor(1000*rand()))
    Random.seed!(s)
    local n = 2
    X1 = randn(n, n)
    X = X1'*X1
    local A = randn(n, n)
    F = eigen(A)
    eigenvalues = diagm(rand(n))
    local A = real(F.vectors*eigenvalues*inv(F.vectors))
    m = 2
    AA = vcat(hcat([ 1.0 -0.5 0; 0 1.0 0], randn(m, n)),
              [1 -1 0 0 0],
              hcat(zeros(n, m+1), A))
    if abs(det(AA))  < 1e-12
        continue
    end
    local B = X - A*X*A'
    BB = vcat(hcat(I(m), zeros(m, n+1)),
              [0 0 0.2 0 0],
              hcat(zeros(n, m+1), B))

    local ws = LyapdWs(n + m + 1)
    Σ = Matrix{Float64}(undef, n+m+1, n+m+1)
    extended_lyapd!(Σ, AA, BB, ws)
    @test Σ[m + 1 .+ (1:n), m + 1 .+ (1:n)] ≈ X
end


text:
include("commit_comment.jl")
event_request = create_event()
event_json = JSON.parse(IOBuffer(HTTP.payload(event_request)))
event = GitHub.event_from_payload!("commit_comment", event_json)

@testset "WebhookEvent" begin
    @test event.repository.name == "BenchmarkTrackers.jl"
    @test event.sender.login == "jrevels"
end # testset

@testset "EventListener" begin
    @test !(GitHub.has_valid_secret(event_request, "wrong"))
    @test GitHub.has_valid_secret(event_request, "secret")
    @test !(GitHub.is_valid_event(event_request, ["wrong"]))
    @test GitHub.is_valid_event(event_request, ["commit_comment"])
    @test !(GitHub.from_valid_repo(event, ["JuliaWeb/GitHub.jl"]))
    @test GitHub.from_valid_repo(event, ["JuliaCI/BenchmarkTrackers.jl"])
    @test GitHub.handle_event_request(event_request, x -> true,
                                    secret = "secret",
                                    events = ["commit_comment"],
                                    repos = ["JuliaCI/BenchmarkTrackers.jl"])

    @test begin
        listener = EventListener(x -> true;
                                secret = "secret",
                                repos = [Repo("JuliaCI/BenchmarkTrackers.jl"), "JuliaWeb/GitHub.jl"],
                                events = ["commit_comment"],
                                forwards = ["http://bob.com", HTTP.URI("http://jim.org")])
        r = listener.handle_request(HTTP.Request())
        r.status == 400
    end
end

@testset "CommentListener" begin
    result = GitHub.handle_comment((e, m) -> m, event, GitHub.AnonymousAuth(), r"`RunBenchmarks\(.*?\)`", false)
    @test result.match == "`RunBenchmarks(\"binary\", \"unary\")`"
    @test begin
        listener = CommentListener((x, y) -> true, r"trigger";
                                secret = "secret",
                                repos = [Repo("JuliaCI/BenchmarkTrackers.jl"), "JuliaWeb/GitHub.jl"],
                                forwards = ["http://bob.com", HTTP.URI("http://jim.org")],
                                check_collab = false)
        r = listener.listener.handle_request(HTTP.Request())
        r.status == 400
    end
end

ext: jl


text:
using Documenter
using Dante

makedocs(
    sitename = "Dante",
    format = Documenter.HTML(),
    modules = [Dante],
    pages = [
        "Home" => "index.md",
        "User Guide" => Any[
            "Getting Started" =>  "getting_started.md",
            "Input Parameters" => "input.md"
            ],
        "Method" => Any[
            "MHD" => "mhd.md",
            "Finite Volume" => "fv.md",
            "Standard Test" => "test.md"
            ],
        "API Reference" => "api.md"
    ]
)

deploydocs(;
    repo="github.com/henry2004y/Dante.jl",
)


text:
#===============================================================================
# Practice using julia modules
#
# Author: Carl Cortright
# Date: 12/31/2016
#
===============================================================================#
# This will precompile the module the first time that it is run
__precompile__()
# __precompile__(false) will cause the module to never be recursively recompliled

module mymodule

export x, y

x() = "x"
p() = "y"
y = p

end

# Modules are mostly independent from files
module mod

include("methods.jl")

using mymodule

println(x())

end

# If we don't want to use base or the default eval, we can use the baremodule
baremodule default

using Base

eval(x) = Core.eval(Mod, x)
eval(m, x) = Core.eval(m,x)

end

# The init function is used when we have code that must be run at runtime
module mymod

function __init__()
  println("In mymod init function")
  # Usually something like a ccall would go here
end

end

max_issues_repo_name: carlcortright/Playground


text:
"""
    split_model(model::EoSModel)

Takes in a model for a multi-component system and returns a vector of model for each pure system.

## Example:
```julia-repl
julia> gerg2 = GERG2008(["propane","pentane"])
GERG008 model with 2 components:
"propane"
"pentane"

julia> split_model(gerg2)
2-element Vector{GERG2008}:
 GERG2008("propane")
 GERG2008("pentane")
```
"""
function split_model end

"""
    is_splittable(model)::Bool

Trait to determine if a `EoSModel` should be splitted by itself or can be simply filled into a vector.
This is useful in the case of models without any parameters, as those models are impossible by definition to split, because they don't have any underlying data.

The Default is `is_splittable(model) = true`.
"""
is_splittable(model) = true

function split_model(param::SingleParam{T},
    splitter =split_model(1:length(param.components))) where T
    function generator(I)
    return SingleParam(
    param.name,
    param.components[I],
    param.values[I],
    param.ismissingvalues[I],
    param.sourcecsvs,
    param.sources
    )
    end
    return [generator(i) for i in splitter]
end

function split_model(param::ClapeyronParam,groups::GroupParam)
    return split_model(param,groups.i_groups)
end
#this conversion is lossy, as interaction between two or more components are lost.

function split_model(param::PairParam{T},
    splitter = split_model(1:length(param.components))) where T
    function generator(I)
        value = param.values[I,I]
        diagvalue = view(value,diagind(value))
        ismissingvalues = param.ismissingvalues[I,I]
        components = param.components[I]
        return PairParam{T}(
                param.name,
                components,
                value,
                diagvalue,
                ismissingvalues,
                param.sourcecsvs,
                param.sources
                )
    end
    return [generator(I) for I in splitter]
end

function split_model(param::AbstractVector)
    return [[xi] for xi in param]
end

#this conversion is lossy, as interaction between two or more components are lost.
#also, this conversion stores the site values for other components. (those are not used)
function split_model(param::AssocParam{T},
    splitter = split_model(1:length(param.components))) where T
    function generator(I)
        _value  = param.values[I,I]
        _ismissingvalues = param.ismissingvalues[I,I]
        return AssocParam{T}(
                param.name,
                param.components[I],
                _value,
                _ismissingvalues,
                param.allcomponentsites[I],
                param.sourcecsvs,
                param.sources
                )
        end
    return [generator(I) for I in splitter]
end

#this param has a defined split form
function split_model(groups::GroupParam)
    len = length(groups.components)
    function generator(i)
        return GroupParam(
        [groups.components[i]],
        [groups.groups[i]],
        [groups.n_groups[i]],
        [collect(1:length(groups.n_groups[i]))],
        groups.flattenedgroups[groups.i_groups[i]],
        [groups.n_groups[i]],
        1:length(groups.n_groups[i]),
        groups.sourcecsvs)
    end
    [generator(i) for i in 1:len]
end

function split_model(param::SiteParam,
    splitter = split_model(param.components))
    function generator(I)
        return SiteParam(
            param.components[I],
            param.sites[I],
            param.n_sites[I],
             param.i_sites[I],
            param.flattenedsites,
            param.n_flattenedsites[I],
            1:length(param.flattenedsites),
        param.sourcecsvs)
    end
    return [generator(i) for i in splitter]
end

function split_model(Base.@nospecialize(params::EoSParam),splitter)
    T = typeof(params)
    split_paramsvals = [split_model(getfield(params,i),splitter) for i  in fieldnames(T)]
    return T.(split_paramsvals...)
end

export SingleParam, SiteParam, PairParam, AssocParam, GroupParam
#

split_model(model::EoSModel) = auto_split_model(model)

function auto_split_model(Base.@nospecialize(model::EoSModel))
    try
        allfields = Dict{Symbol,Any}()

        if has_groups(typeof(model))
            splitter = model.groups.i_groups
        else
            splitter = split_model(1:length(model.components))
        end
        len = length(splitter)

        if hasfield(typeof(model),:groups)
            allfields[:groups] = split_model(model.groups)
        end
        M = typeof(model)

        len_comps = length(model.components)
        allfields[:components] = split_model(model.components)

        if hasfield(typeof(model),:icomponents)
            allfields[:icomponents] = [1:1 for _ in 1:len_comps]
        end

        #process all model fields
        modelfields = filter(x->getproperty(model,x) isa EoSModel,fieldnames(M))
        for modelkey in modelfields
            modelx = getproperty(model,modelkey)
            if is_splittable(modelx)
                allfields[modelkey]= split_model(modelx)
            else
                allfields[modelkey] = fill(modelx,len_comps)
            end
        end

        modelfields = filter(x->getproperty(model,x) isa Vector{<:EoSModel},fieldnames(M))
        for modelkey in modelfields
            modelx = getproperty(model,modelkey)
            allfields[modelkey] = [[tup] for tup in modelx]
        end

        #process all empty (Missing,Nothing) fields
        emptyfields = filter(x->getproperty(model,x) isa Union{Nothing,Missing},fieldnames(M))

        for emptykey in emptyfields
             allfields[emptykey] = fill(model.emptykey,len_comps)
        end

        if hasfield(typeof(model),:params)
            allfields[:params] = split_model(model.params,splitter)
        end

        if hasfield(typeof(model),:references)
            allfields[:references] = fill(model.references,len_comps)
        end

        if hasfield(typeof(model),:absolutetolerance)
            allfields[:absolutetolerance] = fill(model.absolutetolerance,len_comps)
        end
        if hasfield(typeof(model),:sites)
            allfields[:sites] = split_model(model.sites,splitter)
        end

        return [M((allfields[k][i] for k in fieldnames(M))...) for i in 1:len]
    catch e
        @show model
        return simple_split_model(model)
    end
end

##fallback,around 50 times slower if there is any need to read csvs

function simple_split_model(Base.@nospecialize(model::EoSModel))
    MODEL = typeof(model)
    pure = Vector{MODEL}(undef,0)
    for comp ∈ model.components
        push!(pure,MODEL([comp]))
    end
    return pure
end

export split_model


text:
module LightJATS

using PDFExtract

include("tree.jl")
include("jats.jl")
include("sais/StringMatch.jl")
include("align.jl")

end


text:
function LazyLoadedDataContext()
    return LazyLoadedDataContext(())
end

function get_celestial_bodies(obj::LazyLoadedDataContext)
    return jcall(obj, "getCelestialBodies", LazyLoadedCelestialBodies, ())
end

function get_frames(obj::LazyLoadedDataContext)
    return jcall(obj, "getFrames", LazyLoadedFrames, ())
end

function get_geo_magnetic_fields(obj::LazyLoadedDataContext)
    return jcall(obj, "getGeoMagneticFields", LazyLoadedGeoMagneticFields, ())
end

function get_gravity_fields(obj::LazyLoadedDataContext)
    return jcall(obj, "getGravityFields", LazyLoadedGravityFields, ())
end

function get_time_scales(obj::LazyLoadedDataContext)
    return jcall(obj, "getTimeScales", LazyLoadedTimeScales, ())
end


text:
module batch
srcDir  = dirname(@__FILE__)*"\\";

include(srcDir*"Model\\R_sep_Model.jl")
include(srcDir*"Model\\C_sep_Model.jl")
include(srcDir*"Model\\Utils_Model.jl")
include(srcDir*"sep_Compute.jl")
using .R_sep_Model
using .C_sep_Model
using .Utils_Model
using .sep_Compute

using JuMP
using CSV
using DataFrames

export batch_model
default_sdir = dirname(dirname(@__FILE__))*"\\assets\\bounds\\"
""" Produces models for a set of constraints."""
function batch_model(t,ρdic,df;sdir = default_sdir,con_list= ["S1G","S2G","S3G"]) # ,
    for ex in keys(ρdic)
        ρ = ρdic[ex]
        d = filter(:Example => ==(ex),df).Size[1]
        for con ∈  con_list, isR ∈ ['R','C'] # ,
            (!filter(:Example => ==(ex),df).isReal[1] && (isR == 'R')) ? continue : 0
            println("Currently $isR-modeling example $ex with constraints $con")
            sep_mod = (isR == 'R') ? R_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t;con_list=con) : C_sep_Model.Modelξₜˢᵉᵖ(ρ,d,t;con_list=con,noBlock = false)

            Utils_Model.export_model(sep_mod, sdir*"$(ex)_$(con)_$(isR)_$t.dat-s")
        end
    end
end

""" Reads all .dat-s files and computes bounds storing result summary in a .csv file"""
function batch_Computeξₜˢᵉᵖ(boundsDir = default_sdir)
    datsFiles = [file for file in readdir(boundsDir,join = true) if contains(file,".dat-s")]
    file_loc = boundsDir*"Summary.csv"
    touch(file_loc)
    open(file_loc,"a") do io
        write(io, "|Example|Model|Constraint|Primal|Dual|obj_val|\n")
    end
    for file in datsFiles
        name = split(basename(file),".")[1]
        ex   = split(name,"_")[1]
        con  = split(name,"_")[2]
        num  = split(name,"_")[3]
        # try
        mod  = Utils_Model.read_model(file)
        omod = sep_Compute.Computeξₜˢᵉᵖ(mod)
        # catch
        #     open(file_loc,"a") do io; write(io, "|$ex|$con|err.|err.|err.| \n"); end
        #     continue
        # end
        pstat   = JuMP.primal_status(omod)
        dstat   = JuMP.dual_status(omod)
        ov_temp = round(JuMP.objective_value(omod),digits=3)
        ov = (string(pstat) == string(dstat) == "FEASIBLE_POINT") ?
             ov_temp :
             (contains(string(pstat)*string(dstat),"CER") ? "*" : NaN )

        open(file_loc,"a") do io
            write(io, "|$ex|$num|$con|$pstat|$dstat|$ov|\n")
        end
    end
end

"""Postprocesses the output of batch_Computeξₜˢᵉᵖ into a nices tabular structure"""
function unstack_constraints(df,bDir = default_sdir)
    temp_df  = CSV.read("C:\\Users\\andries\\all-my-codes\\sep-rank\\assets\\bounds\\"*"Summary.csv", DataFrame)
    temp2_df = select(innerjoin(df, temp_df, on = :Example),[:Example,:isReal,:isSeparable,:Model,:Size,:Bi_rank,:Constraint,:obj_val])
    temp4_df = unstack(temp2_df, :Model, :obj_val)

    temp4R_df = select(temp4_df,[:Example,:isReal,:isSeparable,:Size,:Bi_rank,:Constraint,:R])
    temp4Rus_df = unstack(temp4R_df, :Constraint, :R)
    rename!(temp4Rus_df,Dict(:S1G=>:RS1G, :S2G => :RS2G, :S3G => :RS3G));

    temp4C_df = select(temp4_df,[:Example,:isReal,:isSeparable,:Size,:Bi_rank,:Constraint,:C])
    temp4Cus_df = unstack(temp4C_df, :Constraint, :C)
    rename!(temp4Cus_df,Dict(:S1G=>:CS1G, :S2G => :CS2G, :S3G => :CS3G));

    temp5_df = innerjoin(temp4Cus_df,select(temp4Rus_df,[:Example ,:RS1G, :RS2G, :RS3G]), on = :Example)

    CSV.write(bDir*"SummaryUnstacked.csv", temp5_df, delim="&")
    final_df = DataFrames.select(temp5_df,[:Example,:Size,:Bi_rank,:CS1G, :CS2G, :CS3G, :RS1G, :RS2G, :RS3G, :isSeparable])
    rename!(final_df,Dict(:Example=>:ρ, :Size => :d₁d₂));
    CSV.write(bDir*"t=__.csv", final_df, delim="&")
    return final_df
end

end

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
init_bounds(bounds, pomdp, sol) = bounds
init_bounds(t::Tuple, pomdp, sol) = (init_bound(first(t), pomdp, sol), init_bound(last(t), pomdp, sol))
bounds(f::Function, pomdp::POMDP, b::ScenarioBelief) = f(pomdp, b)
bounds(t::Tuple, pomdp::POMDP, b::ScenarioBelief) = (lbound(t[1], pomdp, b), ubound(t[2], pomdp, b))

function bounds_sanity_check(pomdp::POMDP, sb::ScenarioBelief, L_0, U_0)
    if L_0 > U_0
        @warn("L_0 ($L_0) > U_0 ($U_0)   |ϕ| = $(length(sb.scenarios))")
    end
    if all(isterminal(pomdp, s) for s in particles(sb))
        if L_0 != 0.0 || U_0 != 0.0
            error(@sprintf("If all states are terminal, lower and upper bounds should be zero (L_0=%-10.2g, U_0=%-10.2g). (try IndependentBounds(l, u, check_terminal=true))", L_0, U_0))
        end
    end
end

# Upper and lower do not depend on each other (most cases)

struct IndependentBounds{L, U}
    lower::L
    upper::U
    check_terminal::Bool
end

IndependentBounds(l, u; check_terminal=false) = IndependentBounds(l, u, check_terminal)

function bounds(bounds::IndependentBounds, pomdp::POMDP, b::ScenarioBelief)
    if bounds.check_terminal && all(isterminal(pomdp, s) for s in particles(b))
        return (0.0, 0.0)
    end
    return (lbound(bounds.lower, pomdp, b), ubound(bounds.upper, pomdp, b))
end

function init_bounds(bounds::IndependentBounds, pomdp::POMDP, sol::DESPOTSolver)
    return IndependentBounds(init_bound(bounds.lower, pomdp, sol),
                             init_bound(bounds.upper, pomdp, sol),
                             bounds.check_terminal
                            )
end
init_bound(bound, pomdp, sol) = bound

ubound(n::Number, pomdp, b) = convert(Float64, n)
lbound(n::Number, pomdp, b) = convert(Float64, n)

ubound(f::Function, pomdp, b) = f(pomdp, b)
lbound(f::Function, pomdp, b) = f(pomdp, b)

# Value if Fully Observed Under a Policy
struct FullyObservableValueUB{P<:Union{Solver, Policy}}
    p::P
end

ubound(ub::FullyObservableValueUB, pomdp::POMDP, b::ScenarioBelief) = mean(value(ub.p, s) for s in particles(b)) # assumes that all are weighted equally

function init_bound(ub::FullyObservableValueUB{S}, pomdp::POMDP, sol::DESPOTSolver) where S <: Solver
    return FullyObservableValueUB(solve(ub.p, pomdp))
end

# Default Policy Lower Bound

struct DefaultPolicyLB{P<:Union{Solver, Policy}, D<:Union{Nothing,Int}}
    policy::P
    max_depth::D
end

function DefaultPolicyLB(policy_or_solver::T; max_depth=nothing) where T <: Union{Solver, Policy}
    return DefaultPolicyLB(policy_or_solver, max_depth)
end

function lbound(lb::DefaultPolicyLB, pomdp::POMDP, b::ScenarioBelief)
    rsum = branching_sim(pomdp, lb.policy, b, lb.max_depth-b.depth)
    return rsum/length(b.scenarios)
end

function init_bound(lb::DefaultPolicyLB{S}, pomdp::POMDP, sol::DESPOTSolver) where S <: Solver
    policy = solve(lb.policy, pomdp)
    return init_bound(DefaultPolicyLB(policy, lb.max_depth), pomdp, sol)
end

function init_bound(lb::DefaultPolicyLB{P}, pomdp::POMDP, sol::DESPOTSolver) where P <: Policy
    max_depth = something(lb.max_depth, sol.D)
    return DefaultPolicyLB(lb.policy, max_depth)
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct DefaultApi <: SwaggerApi
    client::Swagger.Client
end

"""

Checks whether a private IP address is available for use.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: ipAddress::String
Return: IPAddressAvailabilityResult
"""
function virtualNetworksCheckIPAddressAvailability(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, api_version::String, subscriptionId::String; ipAddress=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", IPAddressAvailabilityResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}/CheckIPAddressAvailability", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "ipAddress", ipAddress)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

"""

Lists usage stats.
Param: resourceGroupName::String (required)
Param: virtualNetworkName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: VirtualNetworkListUsageResult
"""
function virtualNetworksListUsage(_api::DefaultApi, resourceGroupName::String, virtualNetworkName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", VirtualNetworkListUsageResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Network/virtualNetworks/{virtualNetworkName}/usages", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "virtualNetworkName", virtualNetworkName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json", "text/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json", "text/json"] : [_mediaType])
    Swagger.exec(_ctx)
end

export virtualNetworksCheckIPAddressAvailability, virtualNetworksListUsage


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# McCormick.jl
# A McCormick operator library in Julia
# See https://github.com/PSORLab/McCormick.jl
#############################################################################
# src/forward_operators/trilinear.jl
# Defines trilinear product-based on the Meyer-Floudas envelope.
#############################################################################

#=
Computes the maximum of six values and returns the position of their occurance.
=#
function max6(a::Float64, b::Float64, c::Float64, d::Float64, e::Float64, f::Float64)
    if a >= b && a >= c && a >= d && a >= e && a >= f
        return a, 1
    elseif b >= a && b >= c && b >= d && b >= e && b >= f
        return b, 2
    elseif c >= a && c >= b && c >= d && c >= e && c >= f
        return c, 3
    elseif d >= a && d >= b && d >= c && d >= e && d >= f
        return d, 4
    elseif e >= a && e >= b && e >= c && e >= d && e >= f
        return e, 5
    end
    return f, 6
end

#=
Computes the maximum of six values and returns the position of their occurance.
=#
function min6(a::Float64, b::Float64, c::Float64, d::Float64, e::Float64, f::Float64)
    if a <= b && a <= c && a <= d && a <= e && a <= f
        return a, 1
    elseif b <= a && b <= c && b <= d && b <= e && b <= f
        return b, 2
    elseif c <= a && c <= b && c <= d && c <= e && c <= f
        return c, 3
    elseif d <= a && d <= b && d <= c && d <= e && d <= f
        return d, 4
    elseif e <= a && e <= b && e <= c && e <= d && e <= f
        return e, 5
    end
    return f, 6
end

#=
Picks the `i + 1` argument provided of six elements.
=#
function coeff6(i::Int64, a::Float64, b::Float64, c::Float64, d::Float64, e::Float64, f::Float64)
    (i === 1) && (return a)
    (i === 2) && (return b)
    (i === 3) && (return c)
    (i === 4) && (return d)
    (i === 5) && (return e)
    return f
end

macro is_neg(x)
    new_sym::Symbol = Symbol(String(x)*"U")
    esc(quote $new_sym <= 0.0 end)
end
macro is_pos(x)
    new_sym::Symbol = Symbol(String(x)*"L")
    esc(quote $new_sym >= 0.0 end)
end
macro is_mix(x)
    new_symU::Symbol = Symbol(String(x)*"U")
    new_symL::Symbol = Symbol(String(x)*"L")
    esc(quote $new_symU >= 0.0 >= $new_symL end)
end

macro unpack_trilinear_bnd()
    esc(quote
            xL = x.Intv.lo;     xU = x.Intv.hi
            yL = y.Intv.lo;     yU = y.Intv.hi
            zL = z.Intv.lo;     zU = z.Intv.hi
            xyLL = xL*yL;       xyLU = xL*yU;
            xyUL = xU*yL;       xyUU = xU*yU
            xzLL = xL*zL;       xzLU = xL*zU;
            xzUL = xU*zL;       xzUU = xU*zU
            yzLL = yL*zL;       yzLU = yL*zU;
            yzUL = yU*zL;       yzUU = yU*zU
            xyzLLL = xyLL*zL;   xyzLLU = xyLL*zU
            xyzLUL = xyLU*zL;   xyzULL = xyUL*zL
            xyzLUU = xyLU*zU;   xyzULU = xyUL*zU
            xyzUUL = xyUU*zL;   xyzUUU = xyUU*zU
        end)
end

macro unpack_trilinear_end()
    esc(quote
        cv, cvind = max6(cv1, cv2, cv3, cv4, cv5, cv6)
        cvsubax = coeff6(cvind, cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6)
        cvsubay = coeff6(cvind, cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6)
        cvsubaz = coeff6(cvind, cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6)
        cv_grad = cvsubax*ifelse(cvsubax > 0.0, x.cv_grad, x.cc_grad) +
                  cvsubay*ifelse(cvsubay > 0.0, y.cv_grad, y.cc_grad) +
                  cvsubaz*ifelse(cvsubaz > 0.0, z.cv_grad, z.cc_grad)

        cc, ccind = min6(cc1, cc2, cc3, cc4, cc5, cc6)
        ccsubax = coeff6(ccind, cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6)
        ccsubay = coeff6(ccind, cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6)
        ccsubaz = coeff6(ccind, cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6)
        cc_grad = ccsubax*ifelse(ccsubax > 0.0, x.cc_grad, x.cv_grad) +
                  ccsubay*ifelse(ccsubay > 0.0, y.cc_grad, y.cv_grad) +
                  ccsubaz*ifelse(ccsubaz > 0.0, z.cc_grad, z.cv_grad)

        MC{N,T}(cv, cc, q, cv_grad, cc_grad, x.cnst && y.cnst && z.cnst)
    end)
end

#=
trilinear_case_1

Case 3.1 + Case 4.1 of Meyer-Floudas 2004
=#
function trilinear_case_1(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()

    delX = xU - xL
    θcv1 = xyzUUL - xyzLUU - xyzULL + xyzULU
    θcv2 = xyzLLU - xyzULL - xyzLUU + xyzLUL

    # define cv and coefficients
    cv_b1 = -2.0*xyzLLL
    cv_b2 = -2.0*xyzUUU
    cv_b3 = -(xyzLLU + xyzULU)
    cv_b4 = -(xyzUUL + xyzLUL)
    cv_b5 = -θcv1*xL/delX - xyzUUL - xyzULU + xyzLUU
    cv_b6 = θcv2*xU/delX - xyzLLU - xyzLUL + xyzULL

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzLL, yzUU, yzLU, yzUL, θcv1/delX, -θcv2/delX
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzLL, xzUU, xzLU, xzUL, xzUL, xzLU
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyLL, xyUU, xyUL, xyLU, xyUL, xyLU

    cv1 = cv_ax1*x.cv + cv_ay1*y.cv + cv_az1*z.cv + cv_b1
    cv2 = cv_ax2*x.cv + cv_ay2*y.cv + cv_az2*z.cv + cv_b2
    cv3 = cv_ax3*x.cv + cv_ay3*y.cv + cv_az3*z.cv + cv_b3
    cv4 = cv_ax4*x.cv + cv_ay4*y.cv + cv_az4*z.cv + cv_b4
    cv5 = cv_ax5*ifelse(cv_ax5 > 0.0, x.cv, -x.cc) + cv_ay5*y.cv + cv_az5*z.cv + cv_b5
    cv6 = cv_ax6*ifelse(cv_ax6 > 0.0, x.cv, -x.cc) + cv_ay6*y.cv + cv_az6*z.cv + cv_b6

    # define cc and coefficients
    cc_b1 = -(xyzUUL + xyzULL)
    cc_b2 = -(xyzUUL + xyzLUL)
    cc_b3 = -(xyzULU + xyzULL)
    cc_b4 = -(xyzLUU + xyzLUL)
    cc_b5 = -(xyzULU + xyzLLU)
    cc_b6 = -(xyzLUU + xyzLLU)

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzLL, yzUL, yzLL, yzUU, yzLU, yzUU
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzUL, xzLL, xzUU, xzLL, xzUU, xzLU
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyUU, xyUU, xyUL, xyLU, xyLL, xyLL

    cc1 = cc_ax1*x.cc + cc_ay1*y.cc + cc_az1*z.cc + cc_b1
    cc2 = cc_ax2*x.cc + cc_ay2*y.cc + cc_az2*z.cc + cc_b2
    cc3 = cc_ax3*x.cc + cc_ay3*y.cc + cc_az3*z.cc + cc_b3
    cc4 = cc_ax4*x.cc + cc_ay4*y.cc + cc_az4*z.cc + cc_b4
    cc5 = cc_ax5*x.cc + cc_ay5*y.cc + cc_az5*z.cc + cc_b5
    cc6 = cc_ax6*x.cc + cc_ay6*y.cc + cc_az6*z.cc + cc_b6

    @unpack_trilinear_end()
end

#=
trilinear_case_2

Case 3.2 + Case 4.2 of Meyer-Floudas 2004
=#
function trilinear_case_2(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()

    delZ = zU - zL
    θcv = xyzLUU - xyzUUL - xyzLLU + xyzULU
    θcc = xyzULL - xyzUUU - xyzLLL + xyzLUL

    # define cv and coefficients
    cv_b1 = -2.0*xyzUUU
    cv_b2 = -(xyzLUL + xyzLUU)
    cv_b3 = -(xyzLUL + xyzLLL)
    cv_b4 = -(xyzULU + xyzULL)
    cv_b5 = -(xyzULL + xyzLLL)
    cv_b6 = -θcv*zL/delZ - xyzLUU - xyzULU + xyzUUL

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzUU, yzUL, yzUL, yzLU, yzLL, yzLU
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzUU, xzLU, xzLL, xzUL, xzUL, xzLU
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyUU, xyLU, xyLL, xyUL, xyLL, θcv/delZ

    # define cc and coefficients
    cc_b1 = -2.0*xyzUUL
    cc_b2 = -(xyzULU + xyzULL)
    cc_b3 = -(xyzLUU + xyzLLU)
    cc_b4 = -(xyzLUU + xyzLUL)
    cc_b5 = -(xyzULU + xyzLLU)
    cc_b6 = θcc*zU/delZ - xyzULL - xyzLUL + xyzUUU

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzUL, yzLL, yzUU, yzUU, yzLU, yzLL
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzUL, xzUU, xzLU, xzLL, xzUU, xzLL
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyUU, xyUL, xyLL, xyLU, xyLL, -θcc/delZ

    cv1 = cv_ax1*x.cv + cv_ay1*y.cv + cv_az1*z.cv + cv_b1
    cv2 = cv_ax2*x.cc + cv_ay2*y.cv + cv_az2*z.cc + cv_b2
    cv3 = cv_ax3*x.cc + cv_ay3*y.cc + cv_az3*z.cc + cv_b3
    cv4 = cv_ax4*x.cv + cv_ay4*y.cc + cv_az4*z.cv + cv_b4
    cv5 = cv_ax5*x.cc + cv_ay5*y.cc + cv_az5*z.cc + cv_b5
    cv6 = cv_ax6*x.cv + cv_ay6*y.cv + cv_az6*ifelse(cv_az6 > 0.0, z.cv, z.cc) + cv_b6

    cc1 = cc_ax1*x.cv + cc_ay1*y.cv + cc_az1*z.cc + cc_b1
    cc2 = cc_ax2*x.cv + cc_ay2*y.cc + cc_az2*z.cv + cc_b2
    cc3 = cc_ax3*x.cc + cc_ay3*y.cc + cc_az3*z.cv + cc_b3
    cc4 = cc_ax4*x.cc + cc_ay4*y.cv + cc_az4*z.cc + cc_b4
    cc5 = cc_ax5*x.cc + cc_ay5*y.cc + cc_az5*z.cv + cc_b5
    cc6 = cc_ax6*x.cv + cc_ay6*y.cv + cc_az6*ifelse(cc_az6 > 0.0, z.cc, z.cv) + cc_b6

    @unpack_trilinear_end()
end

function trilinear_case3_chk_cv(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()
    xyzUUL + xyzLLU <= xyzLUL + xyzULU
end

function trilinear_case3_chk_cc(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()
    xyzLLL + xyzUUU >= xyzULL + xyzLUU
end

function trilinear_case3_cv(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    delY = yU - yL
    delZ = zU - zL

    θcv1 = xyzLLL - xyzUUL - xyzLLU + xyzULU
    θcv2 = xyzULU - xyzUUL - xyzLLU + xyzLUU

    cv_b1 = -2.0*xyzUUU
    cv_b2 = -2.0*xyzULL
    cv_b3 = -(xyzLUL + xyzLUU)
    cv_b4 = -(xyzLUL + xyzLLL)
    cv_b5 = θcv1*yU/delY - xyzLLL - xyzULU + xyzUUL
    cv_b6 = -θcv2*zL/delZ - xyzULU - xyzLUU + xyzUUL

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzUU, yzLL, yzUL, yzUL, yzLU, yzLU
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzUU, xzUL, xzLU, xzLL, -θcv1/delY, xzLU
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyUU, xyUL, xyLU, xyLL, xyLL, θcv2/delZ

    cv1 = cv_ax1*x.cv + cv_ay1*y.cv + cv_az1*z.cv + cv_b1
    cv2 = cv_ax2*x.cv + cv_ay2*y.cc + cv_az2*z.cc + cv_b2
    cv3 = cv_ax3*x.cc + cv_ay3*y.cv + cv_az3*z.cv + cv_b3
    cv4 = cv_ax4*x.cc + cv_ay4*y.cc + cv_az4*z.cc + cv_b4
    cv5 = cv_ax5*x.cc + cv_ay5*ifelse(cv_ay5 > 0.0, y.cv, y.cc) + cv_az5*z.cc + cv_b5
    cv6 = cv_ax6*x.cc + cv_ay6*y.cv + cv_az6*ifelse(cv_az6 > 0.0, z.cv, z.cc) + cv_b6

    cv, cvind = max6(cv1, cv2, cv3, cv4, cv5, cv6)
    cvsubax = coeff6(cvind, cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6)
    cvsubay = coeff6(cvind, cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6)
    cvsubaz = coeff6(cvind, cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6)
    cv_grad = cvsubax*ifelse(cvsubax > 0.0, x.cv_grad, x.cc_grad) +
              cvsubay*ifelse(cvsubay > 0.0, y.cv_grad, y.cc_grad) +
              cvsubaz*ifelse(cvsubaz > 0.0, z.cv_grad, z.cc_grad)

    return cv, cv_grad
end

function trilinear_case3_cc1(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    delY = yU - yL
    delZ = zU - zL

    θcc1 = xyzULL - xyzUUU - xyzLLL + xyzLUL
    θcc2 = xyzULL - xyzUUU - xyzLLL + xyzLLU

    cc_b1 = -2.0*xyzULU
    cc_b2 = -2.0*xyzUUL
    cc_b3 = -(xyzLUU + xyzLLU)
    cc_b4 = -(xyzLUU + xyzLUL)
    cc_b5 = θcc1*zU/delZ - xyzULL - xyzLUL + xyzUUU
    cc_b6 = θcc2*yU/delY - xyzULL - xyzLLU + xyzUUU

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzLU, yzUL, yzUU, yzUU, yzLL, yzLL
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzUU, xzUL, xzLU, xzLL, xzLL, -θcc2/delY
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyUL, xyUU, xyLL, xyLU, -θcc1/delZ, xyLL

    cc1 = cc_ax1*x.cv + cc_ay1*y.cc + cc_az1*z.cv + cc_b1
    cc2 = cc_ax2*x.cv + cc_ay2*y.cv + cc_az2*z.cc + cc_b2
    cc3 = cc_ax3*x.cc + cc_ay3*y.cc + cc_az3*z.cv + cc_b3
    cc4 = cc_ax4*x.cc + cc_ay4*y.cv + cc_az4*z.cc + cc_b4
    cc5 = cc_ax5*x.cc + cc_ay5*y.cv + cc_az5*ifelse(cc_az5 > 0.0, z.cc, z.cv) + cc_b5
    cc6 = cc_ax6*x.cc + cc_ay6*ifelse(cc_ay6 > 0.0, y.cc, y.cv) + cc_az6*z.cv + cc_b6

    cc, ccind = min6(cc1, cc2, cc3, cc4, cc5, cc6)
    ccsubax = coeff6(ccind, cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6)
    ccsubay = coeff6(ccind, cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6)
    ccsubaz = coeff6(ccind, cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6)
    cc_grad = ccsubax*ifelse(ccsubax > 0.0, x.cc_grad, x.cv_grad) +
              ccsubay*ifelse(ccsubay > 0.0, y.cc_grad, y.cv_grad) +
              ccsubaz*ifelse(ccsubaz > 0.0, z.cc_grad, z.cv_grad)

    return cc, cc_grad
end

function trilinear_case3_cc2(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    delY = yU - yL
    delZ = zU - zL

    θcc1 = xyzUUU - xyzULL - xyzLUU + xyzLLU
    θcc2 = xyzUUU - xyzULL - xyzLUU + xyzLUL

    cc_b1 = -2.0*xyzULU
    cc_b2 = -2.0*xyzUUL
    cc_b3 = -(xyzLLL + xyzLUL)
    cc_b4 = -(xyzLLL + xyzLLU)
    cc_b5 = -θcc1*zL/delZ - xyzUUU - xyzLLU + xyzULL
    cc_b6 = -θcc2*yL/delY - xyzUUU - xyzLUL + xyzULL

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzLU, yzUL, yzLL, yzLL, yzUU, yzUU
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzUU, xzUL, xzLL, xzLU, xzLU, θcc1/delZ
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyUL, xyUU, xyLU, xyLL, θcc2/delY, xyLU

    cc1 = cc_ax1*x.cv + cc_ay1*y.cc + cc_az1*z.cv + cc_b1
    cc2 = cc_ax2*x.cv + cc_ay2*y.cv + cc_az2*z.cc + cc_b2
    cc3 = cc_ax3*x.cc + cc_ay3*y.cv + cc_az3*z.cc + cc_b3
    cc4 = cc_ax4*x.cc + cc_ay4*y.cc + cc_az4*z.cv + cc_b4
    cc5 = cc_ax5*x.cc + cc_ay5*y.cc + cc_az5*ifelse(cc_az5 > 0.0, z.cc, z.cv) + cc_b5
    cc6 = cc_ax6*x.cc + cc_ay6*ifelse(cc_ay6 > 0.0, y.cc, y.cv) + cc_az6*z.cc + cc_b6

    cc, ccind = min6(cc1, cc2, cc3, cc4, cc5, cc6)
    ccsubax = coeff6(ccind, cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6)
    ccsubay = coeff6(ccind, cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6)
    ccsubaz = coeff6(ccind, cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6)
    cc_grad = ccsubax*ifelse(ccsubax > 0.0, x.cc_grad, x.cv_grad) +
              ccsubay*ifelse(ccsubay > 0.0, y.cc_grad, y.cv_grad) +
              ccsubaz*ifelse(ccsubaz > 0.0, z.cc_grad, z.cv_grad)

    return cc, cc_grad
end

#=
trilinear_case_3

Case 3.3 + Case 4.3 of Meyer-Floudas 2004
=#
function trilinear_case_3(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N,T<:RelaxTag}
    xL = x.Intv.lo;     xU = x.Intv.hi
    yL = y.Intv.lo;     yU = y.Intv.hi
    zL = z.Intv.lo;     zU = z.Intv.hi

    if @is_pos(x) && @is_mix(y) && @is_mix(z)
        if trilinear_case3_chk_cv(x,y,z)
            cv, cv_grad = trilinear_case3_cv(x,y,z)
        else
            cv, cv_grad = trilinear_case3_cv(x,z,y)
        end
        if trilinear_case3_chk_cc(x,y,z)
            cc, cc_grad = trilinear_case3_cc1(x,y,z)
        elseif trilinear_case3_chk_cc(x,z,y)
            cc, cc_grad = trilinear_case3_cc1(x,z,y)
        else
            cc, cc_grad = trilinear_case3_cc2(x,y,z)
        end
        return MC{N,T}(cv, cc, q, cv_grad, cc_grad, x.cnst && y.cnst && z.cnst)
    elseif @is_mix(x) && @is_pos(y) && @is_mix(z)
        if trilinear_case3_chk_cv(y,x,z)
            cv, cv_grad = trilinear_case3_cv(y,x,z)
        else
            cv, cv_grad = trilinear_case3_cv(y,z,x)
        end
        if trilinear_case3_chk_cc(y,x,z)
            cc, cc_grad = trilinear_case3_cc1(y,x,z)
        elseif trilinear_case3_chk_cc(y,z,x)
            cc, cc_grad = trilinear_case3_cc1(y,z,x)
        else
            cc, cc_grad = trilinear_case3_cc2(y,x,z)
        end
        return MC{N,T}(cv, cc, q, cv_grad, cc_grad, x.cnst && y.cnst && z.cnst)
    end

    if trilinear_case3_chk_cv(z,x,y)
        cv, cv_grad = trilinear_case3_cv(z,x,y)
    else
        cv, cv_grad = trilinear_case3_cv(z,y,x)
    end
    if trilinear_case3_chk_cc(z,x,y)
        cc, cc_grad = trilinear_case3_cc1(z,x,y)
    elseif trilinear_case3_chk_cc(z,y,x)
        cc, cc_grad = trilinear_case3_cc1(z,y,x)
    else
        cc, cc_grad = trilinear_case3_cc2(z,x,y)
    end
    MC{N,T}(cv, cc, q, cv_grad, cc_grad, x.cnst && y.cnst && z.cnst)
end

function trilinear_case4_chk_cv1(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    (xyzUUL + xyzULU + xyzLUU <= xyzLLL + 2*xyzUUU) &&
    (xyzLLL + xyzUUL + xyzULU <= xyzLUU + 2*xyzULL) &&
    (xyzLLL + xyzUUL + xyzLUU <= xyzULU + 2*xyzLUL) &&
    (xyzLLL + xyzULU + xyzLUU <= xyzUUL + 2*xyzLLU)
end
function trilinear_case4_chk_cv2(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    xyzUUL + xyzULU + xyzLUU >= xyzLLL + 2*xyzUUU
end
function trilinear_case4_chk_cv3(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    xyzLLL + xyzUUL + xyzULU >= xyzLUU + 2.0*xyzULL
end

function trilinear_case4_cv1(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    delX = xU - xL
    delY = yU - yL
    delZ = zU - zL

    θcv1 = -0.5*(xyzLUU + xyzLLL - xyzUUL - xyzULU)/delX
    θcv2 = -0.5*(xyzULU + xyzLLL - xyzUUL - xyzLUU)/delY
    θcv3 = -0.5*(xyzUUL + xyzLLL - xyzULU - xyzLUU)/delX
    θcv4 = xyzLLL - θcv1*xL - θcv2*yL - θcv3*zL

    cv_b1 = -2.0*xyzLUL
    cv_b2 = -2.0*xyzLLU
    cv_b3 = -2.0*xyzUUU
    cv_b4 = -2.0*xyzULL
    cv_b5 = θcv4
    cv_b6 = Inf #TODO? eq 6 not in either paper

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzUL, yzLU, yzUU, yzLL, θcv1, Inf  # TODO? eq 6 not in paper
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzLL, xzLU, xzUU, xzUL, θcv2, Inf  # TODO? eq 6 not in paper
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyLU, xyLL, xyUU, xyUL, θcv3, Inf  # TODO? eq 6 not in paper

    cv1 = cv_ax1*x.cc + cv_ay1*y.cv + cv_az1*z.cc + cv_b1
    cv2 = cv_ax2*x.cc + cv_ay2*y.cc + cv_az2*z.cv + cv_b2
    cv3 = cv_ax3*x.cv + cv_ay3*y.cv + cv_az3*z.cv + cv_b3
    cv4 = cv_ax4*x.cv + cv_ay4*y.cc + cv_az4*z.cc + cv_b4
    cv5 = cv_ax5*ifelse(cv_ax5 > 0.0, x.cv, x.cc) +
          cv_ay5*ifelse(cv_ay5 > 0.0, y.cv, y.cc) +
          cv_az5*ifelse(cv_az5 > 0.0, z.cv, z.cc) + cv_b5
    cv6 = -Inf # TODO? eq 6 not in either paper

    cv, cvind = max6(cv1, cv2, cv3, cv4, cv5, cv6)
    cvsubax = coeff6(cvind, cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6)
    cvsubay = coeff6(cvind, cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6)
    cvsubaz = coeff6(cvind, cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6)
    cv_grad = cvsubax*ifelse(cvsubax > 0.0, x.cv_grad, x.cc_grad) +
              cvsubay*ifelse(cvsubay > 0.0, y.cv_grad, y.cc_grad) +
              cvsubaz*ifelse(cvsubaz > 0.0, z.cv_grad, z.cc_grad)

    return cv, cv_grad
end
function trilinear_case4_cv2(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    delX = xU - xL
    delY = yU - yL
    delZ = zU - zL

    θcv1 = xyzUUL - xyzLLL - xyzUUU + xyzULU
    θcv2 = xyzULU - xyzLLL - xyzUUU + xyzLUU
    θcv3 = xyzUUL - xyzLLL - xyzUUU + xyzLUU

    cv_b1 = -2.0*xyzLUL
    cv_b2 = -2.0*xyzLLU
    cv_b3 = -2.0*xyzULL
    cv_b4 = -θcv1*xL/delX - xyzUUL - xyzULU + xyzLLL
    cv_b5 = -θcv2*zL/delZ - xyzULU - xyzLUU + xyzLLL
    cv_b6 = -θcv3*yL/delY - xyzUUL - xyzLUU + xyzLLL

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzUL, yzLU, yzLL, θcv1/delX, yzUU, yzUU
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzLL, xzLU, xzUL, xzUU, xzUU, θcv3/delY
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyLU, xyLL, xyUL, xyUU, θcv2/delZ, xyUU

    cv1 = cv_ax1*x.cc + cv_ay1*y.cv + cv_az1*z.cc + cv_b1
    cv2 = cv_ax2*x.cc + cv_ay2*y.cc + cv_az2*z.cv + cv_b2
    cv3 = cv_ax3*x.cv + cv_ay3*y.cc + cv_az3*z.cc + cv_b3
    cv4 = cv_ax4*ifelse(cv_ax4 > 0.0, x.cv, x.cc) + cv_ay4*y.cv + cv_az4*z.cv + cv_b4
    cv5 = cv_ax5*x.cv + cv_ay5*ifelse(cv_ay5 > 0.0, y.cv, y.cc) + cv_az5*z.cv + cv_b5
    cv6 = cv_ax6*x.cv + cv_ay6*y.cv + cv_az6*ifelse(cv_az6 > 0.0, z.cv, z.cc) + cv_b6

    cv, cvind = max6(cv1, cv2, cv3, cv4, cv5, cv6)
    cvsubax = coeff6(cvind, cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6)
    cvsubay = coeff6(cvind, cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6)
    cvsubaz = coeff6(cvind, cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6)
    cv_grad = cvsubax*ifelse(cvsubax > 0.0, x.cv_grad, x.cc_grad) +
              cvsubay*ifelse(cvsubay > 0.0, y.cv_grad, y.cc_grad) +
              cvsubaz*ifelse(cvsubaz > 0.0, z.cv_grad, z.cc_grad)

    return cv, cv_grad
end
function trilinear_case4_cv3(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    delX = xU - xL
    delY = yU - yL
    delZ = zU - zL

    θcv1 = xyzLLL - xyzLUU - xyzULL + xyzULU
    θcv2 = xyzUUL - xyzLUU - xyzULL + xyzULU
    θcv3 = xyzLLL - xyzLUU - xyzULL + xyzUUL

    cv_b1 = -2.0*xyzLUL
    cv_b2 = -2.0*xyzLLU
    cv_b3 = -2.0*xyzUUU
    cv_b4 = θcv1*yU/delY - xyzLLL - xyzULU + xyzLUU
    cv_b5 = -θcv2*xL/delX - xyzUUL - xyzULU + xyzLUU
    cv_b6 = θcv3*zU/delZ - xyzLLL - xyzUUL + xyzLUU

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzUL, yzLU, yzUU, yzLL, θcv2/delX, yzLL
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzLL, xzLU, xzUU, -θcv1/delY, xzUL, xzUL
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyLU, xyLL, xyUU, xyUL, xyUL, -θcv3/delZ

    cv1 = cv_ax1*x.cc + cv_ay1*y.cv + cv_az1*z.cc + cv_b1
    cv2 = cv_ax2*x.cc + cv_ay2*y.cc + cv_az2*z.cv + cv_b2
    cv3 = cv_ax3*x.cv + cv_ay3*y.cv + cv_az3*z.cv + cv_b3
    cv4 = cv_ax4*x.cv + cv_ay4*ifelse(cv_ay4 > 0.0, y.cv, y.cc) + cv_az4*z.cc + cv_b4
    cv5 = cv_ax5*ifelse(cv_ax5 > 0.0, x.cv, x.cc) + cv_ay5*y.cc + cv_az5*z.cc + cv_b5
    cv6 = cv_ax6*x.cv + cv_ay6*y.cc + cv_az6*ifelse(cv_az6 > 0.0, z.cv, z.cc) + cv_b6

    cv, cvind = max6(cv1, cv2, cv3, cv4, cv5, cv6)
    cvsubax = coeff6(cvind, cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6)
    cvsubay = coeff6(cvind, cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6)
    cvsubaz = coeff6(cvind, cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6)
    cv_grad = cvsubax*ifelse(cvsubax > 0.0, x.cv_grad, x.cc_grad) +
              cvsubay*ifelse(cvsubay > 0.0, y.cv_grad, y.cc_grad) +
              cvsubaz*ifelse(cvsubaz > 0.0, z.cv_grad, z.cc_grad)

    return cv, cv_grad
end

function trilinear_case4_chk_cc1(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    (xyzULL + xyzLUL + xyzLLU) >= (xyzUUU + 2.0*xyzLLL) &&
    (xyzLUL + xyzLLU + xyzUUU) >= (xyzULL + 2.0*xyzLUU) &&
    (xyzULL + xyzLLU + xyzUUU) >= (xyzLUL + 2.0*xyzULU) &&
    (xyzULL + xyzLUL + xyzUUU) >= (xyzLLU + 2.0*xyzUUL)
end
function trilinear_case4_chk_cc2(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    (xyzULL + xyzLUL + xyzLLU) <= (xyzUUU + 2.0*xyzLLL)
end
function trilinear_case4_chk_cc3(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    xyzLUL + xyzLLU + xyzUUU <= xyzULL + 2.0*xyzLUU
end

function trilinear_case4_cc1(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    delX = xU - xL
    delY = yU - yL
    delZ = zU - zL

    θcc1 = 0.5*(xyzULL + xyzUUU - xyzLLU - xyzLUL)/delX
    θcc2 = 0.5*(xyzLUL + xyzUUU - xyzLLU - xyzULL)/delY
    θcc3 = 0.5*(xyzLLU + xyzUUU - xyzLUL - xyzULL)/delZ
    θcc4 = xyzUUU - θcc1*xU - θcc2*yU - θcc3*zU

    cc_b1 = -2.0*xyzLLL
    cc_b2 = -2.0*xyzUUL
    cc_b3 = -2.0*xyzULU
    cc_b4 = -2.0*xyzLUU
    cc_b5 = θcc4
    cc_b6 = Inf #TODO? eq 6 not in paper

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzLL, yzUL, yzLU, yzUU, θcc1, Inf
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzLL, xzUL, xzUU, xzLU, θcc2, Inf
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyLL, xyUU, xyUL, xyLU, θcc3, Inf

    cc1 = cc_ax1*x.cc + cc_ay1*y.cc + cc_az1*z.cc + cc_b1
    cc2 = cc_ax2*x.cv + cc_ay2*y.cv + cc_az2*z.cc + cc_b2
    cc3 = cc_ax3*x.cv + cc_ay3*y.cc + cc_az3*z.cv + cc_b3
    cc4 = cc_ax4*x.cc + cc_ay4*y.cv + cc_az4*z.cv + cc_b4
    cc5 = cc_ax5*ifelse(cc_ax5 > 0.0, x.cc, x.cv) +
          cc_ay5*ifelse(cc_ay5 > 0.0, y.cc, y.cv) +
          cc_az5*ifelse(cc_az5 > 0.0, z.cc, z.cv) + cc_b5
    cc6 = Inf # TODO? eq 6 not in paper

    cc, ccind = min6(cc1, cc2, cc3, cc4, cc5, cc6)
    ccsubax = coeff6(ccind, cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6)
    ccsubay = coeff6(ccind, cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6)
    ccsubaz = coeff6(ccind, cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6)
    cc_grad = ccsubax*ifelse(ccsubax > 0.0, x.cc_grad, x.cv_grad) +
              ccsubay*ifelse(ccsubay > 0.0, y.cc_grad, y.cv_grad) +
              ccsubaz*ifelse(ccsubaz > 0.0, z.cc_grad, z.cv_grad)

    return cc, cc_grad
end
function trilinear_case4_cc2(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    delX = xU - xL
    delY = yU - yL
    delZ = zU - zL

    θcc1 = xyzULL - xyzUUU - xyzLLL + xyzLUL
    θcc2 = xyzLLU - xyzUUU - xyzLLL + xyzLUL
    θcc3 = xyzLLU - xyzUUU - xyzLLL + xyzULL

    cc_b1 = -2.0*xyzUUL
    cc_b2 = -2.0*xyzULU
    cc_b3 = -2.0*xyzLUU
    cc_b4 = θcc1*zU/delZ - xyzULL - xyzLUL + xyzUUU
    cc_b5 = θcc2*xU/delX - xyzLLU - xyzLUL + xyzUUU
    cc_b6 = θcc3*yU/delY - xyzLLU - xyzULL + xyzUUU

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzUL, yzLU, yzUU, yzLL, -θcc2/delX, yzLL
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzUL, xzUU, xzLU, xzLL, xzLL, -θcc3/delY
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyUU, xyUL, xyLU, -θcc1/delZ, xyLL, xyLL

    cc1 = cc_ax1*x.cv + cc_ay1*y.cv + cc_az1*z.cc + cc_b1
    cc2 = cc_ax2*x.cv + cc_ay2*y.cc + cc_az2*z.cv + cc_b2
    cc3 = cc_ax3*x.cc + cc_ay3*y.cv + cc_az3*z.cv + cc_b3
    cc4 = cc_ax4*x.cc + cc_ay4*y.cc + cc_az4*ifelse(cc_az4 > 0.0, z.cc, z.cv) + cc_b4
    cc5 = cc_ax5*ifelse(cc_ax5 > 0.0, x.cc, x.cv) + cc_ay5*y.cc + cc_az5*z.cc + cc_b5
    cc6 = cc_ax6*x.cc + cc_ay6*ifelse(cc_ay6 > 0.0, y.cc, y.cv) + cc_az6*z.cc + cc_b6

    cc, ccind = min6(cc1, cc2, cc3, cc4, cc5, cc6)
    ccsubax = coeff6(ccind, cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6)
    ccsubay = coeff6(ccind, cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6)
    ccsubaz = coeff6(ccind, cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6)
    cc_grad = ccsubax*ifelse(ccsubax > 0.0, x.cc_grad, x.cv_grad) +
              ccsubay*ifelse(ccsubay > 0.0, y.cc_grad, y.cv_grad) +
              ccsubaz*ifelse(ccsubaz > 0.0, z.cc_grad, z.cv_grad)

    return cc, cc_grad
end
function trilinear_case4_cc3(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    delX = xU - xL
    delY = yU - yL
    delZ = zU - zL

    θcc1 = xyzLLU - xyzULL - xyzLUU + xyzLUL
    θcc2 = xyzUUU - xyzULL - xyzLUU + xyzLUL
    θcc3 = xyzLLU - xyzULL - xyzLUU + xyzUUU

    cc_b1 = -2.0*xyzLLL
    cc_b2 = -2.0*xyzUUL
    cc_b3 = -2.0*xyzULU
    cc_b4 = θcc1*xU/delX - xyzLLU - xyzLUL + xyzULL
    cc_b5 = -θcc2*yL/delY - xyzUUU - xyzLUL + xyzULL
    cc_b6 = -θcc3*zL/delZ - xyzLLU - xyzUUU + xyzULL

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzLL, yzUL, yzLU, -θcc1/delX, yzUU, yzUU
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzLL, xzUL, xzUU, xzLU, θcc2/delY, xzLU
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyLL, xyUU, xyUL, xyLU, xyLU, θcc3/delZ

    cc1 = cc_ax1*x.cc + cc_ay1*y.cc + cc_az1*z.cc + cc_b1
    cc2 = cc_ax2*x.cv + cc_ay2*y.cv + cc_az2*z.cc + cc_b2
    cc3 = cc_ax3*x.cv + cc_ay3*y.cc + cc_az3*z.cv + cc_b3
    cc4 = cc_ax4*ifelse(cc_ax4 > 0.0, x.cc, x.cv) + cc_ay4*y.cv + cc_az4*z.cv + cc_b4
    cc5 = cc_ax5*x.cc + cc_ay5*ifelse(cc_ay5 > 0.0, y.cc, y.cv) + cc_az5*z.cv + cc_b5
    cc6 = cc_ax6*x.cc + cc_ay6*y.cv + cc_az6*ifelse(cc_az6 > 0.0, z.cc, z.cv) + cc_b6

    cc, ccind = min6(cc1, cc2, cc3, cc4, cc5, cc6)
    ccsubax = coeff6(ccind, cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6)
    ccsubay = coeff6(ccind, cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6)
    ccsubaz = coeff6(ccind, cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6)
    cc_grad = ccsubax*ifelse(ccsubax > 0.0, x.cc_grad, x.cv_grad) +
              ccsubay*ifelse(ccsubay > 0.0, y.cc_grad, y.cv_grad) +
              ccsubaz*ifelse(ccsubaz > 0.0, z.cc_grad, z.cv_grad)

    return cc, cc_grad
end

#=
trilinear_case_4

Case 3.4 + Case 4.4 of Meyer-Floudas 2004
=#
function trilinear_case_4(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N,T<:RelaxTag}
    if trilinear_case4_chk_cv1(x,y,z)
        cv, cv_grad = trilinear_case4_cv1(x,y,z)
    elseif trilinear_case4_chk_cv1(x,z,y)
        cv, cv_grad = trilinear_case4_cv1(x,z,y)
    elseif trilinear_case4_chk_cv1(y,x,z)
        cv, cv_grad = trilinear_case4_cv1(y,x,z)
    elseif trilinear_case4_chk_cv1(y,z,x)
        cv, cv_grad = trilinear_case4_cv1(y,z,x)
    elseif trilinear_case4_chk_cv1(z,x,y)
        cv, cv_grad = trilinear_case4_cv1(z,x,y)
    elseif trilinear_case4_chk_cv1(z,y,x)
        cv, cv_grad = trilinear_case4_cv1(z,y,x)
    elseif trilinear_case4_chk_cv2(x,y,z)
        cv, cv_grad = trilinear_case4_cv2(x,y,z)
    elseif trilinear_case4_chk_cv2(x,z,y)
        cv, cv_grad = trilinear_case4_cv2(x,z,y)
    elseif trilinear_case4_chk_cv2(y,x,z)
        cv, cv_grad = trilinear_case4_cv2(y,x,z)
    elseif trilinear_case4_chk_cv2(y,z,x)
        cv, cv_grad = trilinear_case4_cv2(y,z,x)
    elseif trilinear_case4_chk_cv2(z,x,y)
        cv, cv_grad = trilinear_case4_cv2(z,x,y)
    elseif trilinear_case4_chk_cv2(z,y,x)
        cv, cv_grad = trilinear_case4_cv2(z,y,x)
    elseif trilinear_case4_chk_cv3(x,y,z)
        cv, cv_grad = trilinear_case4_cv3(x,y,z)
    elseif trilinear_case4_chk_cv3(x,z,y)
        cv, cv_grad = trilinear_case4_cv3(x,z,y)
    elseif trilinear_case4_chk_cv3(y,x,z)
        cv, cv_grad = trilinear_case4_cv3(y,x,z)
    elseif trilinear_case4_chk_cv3(y,z,x)
        cv, cv_grad = trilinear_case4_cv3(y,z,x)
    elseif trilinear_case4_chk_cv3(z,x,y)
        cv, cv_grad = trilinear_case4_cv3(z,x,y)
    else
        cv, cv_grad = trilinear_case4_cv3(z,y,x)
    end

    if trilinear_case4_chk_cc1(x,y,z)
        cc, cc_grad = trilinear_case4_cc1(x,y,z)
    elseif trilinear_case4_chk_cc1(x,z,y)
        cc, cc_grad = trilinear_case4_cc1(x,z,y)
    elseif trilinear_case4_chk_cc1(y,x,z)
        cc, cc_grad = trilinear_case4_cc1(y,x,z)
    elseif trilinear_case4_chk_cc1(y,z,x)
        cc, cc_grad = trilinear_case4_cc1(y,z,x)
    elseif trilinear_case4_chk_cc1(z,x,y)
        cc, cc_grad = trilinear_case4_cc1(z,x,y)
    elseif trilinear_case4_chk_cc1(z,y,x)
        cc, cc_grad = trilinear_case4_cc1(z,y,x)
    elseif trilinear_case4_chk_cc2(x,y,z)
        cc, cc_grad = trilinear_case4_cc2(x,y,z)
    elseif trilinear_case4_chk_cc2(x,z,y)
        cc, cc_grad = trilinear_case4_cc2(x,z,y)
    elseif trilinear_case4_chk_cc2(y,x,z)
        cc, cc_grad = trilinear_case4_cc2(y,x,z)
    elseif trilinear_case4_chk_cc2(y,z,x)
        cc, cc_grad = trilinear_case4_cc2(y,z,x)
    elseif trilinear_case4_chk_cc2(z,x,y)
        cc, cc_grad = trilinear_case4_cc2(z,x,y)
    elseif trilinear_case4_chk_cc2(z,y,x)
        cc, cc_grad = trilinear_case4_cc2(z,y,x)
    elseif trilinear_case4_chk_cc3(x,y,z)
        cc, cc_grad = trilinear_case4_cc3(x,y,z)
    elseif trilinear_case4_chk_cc3(x,z,y)
        cc, cc_grad = trilinear_case4_cc3(x,z,y)
    elseif trilinear_case4_chk_cc3(y,x,z)
        cc, cc_grad = trilinear_case4_cc3(y,x,z)
    elseif trilinear_case4_chk_cc3(y,z,x)
        cc, cc_grad = trilinear_case4_cc3(y,z,x)
    elseif trilinear_case4_chk_cc3(z,x,y)
        cc, cc_grad = trilinear_case4_cc3(z,x,y)
    else
        cc, cc_grad = trilinear_case4_cc3(z,y,x)
    end

    cvc, ccc, cvc_grad, ccc_grad = cut(q.lo, q.hi, cv, cc, cv_grad, cc_grad)

    MC{N,T}(cvc, ccc, q, cvc_grad, ccc_grad, x.cnst && y.cnst && z.cnst)
end

macro trilinear_case5_cv()
    esc(quote
        # define cv and coefficients
        cv_b1 = -(xyzLUL + xyzLLL)
        cv_b2 = -(xyzLUL + xyzLUU)
        cv_b3 = -(xyzULU + xyzULL)
        cv_b4 = -(xyzULU + xyzUUU)
        cv_b5 = -(xyzULL + xyzLLL)
        cv_b6 = -(xyzUUU + xyzLUU)

        cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzUL, yzUL, yzLU, yzLU, yzLL, yzUU
        cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzLL, xzLU, xzUL, xzUU, xzUL, xzLU
        cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyLL, xyLU, xyUL, xyUU, xyLL, xyUU

        cv1 = cv_ax1*x.cc + cv_ay1*y.cc + cv_az1*z.cv + cv_b1
        cv2 = cv_ax2*x.cc + cv_ay2*y.cc + cv_az2*z.cv + cv_b2
        cv3 = cv_ax3*x.cc + cv_ay3*y.cc + cv_az3*z.cv + cv_b3
        cv4 = cv_ax4*x.cc + cv_ay4*y.cc + cv_az4*z.cv + cv_b4
        cv5 = cv_ax5*x.cc + cv_ay5*y.cc + cv_az5*z.cv + cv_b5
        cv6 = cv_ax6*x.cc + cv_ay6*y.cc + cv_az6*z.cv + cv_b6
    end)
end

#=
trilinear_case_5a

Case 3.5 + Case 4.5 of Meyer-Floudas 2004
=#
function trilinear_case_5a(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()
    @trilinear_case5_cv()

    delZ = zU - zL
    θcc1 = xyzULL - xyzUUU - xyzLLL + xyzLUL
    θcc2 = xyzULU - xyzLLL - xyzUUU + xyzLUU

    cc_b1 = -2.0*xyzLLU
    cc_b2 = -2.0*xyzUUL
    cc_b3 = -(xyzULU + xyzULL)
    cc_b4 = -(xyzLUU + xyzLUL)
    cc_b5 = θcc1*zU/delZ - xyzULL - xyzLUL + xyzUUU
    cc_b6 = -θcc2*zL/delZ - xyzULU - xyzLUU + xyzLLL

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzLU, yzUL, yzLL, yzUU, yzLL, yzUU
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzLU, xzUL, xzUU, xzLL, xzLL, xzUU
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyLL, xyUU, xyUL, xyLU, -θcc1/delZ, θcc2/delZ

    cc1 = cc_ax1*x.cv + cc_ay1*y.cv + cc_az1*z.cc + cc_b1
    cc2 = cc_ax2*x.cv + cc_ay2*y.cv + cc_az2*z.cc + cc_b2
    cc3 = cc_ax3*x.cv + cc_ay3*y.cv + cc_az3*z.cc + cc_b3
    cc4 = cc_ax4*x.cv + cc_ay4*y.cv + cc_az4*z.cc + cc_b4
    cc5 = cc_ax5*x.cv + cc_ay5*y.cv + cc_az5*ifelse(cc_az5 > 0.0, z.cc, z.cv) + cc_b5
    cc6 = cc_ax6*x.cv + cc_ay6*y.cv + cc_az6*ifelse(cc_az6 > 0.0, z.cc, z.cv) + cc_b6

    @unpack_trilinear_end()
end

function trilinear_case_5b(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()
    @trilinear_case5_cv()

    delY = yU - yL
    θcc1 = xyzLLL - xyzLUU - xyzULL + xyzULU
    θcc2 = xyzLUL - xyzULL - xyzLUU + xyzUUU

    cc_b1 = -2.0*xyzLLU
    cc_b2 = -2.0*xyzUUL
    cc_b3 = -(xyzLLL + xyzLUL)
    cc_b4 = -(xyzULU + xyzUUU)
    cc_b5 = θcc1*yU/delY - xyzLLL - xyzULU + xyzLUU
    cc_b6 = -θcc2*yL/delY - xyzLUL - xyzUUU + xyzULL

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzLU, yzUL, yzLL, yzUU, yzLL, yzUU
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzLU, xzUL, xzLL, xzUU, -θcc1/delY, θcc2/delY
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyLL, xyUU, xyLU, xyUL, xyUL, xyLU

    cc1 = cc_ax1*x.cv + cc_ay1*y.cv + cc_az1*z.cc + cc_b1
    cc2 = cc_ax2*x.cv + cc_ay2*y.cv + cc_az2*z.cv + cc_b2
    cc3 = cc_ax3*x.cv + cc_ay3*y.cv + cc_az3*z.cv + cc_b3
    cc4 = cc_ax4*x.cv + cc_ay4*y.cv + cc_az4*z.cc + cc_b4
    cc5 = cc_ax5*x.cv + cc_ay5*y.cv + cc_az5*ifelse(cc_az5 > 0.0, z.cc, z.cv) + cc_b5
    cc6 = cc_ax6*x.cv + cc_ay6*y.cv + cc_az6*ifelse(cc_az6 > 0.0, z.cc, z.cv) + cc_b6

    @unpack_trilinear_end()
end

#=
trilinear_case_6

Note: Case 3.6/Case 4.6 of thesis yielded incorrect result. So reformulated...
=#
function trilinear_case_6(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()

    delY = yU - yL
    θcv = xyzULU - xyzUUL - xyzLLU + xyzLLL
    θcc = xyzLUL - xyzULL - xyzLUU + xyzUUU

    # define cv and coefficients
    cv_b1 = -2.0*xyzULL
    cv_b2 = -(xyzLUL + xyzLUU)
    cv_b3 = -(xyzLUL + xyzLLL)
    cv_b4 = -(xyzULU + xyzUUU)
    cv_b5 = -(xyzLUU + xyzUUU)
    cv_b6 = θcv*yU/delY - xyzULU - xyzLLL + xyzUUL

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzLL, yzUL, yzUL, yzLU, yzUU, yzLU
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzUL, xzLU, xzLL, xzUU, xzLU, -θcv/delY
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyUL, xyLU, xyLL, xyUU, xyUU, xyLL

    cv1 = cv_ax1*x.cv + cv_ay1*y.cc + cv_az1*z.cc + cv_b1
    cv2 = cv_ax2*x.cc + cv_ay2*y.cc + cv_az2*z.cv + cv_b2
    cv3 = cv_ax3*x.cc + cv_ay3*y.cc + cv_az3*z.cc + cv_b3
    cv4 = cv_ax4*x.cv + cv_ay4*y.cc + cv_az4*z.cv + cv_b4
    cv5 = cv_ax5*x.cc + cv_ay5*y.cc + cv_az5*z.cv + cv_b5
    cv6 = cv_ax6*x.cv + cv_ay6*ifelse(cv_ay6 > 0, y.cv, y.cc) + cv_az6*z.cc + cv_b6

    # define cc and coefficients
    cc_b1 = -2.0*xyzUUL
    cc_b2 = -(xyzUUU + xyzULU)
    cc_b3 = -(xyzLUL + xyzLLL)
    cc_b4 = -(xyzLLU + xyzLLL)
    cc_b5 = -(xyzLLU + xyzULU)
    cc_b6 = -θcc*yL/delY - xyzLUL - xyzUUU + xyzULL

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzUL, yzUU, yzLL, yzLL, yzLU, yzUU
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzUL, xzUU, xzLL, xzLU, xzLU, θcc/delY
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyUU, xyUL, xyLU, xyLL, xyUL, xyLU

    cc1 = cc_ax1*x.cv + cc_ay1*y.cv + cc_az1*z.cc + cc_b1
    cc2 = cc_ax2*x.cv + cc_ay2*y.cv + cc_az2*z.cv + cc_b2
    cc3 = cc_ax3*x.cc + cc_ay3*y.cv + cc_az3*z.cc + cc_b3
    cc4 = cc_ax4*x.cc + cc_ay4*y.cv + cc_az4*z.cv + cc_b4
    cc5 = cc_ax5*x.cc + cc_ay5*y.cv + cc_az5*z.cv + cc_b5
    cc6 = cc_ax6*x.cv + cc_ay6*ifelse(cc_ay6 > 0, y.cc, y.cv) + cc_az6*z.cc + cc_b6

    @unpack_trilinear_end()
end

function trilinear_case7_chk_cv(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()
    xyzLLL + xyzUUU <= xyzUUL + xyzLLU
end
function trilinear_case7_chk_cc(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()
    xyzULL + xyzLUU >= xyzLUL + xyzULU
end

function trilinear_case7_cv1(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    delX = xU - xL
    delY = yU - yL

    θcv1 = xyzUUL - xyzLLL - xyzUUU + xyzULU
    θcv2 = xyzUUL - xyzLLL - xyzUUU + xyzLUU

    cv_b1 = -2.0*xyzLUL
    cv_b2 = -2.0*xyzULL
    cv_b3 = -(xyzLUU + xyzLLU)
    cv_b4 = -(xyzULU + xyzLLU)
    cv_b5 = -θcv1*xL/delX - xyzUUL - xyzULU + xyzLLL
    cv_b6 = -θcv2*yL/delY - xyzUUL - xyzLUU + xyzLLL

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzUL, yzLL, yzUU, yzLU, θcv1/delX, yzUU
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzLL, xzUL, xzLU, xzUU, xzUU, θcv2/delY
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyLU, xyUL, xyLL, xyLL, xyUU, xyUU

    cv1 = cv_ax1*x.cc + cv_ay1*y.cv + cv_az1*z.cc + cv_b1
    cv2 = cv_ax2*x.cv + cv_ay2*y.cc + cv_az2*z.cc + cv_b2
    cv3 = cv_ax3*x.cc + cv_ay3*y.cv + cv_az3*z.cv + cv_b3
    cv4 = cv_ax4*x.cv + cv_ay4*y.cc + cv_az4*z.cv + cv_b4
    cv5 = cv_ax5*ifelse(cv_ax5 >= 0.0, x.cv, x.cc) + cv_ay5*y.cc + cv_az5*z.cv + cv_b5
    cv6 = cv_ax6*x.cc + cv_ay6*ifelse(cv_ay6 >= 0.0, y.cv, y.cc) + cv_az6*z.cv + cv_b6

    cv, cvind = max6(cv1, cv2, cv3, cv4, cv5, cv6)
    cvsubax = coeff6(cvind, cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6)
    cvsubay = coeff6(cvind, cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6)
    cvsubaz = coeff6(cvind, cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6)
    cv_grad = cvsubax*ifelse(cvsubax > 0.0, x.cv_grad, x.cc_grad) +
              cvsubay*ifelse(cvsubay > 0.0, y.cv_grad, y.cc_grad) +
              cvsubaz*ifelse(cvsubaz > 0.0, z.cv_grad, z.cc_grad)

    return cv, cv_grad
end
function trilinear_case7_cv2(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    delX = xU - xL
    delY = yU - yL

    θcv1 = xyzLLL - xyzUUL - xyzLLU + xyzLUU
    θcv2 = xyzLLL - xyzUUL - xyzLLU + xyzULU

    cv_b1 = -2.0*xyzLUL
    cv_b2 = -2.0*xyzULL
    cv_b3 = -(xyzULU + xyzUUU)
    cv_b4 = -(xyzLUU + xyzUUU)
    cv_b5 = θcv1*xU/delX - xyzLLL - xyzLUU + xyzUUL
    cv_b6 = θcv2*yU/delY - xyzLLL - xyzULU + xyzUUL

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzUL, yzLL, yzLU, yzUU, -θcv1/delX, yzLU
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzLL, xzUL, xzUU, xzLU, xzLU, -θcv2/delY
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyLU, xyUL, xyUU, xyUU, xyLL, xyLL

    cv1 = cv_ax1*x.cc + cv_ay1*y.cv + cv_az1*z.cc + cv_b1
    cv2 = cv_ax2*x.cv + cv_ay2*y.cc + cv_az2*z.cc + cv_b2
    cv3 = cv_ax3*x.cv + cv_ay3*y.cc + cv_az3*z.cv + cv_b3
    cv4 = cv_ax4*x.cc + cv_ay4*y.cv + cv_az4*z.cv + cv_b4
    cv5 = cv_ax5*ifelse(cv_ax6 > 0.0, x.cv, x.cc) + cv_ay5*y.cv + cv_az5*z.cv + cv_b5
    cv6 = cv_ax6*x.cv + cv_ay6*ifelse(cv_ay5 > 0.0, y.cv, y.cc) + cv_az6*z.cv + cv_b6

    cv, cvind = max6(cv1, cv2, cv3, cv4, cv5, cv6)
    cvsubax = coeff6(cvind, cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6)
    cvsubay = coeff6(cvind, cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6)
    cvsubaz = coeff6(cvind, cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6)
    cv_grad = cvsubax*ifelse(cvsubax > 0.0, x.cv_grad, x.cc_grad) +
              cvsubay*ifelse(cvsubay > 0.0, y.cv_grad, y.cc_grad) +
              cvsubaz*ifelse(cvsubaz > 0.0, z.cv_grad, z.cc_grad)

    return cv, cv_grad
end
function trilinear_case7_cc(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N,T<:RelaxTag}

    @unpack_trilinear_bnd()

    delX = xU - xL
    delY = yU - yL

    θcc1 = xyzLUL - xyzULL - xyzLUU + xyzLLU
    θcc2 = xyzUUU - xyzULL - xyzLUU + xyzLUL

    # define cc and coefficients
    cc_b1 = -2.0*xyzLLL
    cc_b2 = -2.0*xyzUUL
    cc_b3 = -(xyzLLU + xyzULU)
    cc_b4 = -(xyzUUU + xyzULU)
    cc_b5 = θcc1*xU/delX - xyzLUL - xyzLLU + xyzULL
    cc_b6 = -θcc2*yL/delY - xyzUUU - xyzLUL + xyzULL

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzLL, yzUL, yzLU, yzUU, -θcc1/delX, yzUU
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzLL, xzUL, xzLU, xzUU, xzLU, θcc2/delY
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyLL, xyUU, xyUL, xyUL, xyLU, xyLU

    cc1 = cc_ax1*x.cc + cc_ay1*y.cc + cc_az1*z.cc + cc_b1
    cc2 = cc_ax2*x.cv + cc_ay2*y.cv + cc_az2*z.cc + cc_b2
    cc3 = cc_ax3*x.cc + cc_ay3*y.cc + cc_az3*z.cv + cc_b3
    cc4 = cc_ax4*x.cv + cc_ay4*y.cv + cc_az4*z.cv + cc_b4
    cc5 = cc_ax5*ifelse(cc_ax5 > 0.0, x.cc, x.cv) + cc_ay5*y.cc + cc_az5*z.cv + cc_b5
    cc6 = cc_ax6*x.cv + cc_ay6*ifelse(cc_ay6 > 0.0, y.cc, y.cv) + cc_az6*z.cv + cc_b6

    cc, ccind = min6(cc1, cc2, cc3, cc4, cc5, cc6)
    ccsubax = coeff6(ccind, cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6)
    ccsubay = coeff6(ccind, cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6)
    ccsubaz = coeff6(ccind, cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6)
    cc_grad = ccsubax*ifelse(ccsubax > 0.0, x.cc_grad, x.cv_grad) +
              ccsubay*ifelse(ccsubay > 0.0, y.cc_grad, y.cv_grad) +
              ccsubaz*ifelse(ccsubaz > 0.0, z.cc_grad, z.cv_grad)

    return cc, cc_grad
end

function trilinear_case_7(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N,T<:RelaxTag}
    if trilinear_case7_chk_cv(x, y, z)
        cv, cv_grad = trilinear_case7_cv1(x, y, z)
    elseif trilinear_case7_chk_cv(y, x, z)
        cv, cv_grad = trilinear_case7_cv1(y, x, z)
    else
        cv, cv_grad = trilinear_case7_cv2(x, y, z)
    end

    if trilinear_case7_chk_cc(x, y, z)
        cc, cc_grad = trilinear_case7_cc(x, y, z)
    else
        cc, cc_grad = trilinear_case7_cc(y, x, z)
    end

    MC{N,T}(cv, cc, q, cv_grad, cc_grad, x.cnst && y.cnst && z.cnst)
end

macro trilinear_case8_cc()
    esc(quote
        # define cc and coefficients
        cc_b1 = -(xyzLLU + xyzLLL)
        cc_b2 = -(xyzLUL + xyzLLL)
        cc_b3 = -(xyzLLU + xyzULU)
        cc_b4 = -(xyzUUU + xyzULU)
        cc_b5 = -(xyzUUU + xyzUUL)
        cc_b6 = -(xyzLUL + xyzUUL)

        cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzLL, yzLL, yzLU, yzUU, yzUU, yzUL
        cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzLU, xzLL, xzLU, xzUU, xzUL, xzUL
        cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyLL, xyLU, xyUL, xyUL, xyUU, xyLU

        cc1 = cc_ax1*x.cc + cc_ay1*y.cv + cc_az1*z.cv + cc_b1
        cc2 = cc_ax2*x.cc + cc_ay2*y.cv + cc_az2*z.cv + cc_b2
        cc3 = cc_ax3*x.cc + cc_ay3*y.cv + cc_az3*z.cv + cc_b3
        cc4 = cc_ax4*x.cc + cc_ay4*y.cv + cc_az4*z.cv + cc_b4
        cc5 = cc_ax5*x.cc + cc_ay5*y.cv + cc_az5*z.cv + cc_b5
        cc6 = cc_ax6*x.cc + cc_ay6*y.cv + cc_az6*z.cv + cc_b6
    end)
end

#=
trilinear_case_8

Case 3.8 + Case 4.8 of Meyer-Floudas 2004
=#
function trilinear_case_8a(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()

    delX = xU - xL
    θcv1 = xyzLLU - xyzUUU - xyzLLL + xyzLUL
    θcv2 = xyzULU - xyzLLL - xyzUUU + xyzUUL

    cv_b1 = -2.0*xyzLUU
    cv_b2 = -2.0*xyzULL
    cv_b3 = -(xyzULU + xyzLLU)
    cv_b4 = -(xyzUUL + xyzLUL)
    cv_b5 = θcv1*xU/delX - xyzLLU - xyzLUL + xyzUUU
    cv_b6 = -θcv2*xL/delX - xyzULU - xyzUUL + xyzLLL

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzUU, yzLL, yzLU, yzUL, -θcv1/delX, θcv2/delX
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzLU, xzUL, xzUU, xzLL, xzLL, xzUU
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyLU, xyUL, xyLL, xyUU, xyLL, xyUU

    cv1 = cv_ax1*x.cv + cv_ay1*y.cc + cv_az1*z.cc + cv_b1
    cv2 = cv_ax2*x.cv + cv_ay2*y.cc + cv_az2*z.cc + cv_b2
    cv3 = cv_ax3*x.cv + cv_ay3*y.cc + cv_az3*z.cc + cv_b3
    cv4 = cv_ax4*x.cv + cv_ay4*y.cc + cv_az4*z.cc + cv_b4
    cv5 = cv_ax5*ifelse(cv_ax5 > 0.0, x.cv, x.cc) + cv_ay5*y.cc + cv_az5*z.cc + cv_b5
    cv6 = cv_ax6*ifelse(cv_ax6 > 0.0, x.cv, x.cc) + cv_ay6*y.cc + cv_az6*z.cc + cv_b6

    @trilinear_case8_cc()
    @unpack_trilinear_end()
end

function trilinear_case_8b(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()

    delY = yU - yL
    θcv1 = xyzUUU - xyzLLU - xyzUUL + xyzLUL
    θcv2 = xyzLLL - xyzUUL - xyzLLU + xyzULU

    cv_b1 = -2.0*xyzLUU
    cv_b2 = -2.0*xyzULL
    cv_b3 = -(xyzULU + xyzUUU)
    cv_b4 = -(xyzLLL + xyzLUL)
    cv_b5 = -θcv1*yL/delY - xyzUUU - xyzLUL + xyzLLU
    cv_b6 = θcv2*yU/delY - xyzLLL - xyzULU + xyzUUL

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzUU, yzLL, yzLU, yzUL, yzUL, yzLU
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzLU, xzUL, xzUU, xzLL, θcv1/delY, -θcv2/delY
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyLU, xyUL, xyUU, xyLL, xyUU, xyLL

    cv1 = cv_ax1*x.cv + cv_ay1*y.cc + cv_az1*z.cc + cv_b1
    cv2 = cv_ax2*x.cv + cv_ay2*y.cc + cv_az2*z.cc + cv_b2
    cv3 = cv_ax3*x.cv + cv_ay3*y.cc + cv_az3*z.cc + cv_b3
    cv4 = cv_ax4*x.cv + cv_ay4*y.cc + cv_az4*z.cc + cv_b4
    cv5 = cv_ax5*x.cv + cv_ay5*ifelse(cv_ay5 > 0.0, y.cv, y.cc) + cv_az5*z.cc + cv_b5
    cv6 = cv_ax6*x.cv + cv_ay6*ifelse(cv_az6 > 0.0, y.cv, y.cc) + cv_az6*z.cc + cv_b6

    @trilinear_case8_cc()
    @unpack_trilinear_end()
end

#=
trilinear_case_9

Case 3.9 + Case 4.9 of Meyer-Floudas 2004
=#
function trilinear_case_9(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()

    delX = xU - xL
    θcv = xyzUUL - xyzLLL - xyzUUU + xyzULU
    θcc = xyzLUL - xyzULL - xyzLUU + xyzLLU

    # define cv and coefficients
    cv_b1 = -2.0*xyzULL
    cv_b2 = -(xyzLUU + xyzLUL)
    cv_b3 = -(xyzLUU + xyzLLU)
    cv_b4 = -(xyzULU + xyzLLU)
    cv_b5 = -(xyzLUL + xyzUUL)
    cv_b6 = -θcv*xL/delX - xyzUUL - xyzULU + xyzLLL

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzLL, yzUU, yzUU, yzLU, yzUL, θcv/delX
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzUL, xzLL, xzLU, xzUU, xzLL, xzUU
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyUL, xyLU, xyLL, xyLL, xyUU, xyUU

    # define cc and coefficients
    cc_b1 = -2.0*xyzLLL
    cc_b2 = -(xyzLLU + xyzULU)
    cc_b3 = -(xyzUUU + xyzULU)
    cc_b4 = -(xyzUUU + xyzUUL)
    cc_b5 = -(xyzLUL + xyzUUL)
    cc_b6 = θcc*xU/delX - xyzLUL - xyzLLU + xyzULL

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzLL, yzLU, yzUU, yzUU, yzUL, -θcc/delX
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzLL, xzLU, xzUU, xzUL, xzUL, xzLU
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyLL, xyUL, xyUL, xyUU, xyLU, xyLU

    cv1 = cv_ax1*x.cv + cv_ay1*y.cc + cv_az1*z.cc + cv_b1
    cv2 = cv_ax2*x.cv + cv_ay2*y.cv + cv_az2*z.cv + cv_b2
    cv3 = cv_ax3*x.cv + cv_ay3*y.cv + cv_az3*z.cv + cv_b3
    cv4 = cv_ax4*x.cv + cv_ay4*y.cc + cv_az4*z.cv + cv_b4
    cv5 = cv_ax5*x.cv + cv_ay5*y.cv + cv_az5*z.cc + cv_b5
    cv6 = cv_ax6*ifelse(cv_ax6 > 0.0, x.cv, x.cc) + cv_ay6*y.cc + cv_az6*z.cc + cv_b6

    cc1 = cc_ax1*x.cc + cc_ay1*y.cc + cc_az1*z.cc + cc_b1
    cc2 = cc_ax2*x.cc + cc_ay2*y.cc + cc_az2*z.cv + cc_b2
    cc3 = cc_ax3*x.cc + cc_ay3*y.cv + cc_az3*z.cv + cc_b3
    cc4 = cc_ax4*x.cc + cc_ay4*y.cv + cc_az4*z.cv + cc_b4
    cc5 = cc_ax5*x.cc + cc_ay5*y.cv + cc_az5*z.cc + cc_b5
    cc6 = cc_ax6*ifelse(cc_ax6 > 0.0, x.cc, x.cv) + cc_ay6*y.cc + cc_az6*z.cc + cc_b6

    @unpack_trilinear_end()
end

#=
trilinear_case_10

Case 3.10 + Case 4.10 of Meyer-Floudas 2004
=#
function trilinear_case_10(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N,T<:RelaxTag}
    @unpack_trilinear_bnd()

    delX = xU - xL

    θcc1 = xyzLUL - xyzULL - xyzLUU + xyzLLU
    θcc2 = xyzUUL - xyzLUU - xyzULL + xyzULU

    # define cv and coefficients
    cv_b1 = -(xyzUUL + xyzULL)
    cv_b2 = -(xyzLUU + xyzLUL)
    cv_b3 = -(xyzUUL + xyzLUL)
    cv_b4 = -(xyzLUU + xyzLLU)
    cv_b5 = -(xyzULU + xyzLLU)
    cv_b6 = -(xyzULU + xyzULL)

    cv_ax1, cv_ax2, cv_ax3, cv_ax4, cv_ax5, cv_ax6 = yzLL, yzUU, yzUL, yzUU, yzLU, yzLL
    cv_ay1, cv_ay2, cv_ay3, cv_ay4, cv_ay5, cv_ay6 = xzUL, xzLL, xzLL, xzLU, xzUU, xzUU
    cv_az1, cv_az2, cv_az3, cv_az4, cv_az5, cv_az6 = xyUU, xyLU, xyUU, xyLL, xyLL, xyUL

    # define cc and coefficients
    cc_b1 = -2.0*xyzLLL
    cc_b2 = -2.0*xyzUUU
    cc_b3 = -(xyzLLU + xyzULU)
    cc_b4 = -(xyzLUL + xyzUUL)
    cc_b5 = θcc1*xU/delX - xyzLUL - xyzLLU + xyzULL
    cc_b6 = -θcc2*xL/delX - xyzUUL - xyzULU + xyzLUU

    cc_ax1, cc_ax2, cc_ax3, cc_ax4, cc_ax5, cc_ax6 = yzLL, yzUU, yzLU, yzUL, -θcc1/delX, θcc2/delX
    cc_ay1, cc_ay2, cc_ay3, cc_ay4, cc_ay5, cc_ay6 = xzLL, xzUU, xzLU, xzUL, xzLU, xzUL
    cc_az1, cc_az2, cc_az3, cc_az4, cc_az5, cc_az6 = xyLL, xyUU, xyUL, xyLU, xyLU, xyUL

    cv1 = cv_ax1*x.cv + cv_ay1*y.cv + cv_az1*z.cv + cv_b1
    cv2 = cv_ax2*x.cv + cv_ay2*y.cv + cv_az2*z.cv + cv_b2
    cv3 = cv_ax3*x.cv + cv_ay3*y.cv + cv_az3*z.cv + cv_b3
    cv4 = cv_ax4*x.cv + cv_ay4*y.cv + cv_az4*z.cv + cv_b4
    cv5 = cv_ax5*x.cv + cv_ay5*y.cv + cv_az5*z.cv + cv_b5
    cv6 = cv_ax6*x.cv + cv_ay6*y.cv + cv_az6*z.cv + cv_b6

    cc1 = cc_ax1*x.cc + cc_ay1*y.cc + cc_az1*z.cc + cc_b1
    cc2 = cc_ax2*x.cc + cc_ay2*y.cc + cc_az2*z.cc + cc_b2
    cc3 = cc_ax3*x.cc + cc_ay3*y.cc + cc_az3*z.cc + cc_b3
    cc4 = cc_ax4*x.cc + cc_ay4*y.cc + cc_az4*z.cc + cc_b4
    cc5 = cc_ax5*x.cc + cc_ay5*y.cc + cc_az5*z.cc + cc_b5
    cc6 = cc_ay6*x.cc + cc_ay6*y.cc + cc_az6*z.cc + cc_b6

    @unpack_trilinear_end()
end

x_mul_y2(x, y) = x*y^2

function trilinear_case1_map_chk(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N, T<:Union{NS,MV}}
    @unpack_trilinear_bnd()
    (xyzULL + xyzLUU <= xyzLUL + xyzULU) && (xyzULL + xyzLUU <= xyzUUL + xyzLLU)
end

function trilinear_case5_map_chk1(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N, T<:Union{NS,MV}}
    @unpack_trilinear_bnd()
    (xyzLLL + xyzUUU >= xyzULL + xyzLUU) &&
    (xyzLLL + xyzUUU >= xyzLUL + xyzULU)
end

function trilinear_case5_map_chk2(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N, T<:Union{NS,MV}}
    @unpack_trilinear_bnd()
    (xyzULL + xyzLUU >= xyzLUL + xyzULU)
end

function trilinear_case8_map_chk(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N, T<:Union{NS,MV}}
    @unpack_trilinear_bnd()
    (xyzLLL + xyzUUU <= xyzLLL + xyzUUU) && (xyzLLL + xyzUUU <= xyzLUL + xyzULU)
end

function trilinear_case10_map_chk(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}) where {N, T<:Union{NS,MV}}
    @unpack_trilinear_bnd()
    (xyzULL + xyzLUU >= xyzLUL + xyzULU) && (xyzULL + xyzLUU >= xyzUUL + xyzLLU)
end

function trilinear_kernel(x::MC{N,T}, y::MC{N,T}, z::MC{N,T}, q::Interval{Float64}) where {N, T<:Union{NS,MV}}
    #=
	if x == y
		if x == z
			return x^3
		else
			return x_mul_y2(z, x)
		end
	elseif y == z
		return x_mul_y2(x, y)
	end
    =#

    @unpack_trilinear_bnd()
	if @is_pos(x) && @is_pos(y) && @is_pos(z)
        if trilinear_case1_map_chk(x, y, z)
            return trilinear_case_1(x, y, z, q)
        elseif trilinear_case1_map_chk(x, z, y)
            return trilinear_case_1(x, z, y, q)
        elseif trilinear_case1_map_chk(y, x, z)
            return trilinear_case_1(y, x, z, q)
        elseif trilinear_case1_map_chk(y, z, x)
            return trilinear_case_1(y, z, x, q)
        elseif trilinear_case1_map_chk(z, x, y)
            return trilinear_case_1(z, x, y, q)
        else
            return trilinear_case_1(z, y, x, q)
        end
    end

    if @is_pos(x) && @is_pos(y) && @is_mix(z)
        return trilinear_case_2(x, y, z, q)
    elseif @is_pos(x) && @is_mix(y) && @is_pos(z)
        return trilinear_case_2(x, z, y, q)
    elseif @is_mix(x) && @is_pos(y) && @is_pos(z)
        return trilinear_case_2(y, z, x, q)
    elseif @is_mix(x) && @is_mix(y) && @is_mix(z)
        return trilinear_case_4(x, y, z, q)
    end

    if (@is_pos(x) && @is_pos(y) && @is_neg(z))
        if trilinear_case5_map_chk1(x, y, z)
            return trilinear_case_5a(x, y, z, q)
        elseif trilinear_case5_map_chk1(y, x, z)
            return trilinear_case_5a(y, x, z, q)
        elseif trilinear_case5_map_chk2(x, y, z)
            return trilinear_case_5b(x, y, z, q)
        elseif trilinear_case5_map_chk2(y, x, z)
            return trilinear_case_5b(y, x, z, q)
        end
    elseif (@is_pos(x) && @is_neg(y) && @is_pos(z))
        if trilinear_case5_map_chk1(x, z, y)
            return trilinear_case_5a(x, z, y, q)
        elseif trilinear_case5_map_chk1(z, x, y)
            return trilinear_case_5a(z, x, y, q)
        elseif trilinear_case5_map_chk2(x, z, y)
            return trilinear_case_5b(x, z, y, q)
        elseif trilinear_case5_map_chk2(z, x, y)
            return trilinear_case_5b(z, x, y, q)
        end
    elseif (@is_neg(x) && @is_pos(y) && @is_pos(z))
        if trilinear_case5_map_chk1(y, z, x)
            return trilinear_case_5a(y, z, x, q)
        elseif trilinear_case5_map_chk1(z, y, x)
            return trilinear_case_5a(z, y, x, q)
        elseif trilinear_case5_map_chk2(y, z, x)
            return trilinear_case_5b(y, z, x, q)
        elseif trilinear_case5_map_chk2(z, y, x)
            return trilinear_case_5b(z, y, x, q)
        end
    end

    if @is_pos(x) && @is_mix(y) && @is_neg(z)
        return trilinear_case_6(x, y, z, q)
    elseif @is_pos(x) && @is_neg(y) && @is_mix(z)
        return trilinear_case_6(x, z, y, q)
    elseif @is_neg(x) && @is_mix(y) && @is_pos(z)
        return trilinear_case_6(z, y, x, q)
    elseif @is_neg(x) && @is_pos(y) && @is_mix(z)
        return trilinear_case_6(y, z, x, q)
    elseif @is_mix(x) && @is_pos(y) && @is_neg(z)
        return trilinear_case_6(y, x, z, q)
    elseif @is_mix(x) && @is_neg(y) && @is_pos(z)
        return trilinear_case_6(z, x, y, q)
    end

    if (@is_mix(x) && @is_mix(y) && @is_neg(z))
        return trilinear_case_7(x, y, z, q)
    elseif (@is_mix(x) && @is_neg(y) && @is_mix(z))
        return trilinear_case_7(x, z, y, q)
    elseif (@is_neg(x) && @is_mix(y) && @is_mix(z))
        return trilinear_case_7(y, z, x, q)
    end

    if @is_pos(x) && @is_neg(y) && @is_neg(z)
        if trilinear_case8_map_chk(x, y, z)
            return trilinear_case_8a(x, y, z, q)
        else
            return trilinear_case_8b(x, z, y, q)
        end
    elseif @is_neg(x) && @is_pos(y) && @is_neg(z)
        if trilinear_case8_map_chk(y, x, z)
            return trilinear_case_8a(y, x, z, q)
        else
            return trilinear_case_8b(y, z, x, q)
        end
    elseif @is_neg(x) && @is_neg(y) && @is_pos(z)
        if trilinear_case8_map_chk(z, x, y)
            return trilinear_case_8a(z, x, y, q)
        else
            return trilinear_case_8b(z, y, x, q)
        end
    elseif @is_mix(x) && @is_neg(y) && @is_neg(z)
        return trilinear_case_9(x, y, z, q)
    elseif @is_neg(x) && @is_mix(y) && @is_neg(z)
        return trilinear_case_9(y, x, z, q)
    elseif @is_neg(x) && @is_neg(y) && @is_mix(z)
        return trilinear_case_9(z, x, y, q)
    end

    if (@is_neg(x) && @is_neg(y) && @is_neg(z))
        if trilinear_case10_map_chk(x, y, z)
            return trilinear_case_10(x, y, z, q)
        elseif trilinear_case10_map_chk(x, z, y)
            return trilinear_case_10(x, z, y, q)
        elseif trilinear_case10_map_chk(y, x, z)
            return trilinear_case_10(y, x, z, q)
        elseif trilinear_case10_map_chk(y, z, x)
            return trilinear_case_10(y, z, x, q)
        elseif trilinear_case10_map_chk(z, x, y)
            return trilinear_case_10(z, x, y, q)
        else
            return trilinear_case_10(z, y, x, q)
        end
    end

    return trilinear_case_3(x,y,z,q)
end

@inline function trilinear(x1::MC{N,T}, x2::MC{N,T}, x3::MC{N,T}) where {N, T<:Union{NS,MV}}
	#if x1 == x2
	#	if x1 == x3
	#		z = x1.Intv^3
	#	else
	#		z = x_mul_y2(x3.Intv, x1.Intv)
	#	end
	#elseif x2 == x3
	#	z = x_mul_y2(x1.Intv, x2.Intv)
	#else
	#	z = x1.Intv*x2.Intv*x3.Intv
	#end
    z = x1.Intv*x2.Intv*x3.Intv
	return trilinear_kernel(x1, x2, x3, z)
end

#=
*(a::MC{N,T}, b::MC{N,T}, c::MC{N,T}, d::MC{N,T}) where {N, T<:Union{NS,MV}} = *(*(a,b,c),d)

*(a::MC{N,T}, b::MC{N,T}, c::MC{N,T}, d::Float64) where {N, T<:Union{NS,MV}} = *(*(a,b,c),d)
*(a::MC{N,T}, b::MC{N,T}, c::Float64, d::MC{N,T}) where {N, T<:Union{NS,MV}} = *(*(a,b,d),c)
*(a::MC{N,T}, b::Float64, c::MC{N,T}, d::MC{N,T}) where {N, T<:Union{NS,MV}} = *(*(a,c,d),b)
*(a::Float64, b::MC{N,T}, c::MC{N,T}, d::MC{N,T}) where {N, T<:Union{NS,MV}} = *(*(b,c,d),a)

*(a::MC{N,T}, b::MC{N,T}, c::MC{N,T}, d::MC{N,T}, e::MC{N,T}) where {N, T<:Union{NS,MV}} = *(*(a,b,c),d,e)

*(a::MC{N,T}, b::MC{N,T}, c::MC{N,T}, d::MC{N,T}, e::Float64) where {N, T<:Union{NS,MV}} = *(a,b,c,d)*e
*(a::MC{N,T}, b::MC{N,T}, c::MC{N,T}, d::Float64, e::MC{N,T}) where {N, T<:Union{NS,MV}} = *(a,b,c,e)*d
*(a::MC{N,T}, b::MC{N,T}, c::Float64, d::MC{N,T}, e::MC{N,T}) where {N, T<:Union{NS,MV}} = *(a,b,d,e)*c
*(a::MC{N,T}, b::Float64, c::MC{N,T}, d::MC{N,T}, e::MC{N,T}) where {N, T<:Union{NS,MV}} = *(a,c,d,e)*b
*(a::Float64, b::MC{N,T}, c::MC{N,T}, d::MC{N,T}, e::MC{N,T}) where {N, T<:Union{NS,MV}} = *(b,c,d,e)*a

*(a::Float64, b::Float64, c::MC{N,T}, d::MC{N,T}, e::MC{N,T}) where {N, T<:Union{NS,MV}} = *(c,d,e)*a*b
*(a::Float64, b::MC{N,T}, c::Float64, d::MC{N,T}, e::MC{N,T}) where {N, T<:Union{NS,MV}} = *(b,d,e)*a*c
*(a::Float64, b::MC{N,T}, c::MC{N,T}, d::Float64, e::MC{N,T}) where {N, T<:Union{NS,MV}} = *(b,c,e)*a*d
*(a::Float64, b::MC{N,T}, c::MC{N,T}, d::MC{N,T}, e::Float64) where {N, T<:Union{NS,MV}} = *(b,c,d)*a*e

*(a::MC{N,T}, b::Float64, c::Float64, d::MC{N,T}, e::MC{N,T}) where {N, T<:Union{NS,MV}} = *(a,d,e)*b*c
*(a::MC{N,T}, b::Float64, c::MC{N,T}, d::Float64, e::MC{N,T}) where {N, T<:Union{NS,MV}} = *(a,c,e)*b*d
*(a::MC{N,T}, b::Float64, c::MC{N,T}, d::MC{N,T}, e::Float64) where {N, T<:Union{NS,MV}} = *(a,c,d)*b*e

*(a::MC{N,T}, b::MC{N,T}, c::Float64, d::Float64, e::MC{N,T}) where {N, T<:Union{NS,MV}} = *(a,b,e)*c*d
*(a::MC{N,T}, b::MC{N,T}, c::Float64, d::MC{N,T}, e::Float64) where {N, T<:Union{NS,MV}} = *(a,b,d)*c*e
*(a::MC{N,T}, b::MC{N,T}, c::MC{N,T}, d::Float64, e::Float64) where {N, T<:Union{NS,MV}} = *(a,b,c)*d*e

function *(a::MC{N,T}, b::MC{N,T}, c::MC{N,T}, d::MC{N,T}, e::MC{N,T}, f::MC{N,T}) where {N, T<:Union{NS,MV}}
    *(*(a,b,c),*(d,e,f))
end
=#

#=
function *(a,b,c,d,e,f,gs...)
    y = *(a,b,c,d,e,f)
    n = length(js)
    t1 = a
    t2 = b
    t3 = c
    for k = 1:n
        y = @inbounds x[k]
        !is_last
        if mod(k)
        end
    end
    for x in js;
        y = op(y,x);
    end
    y
end
=#

#=
function prod(A::Array{MC{N,T},Q}; dims) where {N, Q, T<:Union{NS,MV}}
end
=#


text:
# Use README as the docstring of the module:
function define_docstring()
    path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    doc = replace(read(path, String), r"^```julia"m => "```jldoctest README")
    @eval TaskDAGAnalyzers $Base.@doc $doc TaskDAGAnalyzers
end


text:
##################################################
#  Syntax testing  (uses examples in the doc)
##################################################

#########  rdiff ###########

m.rdiff( :(x^3) , x=2.)             # first order
m.rdiff( :(x^3) , order = 3, x=2.)  # orders up to 3

m.rdiff( :(sin(x)) , order=10, x=2.)  # derivatives up to order 10

res = m.rdiff( :(sin(x)) , order=10, x=2.)
@eval foo(x) = $res
foo(2.)

ex = :( (1 - x[1])^2 + 100(x[2] - x[1]^2)^2 )  # the rosenbrock function
res = m.rdiff(ex, x=zeros(2), order=2)
m.@eval foo(x) = $res
foo([0.5, 2.])

#########  rdiff (function) ###########

rosenbrock(x) = (1 - x[1])^2 + 100(x[2] - x[1]^2)^2   # function to be derived
rosen2 = m.rdiff(rosenbrock, (ones(2),), order=2)       # orders up to 2
rosen2([1,2])
rosen2([0.5,2])

test(x) = exp(x)
m.rdiff(test, (1.,), order=5)

m.rdiff( :(x^3) , x=2.)             # first order
m.rdiff( :(x^3) , order = 3, x=2.)  # orders up to 3

m.rdiff( :(sin(x)) , order=10, x=2.)  # derivatives up to order 10

res = m.rdiff( :(sin(x)) , order=10, x=2.)
@eval foo(x) = $res
foo(2.)

ex = :( (1 - x[1])^2 + 100(x[2] - x[1]^2)^2 )  # the rosenbrock function
res = m.rdiff(ex, x=zeros(2), order=2)
m.@eval foo(x) = $res
foo([0.5, 2.])

#########  @deriv_rule  ###########

m.@deriv_rule *(x::Real         , y::Real )            y     x * ds
m.@deriv_rule *(x::Real         , y::AbstractArray)    y     x .* ds
m.@deriv_rule *(x::AbstractArray, y::Real )            y     sum(x .* ds)
m.@deriv_rule *(x::AbstractArray, y::AbstractArray)    y     x' * ds

foo(x) = log(1+sin(x))

m.@deriv_rule foo(x)   x   cos(x) / ( 1 + sin(x)) * ds

res = m.rdiff( :( 2 ^ foo(x) ) , x=1)
@eval myf(x) = $res

(myf(1.0)[1] - myf(0.999)[1]) * 1000

######### @typeequiv ############

module Sandbox2
    type Bar
        x
        y
    end

    norm(z::Bar) = z.x*z.x + z.y*z.y
end

ex = quote
    z = Sandbox2.Bar(2^a, sin(a))
    Sandbox2.norm(z)
end

m.@deriv_rule  Sandbox2.Bar(x,y)      x  ds[1]   # Derivative accumulator of x is increased by ds[1]
m.@deriv_rule  Sandbox2.Bar(x,y)      y  ds[2]   # Derivative accumulator of y is increased by ds[2]

m.@deriv_rule  Sandbox2.norm(z::Sandbox2.Bar)  z  Any[ 2*z.x*ds , 2*z.y*ds ]  # Note : produces a 2-vector since z is a Bar

res = m.rdiff(ex, a=0.)
@eval df(a) = $res

df(1)

###### internals - tograph example  ######

m.plot( m.tograph( :(a = sin(x+1) ; 2exp(a)) ))

###### higher order   ######

ex = :( (1 - x[1])^2 + 100(x[2] - x[1]^2)^2 )  # the rosenbrock function

res = m.rdiff(ex, x=zeros(2), order=2)
res = m.rdiff(ex, x=zeros(2), order=3)


text:
module CompScienceMeshes

using DelimitedFiles
using LinearAlgebra
using SparseArrays

using StaticArrays
using Compat
using Requires

import Base.getindex

export getcommonedge

# defaults
export index
export euclidianbasis, point

# default mesh creation
export mesh, readmesh, writemesh
export meshsegment, meshrectangle, meshcircle, meshsphere, meshcuboid
export subdMesh
# mesh interface
export dimension, universedimension, vertextype, coordtype
export numvertices, vertices
export numcells, cells
export boundary, skeleton, vertextocellmap
export connectivity, cellpairs # marked for deprecation

# mesh transforms
export translate, translate!, rotate, rotate!
export fliporientation!, fliporientation
export weld

# mesh refinement
export barycentric_refinement, bisecting_refinement

export  Loop_subdivision

# submesh selection and cell retrieval
export GSubdMesh,submesh
export octree, boundingbox
export overlap_gpredicate, interior_tpredicate, inclosure_gpredicate
export interior_vpredicate

export mapper, restriction

# link to the chart concept
export chart

export domain # marked for deprecation
export dimension, universedimension
export volume
export neighborhood
export quadpoints

export isinside, isinclosure, overlap

# specific to simplicial charts
export simplex, center, vertices
export barytocart, carttobary

export intersection

# the neighborhood concept
export barycentric

#export paramtype
export cartesian, parametric
export jacobian, tangents, normal

export trgauss, sqgauss, legendre

# marked for deprecation
export SegmentedAxis
export minmaxdist, rings, ring

using SparseArrays

Pt{N,T} = StaticArrays.SVector{N,T}

include("defaults.jl")
include("utils.jl")
# include("combinatorics.jl")

# quadrature rules for segements, triangles, and squares
include("quadrature/SegmentGauss.jl")
include("quadrature/TriangleGauss.jl")
include("quadrature/SquareGauss.jl")

# mesh component
include("mesh.jl")
include("scomplex.jl")
include("meshes/flippedmesh.jl")
include("meshes/embedding.jl")
include("subdMesh.jl")

# simplices and related algorithms
include("rectangle.jl")
include("charts.jl")
include("subd_chart.jl")
include("sphere.jl")
include("overlap.jl")
include("intersect.jl")
include("isinside.jl")
include("findchart.jl")
include("neighborhood.jl")
include("subd_neighborhood.jl")
include("quadpoints.jl")

include("submesh.jl")
# include("flipped.jl")
include("timeaxis.jl")
include("subd_shape.jl")
include("gaussquarature.jl")
include("fileio/TRI_mesh.jl")
include("fileio/readmesh.jl")
include("fileio/gmsh.jl")
include("fileio/gid.jl")
include("primitives.jl")
#include("../examples/waveguide_with_post.jl")
include("baryref.jl")
include("subdivision.jl")
include("weld.jl")

include("mapper.jl")
include("restrict.jl")

include("plotlyjs_glue.jl")

end # module


text:
using MultiFloats
using Printf

function sqrt_precision(x::MultiFloat{T,N}) where {T,N}
    setprecision(BigFloat, 2200) do
        approx = BigFloat(sqrt(x))
        exact = sqrt(BigFloat(x))
        -Float64(log2(abs((approx - exact) / exact)))
    end
end

function exp_precision(x::MultiFloat{T,N}) where {T,N}
    setprecision(BigFloat, 2200) do
        approx = BigFloat(exp(x))
        exact = exp(BigFloat(x))
        -Float64(log2(abs((approx - exact) / exact)))
    end
end

num_trials = 10^6

println("SQRT (CLEAN)")
use_clean_multifloat_arithmetic()
for N = 2 : 8
    correctbits = floor(Int, minimum(
            sqrt_precision(ldexp(rand(Float64x{N}), rand(-300:300)))
            for _ = 1 : num_trials))
    @printf("Float64x%d: %.3d bits (%2d lost, %.3f%% loss)\n", N, correctbits,
        precision(Float64x{N}) - correctbits,
        100 * abs(precision(Float64x{N}) - correctbits) / precision(Float64x{N}))
end
println()

println("SQRT (STANDARD)")
use_standard_multifloat_arithmetic()
for N = 2 : 8
    correctbits = floor(Int, minimum(
            sqrt_precision(ldexp(rand(Float64x{N}), rand(-300:300)))
            for _ = 1 : num_trials))
    @printf("Float64x%d: %.3d bits (%2d lost, %.3f%% loss)\n", N, correctbits,
        precision(Float64x{N}) - correctbits,
        100 * abs(precision(Float64x{N}) - correctbits) / precision(Float64x{N}))
end
println()

println("SQRT (SLOPPY)")
use_sloppy_multifloat_arithmetic()
for N = 2 : 8
    correctbits = floor(Int, minimum(
            sqrt_precision(ldexp(rand(Float64x{N}), rand(-300:300)))
            for _ = 1 : num_trials))
    @printf("Float64x%d: %.3d bits (%2d lost, %.3f%% loss)\n", N, correctbits,
        precision(Float64x{N}) - correctbits,
        100 * abs(precision(Float64x{N}) - correctbits) / precision(Float64x{N}))
end
println()

println("EXP (CLEAN)")
use_clean_multifloat_arithmetic()
for N = 2 : 8
    correctbits = floor(Int, minimum(
            exp_precision(400 * rand(Float64x{N}) - 200)
            for _ = 1 : num_trials))
    @printf("Float64x%d: %.3d bits (%2d lost, %.3f%% loss)\n", N, correctbits,
        precision(Float64x{N}) - correctbits,
        100 * abs(precision(Float64x{N}) - correctbits) / precision(Float64x{N}))
end
println()

println("EXP (STANDARD)")
use_standard_multifloat_arithmetic()
for N = 2 : 8
    correctbits = floor(Int, minimum(
            exp_precision(400 * rand(Float64x{N}) - 200)
            for _ = 1 : num_trials))
    @printf("Float64x%d: %.3d bits (%2d lost, %.3f%% loss)\n", N, correctbits,
        precision(Float64x{N}) - correctbits,
        100 * abs(precision(Float64x{N}) - correctbits) / precision(Float64x{N}))
end
println()

println("EXP (SLOPPY)")
use_sloppy_multifloat_arithmetic()
for N = 2 : 8
    correctbits = floor(Int, minimum(
            exp_precision(400 * rand(Float64x{N}) - 200)
            for _ = 1 : num_trials))
    @printf("Float64x%d: %.3d bits (%2d lost, %.3f%% loss)\n", N, correctbits,
        precision(Float64x{N}) - correctbits,
        100 * abs(precision(Float64x{N}) - correctbits) / precision(Float64x{N}))
end
println()


text:
@test interwine([1,2], [3,4]) == [1,3,2,4]
@test interwine([1 2; 3 4], [5 6; 7 8]) == [1 5 2 6; 3 7 4 8]

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
# Matrix times row vector
#   cj^T := ai^T * B + cj^T
# with inner ai^T * B computed using dot products:
#   (ai^T * B)_j = ai^T * bj
function mygemm_ijp(C, A, B)
  # C = C + A * B
  (m, n) = size(C)
  (~, k) = size(A)
  @inbounds for i = 1:m
    for j = 1:n
      for p = 1:k
        C[i, j] += A[i, p] * B[p, j]
      end
    end
  end
end

# Matrix times row vector
#   cj^T := ai^T * B + cj^T
# with inner ai^T * B computed using axpy:
#   ai^T * B = [ai1 * b1^T;
#               ai2 * b2^T;
#               ...]
function mygemm_ipj(C, A, B)
  # C = C + A * B
  (m, n) = size(C)
  (~, k) = size(A)
  @inbounds for i = 1:m
    for p = 1:k
      for j = 1:n
        C[i, j] += A[i, p] * B[p, j]
      end
    end
  end
end

# Rank one update
#   C := C + ap * bp^T
# with rank-1 updates computed using axpy with vector bp
#   (ap * bp^T)_i^T = aip * bp^T
function mygemm_pij(C, A, B)
  # C = C + A * B
  (m, n) = size(C)
  (~, k) = size(A)
  @inbounds for p = 1:k
    for i = 1:m
      for j = 1:n
        C[i, j] += A[i, p] * B[p, j]
      end
    end
  end
end

# Rank one update
#   C := C + ap * bp^T
# with rank-1 updates computed using axpy with vector ap
#   (ap * bp^T)_j = bpj * ap
function mygemm_pji(C, A, B)
  # C = C + A * B
  (m, n) = size(C)
  (~, k) = size(A)
  @inbounds for p = 1:k
    for j = 1:n
      for i = 1:m
        C[i, j] += A[i, p] * B[p, j]
      end
    end
  end
end

# Matrix times column vector
#   cj := A * bj + cj
# with inner A * bj computed using axpy:
#   A * bj = [a1 * b1j | (a2 * b2j | ... ]
function mygemm_jpi(C, A, B)
  # C = C + A * B
  (m, n) = size(C)
  (~, k) = size(A)
  @inbounds for j = 1:n
    for p = 1:k
      for i = 1:m
        C[i, j] += A[i, p] * B[p, j]
      end
    end
  end
end

# Matrix times column vector
#   cj := A * bj + cj
# with inner A * bj computed using dot products:
#   (A * bj)_i = ai^T bj
function mygemm_jip(C, A, B)
  # C = C + A * B
  (m, n) = size(C)
  (~, k) = size(A)
  @inbounds for j = 1:n
    for i = 1:m
      for p = 1:k
        C[i, j] += A[i, p] * B[p, j]
      end
    end
  end
end

max_stars_repo_name: jkozdon/Julia-LAFF-On-PfHP

max_issues_repo_name: jkozdon/Julia-LAFF-On-PfHP

max_forks_repo_name: jkozdon/Julia-LAFF-On-PfHP