text:
import ApproxFun: dotu

# This solves via forward substitution
function forwardsubstitution!(ret,B,μ1,μ2,filter=identity)
    n = length(ret)
    if n≥1
        ret[1] = filter(μ1)
    end
    if n≥2
        ret[2] = filter(μ2)
    end
    if n≥3
        @inbounds for k=2:n-1
            ret[k+1] = filter(-(B[k,k-1]*ret[k-1]+B[k,k]*ret[k])/B[k,k+1])
        end
    end
    ret
end

forwardsubstitution(R,n,μ1,μ2) =
    forwardsubstitution!(Array{promote_type(eltype(R),typeof(μ1),typeof(μ2))}(undef,n),R,μ1,μ2)

# This solves as a boundary value provblem
stieltjesbackward(S::Space,z::Number) = JacobiZ(S,z)\[stieltjesmoment(S,0,z)]

stieltjesbackward!(ret,S::Space,z::Number) = ret[:] = stieltjesbackward(S,z)[1:length(ret)]

stieltjesforward(sp::Space,n,z) = forwardsubstitution(JacobiZ(sp,undirected(z)),n,
                                                      stieltjesmoment(sp,0,z),
                                                      stieltjesmoment(sp,1,z))

stieltjesforward!(ret,sp::Space,z,filter=identity) = forwardsubstitution!(ret,JacobiZ(sp,undirected(z)),
                                                    stieltjesmoment(sp,0,z),
                                                    stieltjesmoment(sp,1,z),filter)

hilbertforward(sp::Space,n,z) = forwardsubstitution(JacobiZ(sp,z),n,
                                                    hilbertmoment(sp,0,z),
                                                    hilbertmoment(sp,1,z))

function stieltjesmoment!(ret,S::PolynomialSpace{<:IntervalOrSegment},z,filter=identity)
    if domain(S) == ChebyshevInterval()
        n = length(ret)
        tol = 1/floor(Int,sqrt(n))
        if (abs(real(z)) ≤ 1+tol) && (abs(imag(z)) ≤ tol)
            cfs = stieltjesforward!(ret,S,z,filter)
        else
            cfs = stieltjesbackward!(ret,S,z)
        end

        ret
    else
        stieltjesmoment!(ret,setdomain(S,ChebyshevInterval()),mobius(S,z),filter)
    end
end

function stieltjesintervalrecurrence(S,f::AbstractVector,z)
    tol=1/floor(Int,sqrt(length(f)))
    if isinf(z)
        zero(promote_type(typeof(z), eltype(f)))
    elseif (abs(real(z)) ≤ 1+tol) && (abs(imag(z)) ≤ tol)
        cfs = stieltjesforward(S,length(f),z)
        dotu(cfs,f)
    else
        cfs = stieltjesbackward(S,undirected(z))
        dotu(cfs,f)
    end
end

stieltjesintervalrecurrence(S,f::AbstractVector,z::AbstractArray) =
    reshape(promote_type(eltype(f),eltype(z))[ stieltjesintervalrecurrence(S,f,z[i]) for i in eachindex(z) ], size(z))

function stieltjes(S::PolynomialSpace{<:IntervalOrSegment},f,z::Number)
    if domain(S)==ChebyshevInterval()
        #TODO: check tolerance
        stieltjesintervalrecurrence(Legendre(),coefficients(f,S,Legendre()),z)
    else
        stieltjes(setdomain(S,ChebyshevInterval()),f,mobius(S,z))
    end
end

function hilbert(S::PolynomialSpace{<:IntervalOrSegment},f,z::Number)
    if domain(S) == ChebyshevInterval()
        L = Legendre(domain(S))
        cfs = hilbertforward(L,length(f),z)
        dotu(cfs,coefficients(f, S, L))
    else
        hilbert(setdomain(S,ChebyshevInterval()),f,mobius(S,z))
    end
end

# Sum over all inverses of fromcanonical, see [Olver,2014]
function stieltjes(S::Space{<:Line},f,z)
    if domain(S)==Line()
        # TODO: rename tocanonical
        stieltjes(setcanonicaldomain(S),f,tocanonical(S,z)) +
            stieltjes(setcanonicaldomain(S),f,(-1-sqrt(1+4z^2))/(2z))
    else
        stieltjes(setdomain(S,Line()),f,mappoint(domain(S),Line(),z))
    end
end

## log kernel

function logkernel(S::PolynomialSpace{<:IntervalOrSegment},v,z::Number)
    if domain(S) == ChebyshevInterval()
        DS=JacobiWeight(1,1,Jacobi(1,1))
        D=Derivative(DS)[2:end,:]

        f=Fun(Fun(S,v),Legendre())  # convert to Legendre expansion
        u=D\(f|(2:∞))   # find integral, dropping first coefficient of f

        (f.coefficients[1]*logabslegendremoment(z) +
            real(stieltjes(Fun(u,Legendre()),z+0im)))/π
    else
        Mp=abs(fromcanonicalD(S,0))
        Mp*logkernel(setcanonicaldomain(S),v,mobius(S,z)) +
            linesum(Fun(S,v))*log(Mp)/π
    end
end


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

include("model_CapabilityInformation.jl")
include("model_CheckNameAvailabilityParameters.jl")
include("model_CreateDataLakeStoreAccountParameters.jl")
include("model_CreateDataLakeStoreAccountProperties.jl")
include("model_CreateFirewallRuleWithAccountParameters.jl")
include("model_CreateOrUpdateFirewallRuleParameters.jl")
include("model_CreateOrUpdateFirewallRuleProperties.jl")
include("model_CreateOrUpdateTrustedIdProviderParameters.jl")
include("model_CreateOrUpdateTrustedIdProviderProperties.jl")
include("model_CreateOrUpdateVirtualNetworkRuleParameters.jl")
include("model_CreateOrUpdateVirtualNetworkRuleProperties.jl")
include("model_CreateTrustedIdProviderWithAccountParameters.jl")
include("model_CreateVirtualNetworkRuleWithAccountParameters.jl")
include("model_DataLakeStoreAccount.jl")
include("model_DataLakeStoreAccountBasic.jl")
include("model_DataLakeStoreAccountListResult.jl")
include("model_DataLakeStoreAccountProperties.jl")
include("model_DataLakeStoreAccountPropertiesBasic.jl")
include("model_EncryptionConfig.jl")
include("model_EncryptionIdentity.jl")
include("model_FirewallRule.jl")
include("model_FirewallRuleListResult.jl")
include("model_FirewallRuleProperties.jl")
include("model_KeyVaultMetaInfo.jl")
include("model_NameAvailabilityInformation.jl")
include("model_Operation.jl")
include("model_OperationDisplay.jl")
include("model_OperationListResult.jl")
include("model_Resource.jl")
include("model_SubResource.jl")
include("model_TrustedIdProvider.jl")
include("model_TrustedIdProviderListResult.jl")
include("model_TrustedIdProviderProperties.jl")
include("model_UpdateDataLakeStoreAccountParameters.jl")
include("model_UpdateDataLakeStoreAccountProperties.jl")
include("model_UpdateEncryptionConfig.jl")
include("model_UpdateFirewallRuleParameters.jl")
include("model_UpdateFirewallRuleProperties.jl")
include("model_UpdateFirewallRuleWithAccountParameters.jl")
include("model_UpdateKeyVaultMetaInfo.jl")
include("model_UpdateTrustedIdProviderParameters.jl")
include("model_UpdateTrustedIdProviderProperties.jl")
include("model_UpdateTrustedIdProviderWithAccountParameters.jl")
include("model_UpdateVirtualNetworkRuleParameters.jl")
include("model_UpdateVirtualNetworkRuleProperties.jl")
include("model_UpdateVirtualNetworkRuleWithAccountParameters.jl")
include("model_Usage.jl")
include("model_UsageListResult.jl")
include("model_UsageName.jl")
include("model_VirtualNetworkRule.jl")
include("model_VirtualNetworkRuleListResult.jl")
include("model_VirtualNetworkRuleProperties.jl")

ext: jl


text:
# Test name methods
@testset "Basics" begin
    # initialize model and variable
    m = InfiniteModel()
    num = Float64(0)
    info = VariableInfo(false, num, false, num, false, num, false, num, false, false)
    new_info = VariableInfo(true, 0., true, 0., true, 0., true, 0., true, false)
    bounds = ParameterBounds()
    var = HoldVariable(info, bounds)
    object = VariableData(var, "test")
    idx = HoldVariableIndex(1)
    vref = HoldVariableRef(m, idx)
    gvref = GeneralVariableRef(m, 1, HoldVariableIndex)
    bad_vref = HoldVariableRef(m, HoldVariableIndex(-1))
    # JuMP.owner_model
    @testset "JuMP.owner_model" begin
        @test owner_model(vref) === m
        @test owner_model(gvref) === m
    end
    # JuMP.index
    @testset "JuMP.index" begin
        @test index(vref) == idx
        @test index(gvref) == idx
    end
    # dispatch_variable_ref
    @testset "dispatch_variable_ref" begin
        @test dispatch_variable_ref(m, idx) == vref
        @test dispatch_variable_ref(gvref) == vref
    end
    # _add_data_object
    @testset "_add_data_object" begin
        @test InfiniteOpt._add_data_object(m, object) == idx
    end
    # _data_dictionary
    @testset "_data_dictionary" begin
        @test InfiniteOpt._data_dictionary(m, HoldVariable) === m.hold_vars
        @test InfiniteOpt._data_dictionary(vref) === m.hold_vars
        @test InfiniteOpt._data_dictionary(gvref) === m.hold_vars
    end
    # JuMP.is_valid
    @testset "JuMP.is_valid" begin
        @test is_valid(m, vref)
        @test is_valid(m, gvref)
    end
    # _data_object
    @testset "_data_object" begin
        @test InfiniteOpt._data_object(vref) === object
        @test InfiniteOpt._data_object(gvref) === object
        @test_throws ErrorException InfiniteOpt._data_object(bad_vref)
    end
    # _core_variable_object
    @testset "_core_variable_object" begin
        @test InfiniteOpt._core_variable_object(vref) === var
        @test InfiniteOpt._core_variable_object(gvref) === var
    end
    @testset "_variable_info" begin
        @test InfiniteOpt._variable_info(vref) == info
    end
    # _update_variable_info
    @testset "_update_variable_info" begin
        @test isa(InfiniteOpt._update_variable_info(vref, new_info), Nothing)
        @test InfiniteOpt._variable_info(vref) == new_info
    end
    # _set_core_variable_object
    @testset "_set_core_variable_object" begin
        @test InfiniteOpt._set_core_variable_object(vref, var) isa Nothing
    end
    # JuMP.name
    @testset "JuMP.name" begin
        @test name(vref) == "test"
        @test name(gvref) == "test"
        @test name(bad_vref) == ""
    end
    # JuMP.set_name
    @testset "JuMP.set_name" begin
        @test isa(set_name(vref, "new"), Nothing)
        @test name(vref) == "new"
        @test isa(set_name(gvref, "new2"), Nothing)
        @test name(vref) == "new2"
        @test_throws ErrorException set_name(bad_vref, "")
    end
    # test parameter_refs
    @testset "parameter_refs" begin
        @test parameter_refs(vref) == ()
        @test parameter_refs(gvref) == ()
    end
    # _make_variable_ref
    @testset "_make_variable_ref" begin
        @test InfiniteOpt._make_variable_ref(m, idx) == gvref
    end
    # _var_name_dict
    @testset "_var_name_dict" begin
        @test InfiniteOpt._var_name_dict(m) isa Nothing
    end
    # _update_var_name_dict
    @testset "_update_var_name_dict" begin
        m.name_to_var = Dict{String, ObjectIndex}()
        dict = InfiniteOpt._data_dictionary(vref)
        @test InfiniteOpt._update_var_name_dict(m, dict) isa Nothing
        @test InfiniteOpt._var_name_dict(m) == Dict(name(vref) => idx)
        m.name_to_var = nothing
    end
    # parameter_by_name
    @testset "JuMP.variable_by_name" begin
        # test normal
        @test variable_by_name(m, "new2") == gvref
        @test variable_by_name(m, "bob") isa Nothing
        # prepare variable with same name
        idx2 = HoldVariableIndex(2)
        @test InfiniteOpt._add_data_object(m, object) == idx2
        vref2 = HoldVariableRef(m, idx2)
        @test set_name(vref2, "new2") isa Nothing
        # test multiple name error
        @test_throws ErrorException variable_by_name(m, "new2")
    end
    # _delete_data_object
    @testset "_delete_data_object" begin
        @test InfiniteOpt._delete_data_object(vref) isa Nothing
        @test length(InfiniteOpt._data_dictionary(vref)) == 1
        @test !is_valid(m, vref)
    end
end

# Test variable definition methods
@testset "Definition" begin
    # initialize model and info
    m = InfiniteModel()
    num = Float64(0)
    info = VariableInfo(false, num, false, num, false, num, false, num, false, false)
    info2 = VariableInfo(true, num, true, num, true, num, true, num, true, false)
    info3 = VariableInfo(true, num, true, num, true, num, true, num, false, true)
    @infinite_parameter(m, 0 <= par <= 10)
    @infinite_parameter(m, 0 <= pars[1:2] <= 10)
    @finite_parameter(m, fin, 42)
    # test _check_bounds
    @testset "_check_bounds" begin
        # test normal
        bounds = ParameterBounds((par => IntervalSet(0, 1),))
        @test InfiniteOpt._check_bounds(bounds) isa Nothing
        bounds = ParameterBounds((pars => IntervalSet(0, 1),))
        @test InfiniteOpt._check_bounds(bounds) isa Nothing
        bounds = ParameterBounds((pars => IntervalSet(0, 0),))
        @test InfiniteOpt._check_bounds(bounds) isa Nothing
        bounds = ParameterBounds((pars[1] => IntervalSet(0, 1),))
        @test InfiniteOpt._check_bounds(bounds) isa Nothing
        # test errors
        bounds = ParameterBounds((fin => IntervalSet(0, 1),))
        @test_throws ErrorException InfiniteOpt._check_bounds(bounds)
        bounds = ParameterBounds((par => IntervalSet(-1, 1),))
        @test_throws ErrorException InfiniteOpt._check_bounds(bounds)
        bounds = ParameterBounds((par => IntervalSet(0, 11),))
        @test_throws ErrorException InfiniteOpt._check_bounds(bounds)
        bounds = ParameterBounds((pars[1] => IntervalSet(0, 0),))
        @test_throws ErrorException InfiniteOpt._check_bounds(bounds)
        bounds = ParameterBounds((pars[1] => IntervalSet(0, 0),
                                 pars[2] => IntervalSet(0, 3)))
        @test_throws ErrorException InfiniteOpt._check_bounds(bounds)
    end
    # _make_variable
    @testset "_make_variable" begin
        # test normal
        @test InfiniteOpt._make_variable(error, info, Hold).info == info
        bounds = ParameterBounds((pars => IntervalSet(0, 0),))
        @test InfiniteOpt._make_variable(error, info, Hold,
                                         parameter_bounds = bounds).info == info
        @test InfiniteOpt._make_variable(error, info, Hold,
                            parameter_bounds = bounds).parameter_bounds == bounds
        # test errors
        @test_throws ErrorException InfiniteOpt._make_variable(error, info,
                                                Hold, parameter_values = 3)
        bounds = ParameterBounds((pars[1] => IntervalSet(0, 0),))
        @test_throws ErrorException InfiniteOpt._make_variable(error, info,
                                           Hold, parameter_bounds = bounds)
        bounds = ParameterBounds((par => IntervalSet(-1, -1),))
        @test_throws ErrorException InfiniteOpt._make_variable(error, info,
                                           Hold, parameter_bounds = bounds)
    end
    # build_variable
    @testset "JuMP.build_variable" begin
        # test normal
        @test build_variable(error, info, Hold).info == info
        @test build_variable(error, info, Hold).parameter_bounds == ParameterBounds()
        bounds = ParameterBounds((pars => IntervalSet(0, 5),))
        @test build_variable(error, info, Hold,
                             parameter_bounds = bounds).info == info
        @test build_variable(error, info, Hold,
                             parameter_bounds = bounds).parameter_bounds == bounds
        # test errors
        @test_throws ErrorException build_variable(error, info, Point,
                                                   parameter_bounds = bounds)
        bounds = ParameterBounds((par => IntervalSet(0, 0),
                                  pars => IntervalSet(-1, 1)))
        @test_throws ErrorException build_variable(error, info, Hold,
                                                   parameter_bounds = bounds)
    end
    # _validate_bounds
    @testset "_validate_bounds" begin
        # test normal
        bounds = ParameterBounds((par => IntervalSet(0, 1),))
        @test InfiniteOpt._validate_bounds(m, bounds) isa Nothing
        @test num_supports(par) == 0
        # test error
        @infinite_parameter(InfiniteModel(), par2 in Normal())
        bounds = ParameterBounds((par2 => IntervalSet(0, 1),))
        @test_throws VariableNotOwned{GeneralVariableRef} InfiniteOpt._validate_bounds(m, bounds)
        # test support addition
        bounds = ParameterBounds((par => IntervalSet(0, 0),))
        @test InfiniteOpt._validate_bounds(m, bounds) isa Nothing
        @test supports(par, label = UserDefined) == [0]
        bounds = ParameterBounds((pars => IntervalSet(0, 0),))
        @test InfiniteOpt._validate_bounds(m, bounds) isa Nothing
        @test supports(pars) == zeros(2, 1)
    end
    # _check_and_make_variable_ref
    @testset "_check_and_make_variable_ref" begin
        # test normal
        v = build_variable(error, info, Hold)
        idx = HoldVariableIndex(1)
        vref = HoldVariableRef(m, idx)
        @test InfiniteOpt._check_and_make_variable_ref(m, v, "") == vref
        # test with bounds
        bounds = ParameterBounds((par => IntervalSet(0, 2),))
        v = build_variable(error, info, Hold, parameter_bounds = bounds)
        idx = HoldVariableIndex(2)
        vref = HoldVariableRef(m, idx)
        @test InfiniteOpt._check_and_make_variable_ref(m, v, "") == vref
        @test m.has_hold_bounds
        m.has_hold_bounds = false
        # test bad bounds
        @infinite_parameter(InfiniteModel(), par2 in [0, 2])
        v = build_variable(error, info, Hold,
            parameter_bounds = ParameterBounds((par2 => IntervalSet(0, 1),)))
        @test_throws VariableNotOwned{GeneralVariableRef} InfiniteOpt._check_and_make_variable_ref(m, v, "")
    end
    # add_variable
    @testset "JuMP.add_variable" begin
        idx = HoldVariableIndex(3)
        vref = HoldVariableRef(m, idx)
        gvref = InfiniteOpt._make_variable_ref(m, idx)
        v = build_variable(error, info, Hold)
        @test add_variable(m, v, "name") == gvref
        @test haskey(InfiniteOpt._data_dictionary(vref), idx)
        @test name(vref) == "name"
        # prepare infinite variable with all the possible info additions
        v = build_variable(error, info2, Hold)
        # test info addition functions
        idx = HoldVariableIndex(4)
        vref = HoldVariableRef(m, idx)
        gvref = InfiniteOpt._make_variable_ref(m, idx)
        @test add_variable(m, v, "name") == gvref
        @test !optimizer_model_ready(m)
        @test !m.has_hold_bounds
        # lower bound
        cindex = ConstraintIndex(1)
        cref = InfOptConstraintRef(m, cindex)
        @test has_lower_bound(vref)
        @test InfiniteOpt._lower_bound_index(vref) == cindex
        @test constraint_object(cref) isa ScalarConstraint{GeneralVariableRef,
                                                 MOI.GreaterThan{Float64}}
        @test InfiniteOpt._data_object(cref).is_info_constraint
        # upper bound
        cindex = ConstraintIndex(2)
        cref = InfOptConstraintRef(m, cindex)
        @test has_upper_bound(vref)
        @test InfiniteOpt._upper_bound_index(vref) == cindex
        @test constraint_object(cref) isa ScalarConstraint{GeneralVariableRef,
                                                 MOI.LessThan{Float64}}
        @test InfiniteOpt._data_object(cref).is_info_constraint
        # fix
        cindex = ConstraintIndex(3)
        cref = InfOptConstraintRef(m, cindex)
        @test is_fixed(vref)
        @test InfiniteOpt._fix_index(vref) == cindex
        @test constraint_object(cref) isa ScalarConstraint{GeneralVariableRef,
                                                 MOI.EqualTo{Float64}}
        @test InfiniteOpt._data_object(cref).is_info_constraint
        # binary
        cindex = ConstraintIndex(4)
        cref = InfOptConstraintRef(m, cindex)
        @test is_binary(vref)
        @test InfiniteOpt._binary_index(vref) == cindex
        @test constraint_object(cref) isa ScalarConstraint{GeneralVariableRef,
                                                 MOI.ZeroOne}
        @test InfiniteOpt._data_object(cref).is_info_constraint
        @test InfiniteOpt._constraint_dependencies(vref) == [ConstraintIndex(i)
                                                             for i = 1:4]
        # prepare infinite variable with integer info addition
        bounds = ParameterBounds((par => IntervalSet(0, 2),
                                  pars => IntervalSet(1, 1)))
        v = build_variable(error, info3, Hold, parameter_bounds = bounds)
        # test integer addition functions
        idx = HoldVariableIndex(5)
        vref = HoldVariableRef(m, idx)
        gvref = InfiniteOpt._make_variable_ref(m, idx)
        @test add_variable(m, v, "name") == gvref
        @test !optimizer_model_ready(m)
        @test m.has_hold_bounds
        @test supports(par) == [0]
        @test sortcols(supports(pars)) == [0 1; 0 1]
        cindex = ConstraintIndex(8)
        cref = InfOptConstraintRef(m, cindex)
        @test is_integer(vref)
        @test InfiniteOpt._integer_index(vref) == cindex
        @test constraint_object(cref) isa BoundedScalarConstraint{GeneralVariableRef,
                                                            MOI.Integer}
        @test InfiniteOpt._data_object(cref).is_info_constraint
        @test InfiniteOpt._constraint_dependencies(vref) == [ConstraintIndex(i)
                                                             for i = 5:8]
    end
end

# Test the hold variable macro if no bounds are provided
@testset "Macro" begin
    # initialize model
    m = InfiniteModel()
    @testset "@hold_variable (Not Bounded)" begin
        # test regular
        idx = HoldVariableIndex(1)
        vref = HoldVariableRef(m, idx)
        gvref = InfiniteOpt._make_variable_ref(m, idx)
        @test @hold_variable(m, x >= 1, Bin) == gvref
        @test name(vref) == "x"
        @test lower_bound(vref) == 1
        @test is_binary(vref)
        # test anan
        idx = HoldVariableIndex(2)
        vref = HoldVariableRef(m, idx)
        gvref = InfiniteOpt._make_variable_ref(m, idx)
        @test @hold_variable(m, binary = true, lower_bound = 1,
                               base_name = "x") == gvref
        @test name(vref) == "x"
        @test lower_bound(vref) == 1
        @test is_binary(vref)
        # test array
        idxs = [HoldVariableIndex(3), HoldVariableIndex(4)]
        vrefs = [HoldVariableRef(m, idx) for idx in idxs]
        gvrefs = [InfiniteOpt._make_variable_ref(m, idx) for idx in idxs]
        @test @hold_variable(m, y[1:2] == 2, Int) == gvrefs
        @test name(vrefs[1]) == "y[1]"
        @test fix_value(vrefs[2]) == 2
        @test is_integer(vrefs[1])
        # test errors
        @test_throws AssertionError @hold_variable(Model(), z >= 1, Bin)
        @test_macro_throws ErrorException @hold_variable(m, x >= 1, Bin)
    end
    # test _make_interval_set
    @testset "_make_interval_set" begin
        expected = :(IntervalSet(0, 1))
        @test InfiniteOpt._make_interval_set(error, :([0, 1])) == expected
        @test_throws ErrorException InfiniteOpt._make_interval_set(error, :([0]))
    end
    # test _make_bound_pair (in)
    @testset "_make_bound_pair (in)" begin
        expr = :(t in [0, 1])
        set = :(IntervalSet(0, 1))
        tp = :(t)
        @test InfiniteOpt._make_bound_pair(error, expr,
                                           Val(:in)) == :($(tp) => $(set))
    end
    # test _make_bound_pair (==)
    @testset "_make_bound_pair (==)" begin
        expr = :(t == 0)
        set = :(IntervalSet(0, 0))
        tp = :(t)
        @test InfiniteOpt._make_bound_pair(error, expr,
                                           Val(:(==))) == :($(tp) => $(set))
    end
    # test _make_bound_pair (comparison with <=)
    @testset "_make_bound_pair (<=)" begin
        expr = :(0 <= t <= 1)
        set = :(IntervalSet(0, 1))
        tp = :(t)
        @test InfiniteOpt._make_bound_pair(error, expr, Val(:(<=)),
                                           Val(:(<=))) == :($(tp) => $(set))
    end
    # test _make_bound_pair (comparison with >=)
    @testset "_make_bound_pair (>=)" begin
        expr = :(1 >= t >= 0)
        set = :(IntervalSet(0, 1))
        tp = :(t)
        @test InfiniteOpt._make_bound_pair(error, expr, Val(:(>=)),
                                           Val(:(>=))) == :($(tp) => $(set))
    end
    # test _make_bound_pair (comparison fallback)
    @testset "_make_bound_pair (Fallback)" begin
        expr = :(t <= 0)
        @test_throws ErrorException InfiniteOpt._make_bound_pair(error, expr,
                                                                 Val(:(<=)))
        @test_throws ErrorException InfiniteOpt._make_bound_pair(error, expr,
                                                         Val(:(<=)), Val(:(>=)))
    end
    # test _make_bound_pair (wrapper)
    @testset "_make_bound_pair (Wrapper)" begin
        # test :call expression
        expr = :(t in [0, 1])
        set = :(IntervalSet(0, 1))
        tp = :(t)
        @test InfiniteOpt._make_bound_pair(error, expr,
                                           Val(:in)) == :($(tp) => $(set))
        # test :compariosn expression
        expr = :(0 <= t <= 1)
        @test InfiniteOpt._make_bound_pair(error, expr, Val(:(<=)),
                                           Val(:(<=))) == :($(tp) => $(set))
        # test error
        expr = :([2, 3])
        @test_throws ErrorException InfiniteOpt._make_bound_pair(error, expr)
    end
    # test _parse_parameter_bounds (Vector)
    @testset "_parse_parameter_bounds (Vector)" begin
        args = [:(t in [0, 1]), :(0 <= t <= 1)]
        set = :(IntervalSet(0, 1))
        tp = :(t)
        dict_arg = :($(tp) => $(set))
        dict = :(($(dict_arg), $(dict_arg)))
        bounds = :(ParameterBounds($dict))
        @test InfiniteOpt._parse_parameter_bounds(error, args) == bounds
    end
    # test _parse_parameter_bounds (Expr)
    @testset "_parse_parameter_bounds (Expr)" begin
        expr = :(t in [0, 1])
        set = :(IntervalSet(0, 1))
        tp = :(t)
        dict_arg = :($(tp) => $(set))
        dict = :(($(dict_arg), ))
        bounds = :(ParameterBounds($dict))
        @test InfiniteOpt._parse_parameter_bounds(error, expr) == bounds
    end
    # test _extract_bounds (:call)
    @testset "_extract_bounds (:call)" begin
        # test single anonymous bound
        args = [:in, :t, :([0, 1])]
        set = :(IntervalSet(0, 1))
        tp = :(t)
        dict_arg = :($(tp) => $(set))
        bounds = :(ParameterBounds(($(dict_arg),)))
        @test InfiniteOpt._extract_bounds(error, args,
                                          Val(:call)) == (nothing, bounds)
    end
    # test _extract_bounds (:tuple)
    @testset "_extract_bounds (:tuple)" begin
        args = [:(t in [0, 1]), :(0 <= t <= 1)]
        set = :(IntervalSet(0, 1))
        tp = :(t)
        dict_arg = :($(tp) => $(set))
        bounds = :(ParameterBounds(($(dict_arg), $(dict_arg))))
        @test InfiniteOpt._extract_bounds(error, args,
                                          Val(:tuple)) == (nothing, bounds)
    end
    # test _extract_bounds (:comparison)
    @testset "_extract_bounds (:comparison)" begin
        args = [0, :(<=), :t, :(<=), 1]
        set = :(IntervalSet(0, 1))
        tp = :(t)
        dict_arg = :($(tp) => $(set))
        bounds = :(ParameterBounds(($(dict_arg),)))
        @test InfiniteOpt._extract_bounds(error, args,
                                          Val(:comparison)) == (nothing, bounds)
    end
    # test @hold_variable with bounds
    @testset "@hold_variable (Bounded)" begin
        # prepare parameters
        @infinite_parameter(m, 0 <= par <= 10)
        @infinite_parameter(m, 0 <= pars[1:2] <= 10)
        # test regular call
        idx = HoldVariableIndex(5)
        vref = HoldVariableRef(m, idx)
        gvref = InfiniteOpt._make_variable_ref(m, idx)
        @test @hold_variable(m, xb >= 1, Bin,
                             parameter_bounds = (par == 0)) == gvref
        @test name(vref) == "xb"
        @test lower_bound(vref) == 1
        @test is_binary(vref)
        @test m.has_hold_bounds
        @test parameter_bounds(vref) == ParameterBounds(Dict(par => IntervalSet(0, 0)))
        @test supports(par) == [0]
        # test regular tuple
        idx = HoldVariableIndex(6)
        vref = HoldVariableRef(m, idx)
        gvref = InfiniteOpt._make_variable_ref(m, idx)
        @test @hold_variable(m, yb, parameter_bounds = (pars == 0,
                                                       par in [0, 1])) == gvref
        @test name(vref) == "yb"
        dict = Dict(pars[1] => IntervalSet(0, 0), pars[2] => IntervalSet(0, 0),
                    par => IntervalSet(0, 1))
        @test parameter_bounds(vref) == ParameterBounds(dict)
        @test supports(pars[1]) == [0]
        # test comparison
        idx = HoldVariableIndex(7)
        vref = HoldVariableRef(m, idx)
        gvref = InfiniteOpt._make_variable_ref(m, idx)
        @test @hold_variable(m, zb, parameter_bounds = (0 <= par <= 1),
                             base_name = "bob") == gvref
        @test name(vref) == "bob"
        @test parameter_bounds(vref) == ParameterBounds(Dict(par => IntervalSet(0, 1)))
        # test unrecognized format
        @test_macro_throws ErrorException @hold_variable(m, parameter_bounds = par)
        # test container specification
        idxs = [HoldVariableIndex(8), HoldVariableIndex(9)]
        vrefs = [HoldVariableRef(m, idx) for idx in idxs]
        vals = [0, 1]
        @test @hold_variable(m, [i = 1:2], parameter_bounds = (pars == vals[i]),
                             container = SparseAxisArray) isa JuMPC.SparseAxisArray
        @test name(vrefs[1]) == ""
        @test parameter_bounds(vrefs[1]) == ParameterBounds((pars => IntervalSet(0, 0),))
        # test wrong model type
        @test_macro_throws ErrorException @hold_variable(Model(),
                                                  parameter_bounds = (par == 0))
    end
end

# test usage methods
@testset "Usage" begin
    # initialize model and stuff
    m = InfiniteModel()
    @hold_variable(m, y0)
    vref = dispatch_variable_ref(y0)
    # test used_by_measure
    @testset "used_by_measure" begin
        @test !used_by_measure(vref)
        push!(InfiniteOpt._measure_dependencies(vref), MeasureIndex(1))
        @test used_by_measure(y0)
        @test used_by_measure(vref)
        empty!(InfiniteOpt._measure_dependencies(vref))
    end
    # test used_by_constraint
    @testset "used_by_constraint" begin
        @test !used_by_constraint(vref)
        push!(InfiniteOpt._constraint_dependencies(vref), ConstraintIndex(1))
        @test used_by_constraint(y0)
        @test used_by_constraint(vref)
        empty!(InfiniteOpt._constraint_dependencies(vref))
    end
    # test used_by_objective
    @testset "used_by_objective" begin
        @test !used_by_objective(y0)
        @test !used_by_objective(vref)
        InfiniteOpt._data_object(vref).in_objective = true
        @test used_by_objective(vref)
        InfiniteOpt._data_object(vref).in_objective = false
    end
    # test is_used
    @testset "is_used" begin
        # test not used
        @test !is_used(vref)
        # test used by constraint and/or measure
        push!(InfiniteOpt._constraint_dependencies(vref), ConstraintIndex(1))
        @test is_used(y0)
        empty!(InfiniteOpt._constraint_dependencies(vref))
        # test used by objective
        InfiniteOpt._data_object(vref).in_objective = true
        @test is_used(vref)
    end
end

# Test parameter bound queries for hold variables
@testset "Parameter Bounds" begin
    # initialize the model and other needed information
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 10])
    @infinite_parameter(m, pars[1:2] in [0, 10])
    @infinite_variable(m, inf(par))
    dinf = dispatch_variable_ref(inf)
    @hold_variable(m, x == 0)
    dx = dispatch_variable_ref(x)
    @hold_variable(m, y >= 0, parameter_bounds = par == 0)
    dy = dispatch_variable_ref(y)
    data = TestData(par, 0, 5)
    meas = Measure(x + inf, data, Int[], Int[], true)
    object = MeasureData(meas, "test")
    mindex = MeasureIndex(1)
    @test InfiniteOpt._add_data_object(m, object) == mindex
    dmref = MeasureRef(m, mindex)
    mref = InfiniteOpt._make_variable_ref(m, mindex)
    push!(InfiniteOpt._measure_dependencies(dx), mindex)
    # test parameter_bounds
    @testset "parameter_bounds" begin
        @test parameter_bounds(dx) == ParameterBounds()
        @test parameter_bounds(y) == ParameterBounds((par => IntervalSet(0, 0),))
    end
    # test has_parameter_bounds
    @testset "has_parameter_bounds" begin
        @test !has_parameter_bounds(x)
        @test has_parameter_bounds(dy)
        @test !has_parameter_bounds(dinf)
    end
    # test _update_variable_param_bounds
    @testset "_update_variable_param_bounds" begin
        bounds = ParameterBounds((par => IntervalSet(0, 1),))
        @test isa(InfiniteOpt._update_variable_param_bounds(dy, bounds), Nothing)
        @test parameter_bounds(dy) == bounds
    end
    # test _check_meas_bounds
    @testset "_check_meas_bounds" begin
        # test ok
        bounds = ParameterBounds((par => IntervalSet(0, 5),))
        @test isa(InfiniteOpt._check_meas_bounds(bounds, data), Nothing)
        # test lower bad
        bounds = ParameterBounds((par => IntervalSet(2, 5),))
        @test_throws ErrorException InfiniteOpt._check_meas_bounds(bounds, data)
        # test upper bad
        bounds = ParameterBounds((par => IntervalSet(0, 2),))
        @test_throws ErrorException InfiniteOpt._check_meas_bounds(bounds, data)
    end
    # test _check_meas_bounds (Fallback)
    @testset "_check_meas_bounds (Fallback)" begin
        warn = "Unable to check if hold variables bounds are valid in measure " *
               "with measure data type `BadData`. This can be resolved by " *
               "extending `measure_data_in_hold_bounds`."
        @test_logs (:warn, warn) InfiniteOpt._check_meas_bounds(ParameterBounds(),
                                                                BadData())
    end
    # test _update_bounds
    @testset "_update_bounds" begin
        # test normal
        bounds1 = ParameterBounds((pars[1] => IntervalSet(0, 1),
                                   pars[2] => IntervalSet(-1, 2)))
        bounds2 = ParameterBounds((par => IntervalSet(0, 0),
                                   pars[1] => IntervalSet(0.5, 1.5),
                                   pars[2] => IntervalSet(0, 1)))
        @test isa(InfiniteOpt._update_bounds(bounds1, bounds2), Nothing)
        @test length(bounds1) == 3
        @test bounds1[par] == IntervalSet(0, 0)
        @test bounds1[pars[1]] == IntervalSet(0.5, 1)
        @test bounds1[pars[2]] == IntervalSet(0, 1)
        # test error
        new_bounds = ParameterBounds((par => IntervalSet(1, 1),))
        @test_throws ErrorException InfiniteOpt._update_bounds(bounds1, new_bounds)
        new_bounds = ParameterBounds((par => IntervalSet(-1, -1),))
        @test_throws ErrorException InfiniteOpt._update_bounds(bounds1, new_bounds)
    end
    # test _update_var_bounds (DispatchVariableRef)
    @testset "_update_var_bounds (DispatchVariableRef)" begin
        @test isa(InfiniteOpt._update_var_bounds(dinf, ParameterBounds()),
                  Nothing)
    end
    # test _update_var_bounds (HoldVariableRef)
    @testset "_update_var_bounds (HoldVariableRef)" begin
        bounds = ParameterBounds()
        @test InfiniteOpt._update_variable_param_bounds(dx,
                      ParameterBounds((par => IntervalSet(0, 6),))) isa Nothing
        @test isa(InfiniteOpt._update_var_bounds(x, bounds), Nothing)
        @test bounds[par] == IntervalSet(0, 6)
    end
    # test _update_var_bounds (MeasureRef)
    @testset "_update_var_bounds (MeasureRef)" begin
        bounds = ParameterBounds()
        @test isa(InfiniteOpt._update_var_bounds(dmref, bounds), Nothing)
        @test bounds[par] == IntervalSet(0, 6)
        @test InfiniteOpt._update_variable_param_bounds(dx, ParameterBounds()) isa Nothing
    end
    # test _update_var_bounds (GeneralVariableRef)
    @testset "_update_var_bounds (GeneralVariableRef)" begin
        @test isa(InfiniteOpt._update_var_bounds(inf, ParameterBounds()),
                  Nothing)
    end
    # test _rebuild_constr_bounds (BoundedScalarConstraint)
    @testset "_rebuild_constr_bounds (Bounded)" begin
        # test addition
        bounds = ParameterBounds((par => IntervalSet(0, 2),))
        c = BoundedScalarConstraint(y + inf, MOI.EqualTo(0.), copy(bounds),
                                    copy(bounds))
        expected = ParameterBounds((par => IntervalSet(0, 1),))
        @test InfiniteOpt._rebuild_constr_bounds(c, bounds).bounds == expected
        # test deletion
        bounds = ParameterBounds((par => IntervalSet(0, 0),))
        c = BoundedScalarConstraint(y + inf, MOI.EqualTo(0.), bounds,
                                    ParameterBounds())
        InfiniteOpt._update_variable_param_bounds(dy, ParameterBounds())
        @test isa(InfiniteOpt._rebuild_constr_bounds(c, bounds),
                  ScalarConstraint)
        InfiniteOpt._update_variable_param_bounds(dy, bounds)
    end
    # test _rebuild_constr_bounds (ScalarConstraint)
    @testset "_rebuild_constr_bounds (Scalar)" begin
        bounds = ParameterBounds((par => IntervalSet(0, 0),))
        c = ScalarConstraint(y + inf, MOI.EqualTo(0.))
        @test InfiniteOpt._rebuild_constr_bounds(c, bounds).bounds == bounds
    end
    # test set_parameter_bounds
    @testset "set_parameter_bounds" begin
        bounds1 = ParameterBounds((pars[1] => IntervalSet(0, 1),))
        # test force error
        @test_throws ErrorException set_parameter_bounds(y, bounds1)
        # prepare for main test
        cref1 = FixRef(dx)
        bconstr = BoundedScalarConstraint(mref, MOI.EqualTo(0.0), bounds1, bounds1)
        cref2 = InfOptConstraintRef(m, ConstraintIndex(3))
        @test add_constraint(m, bconstr) == cref2
        # test measure error
        bounds = ParameterBounds((par => IntervalSet(1, 2),))
        @test_throws ErrorException set_parameter_bounds(dx, bounds)
        # test constr error
        bounds = ParameterBounds((pars[1] => IntervalSet(2, 3), ))
        @test_throws ErrorException set_parameter_bounds(dx, bounds)
        @test InfiniteOpt._update_variable_param_bounds(dx, ParameterBounds()) isa Nothing
        # test normal
        bounds = ParameterBounds((par => IntervalSet(0, 5),
                                  pars[1] => IntervalSet(0.5, 2)))
        @test isa(set_parameter_bounds(dx, bounds), Nothing)
        @test parameter_bounds(x) == bounds
        @test m.has_hold_bounds
        @test !optimizer_model_ready(m)
        expected = ParameterBounds((par => IntervalSet(0, 5),
                                    pars[1] => IntervalSet(0.5, 1)))
        @test parameter_bounds(cref2) == expected
        @test parameter_bounds(cref1) == bounds
        # test forced
        bounds = ParameterBounds((par => IntervalSet(1, 1),))
        @test isa(set_parameter_bounds(y, bounds, force = true), Nothing)
        @test parameter_bounds(dy) == bounds
        @test parameter_bounds(LowerBoundRef(dy)) == bounds
    end
    # test _update_constr_bounds (BoundedScalarConstraint)
    @testset "_update_constr_bounds (Bounded)" begin
        dict1 = Dict(pars[1] => IntervalSet(0, 1))
        dict2 = Dict(par => IntervalSet(0, 0), pars[1] => IntervalSet(1.1, 1.5))
        constr = BoundedScalarConstraint(par, MOI.EqualTo(0.0), ParameterBounds(copy(dict1)),
                                         ParameterBounds(copy(dict1)))
        # test error
        @test_throws ErrorException InfiniteOpt._update_constr_bounds(ParameterBounds(dict2), constr)
        @test constr.bounds == ParameterBounds(dict1)
        # test normal
        dict2 = Dict(par => IntervalSet(0, 0), pars[1] => IntervalSet(0.5, 1.5))
        dict = Dict(par => IntervalSet(0, 0), pars[1] => IntervalSet(0.5, 1))
        @test InfiniteOpt._update_constr_bounds(ParameterBounds(dict2), constr).bounds == ParameterBounds(dict)
    end
    # test _update_constr_bounds (ScalarConstraint)
    @testset "_update_constr_bounds (Scalar)" begin
        bounds = ParameterBounds((pars[1] => IntervalSet(0, 1),))
        constr = JuMP.ScalarConstraint(par, MOI.EqualTo(0.0))
        @test InfiniteOpt._update_constr_bounds(bounds, constr).bounds == bounds
    end
    # test add_parameter_bound
    @testset "add_parameter_bounds" begin
        cref1 = FixRef(dx)
        cref2 = InfOptConstraintRef(m, ConstraintIndex(3))
        # test adding normally
        new_bounds = ParameterBounds((pars[2] => IntervalSet(0, 5),))
        @test isa(add_parameter_bounds(x, new_bounds), Nothing)
        bounds = ParameterBounds((par => IntervalSet(0, 5),
                                  pars[1] => IntervalSet(0.5, 2),
                                  pars[2] => IntervalSet(0, 5)))
        @test parameter_bounds(dx) == bounds
        expected = ParameterBounds((par => IntervalSet(0, 5),
                                    pars[1] => IntervalSet(0.5, 1),
                                    pars[2] => IntervalSet(0, 5)))
        @test parameter_bounds(cref2) == expected
        @test parameter_bounds(cref1) == bounds
        # test measure error --> these errors are properly tested in the macros
        new_bounds = ParameterBounds((par => IntervalSet(1, 2),))
        @test_throws ErrorException add_parameter_bounds(dx, new_bounds)
        # test constr error
        new_bounds = ParameterBounds((pars[1] => IntervalSet(2, 3),))
        @test_throws ErrorException add_parameter_bounds(x, new_bounds)
        # test overwritting existing
        new_bounds = ParameterBounds((par => IntervalSet(0, 0),))
        @test isa(add_parameter_bounds(y, new_bounds), Nothing)
        @test parameter_bounds(y) == new_bounds
    end
    # test delete_parameter_bounds
    @testset "delete_parameter_bounds" begin
        cref1 = FixRef(dx)
        cref2 = InfOptConstraintRef(m, ConstraintIndex(3))
        # test normal
        @test isa(delete_parameter_bounds(x), Nothing)
        @test parameter_bounds(dx) == ParameterBounds()
        expected = ParameterBounds((pars[1] => IntervalSet(0, 1),))
        @test parameter_bounds(cref2) == expected
        @test !has_parameter_bounds(cref1)
        # test already gone
        @test isa(delete_parameter_bounds(dx), Nothing)
        @test parameter_bounds(dx) == ParameterBounds()
        @test parameter_bounds(cref2) == expected
        @test !has_parameter_bounds(cref1)
    end
end

# Test Parameter Bound Macros
@testset "Parameter Bound Macros" begin
    # initialize the model and other needed information
    m = InfiniteModel()
    @infinite_parameter(m, par in [0, 10])
    @infinite_parameter(m, pars[1:2] in [0, 10])
    @hold_variable(m, x == 0)
    @hold_variable(m, y, parameter_bounds = par == 0)
    data = TestData(par, 0, 5)
    meas = Measure(x + par, data, Int[], Int[], true)
    object = MeasureData(meas, "test")
    mindex = MeasureIndex(1)
    @test InfiniteOpt._add_data_object(m, object) == mindex
    dmref = MeasureRef(m, mindex)
    mref = InfiniteOpt._make_variable_ref(m, mindex)
    push!(InfiniteOpt._measure_dependencies(x), mindex)
    # test @add_parameter_bounds
    @testset "@add_parameter_bounds" begin
        bounds1 = ParameterBounds(Dict(pars[1] => IntervalSet(0, 1)))
        # prepare for main test
        cref1 = FixRef(x)
        bconstr = BoundedScalarConstraint(mref, MOI.EqualTo(0.0), bounds1, bounds1)
        cref2 = InfOptConstraintRef(m, ConstraintIndex(2))
        @test add_constraint(m, bconstr) == cref2
        # test adding normally
        @test isa(@add_parameter_bounds(x, pars[2] in [0, 5]), Nothing)
        dict = Dict(pars[2] => IntervalSet(0, 5))
        @test parameter_bounds(x) == ParameterBounds(dict)
        expected = Dict(pars[1] => IntervalSet(0, 1),
                        pars[2] => IntervalSet(0, 5))
        @test parameter_bounds(cref2) == ParameterBounds(expected)
        @test parameter_bounds(cref1) == ParameterBounds(dict)
        # test measure error
        @test_macro_throws ErrorException @add_parameter_bounds(x, par in [1, 2])
        # test constr error
        @test_macro_throws ErrorException @add_parameter_bounds(x, pars[1] in [2, 3])
        # test bad input format
        @test_macro_throws ErrorException @add_parameter_bounds(x, pars[1] = 0)
        @test_macro_throws ErrorException @add_parameter_bounds(par, pars[1] == 0)
        # test overwritting existing
        @test isa(@add_parameter_bounds(y, par == 1), Nothing)
        @test parameter_bounds(y) == ParameterBounds((par => IntervalSet(1, 1),))
    end
    # test @set_parameter_bounds
    @testset "@set_parameter_bounds" begin
        @test delete_parameter_bounds(x) isa Nothing
        bounds1 = ParameterBounds((pars[1] => IntervalSet(0, 1),))
        # test force error
        @test_macro_throws ErrorException @set_parameter_bounds(y, pars[1] in [0, 1])
        # test bad arguments
        @test_macro_throws ErrorException @set_parameter_bounds(x, x, 23)
        @test_macro_throws ErrorException @set_parameter_bounds(x, pars[1] == 0,
                                                                bad = 42)
        # test bad input format
        @test_macro_throws ErrorException @set_parameter_bounds(x, pars[1] = 0)
        @test_macro_throws ErrorException @set_parameter_bounds(par, pars[1] == 0)
        # prepare for main test
        cref1 = FixRef(x)
        cref2 = InfOptConstraintRef(m, ConstraintIndex(2))
        # test measure error
        @test_macro_throws ErrorException @set_parameter_bounds(x, par in [1, 2])
        # test constr error
        @test_macro_throws ErrorException @set_parameter_bounds(x, pars[1] in [2, 3])
        # test normal
        bounds = ParameterBounds((par => IntervalSet(0, 5),
                                      pars[1] => IntervalSet(0.5, 2)))
        @test isa(@set_parameter_bounds(x, (par in [0, 5], pars[1] in [0.5, 2])),
                  Nothing)
        @test parameter_bounds(x) == bounds
        @test m.has_hold_bounds
        @test !optimizer_model_ready(m)
        expected = ParameterBounds((par => IntervalSet(0, 5),
                                    pars[1] => IntervalSet(0.5, 1)))
        @test parameter_bounds(cref2) == expected
        @test parameter_bounds(cref1) == bounds
        # test forced
        @test isa(@set_parameter_bounds(y, par == 0, force = true), Nothing)
        @test parameter_bounds(y) == ParameterBounds((par => IntervalSet(0, 0),))
    end
end


text:
function atomlist(xi, pimages, wr, B2R)

dim = size(xi,1);
n = size(xi,2);
m = size(pimages,2);
alist = Int64.(Array(1:n));
k = 0;

if dim==2
    for i = 1:n
        for j = 2:m
            xj = xi[:,i] + pimages[:,j];
            xs = B2R*xj;
            if (wr[1,1] <= xs[1]) && (xs[1] <= wr[1,3]) &&  (wr[2,1] <= xs[2]) && (xs[2] <= wr[2,3])
                k = k + 1;
                xi = [xi xj];
                append!(alist,i);
            end
        end
    end
else
    for i = 1:n
        for j = 2:m
            xj = xi[:,i] + pimages[:,j];
            xs = B2R*xj;
            if (wr[1,1] <= xs[1]) && (xs[1] <= wr[1,7]) &&  (wr[2,1] <= xs[2]) && (xs[2] <= wr[2,7]) &&  (wr[3,1] <= xs[3]) && (xs[3] <= wr[3,7])
                k = k + 1;
                xi = [xi xj];
                append!(alist,i);
            end
        end
    end
end

return xi, alist, n, k

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
#Import .tsv taxbyseq file, create unique sequence identifier, create fasta file for blast
#run blast from shell
#import results, screen and parse, write out fasta file of Syn sequences for next steps

#Kristen R. Hunter-Cevera, Marine Biological Laboratory, Spring 2020

using DataFrames, CSV
using BioSequences
using Gadfly, Compose
##

# ind2 = findall(map(x -> occursin(r"Synecho",x), taxbyseq_df2[:,:Taxonomy])) #make a Syn only database
# taxbyseq_df2=taxbyseq_df2[ind2,:]

#import taxbyseq file:
mvco_df = DataFrame(CSV.File("/home/kristen/Documents/V6V8_analysis/data/vamps_downloads/taxbyseq_tables/mvco_timeseries_taxbyseq.tsv"; delim = '\t',header=2))
mvco_df.uniq_ID=collect(1:size(mvco_df,1)) #just use numbers for now?

#easier handling:
for n in names(mvco_df[!,Between(Symbol("MVCO_2010_2018_timeseries--MV254_run2"),Symbol("MVCO_2010_2018_timeseries--MV400Ster_run7"))])
     rename!(mvco_df, n => replace(string(n),"MVCO_2010_2018_timeseries--" => ""))
end

## setup for blast check:
fastafile="/home/kristen/Documents/V6V8_analysis/analysis_products/mvco_timeseries_uniq.fasta"

w = open(FASTA.Writer, fastafile)
for q=1:size(mvco_df,1)
    rec = FASTA.Record(string(mvco_df[q,:uniq_ID]),mvco_df[q,:Sequence])
    write(w, rec)
end
close(w)

## Run blast:

outfile="/home/kristen/Documents/V6V8_analysis/analysis_products/blast/mvco_timerseries_blast"
run(`bash /home/kristen/Documents/V6V8_analysis/scripts/taxa_processing_scripts/run_blast.sh $fastafile $outfile`)

## once have run blast, reimport results:

#convert output file from blastn search to DataFrame
blast_df = DataFrame(CSV.File("/home/kristen/Documents/V6V8_analysis/analysis_products/blast/mvco_timeseries_blast"; delim = '\t', header = false))

#put headers on output file
blast_col = ["uniq_ID", "sseqid" ,"pident" ,"length", "mismatch","gapopen", "qstart", "qend", "sstart", "send", "evalue", "bitscore"];
names!(blast_df, Symbol.(blast_col))

## merge dataframes:

test_df=join(mvco_df,blast_df,on = :uniq_ID, kind=:left)

test_df.query_length=abs.(test_df[!,:qstart] .- test_df[!,:qend])
test_df.subject_length=abs.(test_df[!,:sstart] .- test_df[!,:send])

## Sometimes blast will come back with two answers for a sequence match -
#not sure if that's a threading issue or something being mismanaged
#this is despite the the fact that only one answer was requested from the program!

#Upon merging, if the dataframe is presented with two options - I think it jsut duplicates the line to add both matches

t=convert(Array,test_df[!,:uniq_ID]) #merged dataframe is longer than mvco_df
g=setdiff(1:size(mvco_df,1),indexin(unique(t), t)) #find the indices of duplicate unqiue IDs
gg=unique(test_df[g,:uniq_ID])

global rec=zeros(size(gg,1),1)

for q=1:length(gg)

    ii=findall(test_df[!,:uniq_ID] .== gg[q]) #find matching indexes
    println(size(ii))
    k=argmin(test_df[ii,:bitscore])
    deleterows!(test_df,ii[k]) #remove the result that had the lower bit score, but keep record of lines removed!
    rec[q]=ii[k]

end

#size of rec should match size of gg!

## great, now check how the syn sequences look:

#find all the Syn based on GAST matches:
synind=findall(occursin.("Synecho",test_df[!,:Taxonomy]))

## prelim plot:

Gadfly.set_default_plot_size(20cm,20cm)
p1=plot(x=test_df[synind,:pident], Geom.histogram,Guide.title("Percent ID")) #most should be > 95
p2=plot(x=test_df[synind,:mismatch], Geom.histogram,Guide.title("Mismatch")) #most should be < 10
p3=plot(x=test_df[synind,:gapopen], Geom.histogram,Guide.title("Gaps Open")) #most should be < 3
p4=plot(x=test_df[synind,:bitscore], Geom.histogram,Guide.title("BitScore")) #should be above 800

gridstack([p1 p2; p3 p4])

##

synind=findall(occursin.(r"Synecho",test_df[!,:Taxonomy]))
noblast_hitsyn = findall(ismissing.(test_df[synind,:pident])) #any syn (classified by gast) that did not match to our blast database
#that should be empty or no hits - all gast syn basically had blast syn hits
ii=setdiff(1:length(synind),noblast_hitsyn)
synind=synind[ii]

#what do min values look like?
jj=findall(test_df[synind,:bitscore] .<= 700)
minimum(test_df[synind,:bitscore])

findall(test_df[synind,:pident] .<= 90)
minimum(test_df[synind,:pident])

findall(test_df[synind,:mismatch] .>= 10)
maximum(test_df[synind,:mismatch])

findall(test_df[synind,:gapopen] .>= 3)
maximum(test_df[synind,:gapopen])

# Some surprisingly low values, but I think for this, we will keep these sequences in and let oligotyping sort it out!

## 🌵 Now, do any sequences need to be added?
# I.e. are there Syn sequences not tagged by GAST but captured by blast?

noblast = findall(ismissing.(test_df[!,:pident])) #no blast results to syn
synblast = setdiff(1:size(test_df,1),noblast) #total blast results that had match to syn
nonsyn = setdiff(synblast,synind) #remove gast matches to syn, what came back as syn?

plot(x=test_df[nonsyn,:bitscore], Geom.histogram)

kk=findall(test_df[nonsyn,:bitscore] .>= 800)
jj=findall(test_df[nonsyn,:bitscore] .>= 700)
unique(test_df[nonsyn[kk],:Taxonomy])
#minimum(test_df[synind,:bitscore])

sum([sum(row) for row = eachrow(test_df[nonsyn[jj],Between(:MV254_run2,:MV400Ster_run7)])])
#So, only a few hundred sequences come back as not-Syn by GAST, but a very good Syn Blast match to our database

#check in general these categories:
ii=findall(test_df[!,:Taxonomy] .== "Bacteria;Cyanobacteria;empty_class;Chroococcales;empty_family")
ii=findall(test_df[!,:Taxonomy] .== "Bacteria;Cyanobacteria;empty_class")
ii=findall(occursin.("Prochloro",test_df[!,:Taxonomy]))
sum([sum(row) for row = eachrow(test_df[ii,Between(:MV254_run2,:MV400Ster_run7)])]) #low pro count! good!

## with the exception of matches to Pro, we're gonna inlcude these extra blast hits :D
#adjusted synind:

ii=findall(occursin.("Prochloro",test_df[!,:Taxonomy]))
ind=setdiff(nonsyn[jj],ii) #jj here refers to bitscore threshold matches
unique(test_df[ind,:Taxonomy])

newind=union(synind,ind)

## Write sequences for Syn to fasta file, in oligotype format

df=test_df[newind,:] #not a deep copy

w = open(FASTA.Writer, "/home/kristen/Documents/V6V8_analysis/analysis_products/syn_sequences/mvco_syn_timeseries_for_oligotyping.fasta")

##
global syncount=DataFrame(event=Union{String,Missing}[],count=Union{Float64,Missing}[])
for n in names(df[!,Between(:MV254_run2,:MV400Ster_run7)])  #names(df[!,Between(:SynA,:SynD)]) # #names(df[!,Between(:MV263_rd2,:MV375)])

    # ## We need to split and clean up dataframe a little bit as there are repeat samples, OOI cruise samples and mock samples...
    # if occursin(r"--\d{3}",string(n)) #mvco sample

        ind=findall(map(x -> x != 0, df[:,n])) #returns the indexes that are not zero
        #df[ind,[n, 68]]
        @show n
        global test=sum(df[ind,n])
        #m=match(r"(?<sample>(M[1,2]{1}|Syn[A,B,C,D]{1}))",string(n))

        m=match(r"(?<sample>MV\d{3}[A,B]{0,1}\w*)",string(n)) #for MVCO time series
        #@show m[:smaple]

        global count=0 #important for oligotyping file - adds number to each sequence
        for j=1:length(ind) #different Syn seqs in this sample
            for i=1:df[ind[j],n] #for each syn seq, duplicate to create fasta file for oligotyping
                #@show
                global count += 1
                rec = FASTA.Record(string("Sample-",m[:sample],"_Read",count),df[ind[j],:Sequence])
                write(w, rec)
            end
        end

        push!(syncount,[m[:sample],count])
        @show test
        @show count
        if count != test
            print("Nooooo - sum does not match sequences printed?")
        end
    #end
end

##
close(w)

#and can check sums here - compare with grep commands on just made fasta file:
sum([sum(row) for row = eachrow(test_df[newind,Between(:MV254_run2,:MV400Ster_run7)])]) #Excellent, numbers should match!

#and if need total sequence counts per sample:
syncount.totals=[sum(cols) for cols = eachcol(test_df[!,Between(:MV254_run2,:MV400Ster_run7)])]
syncount |> CSV.write("/home/kristen/Documents/V6V8_analysis/analysis_products/syn_sequences/mvco_total_and_syn_counts.csv",delim=",")

## OLDER, PREVIOUS CODE:

# m=match(r"(?<event>(\d{3}[A,B]{0,1}_tr|\d{3}[A,B]{0,1}))","MV263A_tr")

#remodel this database to screen out non-mvco timeseries samples: NO LONGER NEEDED SINCE REUPLOAD TO VAMPS!
# redo5=("263","325","336","337","339","357","358")
# redo7=("308","309","310","338")
#
# for name in names(taxbyseq_df)
#         if occursin(r"(--C\d{1,2}N\d{1,2})|(MVBUCKET)|(MVM\d{1})|(Dock)|(SynCommunity)",string(name)) #for now, don't include these samples in analysis
#                 println(name)
#                 delete!(taxbyseq_df,name)
#         elseif occursin(r"MVCO_V6V8_Feb2016--",string(name)) #remove datasetname
#                 rename!(taxbyseq_df, name => replace(string(name),"MVCO_V6V8_Feb2016--" => ""))
#         elseif occursin(r"MVCO_ALLREADS_20190108--",string(name)) #remove dataset name and add MV
#                 m=match(r"MVCO_ALLREADS_20190108--(?<event>\d{3})",string(name))
#                 if m[:event] in redo5
#                         rename!(taxbyseq_df, name => string("MV",m[:event],"_rd2"))
#                 else
#                         rename!(taxbyseq_df, name => string("MV",m[:event]))
#                 end
#         elseif occursin(r"MVCO_TestRun2_August2018--",string(name)) #remove dataset name and add MV
#                 m=match(r"MVCO_TestRun2_August2018--(?<event>MV363A|MV\d{3})",string(name))
#                 rename!(taxbyseq_df, name => string(m[:event],"_tr2"))
#         elseif occursin(r"MVCO_timeseries_run7--",string(name)) #remove dataset name and add MV
#                 m=match(r"MVCO_timeseries_run7--(?<event>\d{3}(Milli|Ster){0,1})",string(name))
#                 if m[:event] in redo7
#                     rename!(taxbyseq_df, name => string("MV",m[:event],"_rd2"))
#                 else
#                     rename!(taxbyseq_df, name => string("MV",m[:event]))
#                 end
#         end
# end

## alter database to only keep mock community samples:
# taxbyseq_df=taxbyseq_df2
# for name in names(taxbyseq_df)
#         if occursin(r"(SynCommunity|MVM\d{1})|Sequence",string(name)) #for now, don't include these samples in analysis
#             println(name)
#             if occursin(r"MVCO_V6V8_Feb2016--",string(name)) #remove datasetname
#                 rename!(taxbyseq_df, name => replace(string(name),"MVCO_V6V8_Feb2016--" => ""))
#             elseif occursin(r"MVCO_timeseries_run7--",string(name)) #remove dataset name
#                 rename!(taxbyseq_df, name => replace(string(name),"MVCO_timeseries_run7--" => ""))
#             end
#         else
#             delete!(taxbyseq_df,name)
#         end
# end

max_stars_repo_name: hsosik/NES-LTER

max_issues_repo_name: hsosik/NES-LTER

max_forks_repo_name: hsosik/NES-LTER


text:
using Distributions
using DataFrames

# Randomly generate covariates, response and errors for a regression problem.
#
#  (A, b, errors) = gen_gaussian_regression_problem(N, P, coefs, probdist, snrDB)
#
# Outputs:
#  A - N*P matrix where each row has P covariates
#  b - N*1 response vector (actual response + errors)
#  errors - N*1 error vector
#
# Inputs:
#  coefs - coefficients for the regression problem
#  probdist - probability ditribution for generating the covariates
#  snrDB - Signal-noise ratio in decibels
#
function gen_regression_problem(N, P, coefs, probdist = Normal(0.0, 1.0), snrDB = 20.0)
  A = zeros(Float64, N, P)
  for i in 1:N
    A[i,:] = rand(probdist, P)
  end
  A

  y = A * coefs

  noisestd = sqrt( noisevar(snrDB, var(y)) )
  noisedist = Normal(0.0, noisestd)
  errors = rand(noisedist, N)

  b = y .+ errors

  return (A, b, errors)
end

noisevar(snrDB, signalvar) = signalvar / (10 ^ (snrDB / 10))

# coefs = rand(5)
# A, b, e = gen_regression_problem(20, 5, coefs)

### Evaluating quality of regression coefficients ###

function confusionmatrix(expected, actual)
  @assert length(expected) == length(actual)
  cm = zeros(2, 2)
  for i in 1:length(expected)
    e = abs(expected[i]) > 0.0 ? 1 : 2
    a = abs(actual[i])   > 0.0 ? 1 : 2
    cm[a, e] += 1.0
  end
  return cm
end

tp(cm) = cm[1,1]
fp(cm) = cm[1,2]
fn(cm) = cm[2,1]
tn(cm) = cm[2,2]
sensitivity(cm) = cm[1,1] / (cm[1,1] + cm[2,1])
accuracy(cm) = (tp(cm) + tn(cm)) / (sum(cm))
mcc(cm) = (tp(cm) * tn(cm) - fp(cm) * fn(cm)) / sqrt( (tp(cm) + fp(cm)) * (tp(cm) + fn(cm)) * (tn(cm) + fp(cm)) * (tn(cm) + fn(cm)))

numactive(coefs) = sum(abs(coefs) .> 0.0)
positivecoefs(p) = 1.0 + rand(p)
normalcoefs(p) = randn(p)

gencoefs(p, k, coefgen = positivecoefs) = shuffle([coefgen(k), zeros(p-k)])

# regressors is a dict from a regressor name to functions from (matrix, vector) to vector of coefficients.
function compare_regressors(regressors, coefgenerator = positivecoefs;
  Ns = [10, 50, 100],
  Pratios = [0.50, 1.0, 2.0],
  Ks = [1, 5, (p) -> sqrt(p)],
  SNRs = [10, 50],
  NumReps = 1,
  CovarProbDists = [Normal(0.0, 1.0)]
  )

  startime = time()

  problemnum = 0

  # For each experiment with a regressor we save its Accuracy, MCC, MSE, CoeffNorm, NumSelectedCoefs
  df = DataFrame()

  for rep in 1:NumReps
    for N in Ns
      for Pratio in Pratios
        P = ifloor(Pratio*N)
        ks = map(Ks) do K
          if typeof(K) <: Function
            K(P)
          else
            K
          end
        end

        for k in ks
          k = ifloor(k)
          # Skip if k is not smaller than P
          if k >= P
            next
          end

          for SNR in SNRs
            for covarprobdist in CovarProbDists

              activecoefs = coefgenerator(k)
              coefs = shuffle([activecoefs, zeros(P-k)])
              A, b, errors = gen_regression_problem(N, P, coefs, covarprobdist, SNR)
              problemnum += 1

              for (name, reg) in regressors

                println("$name: N = $N, P = $P, k = $k, SNR = $SNR")
                tic()
                  betas = reg(A, b)
                t = toq()
                msev = mse(b, A * betas)
                println("  time taken: $t, mse = $msev")

                cm = confusionmatrix(coefs, betas)

                res = DataFrame(
                  # Problem description
                  N = N, P = P, k = k, SNRdb = SNR, CovariateGen = string(covarprobdist), PId = problemnum,

                  Time = t,
                  Accuracy = accuracy(cm),
                  MCC = mcc(cm),
                  nSelected = numactive(betas),
                  AvgAbsNorm = norm(coefs .- betas, 1)/P,
                  Regressor = name,
                  MSE = msev

                )

                df = rbind(df, res)

              end
            end
          end
        end
      end
    end
  end

  filename = strftime("compare_regressors_%Y%m%d_%H%M%S.csv", startime)
  writetable(filename, df)

  return df
end


text:
using ASP

facts("as_nnls") do

    context("Wide and sparse NNLS problem.") do

        m, n = 600, 2560              # No. of rows, cols
        k = 20                        # No. of nonzeros
        p = randperm(n)[1:k]          # Position of nonzeros in x
        x = zeros(n)                  # Generate sparse nonnegative solution
        x[p] = max(randn(k),0)
        A = randn(m,n)                # Gaussian m-by-n ensemble
        b = A*x                       # Compute the RHS vector
        x, r, inform = as_nnls(A,b)   # Solve the nonnegative LS problem

        g = -(A'*r)
        @fact vecnorm( min(x,g) ) --> less_than(1e-8)

    end

    context("Overdetermined NNLS problem.") do

        m, n = 100, 10
        x = max(randn(n), 0)
        A = randn(m,n)
        b = A*x
        x, r, inform = as_nnls(A, b)

        g = -(A'*r)
        @fact vecnorm( min(x,g) ) --> less_than(1e-8)
    end

end


text:
include("Plant.jl")
include("NMPC.jl")
include("Parameters.jl")
using Plots

## Defining Simulation Parameters

## Options for plotting, saving the figure and Eporting the data
#Plotting    = false;
#Fig_Save    = false;
#Export_Vars = false;

## Simulation of the plant
function Simulate_Plant(Data_MPC, Data_Plant)

    #sim_par = Data_MPC;
        T0_sim  = Data_MPC[:T0_sim]
        Tf_sim  = Data_MPC[:Tf_sim]
        dt_sim  = Data_MPC[:dt_sim]
        Tf_MPC  = Data_MPC[:Tf_MPC]
        dt_MPC  = Data_MPC[:dt_MPC]
        NCP     = Data_MPC[:NCP]
        NFE_sim = Data_MPC[:NFE_sim]
        NFE_MPC = Data_MPC[:NFE_MPC]
        uk      = Data_MPC[:uk]
        xₛₚ     = Data_MPC[:xₛₚ]
        xk      = Data_MPC[:xk]
        U_Plot = [];
        X1_Plot = [];
        X2_Plot = [];

        # push the starting point into the final vectors
        push!(U_Plot, uk[1])
        push!(X1_Plot, xk[1])
        push!(X2_Plot, xk[2])
        global dt_sim, uk, xk

    println("Start Simulation!")
    for k in 1:NFE_sim
        global xk, uk, dt_sim
        xₛₚ_MPC = xₛₚ[k:k - 1 + NFE_MPC]
        uk = Solve_MPC(xk, uk, xₛₚ_MPC,NFE_MPC, NCP, Data_Plant)
        xk = Plant(xk, uk, dt_sim, Data_Plant);
        push!(U_Plot, uk[1])
        push!(X1_Plot, xk[1])
        push!(X2_Plot, xk[2])
    end
    println("The objective function is satisfied successfully!")
    myDict = Dict(
        :X1_Plot   => X1_Plot,
        :X2_Plot     => X2_Plot,
        :U_Plot     => U_Plot
    )
    return myDict
end

## Plotting with the option to show/not show the figure. Also different backends can be choose!
function Plot_Plant(MPC_Results, Data_MPC)
    T0_sim  = Data_MPC[:T0_sim];
    dt_sim  = Data_MPC[:dt_sim];
    Tf_sim  = Data_MPC[:Tf_sim];
    NFE_sim = Data_MPC[:NFE_sim];
    xₛₚ     = Data_MPC[:xₛₚ];
    X1_Plot = MPC_Results[:X1_Plot];
    X2_Plot = MPC_Results[:X2_Plot];
    U_Plot  = MPC_Results[:U_Plot];
    t_plot = collect(T0_sim:dt_sim:Tf_sim)
    p11 = plot( t_plot,[X1_Plot[nIter]  for nIter in 1:NFE_sim+1],label="x1")
    p11 = plot!(t_plot,[X2_Plot[nIter]  for nIter in 1:NFE_sim+1],label="x2")
    p11 = plot!(t_plot, xₛₚ[1:NFE_sim + 1],                       label="y_sp", linetype=:steppost, linestyle=:dash)
    p12 = plot(t_plot[1:end], [U_Plot[nIter]  for nIter in 1:NFE_sim+1],    label="u",              linetype=:steppost)
    fig1 = plot(p11, p12, layout=(2, 1));
    plotlyjs()
    #gr()
    println("Start Plotting!")
    display(fig1)
    println("Finish Plotting!")
    return fig1
end

##Checking the Option for saving
function Save_Plant(fig1, Address)
    savefig(fig1, Address)
end
## Exporting Data

function Export_Plant(MPC_Results, Data_MPC, Address)
    X1_Plot = MPC_Results[:X1_Plot];
    X2_Plot = MPC_Results[:X2_Plot];
    U_Plot  = MPC_Results[:U_Plot];
    xₛₚ     = Data_MPC[:xₛₚ];
    NFE_sim = Data_MPC[:NFE_sim];
    T0_sim  = Data_MPC[:T0_sim];
    dt_sim  = Data_MPC[:dt_sim];
    Tf_sim  = Data_MPC[:Tf_sim];
    t_plot  = collect(T0_sim:dt_sim:Tf_sim);
    println("Start Exporting Data!")
    #using CSV
    #using DataFrames
    df = DataFrame(Name = ["x1","x2","U","SetPoint","Time"],
                   Value = [X1_Plot, X2_Plot, U_Plot, xₛₚ[1:NFE_sim+1], t_plot]
                   )
    CSV.write(Address, df)
    println("Finish Exporting Data!")
end


text:
######################################################################
# fovm.jl
#
# like FOV, but allow for more regions?
######################################################################

type FOVm <: Sensor
	region_probs::Vector{NTuple{2,Float64}}

	# for discretized measurements
	num_bins::Int
	bin_range::UnitRange{Int64}

	FOVm(region_probs::Vector{NTuple{2,Float64}}) = new(region_probs, 2, 0:1)

    function FOVm(cone_width::Real, alpha::Float64)
        a1 = cone_width / 2.0
        a2 = 180.0 - a1
        a3 = 180.0

        region_probs = [(a1,1.0-alpha), (a2, 0.5), (a3, alpha)]

        return new(region_probs, 2, 0:1)
    end
end

# 1 means it is in the field of view
# 0 means it is not
function observe(m::SearchDomain, s::FOVm, p::Pose)

	# ensure bearing is reflected across
	rel_bearing = fit_180(p[3] - true_bearing(p, m.theta))
	if rel_bearing < 0.0
		rel_bearing = -1.0 * rel_bearing
	end

	prob_in_view = 0.0
	n = length(s.region_probs)
	for i = 1:n
		temp_angle, temp_prob = s.region_probs[i]
		if rel_bearing  <= temp_angle
			prob_in_view = temp_prob
			break
		end
	end

	o = (rand() < prob_in_view) ? 1.0 : 0.0
	return o
end

# o technically has to be a float, but we only have two obs:
#  1 = in field of view
#  0 = not in field of view
function O(s::FOVm, theta::LocTuple, p::Pose, o::Float64)

	# determine relative bearing and fix it in 0 to 180
	rel_bearing = fit_180(p[3] - true_bearing(p, theta))
	if rel_bearing < 0.0
		rel_bearing = -1.0 * rel_bearing
	end

	prob_in_view = 0.0
	n = length(s.region_probs)
	for i = 1:n
		temp_angle, temp_prob = s.region_probs[i]
		if rel_bearing  <= temp_angle
			prob_in_view = temp_prob
			break
		end
	end

	ret_val = (o == 1.0) ? prob_in_view : (1.0 - prob_in_view)
	return ret_val
end

# I don't really get why this is here...
#function O(x::Vehicle,s::FOV, xp::Pose, theta::LocTuple, o::ObsBin, df::DF)
function O(s::FOVm, theta::LocTuple, xp::Pose, o::ObsBin)
	# determine relative bearing and fix it in 0 to 180
	rel_bearing = fit_180(xp[3] - true_bearing(xp, theta))
	if rel_bearing < 0.0
		rel_bearing = -1.0 * rel_bearing
	end

	prob_in_view = 0.0
	n = length(s.region_probs)
	for i = 1:n
		temp_angle, temp_prob = s.region_probs[i]
		if rel_bearing  <= temp_angle
			prob_in_view = temp_prob
			break
		end
	end

	ret_val = (o == 1.0) ? prob_in_view : (1.0 - prob_in_view)
	return ret_val
end

obs2bin(o::Float64, s::FOVm) = round(Int, o)


text:
MLJBase.is_wrapper(::Type{<:EitherIteratedModel}) = true
MLJBase.caches_data_by_default(::Type{<:EitherIteratedModel}) = false
MLJBase.load_path(::Type{<:DeterministicIteratedModel}) =
    "MLJIteration.DeterministicIteratedModel"
MLJBase.load_path(::Type{<:ProbabilisticIteratedModel}) =
    "MLJIteration.ProbabilisticIteratedModel"
MLJBase.package_name(::Type{<:EitherIteratedModel}) = "MLJIteration"
MLJBase.package_uuid(::Type{<:EitherIteratedModel}) =
    "614be32b-d00c-4edb-bd02-1eb411ab5e55"
MLJBase.package_url(::Type{<:EitherIteratedModel}) =
    "https://github.com/JuliaAI/MLJIteration.jl"
MLJBase.package_license(::Type{<:EitherIteratedModel}) = "MIT"

# inherited traits:
for T in [:DeterministicIteratedModel, :ProbabilisticIteratedModel]
    for trait in [:supports_weights,
                  :supports_class_weights,
                  :is_pure_julia,
                  :input_scitype,
                  :output_scitype,
                  :target_scitype]
        quote
            # try to get trait at level of types ("failure" here just
            # means falling back to `Unknown`):
            MLJBase.$trait(::Type{<:$T{M}}) where M = MLJBase.$trait(M)

            # needed because traits are not always deducable from
            # the type (eg, `target_scitype` and `Pipeline` models):
            MLJBase.$trait(imodel::EitherIteratedModel) = $trait(imodel.model)
        end |> eval
    end
end


text:
#= PLEASE NOTE THAT ALL NCURSES FUNCS ARE ZERO-BASED
=#
libncurses = false
libpanel = false

function initscr()
    global libncurses, libpanel
    if  libncurses == false
        try
            libncurses = Libdl.dlopen("libncurses")
        catch
            libncurses = Libdl.dlopen("libncursesw")
        end
        try
            libpanel = Libdl.dlopen("libpanel")
        catch
            libpanel = Libdl.dlopen("libpanelw")
        end
    end
    ccall( Libdl.dlsym( libncurses, :initscr), Ptr{Void}, () )
end

function endwin()
    ccall( Libdl.dlsym( libncurses, :endwin), Int, () )
end

function isendwin()
    ccall( Libdl.dlsym( libncurses, :isendwin), Bool, () )
end

function newwin( lines::Int, cols::Int, origy::Int, origx::Int )
    ccall( Libdl.dlsym( libncurses, :newwin ), Ptr{Void}, ( Int, Int, Int, Int ),
        lines, cols, origy, origx )
end

function subwin( win::Ptr{Void}, lines::Int, cols::Int, origy::Int, origx::Int )
    ccall( Libdl.dlsym( libncurses, :subwin ), Ptr{Void}, ( Ptr{Void}, Int, Int, Int, Int ),
        win, lines, cols, origy, origx )
end

function derwin( win, lines::Int, cols::Int, origy::Int, origx::Int )
    ccall( Libdl.dlsym( libncurses, :derwin ), Ptr{Void}, ( Ptr{Void}, Int, Int, Int, Int ),
        win, lines, cols, origy, origx )
end

function delwin( win::Ptr{Void} )
    ccall( Libdl.dlsym( libncurses, :delwin ), Void, ( Ptr{Void}, ), win )
end

function mvwaddch( win::Ptr{Void}, y::Int, x::Int, c )
    ccall( Libdl.dlsym( libncurses, :mvwaddch), Void,
        ( Ptr{Void}, Int, Int, Int ), win, y, x, c )
end

function mvwadd_wch( win::Ptr{Void}, y::Int, x::Int, c )
    ccall( Libdl.dlsym( libncurses, :mvwadd_wch), Void,
        ( Ptr{Void}, Int, Int, Int ), win, y, x, Int64(c) )
end

function mvwaddch( w::TwWindow, y::Int, x::Int, c )
    if objtype( w.parent.value ) == :List && typeof( w.parent.value.window ) != TwWindow
        # terminal layer. use its pad
        mvwaddch( w.parent.value.data.pad, y+w.yloc, x+w.xloc, c )
    else
        mvwaddch( w.parent.value.window, y+w.yloc, x+w.xloc, c )
    end
end

function mvwprintw{T<:AbstractString}( win::Ptr{Void}, row::Int, height::Int, fmt::ASCIIString, str::T )
    ccall( Libdl.dlsym( libncurses, :mvwprintw), Void,
        ( Ptr{Void}, Int, Int, Cstring, Cstring ),
        win, row, height, fmt, str )
end

# note that it could in turn call another TwWindow...
function mvwprintw{T<:AbstractString}( w::TwWindow, y::Int, x::Int, fmt::ASCIIString, s::T )
    if objtype( w.parent.value ) == :List && typeof( w.parent.value.window ) != TwWindow
        # terminal layer. use its pad
        mvwprintw( w.parent.value.data.pad, y+w.yloc, x+w.xloc, fmt, s )
    else
        mvwprintw( w.parent.value.window, y+w.yloc, x+w.xloc, fmt, s )
    end
end

function wmove( win::Ptr{Void}, y::Int, x::Int )
    ccall( Libdl.dlsym( libncurses, :wmove), Int, ( Ptr{Void}, Int, Int ), win, y, x )
end

function wrefresh( win::Ptr{Void} )
    ccall( Libdl.dlsym( libncurses, :wrefresh ), Void, ( Ptr{Void}, ), win )
end

function touchwin( win::Ptr{Void} )
    ccall( Libdl.dlsym( libncurses, :touchwin ), Void, ( Ptr{Void}, ), win )
end

function refresh()
    ccall( Libdl.dlsym( libncurses, :refresh ), Void, ( ) )
end

function erase()
    ccall( Libdl.dlsym( libncurses, :erase ), Void, () )
end

function werase( win::Ptr{Void} )
    ccall( Libdl.dlsym( libncurses, :werase ), Void, (Ptr{Void},), win )
end

function werase( win::TwWindow, y::Int=win.yloc, x::Int=win.xloc, h::Int=win.height, w::Int=win.width )
    parwin = win.parent.value.window
    if typeof( parwin ) <: Ptr
        if objtype( win.parent.value ) == :List
            tmpw = win.parent.value.data.pad
        else
            tmpw = parwin
        end
        for r = y:y+h-1
            for c = x:x+w-1
                mvwaddch( tmpw, r, c, ' ' )
            end
        end
    else
        werase( parwin, parwin.yloc+y, parwin.xloc+x, h, w )
    end
end

function wclear( win::Ptr{Void} )
    ccall( Libdl.dlsym( libncurses, :wclear ), Void, (Ptr{Void},), win )
end

function box( win::Ptr{Void}, vchr, hchr )
    ccall( Libdl.dlsym( libncurses, :box ), Void, (Ptr{Void}, Char, Char), win, vchr, hchr )
end

function box( win::TwWindow, vchr::Integer, hchr::Integer, y::Int=win.yloc, x::Int=win.xloc, h::Int=win.height, w::Int=win.width )
    parwin = win.parent.value.window
    if typeof( parwin ) <: Ptr
        if objtype( win.parent.value ) == :List
            tmpw = win.parent.value.data.pad
        else
            tmpw = parwin
        end
        #draw the box myself
        # 4 corners
        mvwaddch( tmpw, y,x, get_acs_val( 'l' ) )
        mvwaddch( tmpw, y+h-1,x, get_acs_val( 'm' ) )
        mvwaddch( tmpw, y,x+w-1, get_acs_val( 'k' ) )
        mvwaddch( tmpw, y+h-1,x+w-1, get_acs_val( 'j' ) )

        if vchr==0
            vchr=get_acs_val('x')
        end
        if hchr==0
            hchr=get_acs_val('q')
        end
        for r = y+1:y+h-2
            mvwaddch( tmpw, r, x, vchr )
            mvwaddch( tmpw, r, x+w-1, vchr )
        end
        for c = x+1:x+w-2
            mvwaddch( tmpw, y, c, hchr )
            mvwaddch( tmpw, y+h-1, c, hchr )
        end
    else
        box( parwin, vchr, hchr, parwin.yloc+win.yloc, parwin.xloc+win.xloc, h, w )
    end
end

function wgetch( win::Ptr{Void} )
    ccall( Libdl.dlsym( libncurses, :wgetch ), UInt32, (Ptr{Void},), win )
end

function keypad( win::Ptr{Void}, bf )
    ccall( Libdl.dlsym( libncurses, :keypad ), Int, (Ptr{Void}, Bool), win, bf )
end

function cbreak()
    ccall( Libdl.dlsym( libncurses, :cbreak), Void, ( ) )
end

function nocbreak()
    ccall( Libdl.dlsym( libncurses, :nocbreak), Void, ( ) )
end

function echo()
    ccall( Libdl.dlsym( libncurses, :echo), Void, ( ) )
end

function noecho()
    ccall( Libdl.dlsym( libncurses, :noecho), Void, ( ) )
end

function nodelay( win, bf )
    ccall( Libdl.dlsym( libncurses, :nodelay ), Int, (Ptr{Void}, Bool), win, bf )
end

function raw()
    ccall( Libdl.dlsym( libncurses, :raw ), Int, () )
end

function noraw()
    ccall( Libdl.dlsym( libncurses, :noraw ), Int, () )
end

function notimeout( win::Ptr{Void}, bf )
    ccall( Libdl.dlsym( libncurses, :notimeout ), Int, (Ptr{Void}, Bool), win, bf )
end

function timeout( delay::Int )
    ccall( Libdl.dlsym( libncurses, :timeout ), Void, (Int,), delay )
end

function wtimeout( win::Ptr{Void}, delay::Int )
    ccall( Libdl.dlsym( libncurses, :wtimeout ), Void, (Ptr{Void}, Int), win, delay )
end

# not standard but convenient
function getwinmaxyx( win::Ptr{Void} )
    maxy = ccall( Libdl.dlsym(libncurses, :getmaxy), Int, ( Ptr{Void}, ), win )
    maxx = ccall( Libdl.dlsym(libncurses, :getmaxx), Int, ( Ptr{Void}, ), win )
    ( maxy, maxx )
end

function getwinbegyx( win::Ptr{Void} )
    maxy = ccall( Libdl.dlsym(libncurses, :getbegy), Int, ( Ptr{Void}, ), win )
    maxx = ccall( Libdl.dlsym(libncurses, :getbegx), Int, ( Ptr{Void}, ), win )
    ( maxy, maxx )
end

function mvwin( win::Ptr{Void}, y::Int, x::Int )
    ccall( Libdl.dlsym( libncurses, :mvwin), Int, ( Ptr{Void}, Int, Int ), win, y, x )
end

function copywin( s::Ptr{Void}, d::Ptr{Void}, sminrow::Int, smincol::Int, dminrow::Int, dmincol::Int, dmaxrow::Int, dmaxcol::Int )
    ccall( Libdl.dlsym( libncurses, :copywin), Int, ( Ptr{Void}, Ptr{Void}, Int, Int, Int, Int, Int, Int ),
       s, d, sminrow,smincol,dminrow, dmincol,dmaxrow,dmaxcol )
end

function beep()
    ccall( Libdl.dlsym( libncurses, :beep ), Void, () )
end

function flash()
    ccall( Libdl.dlsym( libncurses, :flash), Void, () )
end

function is_term_resized( lines::Int, cols::Int )
    ccall( Libdl.dlsym( libncurses, :is_term_resized ), Bool, (Int, Int), lines, cols )
end

function wresize( win::Ptr{Void}, lines::Int, cols::Int )
    ccall( Libdl.dlsym( libncurses, :wresize), Int, (Ptr{Void}, Int, Int), win, lines, cols )
end

function getcuryx( win::Ptr{Void} )
    cury = ccall( Libdl.dlsym(libncurses, :getcury), Cint, ( Ptr{Void}, ), win )
    curx = ccall( Libdl.dlsym(libncurses, :getcurx), Cint, ( Ptr{Void}, ), win )
    ( cury, curx )
end

function start_color()
    ccall( Libdl.dlsym( libncurses, :start_color), Void, () )
end

function init_pair( pair, f, b )
    ccall( Libdl.dlsym( libncurses, :init_pair ), Int, ( Int16, Int16, Int16 ), pair, f, b )
end

function init_color( color, r,g,b )
    ccall( Libdl.dlsym( libncurses, :init_color ), Int, ( Int16, Int16, Int16, Int16 ),
        color, r,g,b )
end

function has_colors()
    ccall( Libdl.dlsym( libncurses, :has_colors), Bool, () )
end

function wattroff( win::Ptr{Void}, attrs )
    ccall( Libdl.dlsym(libncurses, :wattroff), Int, ( Ptr{Void}, UInt32 ), win, attrs )
end

function wattroff( win::TwWindow, attrs )
    parwin = win.parent.value.window
    if objtype( win.parent.value ) == :List && typeof( parwin ) <: Ptr
        wattroff( win.parent.value.data.pad, attrs )
    else
        wattroff( win.parent.value.window, attrs )
    end
end

function wattron( win::Ptr{Void}, attrs )
    ccall( Libdl.dlsym(libncurses, :wattron), Int, ( Ptr{Void}, UInt32 ), win, attrs )
end

function wattron( win::TwWindow, attrs )
    parwin = win.parent.value.window
    if objtype( win.parent.value ) == :List && typeof( parwin ) <: Ptr
        wattron( win.parent.value.data.pad, attrs )
    else
        wattron( win.parent.value.window, attrs )
    end
end

function wattrset( win::Ptr{Void}, attrs )
    ccall( Libdl.dlsym(libncurses, :wattrset), Int, ( Ptr{Void}, UInt32 ), win, attrs )
end

function wbkgdset( win::Ptr{Void}, ch )
    ccall( Libdl.dlsym(libncurses, :wbkgdset ), Void, ( Ptr{Void}, UInt32 ), win, ch )
end

function wbkgd( win::Ptr{Void}, ch )
    ccall( Libdl.dlsym(libncurses, :wbkgd ), Void, ( Ptr{Void}, UInt32 ), win, ch )
end

function curs_set( vis )
    ccall( Libdl.dlsym(libncurses, :curs_set), Int, ( Int, ), vis )
end

function has_mouse()
    ccall( Libdl.dlsym(libncurses, :has_mouse), Bool, () )
end

function mousemask( mask )
    oldmm = Array( UInt64, 1 )
    resultmm = ccall( Libdl.dlsym( libncurses, :mousemask), UInt64, (UInt64, Ptr{UInt64}), mask, oldmm )
    ( resultmm, oldmm[1])
end

function mouseinterval( n::Int )
    ccall( Libdl.dlsym(libncurses, :mouseinterval), Int, (Int, ) , n)
end

#hack!
const mouseByteString = bytestring( Array( UInt8, 64 ) )
function getmouse()
    #=
    type Mouse_Event_t
        id::Int8 # short
        x::Int32 # int
        y::Int32 # int
        z::Int32 # int
        bstate::UInt64 # unsigned long
    end
    =#
    # the 5th byte is x, 9th byte is y
    # 19th byte is x08 if scroll up, 20th byte is x08 if scroll down
    # 17th byte is x02 if button 1 pressed
    # 17th byte is x01 if button 1 released
    # 17-18th is 0xfffd if mousewheel is pressed down
    ccall( Libdl.dlsym( libncurses, :getmouse), Int, (Ptr{UInt8}, ), mouseByteString )
    bs = mouseByteString
    x = @compat UInt8(bs[5])
    y = @compat UInt8(bs[9])
    state=:unknown
    if (@compat UInt(bs[17])) & 0x02 != 0
        state = :button1_pressed
    elseif (@compat UInt(bs[19])) & 0x08 != 0
        state = :scroll_up
    elseif (@compat UInt(bs[20])) & 0x08 != 0
        state = :scroll_down
    end
    ( state, x, y, bs )
end

function screen_to_relative( w::Ptr, y::Integer, x::Integer )
    begy, begx = getwinbegyx( w )
    ( y-begy, x-begx )
end

function screen_to_relative( w::TwWindow, y::Integer, x::Integer )
    # figure out the canvas coordinate, ultimate parent
    xloc = w.xloc
    yloc = w.yloc
    par = w.parent.value
    while( !( typeof( par.window ) <: Ptr ) )
        xloc += par.window.xloc
        yloc += par.window.yloc
        par = par.window.parent.value
    end
    (view_rel_y, view_rel_x) = screen_to_relative( par.window, y, x )

    ( par.data.canvaslocy + view_rel_y, par.data.canvaslocx + view_rel_x )
end

function baudrate()
    ccall( Libdl.dlsym( libncurses, :baudrate), Int, () )
end

function clearok( win::Ptr{Void}, bf )
    ccall( Libdl.dlsym( libncurses, :clearok), Int, ( Ptr{Void}, Bool ), win, bf )
end

function immedok( win::Ptr{Void}, bf )
    ccall( Libdl.dlsym( libncurses, :immedok), Int, ( Ptr{Void}, Bool ), win, bf )
end

function napms( ms::Int )
    ccall( Libdl.dlsym( libncurses, :napms), Int, (Int, ), ms )
end

#===== PANEL library ====#

function update_panels()
    ccall( Libdl.dlsym( libpanel, :update_panels), Void, () )
end

function doupdate()
    ccall( Libdl.dlsym( libpanel, :doupdate ), Void, () )
end

function new_panel( win::Ptr{Void} )
    ccall( Libdl.dlsym( libpanel, :new_panel ), Ptr{Void}, ( Ptr{Void}, ), win )
end

function top_panel( pan::Ptr{Void} )
    ccall( Libdl.dlsym( libpanel, :top_panel ), Int, ( Ptr{Void}, ), pan )
end

function bottom_panel( pan::Ptr{Void} )
    ccall( Libdl.dlsym( libpanel, :bottom_panel ), Int, ( Ptr{Void}, ), pan )
end

function move_panel( pan::Ptr{Void}, starty::Int, startx::Int )
    ccall( Libdl.dlsym( libpanel, :move_panel ), Ptr{Void}, ( Ptr{Void}, Int, Int ), pan, starty, startx )
end

function del_panel( panel::Ptr{Void} )
    ccall(Libdl.dlsym( libpanel, :del_panel ), Void, (Ptr{Void}, ), panel )
end

function show_panel( panel::Ptr{Void} )
    ccall(Libdl.dlsym( libpanel, :show_panel ), Int, (Ptr{Void}, ), panel )
end

function hide_panel( panel::Ptr{Void} )
    ccall(Libdl.dlsym( libpanel, :hide_panel ), Int, (Ptr{Void}, ), panel )
end

function panel_hidden( panel::Ptr{Void } )
    ccall(Libdl.dlsym( libpanel, :panel_hidden ), Int, (Ptr{Void}, ), panel ) != 0
end

function replace_panel( panel::Ptr{Void}, window::Ptr{Void} )
    ccall(Libdl.dlsym( libpanel, :replace_panel ), Void, (Ptr{Void}, Ptr{Void} ), panel, window )
end

function set_panel_userptr( p1::Ptr{Void}, p2::Ptr{Void} )
    ccall(Libdl.dlsym( libpanel, :set_panel_userptr), Void, (Ptr{Void}, Ptr{Void}), p1, p2 )
end

# ------ pad

function newpad( rows::Int, cols::Int )
    ccall(Libdl.dlsym( libncurses, :newpad ), Ptr{Void}, (Int,Int), rows, cols )
end

function subpad( orig::Ptr{Void}, rows::Int, cols::Int, beg_y::Int, beg_x::Int )
    ccall(Libdl.dlsym( libncurses, :subpad ), Ptr{Void}, (Ptr{Void}, Int,Int,Int,Int), orig, rows, cols, beg_y, beg_x )
end

function pnoutrefresh( pad::Ptr{Void}, pminrow::Int, pmincol::Int, sminrow::Int,smincol::Int,smaxrow::Int,smaxcol::Int )
    ccall(Libdl.dlsym( libncurses, :pnoutrefresh), Void, (Ptr{Void}, Int,Int,Int,Int,Int,Int),
        pad, pminrow,pmincol,sminrow,smincol,smaxrow,smaxcol )
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMQuad.jl/blob/master/LICENSE

using FEMQuad, Test
using FEMQuad: integrate_2d, get_rule

@testset "Gauss-Legendre quadratures in triangles" begin
    rules = (:GLTRI1, :GLTRI3, :GLTRI4, :GLTRI6, :GLTRI7, :GLTRI12)
    f(i) = x -> sum(sum(x)^j for j=0:i)
    r(i::Int) = get_rule(i, rules...)

    @test isapprox(integrate_2d(f(0), r(0)), 1/2)
    @test isapprox(integrate_2d(f(1), r(1)), 5/6)
    @test isapprox(integrate_2d(f(2), r(2)), 13/12)
    @test isapprox(integrate_2d(f(3), r(3)), 77/60)
    @test isapprox(integrate_2d(f(4), r(4)), 29/20)
    @test isapprox(integrate_2d(f(5), r(5)), 223/140)
    @test isapprox(integrate_2d(f(6), r(6)), 481/280)
    #@test isapprox(integrate_2d(f(7), r(7)), 4609/2520)
    #@test isapprox(integrate_2d(f(8), r(8)), 4861/2520)
    #@test isapprox(integrate_2d(f(9), r(9)), 55991/27720)
end

@testset "Gauss-Legendre quadratures in triangles B" begin
    rules = (:GLTRI1, :GLTRI3B, :GLTRI4B, :GLTRI6, :GLTRI7, :GLTRI12)
    f(i) = x -> sum(sum(x)^j for j=0:i)
    r(i::Int) = get_rule(i, rules...)

    @test isapprox(integrate_2d(f(2), r(2)), 13/12)
    @test isapprox(integrate_2d(f(3), r(3)), 77/60)
end


text:
#Comment out the next 5 lines if you already have these packages installed.
#using Pkg; #Using the Julia Pkg package....
#packages = ["HTTP","JSON","DataFrames","NCDatasets"] #Add packages to be downloaded here.
#for package in packages
#    Pkg.add(package)  #Download package.
#end

using HTTP, JSON, DataFrames, NCDatasets  #Bring the packages in so that we can use them.

#------CHANGE ME-----#
user = "OOI-API-USER-HERE"
token = "OOI-API-TOKEN-HERE"
directory = "C:/Users/Ian/Desktop/test"
#--------------------#

site = "CE01ISSP"
node = "SP001"
instrument = "09-CTDPFJ000"
method = "recovered_cspp"
stream = "ctdpf_j_cspp_instrument_recovered"
start_date = "2018-06-01"
stop_date = "2019-06-02"

url = ooi_create_url(site,node,instrument,method,stream,start_date = start_date, stop_date = stop_date)
response = ooi_submit_request(url,user,token)
dodsC = ooi_get_location(response)
fileServer = ooi_download_data(dodsC,directory)

max_stars_repo_name: lstoltz/ooi-data-explorations

max_issues_repo_name: lstoltz/ooi-data-explorations

max_forks_repo_name: lstoltz/ooi-data-explorations


text:
export build_SDP_Instance_from_sdpfiles

function build_SDP_Instance_from_sdpfiles(path::String)
  sdp_read = read_SDPPrimal(path)

  # println("VAR_TYPES size:     $(size(sdp_read.VAR_TYPES))")
  # println("BLOCKS size:        $(size(sdp_read.BLOCKS))")
  # println("LINEAR size:        $(size(sdp_read.LINEAR))")
  # println("CONST size:         $(size(sdp_read.CONST))")

  sdp_inst = SDP_Problem()

  # Build numerical descirption of problem, indexed on strings
  set_matrices!(sdp_inst, sdp_read)
  set_linear!(sdp_inst, sdp_read)
  set_const!(sdp_inst, sdp_read)

  # Build structural information and string to id maps
  set_vartypes!(sdp_inst, sdp_read)

  set_blocks!(sdp_inst)
  set_scalvars!(sdp_inst)

  set_constraints!(sdp_inst)

  return sdp_inst
end

function read_SDPPrimal(path::String)
  BLOCKS = readdlm(joinpath(path, "matrix.sdp"), String)
  if isfile(joinpath(path, "lin.sdp")) && length(matchall(r"\n", readstring(joinpath(path, "lin.sdp")))) > 7
    LINEAR = readdlm(joinpath(path, "lin.sdp"), String)
  else
    LINEAR = []
  end
  CONST = readdlm(joinpath(path, "const.sdp"), String)
  VAR_TYPES = readdlm(joinpath(path, "types.sdp"), String)

  SDP_Instance(VAR_TYPES,
               BLOCKS,
               LINEAR,
               CONST)
end

function set_matrices!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)
  for i=1:size(instance.BLOCKS, 1)
    ctr_name = (instance.BLOCKS[i, 1], instance.BLOCKS[i, 2], instance.BLOCKS[i, 3])
    (block_name, var1, var2, coeff) = instance.BLOCKS[i, 4:7]

    # Sort variables for triangular matrix storage
    var1, var2 = min(var1, var2), max(var1, var2)

    # if haskey(sdp.name_to_sdpblock, block_name)
    if !haskey(sdp.matrices, (ctr_name, block_name, var1, var2))
      sdp.matrices[(ctr_name, block_name, var1, var2)] = parse(coeff)
    else
      error("set_matrices!(): sdp.matrices already has key ($ctr_name, $block_name, $var1, $var2) with val $(sdp.matrices[(ctr_name, block_name, var1, var2)]), $(parse(coeff))")
    end

    # else
    #   error("set_matrices!(): Unhandled matrix var $block_name")
    # end
  end

  if debug
  end
end

function set_linear!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)
  if length(instance.LINEAR) != 0
    for i=1:size(instance.LINEAR, 1)
      ctr_name = (instance.LINEAR[i, 1], instance.LINEAR[i, 2], instance.LINEAR[i, 3])
      (var, coeff) = instance.LINEAR[i, 4:5]

      if !haskey(sdp.linear, (ctr_name, var))
        sdp.linear[(ctr_name, var)] = parse(coeff)
      else
        error("set_linear!(): sdp.linear already has key ($ctr_name, $var) with val $(sdp.linear[(ctr_name, var)]). New val is $(parse(coeff))")
      end
    end
  end

  if debug
  end
end

function set_const!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)
  for i=1:size(instance.CONST, 1)
    ctr_name = (instance.CONST[i, 1], instance.CONST[i, 2], instance.CONST[i, 3])
    coeff = instance.CONST[i, 4]

    @assert !haskey(sdp.cst_ctr, ctr_name)
    if coeff != 0
      sdp.cst_ctr[ctr_name] = parse(coeff)
    end
  end

  if debug
  end
end

"""
set_vartypes!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)

Input all matrix varibales structural information regarding name and kind in the appropriate attributes of `SDP_Problem`.
"""
function set_vartypes!(sdp::SDP_Problem, instance::SDP_Instance; debug=false)
  n_sdp = 0
  for i=1:size(instance.VAR_TYPES, 1)
    (block_name, block_type) = instance.VAR_TYPES[i, :]

    if block_type == "SDP"
      n_sdp += 1
      block = SDP_Block(n_sdp, block_name)
      sdp.name_to_sdpblock[block_name] = block
      sdp.id_to_sdpblock[n_sdp] = block

    else
      error("set_vartypes()!: Unknown blockvar type $block_type for variable $block_name")
    end
  end

  if debug
  end
end


text:
module CommunalHelperCassetteZipMover

using ..Ahorn, Maple
using Ahorn.CommunalHelper
using Ahorn.CommunalHelperEntityPresets: CustomCassetteBlockData

const entityData = appendkwargs(CustomCassetteBlockData, :(
    permanent::Bool=false,
    waiting::Bool=false,
    ticking::Bool=false,
    noReturn::Bool=false,
))
@mapdefdata Entity "CommunalHelper/CassetteZipMover" CassetteZipMover entityData

const ropeColors = Dict{Int,Ahorn.colorTupleType}(
    1 => (194, 116, 171, 255) ./ 255,
    2 => (227, 214, 148, 255) ./ 255,
    3 => (128, 224, 141, 255) ./ 255,
)

const defaultRopeColor = (110, 189, 245, 255) ./ 255

const placements = Ahorn.PlacementDict(
    "Cassette Zip Mover ($index - $color) (Communal Helper)" => Ahorn.EntityPlacement(
        CassetteZipMover,
        "rectangle",
        Dict{String,Any}(
            "index" => index,
        ),
        function (entity)
            entity.data["nodes"] = [(
                Int(entity.data["x"]) + Int(entity.data["width"]) + 8,
                Int(entity.data["y"]),
            )]
        end,
    ) for (color, index) in cassetteColorNames
)

Ahorn.editingOptions(entity::CassetteZipMover) = Dict{String,Any}(
    "index" => cassetteColorNames,
)

Ahorn.nodeLimits(entity::CassetteZipMover) = 1, -1

Ahorn.minimumSize(entity::CassetteZipMover) = 16, 16
Ahorn.resizable(entity::CassetteZipMover) = true, true

function Ahorn.selection(entity::CassetteZipMover)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    res = [Ahorn.Rectangle(x, y, width, height)]

    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)

        push!(res, Ahorn.Rectangle(nx + floor(Int, width / 2) - 5, ny + floor(Int, height / 2) - 5, 10, 10))
    end

    return res
end

const cog = "objects/CommunalHelper/cassetteZipMover/cog"
const crossSprite = "objects/CommunalHelper/cassetteMoveBlock/x"

function renderCassetteZipMover(ctx::Ahorn.Cairo.CairoContext, entity::CassetteZipMover)
    x, y = Ahorn.position(entity)
    px, py = x, y

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    index = Int(get(entity.data, "index", 0))

    color = getCassetteColor(index)
    hasCustomColor = false

    hexColor = String(get(entity.data, "customColor", ""))
    if hexColor != "" && length(hexColor) == 6
        color = hexToRGBA(hexColor)
        hasCustomColor = true
    end

    ropeColor = hasCustomColor ? color : get(ropeColors, index, defaultRopeColor)

    # Iteration through all the nodes
    for node in get(entity.data, "nodes", ())
        nx, ny = Int.(node)
        cx, cy = px + width / 2, py + height / 2
        cnx, cny = nx + width / 2, ny + height / 2

        length = sqrt((px - nx)^2 + (py - ny)^2)
        theta = atan(cny - cy, cnx - cx)

        Ahorn.Cairo.save(ctx)

        Ahorn.translate(ctx, cx, cy)
        Ahorn.rotate(ctx, theta)

        Ahorn.setSourceColor(ctx, ropeColor)
        Ahorn.set_antialias(ctx, 1)
        Ahorn.set_line_width(ctx, 1)

        # Offset for rounding errors
        Ahorn.move_to(ctx, 0, 4 + (theta <= 0))
        Ahorn.line_to(ctx, length, 4 + (theta <= 0))

        Ahorn.move_to(ctx, 0, -4 - (theta > 0))
        Ahorn.line_to(ctx, length, -4 - (theta > 0))

        Ahorn.stroke(ctx)

        Ahorn.Cairo.restore(ctx)

        Ahorn.drawSprite(ctx, cog, cnx, cny, tint=color)

        px, py = nx, ny
    end

    renderCassetteBlock(ctx, x, y, width, height, index, color)

    if Bool(get(entity.data, "noReturn", false))
        noReturnSprite = Ahorn.getSprite(crossSprite, "Gameplay")
        Ahorn.drawImage(ctx, noReturnSprite, x + div(width - noReturnSprite.width, 2), y + div(height - noReturnSprite.height, 2), tint=color)
    end
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::CassetteZipMover)
    renderCassetteZipMover(ctx, entity)
end

end


text:
function dataset(dataset_name::String)
    basename = Pkg.dir("EconDatasets", "data")

    ## resolve key names and required package
    fileDict = ["SP500" => "SP500.csv",
                "Sectors" => "sectorAffiliation.csv",
                "UMD" => "UMD.csv",
                "FFF" => "FFF.csv",
                "SP500Ticker" => "SP500TickerSymbols.csv",
                "SP500Industries" => "SP500Industries.csv",
                "Indices" => "Indices.csv",
                "Treasuries" => "Treasuries.csv",
                "DieboldLi" => "DieboldLi.csv"
                ]

    ## get filename
    fname = fileDict[dataset_name]
    filename = joinpath(basename, fname)

    cmdDict = ["SP500" => :(readTimedata($filename)),
               "Sectors" => :(readtable($filename , separator = ' ')),
               "UMD" => :(readTimedata($filename)),
               "FFF" => :(readTimedata($filename)),
               "SP500Ticker" => :(readcsv($filename)),
               "SP500Industries" => :(readtable($filename)),
               "Indices" => :(readTimedata($filename)),
               "Treasuries" => :(readTimedata($filename)),
               "DieboldLi" => :(readTimedata($filename))
               ]

    cmd = cmdDict[dataset_name]

    if !isfile(filename)
        error(string("Unable to locate file $filename - ",
                     "try getDataset(\"$dataset_name\") instead \n"))
    else
        completeCmd = Expr(:(=), :dataVals, cmd)
        eval(completeCmd)
        return dataVals
    end
end


text:
@testset "970.powerful-integers.jl" begin
    @test powerful_integers(2, 3, 10) == Set([2, 3, 4, 5, 7, 9, 10])
    @test powerful_integers(3, 5, 15) == Set([2, 4, 6, 8, 10, 14])
    @test powerful_integers(1, 5, 15) == Set([2, 6])
    @test powerful_integers(1, 1, 5) == Set([2])
end


text:
using TropicalTensors, TropicalTensors.Reversible
using Yao, TropicalYao
using Test
using Viznet
using NiLang, NiLang.AD

@testset "isolve" begin
    Lx, Ly = 3, 4
    lt = ChimeraLattice(Lx, Ly)
    nJ = length(sgbonds(lt))
    sg = Spinglass(lt, ones(Int32, nJ), zeros(Int32, Lx*Ly*8))
    reg = ArrayReg(ones(Tropical{Int32}, 1<<(4*Ly)))
    A = stack4reg(reg, 4*Ly)
    B = stack4reg(reg, Lx-1)
    @test check_inv(isolve, (Int32(0), sg, reg, A, B))
    reg = ArrayReg(ones(Tropical{Int32}, 1<<(4*Ly)))
    A = stack4reg(reg, 4*Ly)
    B = stack4reg(reg, Lx-1)
    @test isolve(Int32(0), sg, reg, A, B)[1] === Int32(nJ)
end

function _solve(Lx, Ly; largemem)
    lt = ChimeraLattice(Lx, Ly)
    nJ = length(sgbonds(lt))
    sg = Spinglass(lt, ones(Int32, nJ), zeros(Int32, Lx*Ly*8))
    reg = ArrayReg(ones(Tropical{Int32}, 1<<(4*Ly)))
    if largemem
        A = stack4reg(reg, Ly*Lx + Ly*(Lx-1)*4 + Lx-1)
        isolve_largemem(Int32(0), sg, reg, A)[1]
    else
        A = stack4reg(reg, Ly*4)
        B = stack4reg(reg, Lx-1)
        isolve(Int32(0), sg, reg, A, B)[1]
    end
end

@testset "isolve largemem" begin
    lt = ChimeraLattice(3, 2)
    _, Lx, Ly = size(lt)
    nJ = length(sgbonds(lt))
    sg = Spinglass(lt, ones(Int32, nJ), zeros(Int32, Lx*Ly*8))
    reg = ArrayReg(ones(Tropical{Int32}, 1<<(4*Ly)))
    A = stack4reg(reg, Ly*Lx + Ly*(Lx-1)*4 + Lx-1)
    @test check_inv(isolve_largemem, (Int32(0), sg, reg, A))
    reg = ArrayReg(ones(Tropical{Int32}, 1<<(4*Ly)))
    A = stack4reg(reg, Ly*Lx + Ly*(Lx-1)*4 + Lx-1)
    @test isolve_largemem(Int32(0), sg, reg, A)[1] === Int32(nJ)
end

@testset "optconfig" begin
    function opt_config_largemem(sg::Spinglass{LT,T}) where {LT,T}
        Ly, Lx = sg.lattice.Ny, sg.lattice.Nx
        reg = ArrayReg(ones(Tropical{T}, 1<<(Ly*4)))
        A = stack4reg(reg, Ly*Lx + Ly*(Lx-1)*4 + Lx-1)
        eng, sg, reg, A = isolve_largemem(T(0.0), sg, reg, A)
        sgg = Spinglass(sg.lattice, GVar.(sg.Js, zero(sg.Js)), GVar.(sg.hs, zero(sg.hs)))
        gres = (~isolve_largemem)(GVar(eng, T(1)), sgg, GVar(reg), GVar(A))
        return TropicalTensors.SpinglassOptConfig(sg, eng, grad.(gres[2].Js), grad.(gres[2].hs))
    end
    sg = rand_spinglass(Int32, ChimeraLattice(3, 4); jt=Randpm(), ht=Zero(), seed=2)
    optc = opt_config(sg)
    optc2 = opt_config_largemem(sg)
    @test optc.grad_Js == optc2.grad_Js
    @test optc.grad_hs == optc2.grad_hs
    @test optc.eng == optc2.eng
end