text:
module RunTests

using UnitJuMP
using Test

import MutableArithmetics
const _MA = MutableArithmetics

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

function test_examples()
    examples = joinpath(@__DIR__, "..", "examples")
    for filename in filter(f -> endswith(f, ".jl"), readdir(examples))
        include(joinpath(examples, filename))
    end
    return
end

function test_jump_variables()
    m = Model()

    @variable(m, x ≥ 0, u"m/s")
    @test x == UnitJuMP.UnitVariableRef(x.variable, u"m/s")
    @test unit(x) == u"m/s"
    @test owner_model(x) === m

    @variable(m, y[1:4], u"km/hr")
    @test y[1] == UnitJuMP.UnitVariableRef(y[1].variable, u"km/hr")
    return
end

function test_jump_constraints()
    m = Model()

    @variable(m, u[1:2] ≥ 0, u"m")
    @variable(m, v[1:4], u"km/hr")
    @variable(m, w, Bin)
    @variable(m, y)

    maxspeed = 4u"ft/s"
    speed = 2.3u"m/s"

    # Various combination of coefficient and variables with and without units
    @constraint(m, 2v[1] ≤ maxspeed)
    @test num_constraints(m, AffExpr, MOI.LessThan{Float64}) == 1

    @constraint(m, 2.3v[1] ≤ maxspeed)
    @test num_constraints(m, AffExpr, MOI.LessThan{Float64}) == 2

    @constraint(m, c1a, speed * w ≤ maxspeed)
    @test num_constraints(m, AffExpr, MOI.LessThan{Float64}) == 3
    @test unit(c1a) == u"m/s"

    @constraint(m, c1b, 40u"km/hr" * w ≤ 15u"m/s")
    @test num_constraints(m, AffExpr, MOI.LessThan{Float64}) == 4
    @test unit(c1b) == u"km/hr"

    @constraint(m, c1c, w * speed ≤ maxspeed)
    @test num_constraints(m, AffExpr, MOI.LessThan{Float64}) == 5
    @test unit(c1c) == u"m/s"

    @constraint(m, c1d, sum(v[i] for i in 1:4) ≤ maxspeed)
    @test num_constraints(m, AffExpr, MOI.LessThan{Float64}) == 6
    @test unit(c1d) == u"km/hr"

    @constraint(m, c1, 2v[1] + 4v[2] ≤ maxspeed)
    @test typeof(c1) <: UnitJuMP.UnitConstraintRef
    @test unit(c1) == u"km/hr"

    @constraint(m, c2, 2v[1] + 4v[2] ≤ maxspeed, u"m/s")
    @test unit(c2) == u"m/s"
    @test normalized_rhs(c2.constraint) ==
          convert(Float64, Unitful.uconvert(u"m/s", maxspeed).val)

    @variable(m, z, Bin)

    maxlength = 1000u"yd"
    period = 1.5u"hr"
    @constraint(m, c3, u[2] + period * v[2] ≤ maxlength * z, u"cm")
    @test unit(c3) == u"cm"

    @constraint(m, c3b, u[2] + 1.5u"hr" * v[2] ≤ 1000u"yd" * z, u"cm")
    @test unit(c3b) == u"cm"
    return
end

function test_mutable_arithmetics()
    m = Model()
    @variable(m, x ≥ 0)
    @variable(m, y)

    xu = UnitJuMP.UnitVariableRef(x, u"m/s")
    yu = UnitJuMP.UnitVariableRef(y, u"km/hr")

    @test _MA.@rewrite(xu) == UnitJuMP.UnitExpression(1x, u"m/s")
    @test _MA.@rewrite(-xu) == UnitJuMP.UnitExpression(-x, u"m/s")
    @test _MA.@rewrite(5xu) == UnitJuMP.UnitExpression(5x, u"m/s")
    @test _MA.@rewrite(xu * 5) == UnitJuMP.UnitExpression(5x, u"m/s")
    @test _MA.@rewrite(xu / 5) == UnitJuMP.UnitExpression(0.2x, u"m/s")

    @test _MA.@rewrite(xu + yu) == UnitJuMP.UnitExpression(x + y / 3.6, u"m/s")
    @test _MA.@rewrite(xu - yu) == UnitJuMP.UnitExpression(x - y / 3.6, u"m/s")
    @test _MA.@rewrite(-xu + yu) ==
          UnitJuMP.UnitExpression(-x + y / 3.6, u"m/s")

    speed = 10u"m/s"
    length = 10u"m"
    freq = 1u"1/s"
    @test _MA.@rewrite(speed * x) == UnitJuMP.UnitExpression(10x, u"m/s")
    @test _MA.@rewrite(x * speed) == UnitJuMP.UnitExpression(10x, u"m/s")
    @test _MA.@rewrite(xu * 10) == UnitJuMP.UnitExpression(10x, u"m/s")
    @test _MA.@rewrite(speed * xu) == UnitJuMP.UnitExpression(10x, u"m^2/s^2")
    @test _MA.@rewrite(xu * speed) == UnitJuMP.UnitExpression(10x, u"m^2/s^2")
    @test _MA.@rewrite(x / speed) == UnitJuMP.UnitExpression(0.1x, u"s/m")
    @test _MA.@rewrite(speed + 2xu) == UnitJuMP.UnitExpression(10 + 2x, u"m/s")
    @test _MA.@rewrite(-speed + xu) == UnitJuMP.UnitExpression(-10 + x, u"m/s")
    @test _MA.@rewrite(2xu - speed) == UnitJuMP.UnitExpression(2x - 10, u"m/s")

    @test _MA.@rewrite(xu + 0.5 * speed) ==
          UnitJuMP.UnitExpression(x + 5, u"m/s")
    @test _MA.@rewrite(xu - 0.5 * speed) ==
          UnitJuMP.UnitExpression(x - 5, u"m/s")
    @test _MA.@rewrite(xu + speed * 0.5) ==
          UnitJuMP.UnitExpression(x + 5, u"m/s")
    @test _MA.@rewrite(xu - speed * 0.5) ==
          UnitJuMP.UnitExpression(x - 5, u"m/s")
    @test _MA.@rewrite(xu + length * freq) ==
          UnitJuMP.UnitExpression(x + 10, u"m/s")
    @test _MA.@rewrite(xu - length * freq) ==
          UnitJuMP.UnitExpression(x - 10, u"m/s")
    @test _MA.@rewrite(0.5 * speed + xu) ==
          UnitJuMP.UnitExpression(x + 5, u"m/s")

    @test _MA.@rewrite(4xu + 2yu - speed) ==
          UnitJuMP.UnitExpression(4x + 2y / 3.6 - 10, u"m/s")
    @test _MA.@rewrite(4xu - 2yu + speed) ==
          UnitJuMP.UnitExpression(4x - 2y / 3.6 + 10, u"m/s")

    ex = xu + yu
    @test _MA.@rewrite(ex + yu + speed) ==
          UnitJuMP.UnitExpression(x + 2y / 3.6 + 10, u"m/s")
    ex = xu + yu
    @test _MA.@rewrite(yu + ex) == UnitJuMP.UnitExpression(3.6x + 2y, u"km/hr")
    ex = xu + yu
    @test _MA.@rewrite(yu - 2ex) == UnitJuMP.UnitExpression(-7.2x - y, u"km/hr")

    @test_throws ErrorException _MA.@rewrite(xu + y)
    @test_throws ErrorException _MA.@rewrite(xu - y)
    @test_throws ErrorException _MA.@rewrite(xu + 5y)
    @test_throws ErrorException _MA.@rewrite(xu - 3y)
    @test_throws ErrorException _MA.@rewrite(y + xu)
    @test_throws ErrorException _MA.@rewrite(y - xu)
    @test_throws ErrorException _MA.@rewrite(y + 5xu)
    @test_throws ErrorException _MA.@rewrite(y - 4xu)
    @test_throws ErrorException _MA.@rewrite(xu + 4)
    @test_throws ErrorException _MA.@rewrite(xu - 4)

    @test_throws ErrorException @constraint(m, x <= xu)
    @test_throws ErrorException @constraint(m, xu <= y)

    return
end

function test_operators()
    m = Model()

    @variable(m, x)
    @variable(m, y)
    @variable(m, z)

    xu = UnitJuMP.UnitVariableRef(x, u"km")
    yu = UnitJuMP.UnitVariableRef(y, u"hr")
    zu = UnitJuMP.UnitVariableRef(z, u"s")

    @test -xu == UnitJuMP.UnitExpression(-x, u"km")
    @test xu + 200u"m" == UnitJuMP.UnitExpression(x + 0.2, u"km")
    @test 200u"m" + xu == UnitJuMP.UnitExpression(1000x + 200, u"m")
    @test xu - 400u"m" == UnitJuMP.UnitExpression(x - 0.4, u"km")
    @test 400u"m" - xu == UnitJuMP.UnitExpression(-1000x + 400, u"m")

    @test 200u"km" * x == UnitJuMP.UnitExpression(200 * x, u"km")
    @test x * 200u"km" == UnitJuMP.UnitExpression(200 * x, u"km")
    @test 200 * xu == UnitJuMP.UnitExpression(200 * x, u"km")
    @test xu * 200 == UnitJuMP.UnitExpression(200 * x, u"km")
    @test 1.5u"s" * xu == UnitJuMP.UnitExpression(1.5 * x, u"s*km")
    @test xu * 1.5u"s" == UnitJuMP.UnitExpression(1.5 * x, u"s*km")

    @test xu / 0.5 == UnitJuMP.UnitExpression(2 * x, u"km")
    @test x / 0.5u"s" == UnitJuMP.UnitExpression(2 * x, u"s^-1")
    @test xu / 0.5u"s" == UnitJuMP.UnitExpression(2 * x, u"km/s")

    @test yu + zu == UnitJuMP.UnitExpression(y + z / 3600, u"hr")
    @test yu - zu == UnitJuMP.UnitExpression(y - z / 3600, u"hr")

    expr = yu + 1800 * zu

    @test expr + 3600u"s" == UnitJuMP.UnitExpression(y + 0.5z + 1, u"hr")
    @test expr - 1800u"s" == UnitJuMP.UnitExpression(y + 0.5z - 0.5, u"hr")
    @test 2u"s" + expr == UnitJuMP.UnitExpression(3600y + 1800z + 2, u"s")
    @test 3600u"s" - expr ==
          UnitJuMP.UnitExpression(-3600y - 1800z + 3600, u"s")

    @test 2 * expr == UnitJuMP.UnitExpression(2y + z, u"hr")
    @test expr * 2 == UnitJuMP.UnitExpression(2y + z, u"hr")
    @test expr / 0.5 == UnitJuMP.UnitExpression(2y + z, u"hr")
    @test 2u"kW" * expr == UnitJuMP.UnitExpression(2y + z, u"kW*hr")
    @test expr * 2u"kW" == UnitJuMP.UnitExpression(2y + z, u"kW*hr")
    @test expr / 0.5u"km" == UnitJuMP.UnitExpression(2y + z, u"hr/km")

    @variable(m, w)
    wu = UnitJuMP.UnitVariableRef(w, u"minute")

    @test expr + wu == UnitJuMP.UnitExpression(y + 0.5z + w / 60, u"hr")
    @test expr - wu == UnitJuMP.UnitExpression(y + 0.5z - w / 60, u"hr")
    @test wu + expr == UnitJuMP.UnitExpression(60y + 30z + w, u"minute")
    @test wu - expr == UnitJuMP.UnitExpression(-60y - 30z + w, u"minute")
    @test expr + yu == UnitJuMP.UnitExpression(2y + 0.5z, u"hr")

    expr2 = 30wu + yu

    @test expr + expr2 == UnitJuMP.UnitExpression(2y + 0.5z + 0.5w, u"hr")
    @test expr2 + expr == UnitJuMP.UnitExpression(120y + 30z + 30w, u"minute")
    @test 2 * expr - expr2 == UnitJuMP.UnitExpression(y + z - 0.5w, u"hr")

    # Combination of unit expression with unitless variables/expressions
    # should throw error
    @test_throws ErrorException w + wu
    @test_throws ErrorException w - wu
    @test_throws ErrorException 4w - wu
    @test_throws ErrorException wu - 4w
    @test_throws ErrorException 4 - wu
    @test_throws ErrorException wu + 4
    @test_throws ErrorException wu - 4
    @test_throws ErrorException x + 2y - wu
    @test_throws ErrorException xu + 2y - wu

    return
end

function test_number_times_quantity_times_variable()
    a = 9.81u"m/s^2"
    model = Model()
    @variable(model, y, u"s")
    expr = @expression(model, 0.5 * a * y)
    @test Unitful.unit(expr) == u"m/s"
    return
end

function test_number_times_quantity_times_variable_to_existing()
    a = 9.81u"m/s^2"
    model = Model()
    @variable(model, x, u"m/s")
    @variable(model, y, u"s")
    expr = @expression(model, x + 0.5 * a * y)
    @test Unitful.unit(expr) == u"m/s"
    return
end

function test_unit_variable()
    model = Model()
    @variable(model, x, u"kW*hr")
    @variable(model, y, u"s")

    # Start values
    set_start_value(x, 1u"MW*hr")
    @test start_value(x) == 1000u"kW*hr"
    @test_throws ErrorException set_start_value(x, 13.2)

    # Lower bounds
    set_lower_bound(x, 10u"MJ")
    @test has_lower_bound(x)
    @test lower_bound(x) == (10 / 3.6)u"kW*hr"
    lref = LowerBoundRef(x)
    @test typeof(lref) <: ConstraintRef
    delete_lower_bound(x)
    @test !has_lower_bound(x)
    @test_throws ErrorException set_lower_bound(x, -4.0)

    # Upper bounds
    set_upper_bound(y, 10u"hr")
    @test has_upper_bound(y)
    @test upper_bound(y) == 36000u"s"
    uref = UpperBoundRef(y)
    @test typeof(uref) <: ConstraintRef
    delete_upper_bound(y)
    @test !has_lower_bound(y)
    @test_throws ErrorException set_upper_bound(y, 42)

    # Fix values
    fix(y, 1u"hr")
    @test is_fixed(y)
    @test fix_value(y) == 3600u"s"
    fref = FixRef(y)
    @test typeof(fref) <: ConstraintRef
    unfix(y)
    @test !is_fixed(y)
    @test_throws ErrorException fix(y, 10)

    # Integer and binary
    set_integer(x)
    @test is_integer(x)
    iref = IntegerRef(x)
    @test typeof(iref) <: ConstraintRef
    unset_integer(x)
    @test !is_integer(x)
    set_binary(y)
    @test is_binary(y)
    bref = BinaryRef(y)
    @test typeof(bref) <: ConstraintRef
    unset_binary(y)
    @test !is_binary(y)

    # Naming
    set_name(x, "xvar")
    @test name(x) == "xvar"

    return
end

end

RunTests.runtests()


text:
# Sioux Falls network data
# http://www.bgu.ac.il/~bargera/tntp/

#Link travel time = free flow time * ( 1 + B * (flow/capacity)^Power ).
#Link generalized cost = Link travel time + toll_factor * toll + distance_factor * distance

# Traffic Assignment Data structure
type TA_Data
    network_name::String

    number_of_zones::Int64
    number_of_nodes::Int64
    first_thru_node::Int64
    number_of_links::Int64

    start_node::Array
    end_node::Array
    capacity::Array
    link_length::Array
    free_flow_time::Array
    B::Array
    power::Array
    speed_limit::Array
    toll::Array
    link_type::Array

    total_od_flow::Float64

    travel_demand::Array
    od_pairs::Array

    toll_factor::Float64
    distance_factor::Float64

    best_objective::Float64
end

function load_ta_network(network_name="Sioux Falls")

    toll_factor = 0
    distance_factor = 0

    if network_name == "Sioux Falls"
        network_data_file = "SiouxFalls_net.txt"
        trip_table_file = "SiouxFalls_trips.txt"
        best_objective = 0
    elseif network_name == "Sioux_simp"
        network_data_file = "SiouxFalls_net_simp.txt"
        trip_table_file = "SiouxFalls_trips_simp.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Jan_AM"
        network_data_file = "East_Massachusetts_net_Jan_AM.txt"
        trip_table_file = "East_Massachusetts_trips_Jan_AM.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Jan_MD"
        network_data_file = "East_Massachusetts_net_Jan_MD.txt"
        trip_table_file = "East_Massachusetts_trips_Jan_MD.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Jan_PM"
        network_data_file = "East_Massachusetts_net_Jan_PM.txt"
        trip_table_file = "East_Massachusetts_trips_Jan_PM.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Jan_NT"
        network_data_file = "East_Massachusetts_net_Jan_NT.txt"
        trip_table_file = "East_Massachusetts_trips_Jan_NT.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Jan_weekend"
        network_data_file = "East_Massachusetts_net_Jan_weekend.txt"
        trip_table_file = "East_Massachusetts_trips_Jan_weekend.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Apr_AM"
        network_data_file = "East_Massachusetts_net_Apr_AM.txt"
        trip_table_file = "East_Massachusetts_trips_Apr_AM.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Apr_MD"
        network_data_file = "East_Massachusetts_net_Apr_MD.txt"
        trip_table_file = "East_Massachusetts_trips_Apr_MD.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Apr_PM"
        network_data_file = "East_Massachusetts_net_Apr_PM.txt"
        trip_table_file = "East_Massachusetts_trips_Apr_PM.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Apr_NT"
        network_data_file = "East_Massachusetts_net_Apr_NT.txt"
        trip_table_file = "East_Massachusetts_trips_Apr_NT.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Apr_weekend"
        network_data_file = "East_Massachusetts_net_Apr_weekend.txt"
        trip_table_file = "East_Massachusetts_trips_Apr_weekend.txt"
	best_objective = 0
    elseif network_name == "East_Massachusetts_Jul_AM"
        network_data_file = "East_Massachusetts_net_Jul_AM.txt"
        trip_table_file = "East_Massachusetts_trips_Jul_AM.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Jul_MD"
        network_data_file = "East_Massachusetts_net_Jul_MD.txt"
        trip_table_file = "East_Massachusetts_trips_Jul_MD.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Jul_PM"
        network_data_file = "East_Massachusetts_net_Jul_PM.txt"
        trip_table_file = "East_Massachusetts_trips_Jul_PM.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Jul_NT"
        network_data_file = "East_Massachusetts_net_Jul_NT.txt"
        trip_table_file = "East_Massachusetts_trips_Jul_NT.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Jul_weekend"
        network_data_file = "East_Massachusetts_net_Jul_weekend.txt"
        trip_table_file = "East_Massachusetts_trips_Jul_weekend.txt"
	best_objective = 0
    elseif network_name == "East_Massachusetts_Oct_AM"
        network_data_file = "East_Massachusetts_net_Oct_AM.txt"
        trip_table_file = "East_Massachusetts_trips_Oct_AM.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Oct_MD"
        network_data_file = "East_Massachusetts_net_Oct_MD.txt"
        trip_table_file = "East_Massachusetts_trips_Oct_MD.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Oct_PM"
        network_data_file = "East_Massachusetts_net_Oct_PM.txt"
        trip_table_file = "East_Massachusetts_trips_Oct_PM.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Oct_NT"
        network_data_file = "East_Massachusetts_net_Oct_NT.txt"
        trip_table_file = "East_Massachusetts_trips_Oct_NT.txt"
        best_objective = 0
    elseif network_name == "East_Massachusetts_Oct_weekend"
        network_data_file = "East_Massachusetts_net_Oct_weekend.txt"
        trip_table_file = "East_Massachusetts_trips_Oct_weekend.txt"
	best_objective = 0
    elseif network_name == "Sioux Falls simplified"
        network_data_file = "SiouxFalls_net_simplified.txt"
        trip_table_file = "SiouxFalls_trips_simplified.txt"
        best_objective = 0
    end

    network_data_file = joinpath("../data_traffic_assignment_uni-class", network_data_file)
    trip_table_file = joinpath("../data_traffic_assignment_uni-class", trip_table_file)

    ##################################################
    # Network Data
    ##################################################

    number_of_zones = 0
    number_of_links = 0
    number_of_nodes = 0
    first_thru_node = 0

    n = open(network_data_file, "r")

    while (line=readline(n)) != ""
        if contains(line, "<NUMBER OF ZONES>")
            number_of_zones = parseint( line[ search(line, '>')+1 : end-1 ] )
        elseif contains(line, "<NUMBER OF NODES>")
            number_of_nodes = parseint( line[ search(line, '>')+1 : end-1 ] )
        elseif contains(line, "<FIRST THRU NODE>")
            first_thru_node = parseint( line[ search(line, '>')+1 : end-1 ] )
        elseif contains(line, "<NUMBER OF LINKS>")
            number_of_links = parseint( line[ search(line, '>')+1 : end-1 ] )
        elseif contains(line, "<END OF METADATA>")
            break
        end
    end

    @assert number_of_links > 0

    start_node = Array(Int, number_of_links)
    end_node = Array(Int, number_of_links)
    capacity = zeros(number_of_links)
    link_length = zeros(number_of_links)
    free_flow_time = zeros(number_of_links)
    B = zeros(number_of_links)
    power = zeros(number_of_links)
    speed_limit = zeros(number_of_links)
    toll = zeros(number_of_links)
    link_type = Array(Int, number_of_links)

    idx = 1
    while (line=readline(n)) != ""
        if contains(line, "~")
            continue
        end

        if contains(line, ";")
            line = strip(line, '\n')
            line = strip(line, ';')

            numbers = split(line)

            start_node[idx] = parseint(numbers[1])
            end_node[idx] = parseint(numbers[2])
            capacity[idx] = parsefloat(numbers[3])
            link_length[idx] = parsefloat(numbers[4])
            free_flow_time[idx] = parsefloat(numbers[5])
            B[idx] = parsefloat(numbers[6])
            power[idx] = parsefloat(numbers[7])
            speed_limit[idx] = parsefloat(numbers[8])
            toll[idx] = parsefloat(numbers[9])
            link_type[idx] = parsefloat(numbers[10])

            idx = idx + 1
        end
    end

    ##################################################
    # Trip Table
    ##################################################

    number_of_zones_trip = 0
    total_od_flow = 0

    f = open(trip_table_file, "r")

    while (line=readline(f)) != ""
        if contains(line, "<NUMBER OF ZONES>")
            number_of_zones_trip = parseint( line[ search(line, '>')+1 : end-1 ] )
        elseif contains(line, "<TOTAL OD FLOW>")
            total_od_flow = parsefloat( line[ search(line, '>')+1 : end-1 ] )
        elseif contains(line, "<END OF METADATA>")
            break
        end
    end

    @assert number_of_zones_trip == number_of_zones # Check if number_of_zone is same in both txt files

    travel_demand = zeros(number_of_zones, number_of_zones)
    od_pairs = Array((Int64,Int64),0)
    while (line=readline(f)) != ""
        if contains(line, "Origin")
            origin = parseint( split(line)[2] )
        elseif contains(line, ";")
            pairs = split(line, ";")
            for i=1:size(pairs)[1]
                if contains(pairs[i], ":")
                    pair = split(pairs[i], ":")
                    destination = parseint( strip(pair[1]) )
                    od_flow = parsefloat( strip(pair[2]) )
                    travel_demand[origin, destination] = od_flow
                    push!(od_pairs, (origin, destination))
                    # println("origin=$origin, destination=$destination, flow=$od_flow")
                end
            end
        end
    end

    # Preparing data to return
    ta_data = TA_Data(
        network_name,
        number_of_zones,
        number_of_nodes,
        first_thru_node,
        number_of_links,
        start_node,
        end_node,
        capacity,
        link_length,
        free_flow_time,
        B,
        power,
        speed_limit,
        toll,
        link_type,
        total_od_flow,
        travel_demand,
        od_pairs,
        toll_factor,
        distance_factor,
        best_objective)

    return ta_data

end # end of load_network function

ext: jl

max_stars_repo_name: jingzbu/InverseVITraffic

max_issues_repo_name: jingzbu/InverseVITraffic

max_forks_repo_name: jingzbu/InverseVITraffic


text:
module MajorizationExtrema

using LinearAlgebra, Random
using UnPack
using LinearAlgebra: dot
using ForwardDiff

export TV, tracedist, ≺, ≻ # distances, majorization
export majmax, majmin, localbound # majorization-extrema
export randprobvec, randunitary, randdm # generate random data

include("random_utilities.jl")

include("classical_case.jl")

include("probvecmults.jl")

# reduce to classical via eigenvalues
include("quantum_case.jl")

include("entropies.jl")

"""
    localbound(f, p, ϵ)

Requires `f` to be Schur convex or Schur concave (but does not check this condition).

* If `p` is a vector, it is a assumed to be a probability vector, and `localbound` returns `δ` such that `|f(p) - f(q)| <= δ` for any probability vector `q` with `TV(p, q) <= ϵ`.
* If `p` is a matrix, it is assumed to be a density matrix, and `localbound` returns `δ` such that `|f(p) - f(q)| <= δ` for any density matrix `q` with `tracedist(p, q) <= ϵ`.
"""
function localbound(f, p, ϵ)
    fp = f(p)
    f_max = f(majmax(p, ϵ))
    f_min = f(majmin(p, ϵ))
    return max(abs(fp - f_max), abs(fp - f_min))
end

end # module


text:
@test string(nd) == "2022-07-28T10:20:08.350789420"
@test string(nd; sep=Underscore) == "2022-07-28T10:20:08.350_789_420"


text:
#=
move types:
 - NNI = nearest neighbor interchange, as defined in
   Gambette, van Iersel, Jones, Lafond, Pardi and Scornavacca (2017).
   Rearrangement moves on rooted phylogenetic networks.
   PLOS Computational Biology 13(8):e1005611.
 - change the direction of a hybrid edge
 - add or remove a hybrid edge
constrained to user-defined topological constraints (like an outgroup)

in `moves.jl`, functions are tailored to level-1 networks;
here, functions apply to semidirected networks of all levels.
=#

"""
Type for various topological constraints, such as:

1. a set of taxa forming a clade in the major tree
2. a set of individuals belonging to the same species
3. a set of taxa forming an outgroup clade or grade in the major tree
4. a set of taxa forming a clade in any one of the displayed trees
5. a set of taxa forming an outgroup clade or grade in any one of the displayed trees
"""
struct TopologyConstraint
    "type of constraint: 1, 2, or 3. Types 4 and 5 not implemented yet"
    type::UInt8
    "names of taxa in the constraint (clade / species / outgroup members)"
    taxonname::Set{String}
    """
    node numbers of taxa in the constraint (clade / species / outgroup members).
    warning: interpretation will be network-dependent.
    """
    taxonnum::Set{Int}
    "number of the stem edge of the constrained group"
    edgenum::Int
end
function TopologyConstraint(type::UInt8, taxonname::Set{String}, net::HybridNetwork)
    taxonnum = Set{Int}() # fixit: use the network instead
    edgenum = -1 # fixit: check that the network satisfies the constraint and return the associated edge
    TopologyConstraint(type, taxonname, taxonnum, edgenum)
end
# fixit: use these constraints in the functions below,
# to check that the proposed NNI move is acceptable

#- see test_moves_semidirected.jl for test file with example uses

#= TODO
- add option to modify branch lengths during NNI:
  `nni!` that take detailed input, and output of same signature
- `unzip` (bool) to constrain branch lengths to 0 below a hybrid node:
  make this part of the optimization process. Or option for `nni!` functions?
- RR move, during optimization: avoid the re-optimization of the likelihood.
  accept the move, just change inheritance values to get same likelihood
  (still update direct / forward likelihoods and tree priors?)
- add options to forbid non-tree-child networks
- add moves that change the root position (and edge directions accordingly)
  without re-calculating the likelihood (accept the move): perhaps recalculate forward / direct likelihoods
  why: the root position affects the feasibility of the NNIs starting from a BR configuration
=#

"""
    nni!(net::HybridNetwork, e::Edge, constraint::Vector{TopologyConstraint},
    no3cycle=true::Bool)

Attempt to perform a nearest neighbor interchange (NNI) around edge `e`,
randomly chosen among all possible NNIs (e.g 3, sometimes more depending on `e`)
satisfying the constraints, and such that the new network is a DAG. The number of
possible NNI moves around an edge depends on whether the edge's parent/child nodes
are tree or hybrid nodes. This is calculated by [`nnimax`](@ref).

Option `no3cycle` forbids moves that would create a 3-cycle in the network.
When `no3cycle` = false, 2-cycle and 3-cycles may be generated by nni!().

Assumptions:
- assumes that the starting network does not have 3-cycles. This function
does not check the input network for the presence of 2- and 3-cycles.
- assumes the edge field `isChild1` is correct in the input network. (This field
will be correct in the output network.)

Output:
information indicating how to undo the move or `nothing` if all NNIs failed.

# examples

```jldoctest
# checks for 3cycle
str_level1 = "(((8,9),(((((1,2,3),4),(5)#H1),(#H1,(6,7))))#H2),(#H2,10));"
net_level1 = readTopology(str_level1);
undoinfo = nni!(net_level1, net_level1.edge[3])
nni!(undoinfo...)
writeTopology(net_level1) == str_level1
true

# does not check for 3cycle
undoinfo = nni!(net_level1, net_level1.edge[3], no3cycle=false)
nni!(undoinfo...)
writeTopology(net_level1) == str_level1
true
```
"""
function nni!(net::HybridNetwork, e::Edge, constraint::Vector{TopologyConstraint},
    no3cycle=true::Bool)
    hybparent = getParent(e).hybrid
    nnirange = 0x01:nnimax(e) # 0x01 = 1 but UInt8 instead of Int
    nnis = Random.shuffle(nnirange)
    for nummove in nnis # iterate through all possible NNIs, but in random order
        moveinfo = nni!(net, e, nummove, no3cycle)
        moveinfo !== nothing || continue # to next possible NNI
        if checknetwork(net, constraint)
            return moveinfo
        else nni!(moveinfo...); end # undo the previous NNI, try again
    end
    return nothing # if we get to this point, all NNIs failed
end

"""
    nnimax(e::Edge)

Return the number of NNI moves around edge `e`,
assuming that the network is semi-directed.
Return 0 if `e` is not internal, and more generally if either node
attached to `e` does not have 3 edges.

Output: UInt8 (e.g. `0x02`)
"""
function nnimax(e::Edge)
    length(e.node[1].edge) == 3 || return 0x00
    length(e.node[2].edge) == 3 || return 0x00
    # e.hybrid and hybrid parent: RR case, 2 possible NNIs only
    # e.hybrid and tree parent:   BR case, 3 or 6 NNIs if e is may contain the root
    # e not hybrid, hyb parent:   RB case, 4 NNIs
    # e not hybrid, tree parent:  BB case, 2 NNIs if directed, 8 if undirected
    hybparent = getParent(e).hybrid
    n = (e.hybrid ? (hybparent ? 0x02 : (e.containRoot ? 0x06 : 0x03)) : # RR & BR
                    (hybparent ? 0x04 : (e.containRoot ? 0x08 : 0x02)))  # RB & BB
    return n
end

"""
    nni!(net::HybridNetwork, uv::Edge, nummove::UInt8, no3cycle=true::Bool)

Modify `net` with a nearest neighbor interchange (NNI) around edge `uv`.
Return the information necessary to undo the NNI, or `nothing` if the move
was not successful (such as if the resulting graph was not acyclic (not a DAG) or if
the focus edge is adjacent to a polytomy). If the move fails, the network is not
modified.
`nummove` specifies which of the available NNIs is performed.

rooted-NNI options according to Gambette et al. (2017), fig. 8:
* BB: 2 moves, both to BB, if directed edges. 8 moves if undirected.
* RR: 2 moves, both to RR.
* BR: 3 moves, 1 RB & 2 BRs, if directed. 6 moves if e is undirected.
* RB: 4 moves, all 4 BRs.
The extra options are due to assuming a semi-directed network, whereas
Gambette et al (2017) describe options for rooted networks.
On a semi-directed network, there might be a choice of how to direct
the edges that may contain the root, e.g. choice of e=uv versus vu, and
choice of labelling adjacent nodes as α/β (BB), or as α/γ (BR).

`no3cycle` prevents moves that would make a 3-cycle by checking for problem 4 cycles.

The edge field `isChild1` is assumed to be correct according to the `net.root`.
"""
function nni!(net::HybridNetwork, uv::Edge, nummove::UInt8, no3cycle=true::Bool)
    nmovemax = nnimax(uv)
    if nmovemax == 0x00 # uv not internal, or polytomy e.g. species represented by polytomy with >2 indiv
        return nothing
    end
    nummove > 0x00 || error("nummove $(nummove) must be >0")
    nummove <= nmovemax || error("nummove $(nummove) must be <= $nmovemax for edge number $(uv.number)")
    u = getParent(uv)
    v = getChild(uv)
    ## TASK 1: grab the edges adjacent to uv: αu, βu, vγ, vδ and adjacent nodes
    # get edges αu & βu connected to u
    # α = u's major parent if it exists, β = u's last child or minor parent
    if u.hybrid
        αu = getMajorParentEdge(u)
        βu = getMinorParentEdge(u)
        α = getParent(αu)
        β = getParent(βu)
    else # u may not have any parent, e.g. if root node
        # pick αu = parent edge if possible, first edge of u (other than uv) otherwise
        labs = [edgerelation(e, u, uv) for e in u.edge]
        pti = findfirst(isequal(:parent), labs) # parent index: there should be 1 or 0
        ci = findall(  isequal(:child), labs)  # vector. there should be 1 or 2
        if pti === nothing
            length(ci) == 2 || error("node $(u.number) should have 2 children other than $(v.number)")
            pti = popfirst!(ci)
            αu = u.edge[pti]
            α = getChild(αu)
        else
            αu = u.edge[pti]
            α = getParent(αu)
        end
        βu = u.edge[ci[1]]
        β = getChild(βu)
    end
    # get edges vδ & vγ connected to v
    # δ = v's last child, γ = other child or v's parent other than u
    labs = [edgerelation(e, v, uv) for e in v.edge]
    if v.hybrid # then v must have another parent edge, other than uv
        vγ = v.edge[findfirst(isequal(:parent), labs)] # γ = getParent(vδ) ?this should be vγ right?
            #on net_hybridladder edge 1 get ERROR: ArgumentError: invalid index: nothing of type Nothing
        vδ = v.edge[findfirst(isequal(:child), labs)]
        γ = getParent(vγ)
        δ = getChild(vδ)
    else
        ci = findall(isequal(:child), labs)
        length(ci) == 2 || error("node $(v.number) should have 2 children")
        vγ = v.edge[ci[1]] # γ = getChild(vγ)
        vδ = v.edge[ci[2]]
        γ = getChild(vγ)
        δ = getChild(vδ)
    end
    ## TASK 2: semi-directed network swaps:
    ## swap u <-> v, α <-> β if undirected and according to move number
    if !u.hybrid && uv.containRoot # case BB or BR, with uv that may contain the root
        nmoves = 0x03 # number of moves in the rooted (directed) case
        if !v.hybrid # BB, uv and all adjacent edges are undirected: 8 moves
            if nummove > 0x04  # switch u & v with prob 1/2
                nummove -= 0x04  # now in 1,2,3,4
                (u,v) = (v,u) #?could this be combined with the line below?
                (α,αu,β,βu,γ,vγ,δ,vδ) = (γ,vγ,δ,vδ,α,αu,β,βu)
            end
            nmoves = 0x02
        end
        # now nmoves should be in 1, ..., 2*nmoves: 1:4 (BB) or 1:6 (BR) #? should this be 2*nummoves?
        # next: switch α & β with probability 1/2
        if nummove > nmoves
            nummove -= nmoves # now nummoves in 1,2 (BB) or 1,2,3 (BR)
            (α,αu,β,βu) = (β,βu,α,αu)
        end
    end
    ## TASK 3: rooted network swaps, then do the NNI
    if u.hybrid # RR or RB cases
        # swap α & β: reduce NNIs from 2 to 1 (RR) or from 4 to 2 (RB)
        if nummove == 0x02 || nummove == 0x04
            (α,αu,β,βu) = (β,βu,α,αu)
        end
        if !v.hybrid && nummove > 0x02 # case RB, moves 3 or 4
            (γ,vγ,vδ,δ) = (δ,vδ,vγ,γ)
        end
        # detach β and graft onto vδ
        if no3cycle && problem4cycle(α,γ, β,δ) return nothing; end
        res = nni!(αu, u, uv, v, vδ)
    else # u = bifurcation: BB or BR cases
        if !v.hybrid # case BB: 2 options
            if nummove == 0x02
                (γ,vγ,vδ,δ) = (δ,vδ,vγ,γ)
            end
            # detach β and graft onto vδ
            if no3cycle && problem4cycle(α,γ, β,δ) return nothing; end
            res = nni!(αu, u, uv, v, vδ)
        else # case BR: 3 options
            # DAG check
            # if α parent of u, check if β -> γ (directed or undirected)
            # if undirected with β parent of u, check α -> γ
            # if undirected and u = root:
            #   moves 1 and 3 will fail if α -> γ
            #   moves 2 and 3 will fail if β -> γ
            # nummove 3 always creates a nonDAG when u is the root
            αparentu = getChild(αu)===u
            βparentu = getChild(βu)===u
            if αparentu
                if isdescendant(γ, β) return nothing; end
            elseif βparentu
                if isdescendant(γ, α) return nothing; end
            else # cases when u is root
                if nummove == 0x01 && isdescendant(γ, α)
                    return nothing
                elseif nummove == 0x02 && isdescendant(γ, β)
                    return nothing
                elseif nummove == 0x03 # fail: not DAG
                    return nothing
                end
            end
            if nummove == 0x01 # graft γ onto α
                if no3cycle && problem4cycle(α,γ, β,δ) return nothing; end
                res = nni!(vδ, v, uv, u, αu)
            elseif nummove == 0x02 # graft γ onto β
                if no3cycle && problem4cycle(β,γ, α,δ) return nothing; end
                res = nni!(vδ, v, uv, u, βu)
            else # nummove == 0x03
                if αparentu # if α->u: graft δ onto α
                    if no3cycle && problem4cycle(α,δ, β,γ) return nothing; end
                    res = nni!(vγ, v, uv, u, αu)
                else # if β->u: graft δ onto β
                    if no3cycle && problem4cycle(α,γ, β,δ) return nothing; end
                    res = nni!(vγ, v, uv, u, βu)
                # case when u is root already excluded (not DAG)
                end
            end
        end
    end
    return res
end

"""
    nni!(αu, u, uv::Edge, v, vδ)
    nni!(αu,u,uv,v,vδ, flip::Bool, inner::Bool, indices)

Transform a network locally around the focus edge `uv` with
the following NNI, that detaches u-β and grafts it onto vδ:

    α - u -- v ------ δ
        |    |
        β    γ

    α ------ v -- u - δ
             |    |
             γ    β

`flip` boolean indicates if the uv edge was flipped
`inner` boolean indicates if edges αu and uv both point toward node u,
i.e. α->u<-v<-δ. If this is true, we flip the hybrid status of αu and vδ.

`indices` give indices for nodes and edges u_in_αu, αu_in_u, vδ_in_v, and v_in_vδ.
These are interpreted as:
    u_in_αu: the index for u in the edge αu
    αu_in_u: the index for αu in node u
    vδ_in_v: the index for vδ in node v
    v_in_vδ: the index for v in edge vδ

**Warnings**:
- *No* check of assumed adjacencies
- Not implemented for cases that are not necessary thanks to symmetry,
  such as cases covered by `nni!(vδ, v, uv, u, αu)` or `nni!(βu, u, v, vγ)`.
  More specifically, these cases are not implemented (and not checked):
  * u not hybrid & v hybrid
  * u hybrid, v not hybrid, α -> u <- v -> δ
- Because of this, `nni(αu,u,uv,v,vδ, ...)` should not be used directly;
  use instead `nni!(net, uv, move_number)`.

Node numbers and edge numbers are not modified.
Edge `uv` keeps its direction unchanged *unless* the directions were
`α -> u -> v -> δ` or `α <- u <- v <- δ`,
in which case the direction of `uv` is flipped.

The second version's input has the same signature as the output, but
will undo the NNI more easily. This means that if `output = nni!(input)`,
then `nni!(output...)` is valid and undoes the first operation.

Right now, branch lengths are not modified. Future versions might implement
options to modify branch lengths.
"""
function nni!(αu::Edge, u::Node, uv::Edge, v::Node, vδ::Edge)
    # find indices to detach αu from u, and to detach v from vδ
    u_in_αu = findfirst(n->n===u, αu.node)
    αu_in_u = findfirst(e->e===αu, u.edge)
    vδ_in_v = findfirst(e->e===vδ, v.edge)
    v_in_vδ = findfirst(n->n===v, vδ.node)
    # none of them should be 'nothing' --not checked
    αu_child = getChild(αu)
    uv_child = getChild(uv)
    vδ_child = getChild(vδ)
    # flip the direction of uv if α->u->v->δ or α<-u<-v<-δ
    flip = (αu_child === u && uv_child === v && vδ_child !== v) ||
             (αu_child !== u && uv_child === u && vδ_child === v)
    inner = (flip ? false : αu_child === u && vδ_child === v )
    # fixit: calculate new edge lengths here; give them as extra arguments to nni! below
    res = nni!(αu,u,uv,v,vδ, flip,inner, u_in_αu,αu_in_u,vδ_in_v,v_in_vδ)
    return res
end
function nni!(αu::Edge, u::Node, uv::Edge, v::Node, vδ::Edge,
              flip::Bool, inner::Bool,
              u_in_αu::Int, αu_in_u::Int, vδ_in_v::Int, v_in_vδ::Int)
              # extra arguments: new edge lengths
    # fixit: extract old edge lengths; update to new lengths; return old lengths as extra items
    ## TASK 1: do the NNI swap: to detach & re-attach at once
    (αu.node[u_in_αu], u.edge[αu_in_u], v.edge[vδ_in_v], vδ.node[v_in_vδ]) =
        (v, vδ, αu, u)
    ## TASK 2: update edges' directions
    if flip
        uv.isChild1 = !uv.isChild1
    end
    ## TASK 3: update hybrid status of nodes & edges, and edges containRoot field
    if u.hybrid && !v.hybrid
        if flip # RB, BR 1 or BR 2': flip hybrid status of uv and (αu or vδ)
            if uv.hybrid # BR 1 or BR 2'
                uv.hybrid = false
                vδ.hybrid = true
                norootbelow!(uv)
            else # RB
                uv.hybrid = true
                αu.hybrid = false
                if αu.containRoot
                    allowrootbelow!(v, αu)
                end
            end
        else # assumes α<-u or v<-δ --- but not checked
            # not implemented: α->u<-v->δ: do uv.hybrid=false and γv.hybrid=true
            if inner # BR 3 or 4': α->u<-v<-δ: flip hybrid status of αu and vδ
                αu.hybrid = false
                vδ.hybrid = true
                # containRoot: nothing to update
            else # BR 1' or 2: α<-u<-v->δ : hybrid status just fine
                norootbelow!(vδ)
                if uv.containRoot
                    allowrootbelow!(αu)
                end
            end
        end
    end
    # throw error if u not hybrid & v hybrid? (does not need to be implemented) #?
    return vδ, u, uv, v, αu, flip, inner, v_in_vδ, αu_in_u, vδ_in_v, u_in_αu
end

"""
    problem4cycle(β::Node, δ::Node, α::Node, γ::Node)

Checks if the focus edge uv has a 4 cycle that could lead to a 3 cycle
    after an chosen NNI. Return true if there is a problem 4 cycle, false if none.
"""
function problem4cycle(β::Node, δ::Node, α::Node, γ::Node)
    isconnected(β, δ) || isconnected(α, γ)
end

"""
    cladesviolated(network::HybridNetwork, cladeconstraints)

Check if network violates user-given clade constraints. Return false if passes,
return true if violates clades.
"""
function cladesviolated(net::HybridNetwork, cladeconstraints::Dict)
    return false
end


text:
include("LD.jl")
using MATLAB
using ForwardDiff
using ProgressMeter
x=range(.25,stop=1.2,length=Int(1e4))
y1=@showprogress map(Dimer_Lagrangian_Descriptor_Function,x);
logz1=log1p.(abs.(y1))

mat"figure()"
mat"plot($x,$y1)"
# mat"plot($x,$logz1);hold on;"
delta=1e-12
g(x)=( Dimer_Lagrangian_Descriptor_Function(x+delta)-Dimer_Lagrangian_Descriptor_Function(x-delta))/(2*delta)
f(x)=abs.(( Dimer_Lagrangian_Descriptor_Function(x-delta)-Dimer_Lagrangian_Descriptor_Function(x)+Dimer_Lagrangian_Descriptor_Function(x+delta))/(delta^2))
# g(x)=( Dimer_Lagrangian_Descriptor_Function(x-.001)-Dimer_Lagrangian_Descriptor_Function(x+.001))/(.002)
y2=@showprogress map(g,x);
logz2=log10.(abs.(y2))
mat"figure()"
mat"plot($x,$y2,'k')"

mat"figure()"
mat"plot($x,$logz2,'k.')"

# y3=@showprogress map(f,x)
# # logz=log1p.(abs.(y2))
# # logz3=log1p.(abs.(y3))

# mat"figure()"
# mat"plot($x,$y3)"
# # mat"plot($x,$logz3)"


text:
using JuliaFormatter
format(
    ".",
    verbose = true,
    format_markdown = true,
    conditional_to_if = true,
    format_docstrings = false,
    always_use_return = true,
    import_to_using = true,
    remove_extra_newlines = true,
    whitespace_ops_in_indices = true,
    whitespace_typedefs = true,
    margin = 80,
    always_for_in = true,
    short_to_long_function_def = true,
)
nothing


text:
using PkgBenchmark
result = PkgBenchmark.readresults(joinpath(@__DIR__, "result.json"))

using Markdown
# pretty-print to screen
# let io = IOBuffer()
#     PkgBenchmark.export_markdown(io, result)
#     seekstart(io)
#     md = Markdown.parse(io)
#     display(md)
# end

# pretty-print to file
@info "Pretty-printing markdown to file"
PkgBenchmark.export_markdown("result.md", result)

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

struct Nullable{T}
    hasvalue::Bool
    value::T

    Nullable{T}() where {T} = new(false)
    Nullable{T}(value::T, hasvalue::Bool=true) where {T} = new(hasvalue, value)
end

max_stars_repo_name: xorJane/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
@testset "colorfulness" begin

    @info "test: colorfulness"

    # Test against a simple result calculated by hand

    x = [RGB(1,0,0), RGB(0,1,0), RGB(0,0,1)]
    @test trunc(colorfulness(x)) == 337

    # Black is not colorful

    x = [RGB(0,0,0), RGB(0,0,0), RGB(0,0,0)]
    @test trunc(colorfulness(x)) == 0

    # White is not colorful

    x = [RGB(1,1,1), RGB(1,1,1), RGB(1,1,1)]
    @test trunc(colorfulness(x)) == 0

    # A grayscale image is not colorful

    cameraman = testimage("cameraman")
    @test colorfulness(cameraman) == 0

    # A color image with only grays is not colorful

    x = convert(Array{Float64}, cameraman)
    img = RGB.(x, x, x)
    @test colorfulness(img) == 0

    # Lena 256 is a reduced color image and so should be less colorful
    # than the original

    imga = testimage("lena_color_256")
    imgb = testimage("lena_color_512")

    @test  colorfulness(imga) < colorfulness(imgb)

    # Test all invocation styles are equivalent

    c1 = colorfulness(imga)
    c2 = hasler_and_susstrunk_m3(imga)
    c3 = colorfulness(HASLER_AND_SUSSTRUNK_M3(), imga)
    c4 = assess(HASLER_AND_SUSSTRUNK_M3(), imga)

    @test c1 == c2 == c3 == c4

end


text:
# SnoopComplie

using SnoopCompile

### Log the compiles
# This only needs to be run once (to generate "/tmp/images_compiles.csv")
SnoopCompile.@snoop "C:\\temp\\__snoop_dtm_compiles.csv" begin
    include(Pkg.dir("DecisionTrees", "test", "runtests.jl"))
end

### Parse the compiles and generate precompilation scripts
# This can be run repeatedly to tweak the scripts
data = SnoopCompile.read("C:\\temp\\__snoop_dtm_compiles.csv")

pc = SnoopCompile.parcel(reverse!(data[2]))
SnoopCompile.write("c:\\temp\\precompile", pc)


text:
module Zeiss

using ZeissCAN29
import ..Gamepad.GamepadState
export update_velocity

const V_SCALE = Int32(100)
const VELOCITY = Int32[50, 100, 250, 500, 1000, 4000] .* V_SCALE # zeiss focus speed ≈ 10x UMP speed

function update_velocity(gamepad::GamepadState, speed_step::Integer)::Int32
    LT, RT = gamepad.L.trigger, gamepad.R.trigger
    return LT > 0.4 ? (RT > 0.4 ? 0 : -VELOCITY[speed_step]) :
                      (RT > 0.4 ? VELOCITY[speed_step] : 0)
end

end # module


text:
using MoreMatrices

trim{T}(coords::Vector{Matrix{T}}, n::Int, m::Int) = [coords[i][:, 1:n] for i in 1:m]

function random_perturb!{T}(A::Vector{Matrix{T}}, σ::T)
    m = length(A)
    d, n = size(A[1])
    N = Vector{Matrix{T}}(m)

    for i in 1:m
        N[i] = σ * randn(d, n)
        copy!(A[i], A[i] + N[i])
    end

    return -N
end

function random_permute!{T}(A::Vector{Matrix{T}})
    m = length(A)
    d, n = size(A[1])
    P = Vector{Matrix{T}}(m)

    for i in 1:m
        P[i] = eye(n)[:, randperm(n)]
        copy!(A[i], A[i] * P[i]')
    end

    return P
end

function random_rotate!{T}(A::Vector{Matrix{T}})
    m = length(A)
    d, n = size(A[1])
    R = Vector{Matrix{T}}(m)

    for i in 1:m
        R[i] = random_rotation()
        copy!(A[i], R[i]' * A[i])
    end

    return R
end

function random_rotate_reflect!{T}(A::Vector{Matrix{T}})
    m = length(A)
    d, n = size(A[1])
    R = Vector{Matrix{T}}(m)

    for i in 1:m
        R[i] = random_orthogonal(3)
        copy!(A[i], R[i]' * A[i])
    end

    return R
end

function pointsets_transform{T}(A::Vector{Matrix{T}}, R::Vector{Matrix{T}}, P::Vector{Matrix{T}})
    d, n = size(A[1])
    m = min(length(R), length(A), length(P))
    A_hat = [R[i] * A[i] * P[i] for i in 1:m]
end

function pointsets_transform{T}(A::Vector{Matrix{T}}, R::Vector{Matrix{T}})
    d, n = size(A[1])
    m = min(length(R), length(A))
    A_hat = [R[i] * A[i] for i in 1:m]
end

function pointsets_rmsd{T}(A::Vector{Matrix{T}})
    d, n = size(A[1])
    m = length(A)
    X = mean(A)

    return vecnorm(map(A) do B vecnorm(X - B) end) / sqrt(d*n*m)
end

pointsets_rmsd{T}(A::Vector{Matrix{T}}, R::Vector{Matrix{T}}, P::Vector{Matrix{T}}) = pointsets_rmsd(pointsets_transform(A, R, P))


text:
@testset "Long read datastores" begin
    function get_fastq_seqs(file)
        seqs = map(open(FASTQ.Reader, file) do rdr
            collect(rdr)
        end) do rec
            FASTQ.sequence(LongDNASeq, rec)
        end
        return seqs
    end

    function check_round_trip(FQ)
        seqs = get_fastq_seqs(FQ)
        fq = open(FASTQ.Reader, FQ)
        ds = LongReads(fq, "human-nanopore.lrds", "human-nanopore", UInt64(0))
        ds2 = open(LongReads, "human-nanopore.lrds")
        ds_seqs = collect(ds)
        ds2_seqs = collect(ds2)
        return ds_seqs == seqs == ds2_seqs
    end

    @testset "Human oxford nanopore 2D consensus reads tester" begin
        @test check_round_trip("human_nanopore_tester_2D.fastq")

        ds = open(LongReads, "human-nanopore.lrds")
        ds2 = open(LongReads, "human-nanopore.lrds")

        @test ReadDatastores.index(ds) == ReadDatastores.index(ds2)
        @test ReadDatastores.index(ds)[1] == ReadDatastores.index(ds2)[1]
        @test firstindex(ds) == firstindex(ds2) == 1
        @test lastindex(ds) == lastindex(ds2) == 10

        @test Base.IteratorSize(ds) == Base.IteratorSize(ds2) == Base.HasLength()
        @test Base.IteratorEltype(ds) == Base.IteratorEltype(ds2) == Base.HasEltype()
        @test Base.eltype(ds) == Base.eltype(ds2) == LongSequence{DNAAlphabet{4}}

    end
end


text:
module CategoricalArrays
    export CategoricalPool, CategoricalValue
    export AbstractCategoricalArray, AbstractCategoricalVector, AbstractCategoricalMatrix,
           CategoricalArray, CategoricalVector, CategoricalMatrix
    export LevelsException, OrderedLevelsException

    export categorical, compress, decompress, droplevels!, levels, levels!, levelcode,
           isordered, ordered!
    export cut, recode, recode!

    using JSON
    using DataAPI
    using Missings
    using Printf

    include("typedefs.jl")

    include("pool.jl")
    include("value.jl")

    include("array.jl")
    include("missingarray.jl")
    include("subarray.jl")

    include("extras.jl")
    include("recode.jl")

    include("deprecated.jl")
end


text:
using NCDatasets
using Formatting
using ArgParse

function parse_commandline()

    s = ArgParseSettings()

    @add_arg_table s begin

        "--domain-file"
            help = "Domain file."
            arg_type = String
            required = true
    end

    return parse_args(s)
end

parsed = parse_commandline()

function maketime(years)

    dom = [31.0, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
    sum_dom = sum(dom)
    sum_dom == 365 || throw(ErrorException("Sum of dom is $(sum_dom) rather than 365."))

    _t    = zeros(Float64, length(dom))
    _bnds = zeros(Float64, 2, length(dom))

    for m=1:length(dom)
        #bnds[m, 1] = beg of month  m
        #bnds[m, 2] = end of month  m
        if m==1
            _bnds[1, m] = 0.0
        else
            _bnds[1, m] = _bnds[2, m-1]
        end

        _bnds[2, m] = _bnds[1, m] + dom[m]

        _t[m] = (_bnds[1, m] + _bnds[2, m]) / 2.0
    end

    t    = zeros(Float64,    12*years)
    bnds = zeros(Float64, 2, 12*years)

    for y = 1:years
        i_offset = (y-1)*12
        t_offset = (y-1)*sum_dom
        t[i_offset+1:i_offset+12]       .+= _t    .+ t_offset
        bnds[:, i_offset+1:i_offset+12] .+= _bnds .+ t_offset
    end

    return t, bnds
end

Dataset(parsed["domain-file"], "r") do ds

    global Nx = ds.dim["ni"]
    global Ny = ds.dim["nj"]

    global lon = nomissing(ds["xc"][:], NaN) .|> deg2rad
    global lat = nomissing(ds["yc"][:], NaN) .|> deg2rad

end

z_w = collect( Float64, range(0.0, -600.0, length=11) )
Nz = length(z_w)-1

t, bnds = maketime(5)

tt = reshape(t, 1, 1, :)
llat = reshape(lat, size(lat)..., 1)
llon = reshape(lon, size(lon)..., 1)

TAUX  = sin.((tt.-0) / 365.0 * 2π .+ llon) .* cos.(llat) .* sin.(llon       ) * 0.1
TAUY = sin.((tt.-0) / 365.0 * 2π .+ llon) .* cos.(llat) .* sin.(llon + llat) * 0.1
SWFLX      = (1.0 .+ sin.((tt.-0) / 365.0 * 2π .+ llon * 2))/2.0 .* cos.(llat) * (-1000.0)
NSWFLX     = (1.0 .+ cos.((tt.-0) / 365.0 * 2π .+ llon * 2))/2.0 .* cos.(llat) * (  500.0)
VSFLX      = (1.0 .+ cos.((tt.-0) / 365.0 * 2π .+ llon * 2))/2.0 .* cos.(llat) * ( 1.0 / 86400.0 )

println("Outputting atmospheric forcing.nc")
Dataset("atm_forcing.nc", "c") do ds

    defDim(ds, "time", Inf)
    defDim(ds, "d2", 2)
    defDim(ds, "nlon", Nx)
    defDim(ds, "nlat", Ny)

    defVar(ds, "time", t, ("time", ), ; attrib = Dict(
        "long_name" => "time",
        "bounds"    => "time_bound",
        "calendar"  =>  "noleap",
        "units"     => "days since 0001-01-01 00:00:00",
    ))

    defVar(ds, "time_bound", bnds, ("d2", "time"), ; attrib = Dict(
        "long_name" => "boundaries for time-averaging interval",
        "units"     => "days since 0001-01-01 00:00:00",
    ))

    defVar(ds, "TAUX", TAUX, ("nlon", "nlat", "time",), ; attrib = Dict(
        "long_name" => "TAUX east",
        "units"     => "N/m^2",
    ))

    defVar(ds, "TAUY", TAUY, ("nlon", "nlat", "time", ), ; attrib = Dict(
        "long_name" => "TAUY north",
        "units"     => "N/m^2",
    ))

    defVar(ds, "SWFLX", SWFLX, ("nlon", "nlat", "time", ), ; attrib = Dict(
        "long_name" => "SWFLX",
        "units"     => "W/m^2",
    ))

    defVar(ds, "NSWFLX", NSWFLX, ("nlon", "nlat", "time", ), ; attrib = Dict(
        "long_name" => "NSWFLX",
        "units"     => "W/m^2",
    ))

    defVar(ds, "VSFLX", VSFLX, ("nlon", "nlat", "time", ), ; attrib = Dict(
        "long_name" => "VSFLX",
        "units"     => "kg/s/m^2",
    ))

    defVar(ds, "z_w_top", z_w[1:end-1], ("z_t", ), ; attrib = Dict(
        "long_name" => "Z coordinate at grid top.",
        "units"     => "m",
    ))

end

t, bnds = maketime(1)

tt = reshape(t, 1, 1, :)
llat = reshape(lat, size(lat)..., 1)
llon = reshape(lon, size(lon)..., 1)

z_w = collect( Float64, range(0.0, -300.0, length=11) )
z_t = (z_w[1:end-1] + z_w[2:end])/2.0
Nz = length(z_w)-1

HMXL = (1.0 .+ sin.((tt.-0) / 365.0 * 2π .+ llon * 2))/2.0 .* cos.(llat) * 100 .+ 30.0

tt = reshape(t, 1, 1, 1, :)
llat = reshape(lat, size(lat)..., 1, 1)
llon = reshape(lon, size(lon)..., 1, 1)
zz_t = reshape(z_t, 1, 1, :, 1)

TEMP  = exp.(zz_t / 50.0) .* (1.0 .+ sin.((tt.-0) / 365.0 * 2π .+ llon * 2))/2.0 .* cos.(llat) * 30.0
SALT  = exp.(zz_t / 50.0) .* (1.0 .+ sin.((tt.-0) / 365.0 * 2π .+ llon * 2))/2.0 .* cos.(llat) * 30.0

QFLX_TEMP  = exp.(zz_t / 50.0) .* cos.((tt.-0) / 365.0 * 2π .+ llon * 2) .* cos.(llat) * 1e-6 * 3996 * 1026
QFLX_SALT  = exp.(zz_t / 50.0) .* sin.((tt.-0) / 365.0 * 2π .+ llon * 2) .* cos.(llat) * 1e-6

println("Outputting ocean forcing.")
Dataset("ocn_forcing.nc", "c") do ds

    defDim(ds, "time", Inf)
    defDim(ds, "d2", 2)
    defDim(ds, "nlon", Nx)
    defDim(ds, "nlat", Ny)
    defDim(ds, "z_t", Nz)

    defVar(ds, "time", t, ("time", ), ; attrib = Dict(
        "long_name" => "time",
        "bounds"    => "time_bound",
        "calendar"  =>  "noleap",
        "units"     => "days since 0001-01-01 00:00:00",
    ))

    defVar(ds, "time_bound", bnds, ("d2", "time"), ; attrib = Dict(
        "long_name" => "boundaries for time-averaging interval",
        "units"     => "days since 0001-01-01 00:00:00",
    ))

    defVar(ds, "TEMP", TEMP, ("nlon", "nlat", "z_t", "time",), ; attrib = Dict(
        "long_name" => "TEMP",
        "units"     => "degC",
    ))

    defVar(ds, "SALT", SALT, ("nlon", "nlat", "z_t", "time",), ; attrib = Dict(
        "long_name" => "SALT",
        "units"     => "kg/m^3",
    ))

    defVar(ds, "QFLX_TEMP", QFLX_TEMP, ("nlon", "nlat", "z_t", "time", ), ; attrib = Dict(
        "long_name" => "QFLX_TEMP",
        "units"     => "W/m^3",
    ))

    defVar(ds, "QFLX_SALT", QFLX_SALT, ("nlon", "nlat", "z_t", "time", ), ; attrib = Dict(
        "long_name" => "QFLX_SALT",
        "units"     => "kg/m^3/s/m^3",
    ))

    defVar(ds, "HMXL", HMXL, ("nlon", "nlat", "time", ), ; attrib = Dict(
        "long_name" => "HMXL",
        "units"     => "m",
    ))

    defVar(ds, "z_w_top", z_w[1:end-1], ("z_t", ), ; attrib = Dict(
        "long_name" => "Z coordinate at grid top.",
        "units"     => "m",
    ))

    defVar(ds, "z_w_bot", z_w[2:end], ("z_t", ), ; attrib = Dict(
        "long_name" => "Z coordinate at grid bottom.",
        "units"     => "m",
    ))

end

max_stars_repo_name: meteorologytoday/IOM

max_issues_repo_name: meteorologytoday/IOM

max_forks_repo_name: meteorologytoday/IOM


text:
using Mocking

batch_get_item_1_patch = @patch function batch_get_item(a...; b...)
    return Dict{String, Any}("Responses" => Dict{String, Any}("Music" => Any[Dict{String, Any}("AlbumTitle" => Dict{String, Any}("S" => "Somewhat Famous")), Dict{String, Any}("AlbumTitle" => Dict{String, Any}("S" => "Blue Sky Blues")), Dict{String, Any}("AlbumTitle" => Dict{String, Any}("S" => "Louder Than Ever"))]))
end

batch_write_item_1_patch = @patch function batch_write_item(a...; b...)
    return Dict{String, Any}()
end

create_table_1_patch = @patch function create_table(a...; b...)
    return Dict{String, Any}("TableDescription" => Dict{String, Any}("TableSizeBytes" => 0, "AttributeDefinitions" => Any[Dict{String, Any}("AttributeName" => "Artist", "AttributeType" => "S"), Dict{String, Any}("AttributeName" => "SongTitle", "AttributeType" => "S")], "ItemCount" => 0, "TableName" => "Music", "ProvisionedThroughput" => Dict{String, Any}("ReadCapacityUnits" => 5, "WriteCapacityUnits" => 5), "CreationDateTime" => "1421866952.062", "KeySchema" => Any[Dict{String, Any}("AttributeName" => "Artist", "KeyType" => "HASH"), Dict{String, Any}("AttributeName" => "SongTitle", "KeyType" => "RANGE")], "TableStatus" => "CREATING"))
end

delete_item_1_patch = @patch function delete_item(a...; b...)
    return Dict{String, Any}("ConsumedCapacity" => Dict{String, Any}("CapacityUnits" => 1, "TableName" => "Music"))
end

delete_table_1_patch = @patch function delete_table(a...; b...)
    return Dict{String, Any}("TableDescription" => Dict{String, Any}("TableSizeBytes" => 0, "ItemCount" => 0, "TableName" => "Music", "ProvisionedThroughput" => Dict{String, Any}("ReadCapacityUnits" => 5, "NumberOfDecreasesToday" => 1, "WriteCapacityUnits" => 5), "TableStatus" => "DELETING"))
end

describe_limits_1_patch = @patch function describe_limits(a...; b...)
    return Dict{String, Any}("AccountMaxWriteCapacityUnits" => 20000, "TableMaxWriteCapacityUnits" => 10000, "AccountMaxReadCapacityUnits" => 20000, "TableMaxReadCapacityUnits" => 10000)
end

describe_table_1_patch = @patch function describe_table(a...; b...)
    return Dict{String, Any}("Table" => Dict{String, Any}("TableSizeBytes" => 0, "AttributeDefinitions" => Any[Dict{String, Any}("AttributeName" => "Artist", "AttributeType" => "S"), Dict{String, Any}("AttributeName" => "SongTitle", "AttributeType" => "S")], "ItemCount" => 0, "TableName" => "Music", "ProvisionedThroughput" => Dict{String, Any}("ReadCapacityUnits" => 5, "NumberOfDecreasesToday" => 1, "WriteCapacityUnits" => 5), "CreationDateTime" => "1421866952.062", "KeySchema" => Any[Dict{String, Any}("AttributeName" => "Artist", "KeyType" => "HASH"), Dict{String, Any}("AttributeName" => "SongTitle", "KeyType" => "RANGE")], "TableStatus" => "ACTIVE"))
end

get_item_1_patch = @patch function get_item(a...; b...)
    return Dict{String, Any}("Item" => Dict{String, Any}("SongTitle" => Dict{String, Any}("S" => "Happy Day"), "AlbumTitle" => Dict{String, Any}("S" => "Songs About Life"), "Artist" => Dict{String, Any}("S" => "Acme Band")))
end

list_tables_1_patch = @patch function list_tables(a...; b...)
    return Dict{String, Any}("TableNames" => Any["Forum", "ProductCatalog", "Reply", "Thread"])
end

put_item_1_patch = @patch function put_item(a...; b...)
    return Dict{String, Any}("ConsumedCapacity" => Dict{String, Any}("CapacityUnits" => 1, "TableName" => "Music"))
end

query_1_patch = @patch function query(a...; b...)
    return Dict{String, Any}("ScannedCount" => 2, "Items" => Any[Dict{String, Any}("SongTitle" => Dict{String, Any}("S" => "Call Me Today"))], "ConsumedCapacity" => Dict{String, Any}(), "Count" => 2)
end

scan_1_patch = @patch function scan(a...; b...)
    return Dict{String, Any}("ScannedCount" => 3, "Items" => Any[Dict{String, Any}("SongTitle" => Dict{String, Any}("S" => "Call Me Today"), "AlbumTitle" => Dict{String, Any}("S" => "Somewhat Famous")), Dict{String, Any}("SongTitle" => Dict{String, Any}("S" => "Scared of My Shadow"), "AlbumTitle" => Dict{String, Any}("S" => "Blue Sky Blues"))], "ConsumedCapacity" => Dict{String, Any}(), "Count" => 2)
end

update_item_1_patch = @patch function update_item(a...; b...)
    return Dict{String, Any}("Attributes" => Dict{String, Any}("SongTitle" => Dict{String, Any}("S" => "Happy Day"), "AlbumTitle" => Dict{String, Any}("S" => "Louder Than Ever"), "Artist" => Dict{String, Any}("S" => "Acme Band"), "Year" => Dict{String, Any}("N" => "2015")))
end

update_table_1_patch = @patch function update_table(a...; b...)
    return Dict{String, Any}("TableDescription" => Dict{String, Any}("TableSizeBytes" => 0, "AttributeDefinitions" => Any[Dict{String, Any}("AttributeName" => "Artist", "AttributeType" => "S"), Dict{String, Any}("AttributeName" => "SongTitle", "AttributeType" => "S")], "ItemCount" => 0, "TableName" => "MusicCollection", "ProvisionedThroughput" => Dict{String, Any}("ReadCapacityUnits" => 1, "LastIncreaseDateTime" => "1421874759.194", "NumberOfDecreasesToday" => 1, "WriteCapacityUnits" => 1), "CreationDateTime" => "1421866952.062", "KeySchema" => Any[Dict{String, Any}("AttributeName" => "Artist", "KeyType" => "HASH"), Dict{String, Any}("AttributeName" => "SongTitle", "KeyType" => "RANGE")], "TableStatus" => "UPDATING"))
end


text:
include(joinpath("src", "EpidemicModel.jl"))
begin
    configs = Dict()
    for (i, behavior_reduction_factor) in enumerate(0:0.25:1)
        configs[i] = EpidemicModel.Config(
            n_nodes=100,
            n0=3,
            k=3,
            agentcount=100,
            transmission_prob=0.05,
            test_prob=0.8,
            quarantine_duration=140,
            behavior_reduction_factor=behavior_reduction_factor,
            behavior_activated=true,
            tick_authority_recommendation=100,
            tick_authority_policy=200,
            n_public_spaces=10,
            initial_infectives=1,
            n_iter=1000
        )
    end
    EpidemicModel.run_batch(configs=configs, replicates=2, batchname="test", compress=true)
end


text:
using MKL

using Logging: global_logger
using TerminalLoggers: TerminalLogger
global_logger(TerminalLogger())

using OrdinaryDiffEq, LinearAlgebra, Plots, DiffEqOperators, SparseArrays, Distributions, CUDA
using Symbolics, SparseDiffTools
using PreallocationTools
using PatternFormation
using LinearSolve
using IncompleteLU
using BenchmarkTools
using LazyArrays

BLAS.set_num_threads(16)
CUDA.allowscalar(false)

#=============================================================
Preconditioners
==============================================================#
function incompletelu(W,du,u,p,t,newW,Plprev,Prprev,solverdata)
  if newW === nothing || newW
    Pl = ilu(convert(AbstractMatrix,W), τ = 150.0)
  else
    Pl = Plprev
  end
  Pl,nothing
end
Base.eltype(::IncompleteLU.ILUFactorization{Tv,Ti}) where {Tv,Ti} = Tv
#==============================================================
Problem parameters
===============================================================#
println("Setting up problem parameters")
ū = 0.88
v̄ = 0.37
D = 0.5
β = -1.0 #FIXE
r₁ = 3.0
r₂ = 10.0
α = -(v̄ - ū*v̄*r₂)/(ū - ū*v̄^2*r₁)
γ = -α
δ = 1.0; #1/(5*pi)^2*(α + β*D)/(2D)
γᵤ =  10.0
γᵥ =   1.0
#==============================================================
Grid parameters
===============================================================#
dx = Float64(1/5)
dy = Float64(1/5)
N = 500
tspan = (0.0, Float64(1e2))
#=============================================================
Initial Conditions
=============================================================#
println("Initial Conditions")
r01 = Float64.(ū .+ rand(Uniform(-ū/100, +ū/100), N, N))
r02 = Float64.(v̄ .+ rand(Uniform(-v̄/100, +v̄/100), N, N))
x = range(0, N*dx, length=N)
y = range(0, N*dy, length=N)
#r01, x, y, dx, dy = init_cond(:NoisePatches, N, c01 = ū, dx = dx)
#r02, x, y, dx, dy = init_cond(:NoisePatches, N, c01 = v̄, dx = dx)
#r01, x, y, dx, dy = init_cond(:ErMnO3_MnOnly, N, M = 15, α=3000, BC=:Neumann, c01 =ū, f = +0.01, dx = dx)
#r02, x, y, dx, dy = init_cond(:ErMnO3_ErOnly, N, M = 15, α=500, BC=:Neumann, c01 =v̄, f = -0.01, dx = dx)
#r01 .= 1 .- r01
r0 = cat(CuArray(r01), CuArray(r02), dims=3)
r0_cpu = cat(r01, r02, dims=3)
#=============================================================
Finite differences matrices
==============================================================#
println("Dealing with finite differencing.")
Q = Neumann0BC(1.0, 1)
#Q = PeriodicBC(Float64)
D2 = CenteredDifference(2, 2, 1.0, N)
L2 = D2*Q

A2 = CuArray(sparse(L2)[1])
A2_cpu = sparse(L2)[1]
#============================================================
Caches
=============================================================#
println("Dealing with caches.")
A2u = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
uA2 = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
D2u = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
A2v = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
vA2 = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
D2v = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
A2u_i = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
u_iA2 = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
u_i = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
A2v_i = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
v_iA2 = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
v_i = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
u_ni = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))
v_ni = DiffEqBase.dualcache(CuArray(zeros(Float64, N,N)))

cache =(A2u, uA2, D2u, A2v, vA2, D2v, A2u_i, u_iA2, u_i, A2v_i, v_iA2, v_i)
p = (D, δ, γᵤ, γᵥ, α, β, γ, r₁, r₂, dx, dy, A2, cache)
p_cpu = (D, δ, γᵤ, γᵥ, α, β, γ, r₁, r₂, dx, dy, A2_cpu, cache)
#==========================================================
Advanced versoin of the PDE, for solving
===========================================================#
function advanced_version!(dr,r,p,t)
  D, δ, γᵤ, γᵥ, α, β, γ, r₁, r₂, dx, dy, A2, cache = p
  A2u, uA2, D2u, A2v, vA2, D2v, A2u_i, u_iA2, u_i, A2v_i, v_iA2, v_i =  cache

  u = @view r[:,:,1]
  v = @view r[:,:,2]
  du = @view dr[:,:,1]
  dv = @view dr[:,:,2]

  A2u = get_tmp(A2u, 0.0)
  uA2 = get_tmp(uA2, 0.0)
  D2u = get_tmp(D2u, 0.0)
  A2v = get_tmp(A2v, 0.0)
  vA2 = get_tmp(vA2, 0.0)
  D2v = get_tmp(D2v, 0.0)
  A2u_i = get_tmp(A2u_i, 0.0)
  u_iA2 = get_tmp(u_iA2, 0.0)
  u_i = get_tmp(u_i, 0.0)
  A2v_i = get_tmp(A2v_i, 0.0)
  v_iA2 = get_tmp(v_iA2, 0.0)
  v_i = get_tmp(v_i, 0.0)

  mul!(A2u,A2,u)
  mul!(uA2,u,CuArray(A2'))
  mul!(A2v,A2,v)
  mul!(vA2,v,CuArray(A2'))
  @. u_i = 2.0*u*(u-1)*(2.0*u-1) - γᵤ*(1.0/dy^2*A2u + 1.0/dx^2*uA2) # 4 allocations
  @. v_i = 2.0*v*(v-1)*(2.0*v-1) - γᵥ*(1.0/dy^2*A2v + 1.0/dx^2*vA2) # 4 allocatoins
  mul!(A2u_i,A2,u_i)
  mul!(u_iA2,u_i,CuArray(A2'))
  mul!(A2v_i,A2,v_i)
  mul!(v_iA2,v_i,CuArray(A2'))

  @. D2u = δ*D*(1.0/dy^2*A2u_i + 1.0/dx^2*u_iA2)
  @. D2v = δ*(1.0/dy^2*A2v_i + 1.0/dx^2*v_iA2)
  @. du = D2u + α*u*(1.0-r₁*v^2) + v*(1.0 - r₂*u)
  @. dv = D2v + β*v*(1.0+α*r₁/β*u*v) + u*(γ + r₂*v)
  nothing
end
#=============================================================
Basic Version for computing Jacobian sparsity
==============================================================#
function basic_version!(dr,r,p,t)
  D, δ, γᵤ, γᵥ, α, β, γ, r₁, r₂, dx, dy, A2, cache = p

  u = @view r[:,:,1]
  v = @view r[:,:,2]

  u_n = 2u.*(u.-1).*(2u.-1) .- γᵤ*(1/dy^2*A2*u + 1/dx^2*u*A2')
  v_n = 2v.*(v.-1).*(2v.-1) .- γᵥ*(1/dy^2*A2*v + 1/dx^2*v*A2')
  D2u = D*δ*(1/dy^2*A2*u_n + 1/dx^2*u_n*A2')
  D2v = δ*(1/dy^2*A2*v_n + 1/dx^2*v_n*A2')
  dr[:,:,1] = D2u .+ α*u.*(1 .- r₁*v.^2) .+ v.*(1 .- r₂*u)
  dr[:,:,2] = D2v .+ β*v.*(1 .+ α*r₁/β*u.*v) .+ u.*(γ .+ r₂*v)
  nothing
end
#=================================================================
Compute the Jacobian sparsity pattern
==================================================================#
#println("Jacobian sparsity")
#dr0 = copy(r0_cpu)
#@time jac_sparsity = Symbolics.jacobian_sparsity((dr,r)->basic_version!(dr,r,p_cpu,0.),dr0,r0_cpu)
#f = ODEFunction(advanced_version!;jac_prototype=Float64.(jac_sparsity))

#==================================================================
Solve the system
==================================================================#
println("Solving")
prob = ODEProblem(advanced_version!,r0,tspan,p)
CUDA.allowscalar(false)
#@time sol = solve(prob, TRBDF2(linsolve=KrylovJL_GMRES(), precs=incompletelu, concrete_jac=true), saveat=range(0.0, stop=tspan[2], length=101), progress=true, progress_steps=1);
@time sol = solve(prob, ROCK2(), saveat=range(0.0, stop=tspan[2], length=101), progress=true, progress_steps=1);
#@time sol = solve(prob, ROCK2(), saveat=[tspan[2]], progress=true, progress_steps=1);

#==================================================================
Plot
==================================================================#
#println("Plotting")
anim = @animate for i in 1:length(sol.t)
  title = "t = $(sol.t[i])"
  #u = kron(ones(3,3), sol.u[i][:,:,1])
  u = Array(sol.u[i])[:,:,2]
  heatmap(x, y, u, c=:roma, aspect_ratio=1, axis=([], false), title=title, colorbar=false)
  #heatmap(u, c=:roma, aspect_ratio=dy/dx, axis=([], false), title=title, colorbar=false)
end

gif(anim, fps=10)
#title = "t = $(tspan[2])"
#u = Array(sol.u[end])[:,:,2]
#heatmap(u, c=:roma, aspect_ratio=dy/dx, axis=([], false), title=title, colorbar=true)


text:
function getij(x,y,paramsx0,paramsy0,paramsx1,paramsy1,paramsdx,paramsdy,paramsnx,paramsny)

i=-1; j=-1;
valid = true;
# i / y
y0 = paramsy0;
for n=1:paramsny

    if y0 <= y && y < y0+paramsdy
        i=n;
        break
    end
    y0 = y0+paramsdy;
end

# j / x
x0 = paramsx0;
for n=1:paramsnx
    if x0 <= x && x < x0+paramsdx
        j=n;
        break
    end
    x0 = x0+paramsdx;
end

if i==-1 || j==-1
    valid = false;
end

return i,j,valid;

end


text:
@testset "minimal_overlap_test" begin
    # overlap [empty] [empty] = bothEmpty
    # overlap [empty] [1.0,2.0] = firstEmpty
    # overlap [1.0,2.0] [empty] = secondEmpty
    # overlap [-infinity,2.0] [3.0,infinity] = before
    # overlap [-infinity,2.0] [3.0,4.0] = before
    # overlap [2.0,2.0] [3.0,4.0] = before
    # overlap [1.0,2.0] [3.0,4.0] = before
    # overlap [1.0,2.0] [3.0,3.0] = before
    # overlap [2.0,2.0] [3.0,3.0] = before
    # overlap [2.0,2.0] [3.0,infinity] = before
    # overlap [-infinity,2.0] [2.0,3.0] = meets
    # overlap [1.0,2.0] [2.0,3.0] = meets
    # overlap [1.0,2.0] [2.0,infinity] = meets
    # overlap [1.0,2.0] [1.5,2.5] = overlaps
    # overlap [1.0,2.0] [1.0,infinity] = starts
    # overlap [1.0,2.0] [1.0,3.0] = starts
    # overlap [1.0,1.0] [1.0,3.0] = starts
    # overlap [1.0,2.0] [entire] = containedBy
    # overlap [1.0,2.0] [-infinity,3.0] = containedBy
    # overlap [1.0,2.0] [0.0,3.0] = containedBy
    # overlap [2.0,2.0] [0.0,3.0] = containedBy
    # overlap [2.0,2.0] [0.0,infinity] = containedBy
    # overlap [1.0,2.0] [-infinity,2.0] = finishes
    # overlap [1.0,2.0] [0.0,2.0] = finishes
    # overlap [2.0,2.0] [0.0,2.0] = finishes
    # overlap [1.0,2.0] [1.0,2.0] = equals
    # overlap [1.0,1.0] [1.0,1.0] = equals
    # overlap [-infinity,1.0] [-infinity,1.0] = equals
    # overlap [entire] [entire] = equals
    # overlap [3.0,4.0] [2.0,2.0] = after
    # overlap [3.0,4.0] [1.0,2.0] = after
    # overlap [3.0,3.0] [1.0,2.0] = after
    # overlap [3.0,3.0] [2.0,2.0] = after
    # overlap [3.0,infinity] [2.0,2.0] = after
    # overlap [2.0,3.0] [1.0,2.0] = metBy
    # overlap [2.0,3.0] [-infinity,2.0] = metBy
    # overlap [1.5,2.5] [1.0,2.0] = overlappedBy
    # overlap [1.5,2.5] [-infinity,2.0] = overlappedBy
    # overlap [1.0,infinity] [1.0,2.0] = startedBy
    # overlap [1.0,3.0] [1.0,2.0] = startedBy
    # overlap [1.0,3.0] [1.0,1.0] = startedBy
    # overlap [-infinity,3.0] [1.0,2.0] = contains
    # overlap [entire] [1.0,2.0] = contains
    # overlap [0.0,3.0] [1.0,2.0] = contains
    # overlap [0.0,3.0] [2.0,2.0] = contains
    # overlap [-infinity,2.0] [1.0,2.0] = finishedBy
    # overlap [0.0,2.0] [1.0,2.0] = finishedBy
    # overlap [0.0,2.0] [2.0,2.0] = finishedBy
end


text:
function get_filter(::Type{CartesianDerivativesFilter}, arg0::jint)
    return jcall(CartesianDerivativesFilter, "getFilter", CartesianDerivativesFilter, (jint,), arg0)
end

function get_max_order(obj::CartesianDerivativesFilter)
    return jcall(obj, "getMaxOrder", jint, ())
end

function value_of(::Type{CartesianDerivativesFilter}, arg0::JString)
    return jcall(CartesianDerivativesFilter, "valueOf", CartesianDerivativesFilter, (JString,), arg0)
end

function values(::Type{CartesianDerivativesFilter})
    return jcall(CartesianDerivativesFilter, "values", Vector{CartesianDerivativesFilter}, ())
end


text:
struct ComputePipelineCreateInfo
    mHandleRef::Ref{vk.VkComputePipelineCreateInfo}
    mLayout

    function ComputePipelineCreateInfo(
        stage::PipelineShaderStageCreateInfo, # required
        layout; # ::PipelineLayout
        pNext::Ptr{Cvoid} = C_NULL,
        flags::vk.VkPipelineCreateFlags = vk.VkFlags(0),
        basePipelineHandle::vk.VkPipeline = vk.VkPipeline(vk.VK_NULL_HANDLE), # ::VkPipeline
        basePipelineIndex::Int32 = Int32(0)) # ::Int32

            new(Ref(vk.VkComputePipelineCreateInfo(
                vk.VK_STRUCTURE_TYPE_COMPUTE_PIPELINE_CREATE_INFO, # sType::VkStructureType
                pNext, # ::Ptr{Cvoid}
                flags, # ::VkPipelineCreateFlags
                stage.handleRef()[], # ::VkPipelineShaderStageCreateInfo
                layout.handleRef()[], # ::VkPipelineLayout
                basePipelineHandle, # ::VkPipeline
                basePipelineIndex)), # ::Int32)
                layout
            )
    end
end
@class ComputePipelineCreateInfo

function handleRef(this::ComputePipelineCreateInfo)::Ref{vk.VkComputePipelineCreateInfo}
    return this.mHandleRef
end


text:
using EconPDEs, Distributions

Base.@kwdef struct AchdouHanLasryLionsMoll_TwoAssetsModel
    # income process parameters
    κy::Float64 = 0.1
    ybar::Float64 = 1.0
    σy::Float64 = 0.07

    r::Float64 = 0.03
    μR::Float64 = 0.04
    σR::Float64 = 0.1

    # utility parameters
    ρ::Float64 = 0.05
    γ::Float64 = 2.0

    amin::Float64 = 0.0
    amax::Float64 = 1000.0
end

function (m::AchdouHanLasryLionsMoll_TwoAssetsModel)(state::NamedTuple, value::NamedTuple)
    (; κy, σy, ybar, r, μR, σR, ρ, γ, amin, amax) = m
    (; y, a) = state
    (; v, vy_up, vy_down, va_up, va_down, vyy, vya, vaa) = value
    μy = κy * (ybar - y)
    vy = (μy >= 0) ? vy_up : vy_down
    va = va_up
    iter = 0
    @label start
    va = max(va, eps())
    c = va^(-1 / γ)
    k = (μR - r) / σR^2 * (- va / vaa)
    k = clamp(k, 0.0, a - amin)
    μa = y + r * a + (μR - r) * k - c
    if (iter == 0) & (μa <= 0)
        iter += 1
        va = va_down
        @goto start
    end
    σa = k * σR
    # There is no second derivative at 0 so just specify first order derivative
    if (a ≈ amin) && (μa <= 0.0)
        va = (y + r * amin)^(-γ)
        k = 0.0
        c = y + r * amin
        μa = 0.0
    end
    # this branch is unnecessary when individuals dissave at the top (default)
    if (a ≈ amax) && (μa >= 0.0)
        va = ((ρ - r) / γ + r - (1-γ) / (2 * γ) * (μR - r)^2 / (γ * σR^2))^(-γ) * a^(-γ)
    end
    # Since first derivative is upwinded, I can directly impose value of second derivative
    if (a ≈ amax)
        vaa = - m.γ * va / a
    end
    c = va^(-1 / γ)
    k = (μR - r) / σR^2 * (- va / vaa)
    k = clamp(k, 0.0, a - amin)
    μa = y + r * a + (μR - r) * k - c
    σa = k * σR
    vt = - (c^(1 - γ) / (1 - γ) + va * μa + 0.5 * vaa * σa^2 + vy * μy + 0.5 * vyy * σy^2 - ρ * v)
    return (; vt), (; v, c, k, va, vaa, vy, y, a, μa)
end

m = AchdouHanLasryLionsMoll_TwoAssetsModel()
distribution = Gamma(2 * m.κy * m.ybar / m.σy^2, m.σy^2 / (2 * m.κy))
ys = range(quantile(distribution, 0.001), quantile(distribution, 0.999), length = 5)
as = range(m.amin, m.amax, length = 100)
stategrid = OrderedDict(:y => ys, :a => as)
yend = OrderedDict(:v => [log(y + a) for y in stategrid[:y], a in stategrid[:a]])
y, residual_norm = pdesolve(m, stategrid, yend)
#
# # Important: check marginal value of wealth converges to 1.0
# # This happens ONLY if a >= 1000.0. Otherwise with 300 it does not work. This is interesting. Maybe it means there should be a better way to have bordering condition at top
# b = ((m.r + (m.ρ - m.r)/m.γ - (1-m.γ) / (2 * m.γ) * (m.μR - m.r)^2 / (m.γ * m.σR^2)))^(1/(1 - 1/m.γ))
# pw = (result[:v] * (1-m.γ)).^(1/(1-m.γ)-1) .* result[:va] ./ b


text:
# Translated from the Golang Tour's select intro: https://tour.golang.org/concurrency/5
module Fibonacci

using Select
import Base.Threads: @spawn

function fibonacci(c::Channel{>:Int}, quit::Channel{>:Int})
    x, y = 0, 1
    while true
        @select begin
            c <| x    => begin
                x, y = y, x+y
            end
            quit      => begin
                println("quit")
                return
            end
        end
    end
end

function main()
    @sync begin
        c = Channel{Int}()
        quit = Channel{Int}()
        @spawn begin
            for i in 0:9
                println(take!(c))
            end
            put!(quit, 0)
        end
        fibonacci(c, quit)
    end
end

main()

end


text:
# GMM estimation for a sample from Chi^2(theta)
# compare to two method of moments estimators (see chi2mm.m)
using Distributions
n = 30
theta = 3.0
y = rand(Chisq(theta), n)
thetahat = mean(y)
println("mm estimator based on mean: ", thetahat)

thetahat = 0.5*var(y)
println("mm estimator based on variance: ", thetahat)

# GMM
W = eye(2)
obj = theta -> ([theta.-mean(y); theta.-0.5*var(y)]'W*[theta.-mean(y); theta.-0.5*var(y)])
thetahat, junk, junk = fminunc(obj, [2.0])
println("GMM estimator based on mean and variance: ", thetahat)

max_stars_repo_name: nilshg/Econometrics

max_issues_repo_name: nilshg/Econometrics

max_forks_repo_name: nilshg/Econometrics


text:
# Use Welch's t-test to cope with unequal sample sizes and unequal variance.
# Could (should?) use other methods in other cases.

immutable TTest <: HypothesisTest
    test_name::ASCIIString
    # t statistic
    statistic::Float64
    parameter::Float64
    p_value::Float64
    confidence_interval::Vector{Float64}
    estimate::Float64
    null_value::Float64
    alternative_hypothesis::ASCIIString
    method::ASCIIString
    degrees_of_freedom::Float64
end

function test_statistic{S <: Real,
                        T <: Real}(::Type{TTest},
                                   x::Vector{S},
                                   y::Vector{T})
    n_x, n_y = length(x), length(y)
    x_bar, y_bar = mean(x), mean(y)
    s_x, s_y = var(x), var(y)
    s_welch = s_x / n_x + s_y / n_y
    return (x_bar - y_bar) / sqrt(s_welch)
end

test_statistic(t::TTest) = t.statistic

function ci{S <: Real,
            T <: Real}(::Type{TTest},
                       x::Vector{S},
                       y::Vector{T})
    n_x, n_y = length(x), length(y)
    x_bar, y_bar = mean(x), mean(y)
    s_x, s_y = var(x), var(y)
    s_welch = s_x / n_x + s_y / n_y
    t = (x_bar - y_bar) / sqrt(s_welch)
    df_num = (s_x / n_x + s_y / n_y)^2
    df_denom = ((s_x / n_x)^2) / (n_x - 1) + ((s_y / n_y)^2) / (n_y - 1)
    df = df_num / df_denom
    tau = abs(quantile(TDist(df), 0.05 / 2))
    return [(x_bar - y_bar) - tau * sqrt(s_welch), (x_bar - y_bar) + tau * sqrt(s_welch)]
end

ci(t::TTest) = t.confidence_interval

function p_value{S <: Real,
                 T <: Real}(::Type{TTest},
                            x::Vector{S},
                            y::Vector{T})
    n_x, n_y = length(x), length(y)
    x_bar, y_bar = mean(x), mean(y)
    s_x, s_y = var(x), var(y)
    s_welch = s_x / n_x + s_y / n_y
    t = (x_bar - y_bar) / sqrt(s_welch)
    df_num = (s_x / n_x + s_y / n_y)^2
    df_denom = ((s_x / n_x)^2) / (n_x - 1) + ((s_y / n_y)^2) / (n_y - 1)
    df = df_num / df_denom
    return 2.0 * (1.0 - cdf(TDist(df), abs(t)))
end

p_value(t::TTest) = t.p_value

function Base.run{S <: Real,
                  T <: Real}(::Type{TTest},
                             x::Vector{S},
                             y::Vector{T})
    n_x, n_y = length(x), length(y)
    x_bar, y_bar = mean(x), mean(y)
    s_x, s_y = var(x), var(y)
    s_welch = s_x / n_x + s_y / n_y
    t = (x_bar - y_bar) / sqrt(s_welch)
    df_num = (s_x / n_x + s_y / n_y)^2
    df_denom = ((s_x / n_x)^2) / (n_x - 1) + ((s_y / n_y)^2) / (n_y - 1)
    df = df_num / df_denom
    tau = abs(quantile(TDist(df), 0.05 / 2))
    ci = [(x_bar - y_bar) - tau * sqrt(s_welch), (x_bar - y_bar) + tau * sqrt(s_welch)]
    p = 2.0 * (1.0 - cdf(TDist(df), abs(t)))
    return TTest("t-Test",
                 t,
                 x_bar - y_bar,
                 p,
                 ci,
                 x_bar - y_bar,
                 0.0,
                 "Two-sided",
                 "Welch Two Sample t-test",
                 df)
end


text:
include("SupportFunctions.jl")

function updateS(Y,A::Array{Float64,2},lambda,S::Array{Float64,2};maxIter=10000,threshold=1e-5)

    H=A'*A
    AtY=A'*Y

    R=deepcopy(S)

    Sp=S
    L=0
    try
        L=eigmax(H)
    catch
        return (S,0,false)
        error("Error in eigenvalue calculation")
    end

    At=A'
    t=1
    tp=t

    i=1

    converged=false
    while i<maxIter
        tp=t
        t=(1+(1+4*tp^2)^0.5)/2

        w=(tp-1)/t

        R=(1+w)*S-w*Sp
        Sp=S
        try
            #S=plusOp(soft(R-(H*R-AtY)/L,lambda/L))
            S=max(R-(H*R-AtY)/L-lambda/L,0)
        catch
            error("Error in S Iteration calculation")
        end
        if sum(abs(Sp-S))<(threshold*sum(S))
            converged=true
            break
        end

        i+=1
    end
    return (S,i,converged)
end

function updateA(Y,S::Array{Float64,2},A::Array{Float64,2};maxIter=10000,threshold=1e-5)

    Ap=A
    B=zeros(size(A))

    H=S*S'
    L=eigmax(H)
    YSt= Y*S'

    t=1
    tp=t

    i=1
    converged=false
    while i<maxIter

        tp=t
        t=(1+(1+4*tp^2)^0.5)/2

        B = (1+(tp-1)/(t)) * A - (tp-1)/(t) * Ap
        Ap=A
        A = max(B-(B*H-YSt)/L,0)

        if sum(abs(Ap-A))<(threshold*sum(A))
            converged=true
            break
        end
        i+=1
    end
    return (A,i,converged)
end

#this definition was taken from thier matlab code with some modifications
function updateLambda!(A::Array{Float64,2},S::Array{Float64,2},Y,iteration,maxIter, phaseRatio,lambda)
    refinement_beginning = floor(phaseRatio * maxIter)

    if refinement_beginning > iteration

        sigma_residue = 1.4826 * median(abs((Y - A*S)[:].-median((Y - A*S)[:])))
        #linear decrease to tau_MAD*sigma_residue when reaching the refinement steps
        lambda = maximum([sigma_residue, lambda - 1/(refinement_beginning - iteration) * (lambda - sigma_residue)])

    else #during refinement, do not modify lambda
        lambda = lambda
    end

    return lambda
end

function nGMCA(Y::Array{Float64},r;verbose=false,maxIter=5000,threshold=1e-6,phaseRatio=0.30,kickstart=true,submaxiter=200)

    #Randomly build one of the two output matricies
    S=rand(size(Y,2),r)'

    #Fit the first as an approximation of the other
    A=Y/S

    #Estimate an L0
    l=maximum(abs(A'*(A*S-Y)))

    if verbose && kickstart
        print("\nKickstarting BSS...")
    end

    if kickstart
        #Kick Start this show!
        A=plusOp(Y/S)
        S=plusOp(A\Y)
        S=updateS(Y,A,0,S;maxIter=submaxiter)[1]
        (A,S)=m_reinitializeS(A,S, Y, verbose)
        A=updateA(Y,S,A;maxIter=submaxiter)[1]
        if verbose
          print("\nKickstart complete.")
        end
    end

    if verbose
        print("\nL0=",""*string(l,"\nStart Full BSS:\n"))
    end

    #Run the algorithm
    j=1
    for j=1:maxIter
        Ap=A
        Sp=S

        #normalize columns of A
        for i=1:size(A,2)
            A[:,i]= A[:,i]./(sum(A.^2.,1).^.5)[i]
        end

        (A,S)=m_reinitializeS(A,S, Y, verbose)
        (S,S_iter,S_converge)=updateS(Y,A,l,S;maxIter=submaxiter)
        (A,S)=m_reinitializeS(A,S, Y, verbose)
        (A,A_iter,A_converge)=updateA(Y,S,A;maxIter=submaxiter)
        l=updateLambda!(A,S,Y,j,submaxiter,phaseRatio,l)

        if verbose && mod(j,max(div(maxIter,100),1))==0
            print("\r"*string("GI(S) ",round(GI(S),3)," GI(A) ",
            round(GI(A),3),"\t Iterations Total:",j,"  S:",S_iter," A: ",A_iter,
            "\tL:",round(l,5),"\t Dist to Y: ",round(sum((A*S-Y).^2.).^.5),5))
        end

        if mod(j,10)==0 && sum(abs(A-Ap))/length(A)<threshold && sum(abs(S-Sp))/length(S)<threshold
            break
        end
    end
    if verbose
        print("\nDone.\n")
        if j>=maxIter
            print("\nFailed to converge!\n")
        end
    end
    return {A,S}
end

# reinitialization when a source was set to 0
function m_reinitializeS(A,S, Y, verbose)
    #straightforward reinitialization of a line of S and a column of A
    #by picking one column in the residue
    attempt=0
    while sum(sum(S, 2) .== 0) > 0 && attempt<10
        if verbose
            print("\r"*string("Reinitilizing blank rows/columns: ",sum(sum(S,2).==0)))
        end
        indices = [sum(S, 2) .== 0][:]

        (A[:, indices], S[indices, :]) = m_fastExtractNMF(Y - A * S, sum(indices))
        attempt+=1
    end
    (A,S)
end

function m_fastExtractNMF(residual, r)
    if r > 0
        (m, n) = size(residual);
        A = zeros(m, r);
        S = zeros(r, n);
        for i = 1 : r
            residual = max(residual,0)
            if sum(residual[:]) != 0
                #compute square norm of residual to select maximum one
                res2 = sum(residual.^2,1);
                j = findin(res2,maximum(res2));
                if !isempty(j)
                    j = j[1];
                    if res2[j] > 0
                        #normalize maximum residual
                        A[:, i] = residual[:, j] / res2[j].^.5;
                        #compute scalar product with the rest of the residual and keep only
                        #positive coefficients
                        S[i, :] = A[:, i]' * residual;
                        S[i, :] = max(S[i, :], zeros(size(S[i,:])) );
                        #compute new residual
                        residual = residual - A[:, i] * S[i, :];
                    end
                end
            end
        end
    else
        A = [];
        S = [];
    end

    return (A,S)
end

max_stars_repo_name: dahlend/CompressedSensing


text:
# Code for Video 1.3:

year = 2015
condition = true
letter = 'J'
mass_sun = 1.98855e30 # in kg
flower = "iris"

year == 2015 # true
year != 2016 # true

# getting the types of data:
typeof(year)  #Int64
typeof(condition) # Bool
typeof(letter) # Char
typeof(mass_sun) # Float64
typeof(flower) # ASCIIString

# some errors:
#  typeof(mass_sum)
# ERROR: UndefVarError: mass_sum not defined

#  flower + 28
# ERROR: MethodError: `+` has no method matching +(::ASCIIString, ::Int64)
# Closest candidates are:
#   +(::Any, ::Any, ::Any, ::Any...)
#   +(::Int64, ::Int64)
#   +(::Complex{Bool}, ::Real)
#   ...

isa(flower, ASCIIString) # true

convert(Float32, year) # 2015.0f0
Float32(year)
# convert(Int64, mass_sun) # ERROR: InexactError() in convert at int.jl:209

1:24
1.0:0.5:20.0
typeof(1:24) # UnitRange{Int64}
typeof(1.0:0.5:20.0) # FloatRange{Float64}
range(1, 10) # 1:10

char1 = 'I'
char2 = 'β'
Int(char2) # 946
char3 = '\u4E80' # '亀'

species1 = "I. setosa"     # "I. setosa"
typeof(species1) # ASCIIString

arab1 = "الله" # "الله"
typeof(arab1) # UTF8String

species2 = "I. versicolor" # "I. versicolor"
species3 = "I. virginica"  # "I. virginica"
species1[1] # 'I'
species1[1:3] # "I. "

julia> species1[0]
# ERROR: BoundsError: attempt to access 9-element Array{UInt8,1}:
#  0x49
#  0x2e
#  0x20
#  0x73
#  0x65
#  0x74
#  0x6f
#  0x73
#  0x61
#   at index [0]
#  in getindex at ascii.jl:13

# species1[5] = '3' # ERROR: MethodError: `setindex!` has no method matching setindex!(::ASCIIString,
                  # ::Char, ::Int64)

println("The first iris species is: $species1") # The first iris species is: I. setosa
@printf("The 3rd iris species is: %s", species3) # The 3rd iris species is: I. virginica

s1 = :male
typeof(s1) # Symbol

println("")
for c in species2
    print(c, " - ")
end
# I - . -   - v - e - r - s - i - c - o - l - o - r -

search(species1, '.') # 2

email_pattern = r"(.+)@(.+)"
input = "john.doe@mit.edu"
ismatch(email_pattern, input) # true

split(species1) # 2-element Array{SubString{ASCIIString},1}:
                # "I."
                # "setosa"

temperature = [98.1, 98.6, 101.4] # 3-element Array{Float64,1}:
flu_status = [false, false, true] # 3-element Array{Bool,1}:
iris_species = [species1, species2, species3]
# 3-element Array{ASCIIString,1}:
#  "I. setosa"
#  "I. versicolor"
#  "I. virginica"

iris_species[3] # "I. virginica"
iris_species[2:3] # 2-element Array{ASCIIString,1}:
 # "I. versicolor"
 # "I. virginica"
iris_species[3] = "I. albicans"
iris_species
# 3-element Array{ASCIIString,1}:
#  "I. setosa"
#  "I. versicolor"
#  "I. albicans"
length(iris_species) # 3

arr1 = collect(1:5)
# 5-element Array{Int64,1}:
#  1
#  2
#  3
#  4
#  5

arr2 = Float64[] # 0-element Array{Float64,1}
push!(arr2, 5.1) # 1-element Array{Float64,1}: 5.1
push!(arr2, 3.5) # 2-element Array{Float64,1}: 5.1 3.5
println("")
show(arr2) # [5.1,3.5]
println("")

randn(5)
5-element Array{Float64,1}:
 -1.01086
 -0.154656
 -1.8495
  1.01637
  0.184996

for s in iris_species
    print(s, " / ")
end
# I. setosa / I. versicolor / I. albicans /
println("")

arr3 = [11, "Julia", 3.14, 'α']
# 4-element Array{Any,1}:
#  11
#  "Julia"
#  3.14
#  'α'

iris_species
# 3-element Array{ASCIIString,1}:
#  "I. setosa"
#  "I. versicolor"
#  "I. albicans"

push!(iris_species, "I. virginica")
# 4-element Array{ASCIIString,1}:
#  "I. setosa"
#  "I. versicolor"
#  "I. albicans"
#  "I. virginica"

pop!(iris_species)
# "I. virginica"

sort!(iris_species)
# 3-element Array{ASCIIString,1}:
#  "I. albicans"
#  "I. setosa"
#  "I. versicolor"

max_stars_repo_path: Section 1/1_3.jl

max_stars_repo_name: lytemar/Julia-for-Data-Science-Video

max_issues_repo_path: Section 1/1_3.jl

max_issues_repo_name: lytemar/Julia-for-Data-Science-Video

max_forks_repo_path: Section 1/1_3.jl

max_forks_repo_name: lytemar/Julia-for-Data-Science-Video


text:
struct marginal{Ω, T <: ApproxCDF}
  wv::Ω
  μ::Float64
  σ::Float64
  itp::T
end
function update_MarginalBuffer!(jp::JointPosterior{P}, f::Function) where P
    m = get!(() -> MarginalBuffer(length(jp.density)), jp.M.MarginalBuffers, f )
    @inbounds for (i,j) ∈ enumerate(jp.Θ)
        m.v[i] = f(j.Θ...)
    end
    Vandermonde!(m, jp.density, calc_μ_σ(m.v, jp.density)...)
end

function Vandermonde(m::MarginalBuffer, density::Vector{Float64})
    sortperm!(m.ind, m.v)
    cumulative_w = 0.0
    μ = 0.0
    EX² = 0.0
    k = 2
    @inbounds for (i, j) ∈ enumerate(m.ind)
        dᵢ = density[j]
        cumulative_w += dᵢ
        m.w[i] = cumulative_w
        vᵢ = m.v[j]
        μ += vᵢ*dᵢ
        vᵢ² = vᵢ^2
        EX² += vᵢ²*dᵢ
        m.V[k] = vᵢ
        m.V[k+1] = vᵢ²  * exp(-2)
        m.V[k+2] = vᵢ^3 * exp(-4)
        m.V[k+3] = vᵢ^4 * exp(-6)
        m.V[k+4] = vᵢ^5 * exp(-8)
        m.V[k+5] = vᵢ^6 * exp(-10)
        m.V[k+6] = vᵢ^7 * exp(-12)
        m.V[k+7] = vᵢ^8 * exp(-14)
        m.V[k+8] = vᵢ^9 * exp(-16)
        k += 10
    end
    m, μ, sqrt(EX² - μ^2)
end
function Vandermonde!(m::MarginalBuffer, density::Vector{T}, μ::Number, σ::Number) where T <:Number
    sortperm!(m.ind, m.v)
    cumulative_w = zero(T)
    k = 2
    @inbounds for (i, j) ∈ enumerate(m.ind)
        dᵢ = density[j]
        cumulative_w += dᵢ
        m.w[i] = cumulative_w
        vᵢ = (m.v[j] - μ)/σ
        vᵢ² = vᵢ^2
        vᵢ4 = vᵢ²^2
        m.V[k] = vᵢ
        m.V[k+1] = vᵢ²
        m.V[k+2] = vᵢ*vᵢ²
        m.V[k+3] = vᵢ4
        m.V[k+4] = vᵢ4*vᵢ
        m.V[k+5] = vᵢ4*vᵢ²
        m.V[k+6] = vᵢ4*vᵢ²*vᵢ
        m.V[k+7] = vᵢ4^2
        m.V[k+8] = vᵢ4^2*vᵢ
        k += 10
    end
    m, μ, σ
end
@inline update_MarginalBuffer!(jp::JointPosteriorRaw, f::Function) = update_MarginalBuffer!(jp.M, jp.grid.cache, jp.grid.density, f::Function)
@inline update_MarginalBuffer!(M::Model, cache::MatrixVecSVec, density::Vector, f::Function ) = update_MarginalBuffer!(M, cache.M, density, f::Function)
function update_MarginalBuffer!(M::Model, cache::Matrix, density::Vector, f::Function)
    m = get!(() -> MarginalBuffer(length(density)), M.MarginalBuffers, f )
    svf = SlidingVecFun(M.Θ, f, cache)
    @inbounds for i ∈ eachindex(m.v)
        m.v[i] = svf()
    end
    Vandermonde!(m, density, calc_μ_σ(m.v, density)...)
end

function calc_μ_σ(x::AbstractVector{T}, w::AbstractVector{T}) where T
  μ = dot(x, w)
  EX² = zero(T)
  @inbounds for i ∈ eachindex(x)
    EX² += x[i]^2 * w[i]
  end
  μ, √(EX² - μ^2)
end

function SlidingVecFun( Θ::ModelParam, f::Function, M::Matrix )
    g = function ()
        update!(Θ)
        f(Θ.Θ...)
    end
    svf = SlidingVecFun(g, Θ.v)
    set!(svf, M)
    svf
end

function weights_values(jp::JointPosterior{P}, f::Function) where {P}
    values = similar(jp.density)
    weights = similar(jp.density)
    @inbounds for (i,j) in enumerate(jp.Θ)
        values[i] = f(j.Θ...)
    end
    weights_values(copy(jp.density), values)
end
@inline weights_values(jp::JointPosteriorRaw, f::Function) = weights_values(jp.M.Θ, jp.grid.cache, copy(jp.grid.density), f::Function)
@inline weights_values(Θ::ModelParam, cache::MatrixVecSVec, density::Vector, f::Function ) = weights_values(Θ, cache.M, density, f::Function)
function weights_values(Θ::ModelParam, cache::Matrix, density::Vector, f::Function)
  values = similar(density) #May add dict to model that will allow for saving these.
  svf = SlidingVecFun(Θ, f, cache)
  @inbounds for i ∈ eachindex(values)
    values[i] = svf()
  end
  weights_values(density, values)
end

@inline marginal(jp::JointDist, f::Function, ::Type{Grid}) = marginal(jp::JointDist, f::Function)
function marginal(jp::JointDist, f::Function)
  wv = weights_values(jp, f)
  μ = dot(wv.weights, wv.values)
  σ = sqrt(dot(wv.weights, wv.values .^ 2) - μ^2)
  marginal(wv, μ, σ, Grid(wv))
end
function marginal(jp::JointDist, f::Function, ::Type{T}) where {T <: ContinuousUnivariateDistribution}
    marginal(update_MarginalBuffer!(jp, f)..., T)
end
function marginal(m::MarginalBuffer, μ::Real, σ::Real, ::Type{Normal})
    marginal(m, μ, σ, NestedPolyGLM(m, Normal(μ, σ)))
end
function marginal(jp::MarginalBuffer, μ::Real, σ::Real, ::Type{Gamma})
  throw("Currently unsupported.")
  marginal(m, μ, σ, NestedPolyGLM(m, Gamma((μ/σ)^2, σ^2/μ)))
end

function marginal(wv::weights_values, μ::Float64, σ::Float64, ::Type{Grid})
  marginal(wv, μ, σ, Grid(wv))
end

function Distributions.cdf(m::marginal, x)
  Distributions.cdf(m.itp, x)
end
function Distributions.pdf(m::marginal, x)
  Distributions.pdf(m.itp, x)
end
function Base.quantile(m::marginal, x)
  Base.quantile(m.itp, x)
end

function Base.show(io::IO, m::marginal)
  println("Marginal parameter")
  println("μ: ", m.μ)
  println("σ: ", m.σ)
  println("Quantiles: ", quantile.(m, [.025 .25 .5 .75 .975]))
end


text:
function empty_xcom_table()
    (
     E                                        = empty_col(uenergy ),
     CoherentScattering                       = empty_col(umassatt),
    IncoherentScattering                      = empty_col(umassatt),
    PhotoelectricAbsorption                   = empty_col(umassatt),
    PairProductionNuclearField                = empty_col(umassatt),
    PairProductionElectronField               = empty_col(umassatt),
    TotalAttenuation                          = empty_col(umassatt),
    TotalAttenuationWithoutCoherentScattering = empty_col(umassatt),
    )
end

const XCOMTable = typeof(empty_xcom_table())
struct XCOMData  <: DataSource
    element_tables::Vector{XCOMTable}
    compound_tables::Dict{Symbol, XCOMTable}
end

emptytable(::Type{XCOMData}) = empty_xcom_table()

const XCOM  = load(XCOMData, Zs=1:100, compounds=[:H2O],
                   dir=datapath("XCOM"))

function get_table(s::XCOMData, key::Symbol)::XCOMTable
    res = get(s.compound_tables, key, nothing)
    if res === nothing
        get_element_table(s, key)
    else
        res
    end
end

for P in columnnames(XCOMData)
    P == :E && continue
    @eval default_data_source(mat, E, ::$P) = XCOM
end


text:
using RandomizedPropertyTest, NOVAS, BenchmarkTools, Random, DataFrames

rng = MersenneTwister(0)

# Include un-exported c wrappers
include("wrapper.jl")

# Add custom isapprox to test tuples
import Base.isapprox
Base.isapprox(x::Tuple, y::Tuple; kws...) = isapprox(collect(x), collect(y); kws...)

# Setup record of benchmark results
results = DataFrame(; fname=String[], c_time=Float64[], julia_time=Float64[],
                    speedup=Float64[])

# Benchmark using curried f from RPT
function quickbench(f::Function, types)
    genargs = () -> RandomizedPropertyTest.generate(rng, types)
    return median(@benchmark $f($genargs()...)).time
end

# Stolen from RandomizedPropertyTest
macro testbench(args...)
    names = Symbol[]
    types = []
    nexpr = 10^3
    expr = args[1]
    vartypes = args[2:length(args)]
    for e in vartypes
        if e.args[1] isa Symbol
            newsymbols = Symbol[e.args[1]]
        elseif e.args[1].head == :tuple && all(x -> x isa Symbol, e.args[1].args)
            newsymbols = e.args[1].args
        end
        for symb in newsymbols
            push!(names, symb)
            push!(types, e.args[2])
        end
    end
    nametuple = Expr(:tuple, names...)
    typetuple = esc(Expr(:tuple, types...))
    exprstr = let io = IOBuffer()
        print(io, expr)
        seek(io, 0)
        read(io, String)
    end
    namestrs = [String(n) for n in names]
    # Build equivalence test
    NOVAS_expr = :(NOVAS.$(expr.args[1])($(expr.args[2:end]...)))
    test_expr = :($expr ≈ $NOVAS_expr)
    test_fexpr = esc(Expr(:(->), nametuple, test_expr))
    c_fexpr = esc(Expr(:(->), nametuple, expr))
    j_fexpr = esc(Expr(:(->), nametuple, NOVAS_expr))
    return quote
        test_result = RandomizedPropertyTest.do_quickcheck($test_fexpr, $exprstr, $namestrs,
                                                           $typetuple, $nexpr)
        if test_result && !haskey(ENV, "JULIA_DONT_BENCH")
            ctime = quickbench($c_fexpr, $typetuple)
            jtime = quickbench($j_fexpr, $typetuple)
            speedup = ctime / jtime
            push!(results, ($exprstr, ctime, jtime, speedup))
        end
        test_result
    end
end

# Custom types for common arguments
struct Position end
RandomizedPropertyTest.generate(rng::AbstractRNG, _::Type{Position}) = rand(Float64, 3)

struct Location end

function RandomizedPropertyTest.generate(rng::AbstractRNG, _::Type{Location})
    lat = rand(Cdouble) * 180 - 90
    lon = rand(Cdouble) * 180 - 90
    alt = rand(Cdouble) * 8000
    temp = rand(Cdouble) * 50
    pressure = rand(Cdouble) * 1000
    return NOVAS.OnSurface(lat, lon, alt, temp, pressure)
end


text:
using Test, bad

@info "testing probability mass function..."
@time include("test_pmf.jl")

@info "testing priors..."
@time include("test_Beta.jl")

@info "testing scores..."
@time include("test_scores.jl")

@info "testing point estimators..."
@time include("test_point-estimators.jl")

@info "testing p-values..."
@time include("test_p-values.jl")

@info "testing interval estimators..."
@time include("test_intervals.jl")

@info "testing one stage designs..."
@time include("test_one_stage.jl")

@info "testing vs simons designs..."
@time include("test_simons.jl")

@info "testing minimmax sample size..."
@time include("test_MiniMaxSampleSize.jl")

@info "testing expected utility..."
@time include("test_ExpectedUtility.jl")

@info "testing design adaptations..."
@time include("test_adapt.jl")

@info "... done - stay vigilent!"