text:
###################################################
## scenarios.jl
##      Compute routing solutions for a single school
## Authors: Arthur Delarue, Sébastien Martin, 2018
###################################################

"""
    The main euclidian traveltime function.
    Select the time that is the closest to `time` and return the euclidian travel-time
    for the corresponding speed.
    o:: Origin
    d:: Destination
"""
 # using JuMP, Gurobi
using JuMP, Cbc
include("problem.jl")
include("load.jl")
#Sinclude("SchoolBusRouting.jl")
import Statistics: mean
using Base:GC
using ProgressMeter
using Random

data = loadSyntheticBenchmark("data/input/CSCB01/Schools.txt",
                                  "data/input/CSCB01/Stops.txt")

# data = loadSyntheticBenchmark("../data/input/CSCB01/Schools.txt",
                                  # "../data/input/CSCB01/Stops.txt")
data.stops[1]

function traveltime(data::SchoolBusData, o::Point, d::Point)
    if data.params.metric == MANHATTAN
        return manhattandistance(o,d) / data.params.velocity
    else
        return euclideandistance(o,d) / data.params.velocity
    end
end

"""
    Manhattan Distance between two points, in meters
"""
function manhattandistance(o::Point, d::Point)
    return abs(o.x - d.x) + abs(o.y - d.y)
end

"""
    Euclidean distance between two points, in meters
"""
function euclideandistance(o::Point, d::Point)
    return sqrt((o.x - d.x) ^ 2 + (o.y - d.y) ^ 2)
end

"""
School, Yard, Stop  2^3
"""
traveltime(data::SchoolBusData, o::Yard,   d:: Stop)   = traveltime(data, o.position, d.position)
traveltime(data::SchoolBusData, o::School, d:: Stop)   = traveltime(data, o.position, d.position)
traveltime(data::SchoolBusData, o::Stop,   d:: Stop)   = traveltime(data, o.position, d.position)
traveltime(data::SchoolBusData, o::Stop,   d:: School) = traveltime(data, o.position, d.position)
traveltime(data::SchoolBusData, o::School, d:: Yard)   = traveltime(data, o.position, d.position)
traveltime(data::SchoolBusData, o::Yard,   d:: School) = traveltime(data, o.position, d.position)
traveltime(data::SchoolBusData, o::Stop,   d:: Yard)   = traveltime(data, o.position, d.position)

"""
    Get number of students at the stop
"""
nStudents(data::SchoolBusData, stop::Stop) = stop.nStudents
nStudents(data::SchoolBusData, school::Int, stop::Int) = data.stops[school][stop].nStudents
#
# bus_data = loadSchoolsReduced("../data/input/CSCB01/Schools.txt")
#
# bus_data[1].position
data.schools[1].position
"""
    Get time that bus must remain at stop
"""
function stopTime(data::SchoolBusData, stop::Stop)
    return data.params.constant_stop_time + data.params.stop_time_per_student *
    stop.nStudents
end

# stop_data = CSV.read("../data/input/CSCB01/Stops.txt")
# stops = stop_data[1,:]
# stops = data.stops[1][1]
# stopTime(data, stops)
"""
    Get the maximum allowed travel time for a given bus stop
"""
function maxTravelTime(data::SchoolBusData, stop::Stop)
    return data.params.max_time_on_bus
end

maxTravelTime(data, data.stops[2][3]) # constant time 2700s
"""
    The current state of the greedy algorithm
    Contains:
        - the current route
        - an occupancy object that tells us how many students are currently on the bus
        - for each stop, the amount of extra time we can leave the students on the bus
        - for each stop, the current time spent on the bus
        - the total service time of the route (including time of first stop)
"""
struct GreedyState
    "the current route"
    route::Route
    "the current occupancy"
    nStudents::Int
    "amount of extra time students at each stop can spend on bus"
    slackTimes::Vector{Float64}
    "amount of time students at each stop are already spending on bus"
    stopTimes::Vector{Float64}
    "total time of current route"
    routeTime::Float64
end

"""
    Greedy single-school routing heuristic
    Args:
        - the data
        - the ID of the school
        - the maximum time allowed on a route
        - the penalty (in seconds) for adding D2D students to routes
    Returns:
        - all routes for the school as a Vector{Route}
"""

maxRouteTime = 2700.0

Random.seed!(1234)

function greedy(data::SchoolBusData, schoolID::Int, maxRouteTime::Float64)
    routes = Route[]
    availableStops = trues(length(data.stops[schoolID]))
    while sum(availableStops) > 0
        # randomly select a starting stop
        startingStopID = rand((1:length(availableStops))[availableStops])
        # create initial route (select initial bus)
        currentState = initialRoute(data, schoolID, startingStopID,
                                    length(routes)+1)
        availableStops[startingStopID] = false
        while true
            bestStopID = 0
            bestInsertId = -1
            bestTimeDiff = Inf
            for stopID in collect(1:length(availableStops))[availableStops]
                # check the students number when adding a stop
                if (data.stops[schoolID][stopID].nStudents + currentState.nStudents <=
                    data.params.bus_capacity)

                    insertId, timeDiff = bestInsertion(data, schoolID, stopID,
                                                       currentState, maxRouteTime)
                    if timeDiff < bestTimeDiff
                        bestTimeDiff = timeDiff
                        bestStopID = stopID
                        bestInsertId = insertId
                    end
                end
            end

            if bestTimeDiff < Inf
                currentState = buildRoute(data, currentState, schoolID,
                                          bestStopID, bestInsertId)
                availableStops[bestStopID] = false
            else
                push!(routes, currentState.route)
                break
            end
        end
    end
    return routes
end

schoolID = 1
stopID = 1

# maxRouteTime = 2700.0
# all_routes = greedy(data, 1, maxRouteTime)
# all_routes[1]

"""
    Get total travel time difference when adding stop to the route, together with
        a "cost" of how good this insertion is.
    The cost is equal to the total amount of time by which this insertion
        increases the length of the route, not including stop time
    Args:
        - data          : the data object, assumes it has the stops computed
        - schoolID      : the school we are routing
        - newStopID     : the stop we are trying to insert
        - cs            : the current GreedyState
        - maxRouteTime  : the maximal allowed length of the route
    Returns:
        - insertID      : where we would optimally insert it
        - bestTimeDiff  : the optimal time difference imposed on the route
   Comment:
       when adding a new stop, if we should add in the 1st,
       or in the middle or the last stop before school.
"""
function bestInsertion(data::SchoolBusData,
                       schoolID::Int,
                       newStopID::Int,
                       cs::GreedyState,
                       maxRouteTime::Float64)
    newStop      = data.stops[schoolID][newStopID]
    school       = data.schools[schoolID]
    bestTimeDiff = Inf
    insertId     = -1
    # before first stop
    timeDiff = traveltime(data, newStop, data.stops[schoolID][cs.route.stops[1]])
    totTime = timeDiff + cs.routeTime + stopTime(data, newStop)
    if (totTime - stopTime(data, newStop) <= maxTravelTime(data, newStop) &&
                totTime <= maxRouteTime)
        bestTimeDiff = timeDiff
        insertId = 0
    end
    # between stops
    # loop over all the stops in the middle
    # local search to get the best feasible route.

    for i = 1:(length(cs.route.stops)-1)
        timeToNextStop = traveltime(data, newStop, data.stops[schoolID][cs.route.stops[i+1]])
        timeDiff = traveltime(data, data.stops[schoolID][cs.route.stops[i]], newStop) +
                   timeToNextStop -
                   traveltime(data, data.stops[schoolID][cs.route.stops[i]],
                      data.stops[schoolID][cs.route.stops[i+1]])
        if timeDiff < bestTimeDiff
            totTime = timeDiff + cs.routeTime + stopTime(data, newStop)
            # check feasibility for all stops preceding this one
            isFeasible = (timeDiff + stopTime(data, newStop) <= cs.slackTimes[i])
            # check feasibility for the potential new stop
            isFeasible = isFeasible && (cs.stopTimes[i+1] + timeToNextStop +
                                        stopTime(data,
                                             data.stops[schoolID][cs.route.stops[i+1]]) <=
                                        maxTravelTime(data, newStop))
            isFeasible = isFeasible && totTime <= maxRouteTime
            if isFeasible
                bestTimeDiff = timeDiff
                insertId = i
            end
        end
    end
    # after last stop
    timeDiff =  traveltime(data, data.stops[schoolID][cs.route.stops[end]], newStop) +
                traveltime(data, newStop, school) -
                traveltime(data, data.stops[schoolID][cs.route.stops[end]], school)
    if timeDiff < bestTimeDiff
        totTime = timeDiff + cs.routeTime + stopTime(data, newStop)
        isFeasible = timeDiff + stopTime(data, newStop) <= cs.slackTimes[end] &&
                        (traveltime(data, newStop, school) <= maxTravelTime(data, newStop))
        isFeasible = isFeasible && totTime <= maxRouteTime
        if isFeasible
            bestTimeDiff = timeDiff
            insertId = length(cs.route.stops)
        end
    end

    return insertId, bestTimeDiff
end

"""
    Given new stop ID and current route/GreedyState, update route
    Args:
        - data          : the data object, assumes it has the stops computed
        - cs            : the current GreedyState
        - schoolID      : the school we are routing
        - newStopID     : the stop we are inserting
        - insertID      : where we insert it
    Returns:
        - a new GreedyState with the updated route
"""
function buildRoute(data::SchoolBusData,
                    cs::GreedyState,
                    schoolID::Int,
                    newStopID::Int,
                    insertID::Int)
    newStop = data.stops[schoolID][newStopID]
    school = data.schools[schoolID]
    if insertID == 0 # at the beginning
        nextStop = data.stops[schoolID][cs.route.stops[1]]
        "amount of time students at each stop are already spending on bus"
        newStopTimeOnBus = traveltime(data, newStop, nextStop) +
                           cs.stopTimes[1] + stopTime(data,nextStop)
        timeDiff = 0.
    elseif insertID == length(cs.route.stops) # at the end
        previousStop = data.stops[schoolID][cs.route.stops[end]] # last stop on route
        newStopTimeOnBus = traveltime(data, newStop, school) # time from new stop to school
        timeDiff = newStopTimeOnBus +
                   traveltime(data, previousStop, newStop) -
                   traveltime(data, previousStop, school)
    else # in the middle
        previousStop = data.stops[schoolID][cs.route.stops[insertID]]
        nextStop = data.stops[schoolID][cs.route.stops[insertID+1]]
        timeToNextStop = traveltime(data, newStop, nextStop)
        timeDiff = traveltime(data, previousStop, newStop) +
                   timeToNextStop -
                   traveltime(data, previousStop, nextStop)
        newStopTimeOnBus = timeToNextStop + cs.stopTimes[insertID+1] + stopTime(data, nextStop)
    end

    # it inclueds all stops' stoptime and slackTimes
    # if insertId == 0, then 1:insertID will aalwyas return 0.
    newStopTimes = vcat(cs.stopTimes[1:insertID] .+ stopTime(data, newStop) .+ timeDiff,
                        [newStopTimeOnBus],
                        cs.stopTimes[(insertID+1):end])
    newStops = vcat(cs.route.stops[1:insertID], [newStopID], cs.route.stops[(insertID+1):end])
    newSlackTimes = vcat(cs.slackTimes[1:insertID] .- stopTime(data, newStop) .- timeDiff,
                         [maxTravelTime(data, newStop)] .- newStopTimeOnBus,
                         cs.slackTimes[(insertID+1):end])
    # fix slack time property
    for i=eachindex(newSlackTimes)
        if i > 1
            newSlackTimes[i] = min(newSlackTimes[i-1], newSlackTimes[i])
        end
    end
    routeTime = newStopTimes[1] + stopTime(data, data.stops[schoolID][newStops[1]])
    return GreedyState(Route(cs.route.id, newStops),
                       cs.nStudents + newStop.nStudents,
                       newSlackTimes, newStopTimes, routeTime)
end

# cs = initial_route_test
# println(cs)

# newStopID = 2
# insertID = 0
# schoolID = 1
# route_example = buildRoute(data, cs,
#  schoolID, newStopID, insertID )
# println(route_example)
# cs = route_example
#
# newStopID = 2
# insertID = 1
# schoolID = 1
#
# route_example1 = buildRoute(data, route_example,
#  schoolID, newStopID, insertID )
#  println(route_example1)

"""
    Create initial route for greedy
    Args:
        - data          : the data object, assumes it has the stops computed
        - schoolID      : the school we are routing
        - stopID        : the stop we start with
        - routeID       : the ID of the route we are creating
    Returns a GreedyState object
"""
schoolID = 1
stopID = 1
function initialRoute(data::SchoolBusData,
                      schoolID::Int,
                      stopID::Int,
                      routeID::Int)
    timeOnBus = traveltime(data, data.stops[schoolID][stopID], data.schools[schoolID])
    nStudents = data.stops[schoolID][stopID].nStudents
    slackTimes = [maxTravelTime(data, data.stops[schoolID][stopID]) - timeOnBus]
    "amount of time students at each stop are already spending on bus"
    stopTimes = [timeOnBus]
        routeTime = timeOnBus + stopTime(data, data.stops[schoolID][stopID])
    return GreedyState(Route(routeID, [stopID]), nStudents, slackTimes, stopTimes, routeTime)
end

initial_route_test = initialRoute(data,data.stops[1][1].id, data.stops[1][1].id, 1)
println(initial_route_test)
"""
    Returns the sum of each student's travel times on this route
"""

function sumIndividualTravelTimes(data::SchoolBusData, schoolID::Int, r::Route)
    allStops = data.stops[schoolID]
    numStops = length(r.stops)
    t = traveltime(data, allStops[r.stops[numStops]], data.schools[schoolID]) * numStops
    while numStops > 1
        t += (traveltime(data, allStops[r.stops[numStops-1]], allStops[r.stops[numStops]]) +
              stopTime(data, allStops[r.stops[numStops]])) * (numStops - 1)
        numStops -= 1
    end
    return t
end
#
# r = all_routes[2]
# sumIndividualTravelTimes(data, 1, r)
function sumIndividualTravelTimes(data::SchoolBusData, schoolID::Int,
                                  routes::Vector{Route})
    return sum(sumIndividualTravelTimes(data, schoolID, route) for route in routes)
end

# sumIndividualTravelTimes(data, 1, all_routes)

"""
    Returns the total travel time from the time the first student enters the bus to the time of pickup/dropoff in school
"""
function serviceTime(data::SchoolBusData, schoolID::Int, stoplist::Vector{Int})
    allStops = data.stops[schoolID]
    numStops = length(stoplist)
    t = traveltime(data, allStops[stoplist[numStops]], data.schools[schoolID])
    while numStops > 1
        t += (traveltime(data, allStops[stoplist[numStops-1]], allStops[stoplist[numStops]]) +
              stopTime(data, allStops[stoplist[numStops]]))
        numStops -= 1
    end
    t += stopTime(data, allStops[stoplist[1]])
    return t
end
serviceTime(data::SchoolBusData, schoolID::Int, r::Route) = serviceTime(data, schoolID, r.stops)

# serviceTime(data, 1, r)

"""
    Simple route representation : just list of IDs and associated cost
"""
struct FeasibleRoute
    "The list of stop ids"
    stopIds::Vector{Int}
    "The cost associated with the route"
    cost::Float64
end

function FeasibleRoute(data::SchoolBusData, schoolID::Int, r::Route)
    return FeasibleRoute(r.stops, sumIndividualTravelTimes(data, schoolID, r))
end

# a_feasible_route = FeasibleRoute(data, 1, r)

"""
    Stores a list of routes in a way that makes column generation easy
"""
struct FeasibleRouteSet
    "The list of FeasibleRoutes available"
    list::Vector{FeasibleRoute}
    "The set of stopIds list (one for each route)"
    set::Set{Vector{Int}}
    "For each stop, the index of the routes in the list that go through this stop"
    atStop::Vector{Vector{Int}}

    FeasibleRouteSet(data::SchoolBusData, schoolID::Int) =
        new(FeasibleRoute[], Set{Vector{Int}}(), [Int[] for i = 1:length(data.stops[schoolID])])
end

# a_feasible_set = FeasibleRouteSet(data, 1)
"""
    Generate N random greedy solutions, and combines them smartly to get the best
    possible solution.
"""
function greedyCombined(data::SchoolBusData,
                        schoolID::Int,
                        N::Int,
                        maxRouteTimeLower::Float64,
                        maxRouteTimeUpper::Float64,
                        λ::Float64;
                        # ,
                        # env::Gurobi.Env=Gurobi.Env();
                        args...)
    routeList = generateRoutes(data, schoolID, N, maxRouteTimeLower, maxRouteTimeUpper)
    routes = FeasibleRouteSet(data, schoolID)
    addRoute!(routes, routeList)
    # selectedRoutes = bestRoutes(data, schoolID, routes, λ, env; args...)
    selectedRoutes = bestRoutes(data, schoolID, routes, λ;args... )

    return buildSolution(data, schoolID, routes, selectedRoutes)
end

# greedyCombined(data, schoolID, N, maxRouteTimeLower, maxRouteTimeUpper, λ  )

function greedyCombined(data::SchoolBusData,
                        schoolID::Int,
                        startRoutes::Vector{Route},
                        N::Int,
                        maxRouteTimeLower::Float64,
                        maxRouteTimeUpper::Float64,
                        λ::Float64;
                        # change it back if gurobi is available
                        # ,
                        # env::Gurobi.Env=Gurobi.Env();
                        args...
                        )
    routeList = generateRoutes(data, schoolID, N, maxRouteTimeLower, maxRouteTimeUpper)
    routes = FeasibleRouteSet(data, schoolID)
    addRoute!(routes, routeList)
    addRoute!(routes, collect(FeasibleRoute(data, schoolID, r) for r in startRoutes))
    # selectedRoutes = bestRoutes(data, schoolID, routes, λ, env; args...)
    selectedRoutes = bestRoutes(data, schoolID, routes, λ;args... )

    return buildSolution(data, schoolID, routes, selectedRoutes)
end

# startRoutes =  greedy(data, schoolID, Inf)

# greedyCombined(data, schoolID, startRoutes, N, maxRouteTimeLower, maxRouteTimeUpper, λ  )

"""
    Same as greedy combined, but iterates it to keep improving the best solution
"""
function greedyCombinedIterated(data::SchoolBusData,
                                schoolID::Int,
                                maxRouteTimeLower::Float64,
                                maxRouteTimeUpper::Float64,
                                nGreedy::Int,
                                nIteration::Int,
                                λ::Float64;
                                verbose::Bool=false,
                                # ,
                                args...
                                )
    # env = Gurobi.Env()
    routes = greedy(data, schoolID, Inf)
    verbose && @printf("Iteration 0: %d buses, %2.fs\n", length(routes),
                        sumIndividualTravelTimes(data, schoolID, routes))
    for i=1:nIteration
        routes = greedyCombined(data, schoolID, routes, nGreedy,
                                maxRouteTimeLower, maxRouteTimeUpper, λ;
                                # , env;
                                args...
                                )

        verbose && @printf("Iteration %d: %d buses, %2.fs\n", i, length(routes),
                           sumIndividualTravelTimes(data, schoolID, routes))
    end
    # gc()
    GC.gc()
    return routes
 end

a = greedyCombinedIterated(data, schoolID, 1700.0, 2700.0,10, 500, 100.0; verbose = false, loglevel = 1,seconds = 10, threads = 8, allowableGap = 0.5 )
"""
    Add feasible route to a feasible set
"""
function addRoute!(routes::FeasibleRouteSet, newRoute::FeasibleRoute)
    if ! (newRoute.stopIds in routes.set)
        push!(routes.set, newRoute.stopIds)
        push!(routes.list, newRoute)
        newRouteId = length(routes.list)
        for stopId in newRoute.stopIds
            push!(routes.atStop[stopId], newRouteId)
        end
    end
end
# r = all_routes[2]
#
# newRoute = a_feasible_route
#
# routes = a_feasible_set
# routes.atStop[20:60]
#
# routeList = generateRoutes(data, schoolID, N, maxRouteTimeLower, maxRouteTimeUpper)
# routes_example = FeasibleRouteSet(data, schoolID)
# routes.atStop[20:60]

"""
    Add list of feasible routes to a feasible set
"""
function addRoute!(routes::FeasibleRouteSet, newRoutes::Vector{FeasibleRoute})
    for r in newRoutes
        addRoute!(routes, r)
    end
end

"""
    Generate N sets of Routes using the greedy heuristic.
"""
#
# maxRouteTimeLower = 2400.0
# maxRouteTimeUpper = 3600.0
# N = 3
function generateRoutes(data::SchoolBusData,
                        schoolID::Int,
                        N::Int,
                        maxRouteTimeLower::Float64=Inf,
                        maxRouteTimeUpper::Float64=Inf)
    routes = FeasibleRoute[]
    for i in 1:N
        if maxRouteTimeLower < Inf
            maxTime = (maxRouteTimeUpper - maxRouteTimeLower) * rand() + maxRouteTimeLower
        else
            maxTime = Inf
        end
        singleRoutes = greedy(data, schoolID, maxTime)
        append!(routes, [FeasibleRoute(data, schoolID,route) for route in singleRoutes])
    end
    return routes
end

# feasible_set = generateRoutes(data, 1, 3, maxRouteTimeLower, maxRouteTimeUpper )

"""
    Solves the routing problem given a set of routes. (MIP)
"""

routeList = generateRoutes(data, schoolID, 3, 1700.0, 2700.0)
routes = FeasibleRouteSet(data, schoolID)
addRoute!(routes, routeList)
function bestRoutes(data::SchoolBusData, schoolID::Int, routes::FeasibleRouteSet,
                    λ::Float64;
                    # , env::Gurobi.Env;
                     args...
                    )
    # model = Model(solver=GurobiSolver(env, Threads=getthreads(); args...))

    # model =  Model(with_optimizer(Gurobi.Optimizer; args...))
    model =  Model(with_optimizer(Cbc.Optimizer; args...))
    # The binaries, whether we choose the routes
    @variable(model, r[k in 1:length(routes.list)], Bin)
    # Minimize the number of buses first, then the travel time
    @objective(model, Min, sum(r[k] * (λ + routes.list[k].cost) for k in 1:length(routes.list)))
    # At least one route per stop
    @constraint(model, stopServed[i in 1:length(data.stops[schoolID])],
        sum(r[k] for k in routes.atStop[i]) >= 1)
    # solve(model)
     JuMP.optimize!(model)
    # selectedRoutes = [k for k in 1:length(routes.list) if getvalue(r[k]) >= 0.5]
    selectedRoutes = [k for k in 1:length(routes.list) if JuMP.value(r[k]) == 1]
    return selectedRoutes
end

selectedRoutes = bestRoutes(data, 1, routes, 100.0; seconds = 1)
"""
    Given a covering list of FeasibleRoutes, create correct route object
"""
function buildSolution(data::SchoolBusData, schoolID::Int,
                       routeSet::FeasibleRouteSet, selectedRoutes::Vector{Int})
    # retrieve the route stopid and cost for the covering model route
    routes = copy(routeSet.list[selectedRoutes])
    routesAtStop = Vector{Int}[Int[] for s in data.stops[schoolID]]
    for (routeId,r) in enumerate(routes)
        for stopId in r.stopIds
            push!(routesAtStop[stopId], routeId)
        end
    end
    # which trip covered the stop[1:END]
    for (stopId,intersectingRoutes) in enumerate(routesAtStop)
        if length(intersectingRoutes) > 1
            newRoutes = splitRoutes(data, schoolID, routes[intersectingRoutes], stopId)
            for (i, routeId) in enumerate(intersectingRoutes)
                routes[routeId] = newRoutes[i]
            end
        end
    end
    # remove the feasieble route with 0 stop because the previous deleting stop
    return [Route(i, fr.stopIds) for (i, fr) in enumerate(routes) if length(fr.stopIds) > 0]
end
#
# routeSet = routes
# routeSet.list[21]
# routes[21]
"""
    When several routes intersect in stopId, choose the best one to serve the stopId
    and remove the others
"""
function splitRoutes(data::SchoolBusData, schoolID::Int,
                     routes::Vector{FeasibleRoute}, stopId::Int)
    costs = collect(deletionCost(data, schoolID, route, stopId) for route in routes)
    # selectedRoute = indmax(costs)
    selectedRoute = argmax(costs) # max(-20, -1), lower means by reomoveing the stop, it saves more time -20
    for (routeId, route) in enumerate(routes)
        if selectedRoute != routeId
            newRouteIds = collect(s for s in route.stopIds if s != stopId)
            routes[routeId] = FeasibleRoute(newRouteIds, route.cost + costs[routeId])
        end
    end
    return routes
end
#
# stopId = 1
# intersectingRoutes = routesAtStop[stopId]
# routes =  routeSet[intersectingRoutes]
# # route = routes[intersectingRoutes][1]
#
# newRoutes =  splitRoutes(data, 1, routes, stopId)

"""
    Cost of removing a stop from a route (usually negative)
    The cost is just the difference between the old travel time and the newStop
"""

function deletionCost(data::SchoolBusData, schoolID::Int,
                      route::FeasibleRoute, stopId::Int)
    # stop = findfirst(route.stopIds, stopId)
    stop = findfirst(isequal(stopId), route.stopIds)
    length(route.stopIds) <= 1 && return (-Inf)
    stops = data.stops[schoolID]
    school = data.schools[schoolID]
    # remove travel time of that particular stop
    cost = 0.
    for id = stop:(length(route.stopIds)-1)
        # cost -= traveltime(data, stops[route.stopIds[id]], stops[route.stopIds[id+1]])
        # cost -= stopTime(data, stops[route.stopIds[id+1]])
        global cost -= traveltime(data, stops[route.stopIds[id]], stops[route.stopIds[id+1]])
        global cost -= stopTime(data, stops[route.stopIds[id+1]])
    end
    cost -= traveltime(data, stops[route.stopIds[end]], school)
    # remove travel time effect on other stops
    if stop == length(route.stopIds)
        global cost -= (traveltime(data, stops[route.stopIds[stop-1]], stops[stopId]) +
                 traveltime(data, stops[stopId], school) +
                 stopTime(data, stops[stopId]) -
                 traveltime(data, stops[route.stopIds[stop-1]], school)) *
                (length(route.stopIds) - 1)
    elseif stop > 1
        global cost -= (traveltime(data, stops[route.stopIds[stop-1]], stops[stopId]) +
                 traveltime(data, stops[stopId], stops[route.stopIds[stop+1]]) +
                 stopTime(data, stops[stopId]) -
                 traveltime(data, stops[route.stopIds[stop-1]], stops[route.stopIds[stop+1]]))*
                (stop - 1)
    end
    return cost
end

# deletionCost(data, 1, routes[1], 1 )

"""
    Contains the parameters that were used to compute a particular scenario
"""
struct ScenarioParameters
    "Maximum time of routes - lower end of interval"
    maxRouteTimeLower::Float64
    "Maximum time of routes - upper end of interval"
    maxRouteTimeUpper::Float64
    "Number of greedy solutions optimized over"
    nGreedy::Int
    "Tradeoff parameter in optimization"
    λ::Float64
    "Number of iterations for greedy combined iterated"
    nIterations::Int
end
"""
    Constructor for ScenarioParameters, with keyword parameters for readability
"""
function ScenarioParameters(;maxRouteTimeLower=Inf,
                            maxRouteTimeUpper=Inf,
                            nGreedy::Int=10,
                            λ=5e3,
                            nIterations::Int=10)
    return ScenarioParameters(maxRouteTimeLower, maxRouteTimeUpper, nGreedy, λ, nIterations)
end

λs = [1e2, 5e2, 1e3, 2e3, 5e3];

maxtime = 2700.0
nGreedy  = 10
params = [ScenarioParameters(maxRouteTimeLower=maxtime-2000,
                                 maxRouteTimeUpper=maxtime+2000,
                                 nGreedy=nGreedy, λ=λ,
                                 nIterations=10) for λ in λs]
scenariolist = computescenarios(data, params;
                            ## Gurobi arguments
                                 # OutputFlag=0,
                                 # MIPGap=ifelse(maxtime < 4000, 1e-4, 0.05),
                                 # TimeLimit=ifelse(maxtime < 4000, 90, 30)
                            ## Cbc arguments
                                 logLevel=1,
                                 allowableGap=ifelse(maxtime < 4000, 1e-4, 0.05),
                                 seconds=ifelse(maxtime < 4000, 90, 30),
                                 threads = 8
                                 );
"""
    Compute scenario
"""
function getscenario(data, scenarioinfo; args...)
    school, scenarioid, params = scenarioinfo
    routes = greedyCombinedIterated(data, school.id,
                                    params.maxRouteTimeLower, params.maxRouteTimeUpper,
                                    params.nGreedy, params.nIterations, params.λ;
                                    verbose=false, args...)
    return Scenario(school.id, scenarioid, collect(eachindex(routes))), routes
end

"""
    Compute multiple scenarios
    tocompute:
30-element Array{Tuple{School,Int64,ScenarioParameters},1}:
 (School 5 - 200005, 4, ScenarioParameters(700.0, 4700.0, 10, 2000.0, 10))
 (School 4 - 200004, 3, ScenarioParameters(700.0, 4700.0, 10, 1000.0, 10))
 (School 4 - 200004, 2, ScenarioParameters(700.0, 4700.0, 10, 500.0, 10))
 (School 5 - 200005, 2, ScenarioParameters(700.0, 4700.0, 10, 500.0, 10))
"""
function computescenarios(data, params; args...)
    tocompute = shuffle!(vec([(school,paramid,param) for school in data.schools, (paramid,param) in enumerate(params)]))
    results = Tuple{Scenario,Vector{Route}}[]
    @showprogress for scenarioinfo in tocompute
        push!(results, getscenario(data, scenarioinfo; args...))
    end
    return results
end
function computescenariosparallel(data, params; args...)
    tocompute = shuffle!(vec([(school,paramid,param) for school in data.schools, (paramid,param) in enumerate(params)]))
    results = Tuple{Scenario,Vector{Route}}[]
    results = pmap(x->getscenario(data,x;args...), tocompute)
    return results
end

"""
    Put the scenarios together
"""
function loadroutingscenarios!(data, scenariolist)
    scenarios = [Scenario[] for school in data.schools]
    routes = [Route[] for school in data.schools] # each school has a set of routes
    ids = vec([(scenario.school, scenario.id) for (scenario, routelist) in scenariolist])

    for k in sortperm(ids)
        (scenario, routelist) = scenariolist[k]
        for (i,route) in enumerate(routelist)
            idx = findfirst(x -> x.stops == route.stops, routes[scenario.school])
            # if idx == 0 # need to add new route and update its id
            if idx == nothing # need to add new route and update its id
                scenario.routeIDs[i] = length(routes[scenario.school]) + 1
                push!(routes[scenario.school],
                      Route(length(routes[scenario.school]) + 1, route.stops))
            else # need to update reference in the scenario
                scenario.routeIDs[i] = idx
            end
        end
        push!(scenarios[scenario.school], scenario)
    end
    data.scenarios = scenarios
    data.routes = routes
    data.withRoutingScenarios = true
    return data
end

k = 11
i = 1
route = Route(1, [67])
data = loadroutingscenarios!(data, scenariolist);

data.scenarios[1]
data.routes[1][9]


text:
# LICENSE:
#-------------------------------------------------------------------
# The code here is essentially copy-pasted from the MIT licensed
#
# https://github.com/mclements/Splines2.jl
#
# Only keeping it here temporarily until Splines2.jl is officially
# registered.
#-------------------------------------------------------------------

# utility functions
zeroArray(a::AbstractArray{T}) where {T<:Real} = OffsetArray(a, (size(a) .* 0) .- 1)

# type definitions
abstract type AbstractSplineBasis{T<:Real} end
# Note: we have used OffsetArray for converting from C code
mutable struct SplineBasis{T<:Real} <: AbstractSplineBasis{T}
    order::Int # order of the spline
    nknots::Int # number of knots
    ncoef::Int # number of coefficients
    knots::OffsetArray{T,1} # knot vector
end
mutable struct BSplineBasis{T<:Real} <: AbstractSplineBasis{T}
    spline_basis::SplineBasis{T}
    boundary_knots::Tuple{T,T}
    interior_knots::Union{Array{T,1},Nothing}
    intercept::Bool
    df::Int
end
mutable struct NSplineBasis{T<:Real} <: AbstractSplineBasis{T}
    b_spline_basis::BSplineBasis{T}
    qmat::Array{T,2}
    tl0::Array{T,1}
    tl1::Array{T,1}
    tr0::Array{T,1}
    tr1::Array{T,1}
end

# constructors
function SplineBasis(knots::Array{T,1}, order::Int=4) where {T<:Real}
    return SplineBasis(order, length(knots), length(knots) - order, zeroArray(knots))
end
function BSplineBasis(
    boundary_knots::Tuple{T,T},
    interior_knots::Union{Array{T,1},Nothing}=nothing,
    order::Int=4,
    intercept::Bool=false,
) where {T<:Real}
    l_interior_knots = interior_knots === nothing ? 0 : length(interior_knots)
    df = Int(intercept) + order - 1 + l_interior_knots
    nknots = l_interior_knots + 2 * order
    ncoef = nknots - order
    knots = zeros(T, nknots)
    for i in 1:order
        knots[i] = boundary_knots[1]
        knots[nknots - i + 1] = boundary_knots[2]
    end
    if (l_interior_knots > 0)
        for i in 1:l_interior_knots
            knots[i + order] = interior_knots[i]
        end
    end
    return BSplineBasis(
        SplineBasis(knots, order), boundary_knots, interior_knots, intercept, df
    )
end
function NSplineBasis(
    boundary_knots::Tuple{T,T},
    interior_knots::Union{Array{T,1},Nothing}=nothing,
    order::Int=4,
    intercept::Bool=false,
) where {T<:Real}
    bs = BSplineBasis(boundary_knots, interior_knots, order, intercept)
    co = basis(bs, [bs.boundary_knots[1], bs.boundary_knots[2]], 2)
    qmat_ = LinearAlgebra.qr(transpose(co)).Q
    qmat = collect(transpose(qmat_[:, 3:size(qmat_, 2)]))
    tl0 = qmat * basis(bs, boundary_knots[1])
    tl1 = qmat * basis(bs, boundary_knots[1], 1)
    tr0 = qmat * basis(bs, boundary_knots[2])
    tr1 = qmat * basis(bs, boundary_knots[2], 1)
    return NSplineBasis(bs, qmat, tl0, tl1, tr0, tr1)
end

function basis(bs::SplineBasis{T}, x::T, ders::Int=0) where {T<:Real}
    function find_interval(x::T)
        k = bs.order - 1
        n = bs.nknots - k - 1
        l = k
        while (x < bs.knots[l] && l != k)
            l = l - 1
        end
        l = l + 1
        while (x >= bs.knots[l] && l != n)
            l = l + 1
        end
        return l - 1
    end
    t = bs.knots
    k = bs.order - 1
    m = ders
    hh = bs.order
    ell = find_interval(x)
    result = zeroArray(zeros(T, 2 * k + 2))
    one = T(1)
    zero = T(0)
    result[0] = one
    for j in 1:(k - m)
        for n in 0:(j - 1)
            result[hh + n] = result[n]
        end
        result[0] = zero
        for n in 1:j
            ind = ell + n
            xb = t[ind]
            xa = t[ind - j]
            if (xb == xa)
                result[n] = zero
                continue
            end
            w = result[hh + n - 1] / (xb - xa)
            result[n - 1] = result[n - 1] + w * (xb - x)
            result[n] = w * (x - xa)
        end
    end
    for j in (k - m + 1):k
        for n in 0:(j - 1)
            result[hh + n] = result[n]
        end
        result[0] = zero
        for n in 1:j
            ind = ell + n
            xb = t[ind]
            xa = t[ind - j]
            if (xb == xa)
                result[m] = zero
                continue
            end
            w = j * result[hh + n - 1] / (xb - xa)
            result[n - 1] = result[n - 1] - w
            result[n] = w
        end
    end
    offset = ell - k
    v = zeros(T, bs.ncoef)
    v[(1 + offset):(k + 1 + offset)] = result[0:k]
    return v
end

function basis(bs::BSplineBasis{T}, x::T, ders::Int=0) where {T<:Real}
    if (x < bs.boundary_knots[1] || x > bs.boundary_knots[2])
        if (x < bs.boundary_knots[1])
            k_pivot =
                T(0.75) * bs.boundary_knots[1] + T(0.25) * bs.spline_basis.knots[5 - 1] # 0-based
        else
            k_pivot =
                T(0.75) * bs.boundary_knots[2] +
                T(0.25) * bs.spline_basis.knots[length(bs.spline_basis.knots) - 4 - 1] # 0-based
        end
        delta = x - k_pivot
        if (ders == 0)
            vec =
                basis(bs.spline_basis, k_pivot, 0) +
                basis(bs.spline_basis, k_pivot, 1) * delta +
                basis(bs.spline_basis, k_pivot, 2) * delta * delta / T(2.0) +
                basis(bs.spline_basis, k_pivot, 3) * delta * delta * delta / T(6.0)
        elseif (ders == 1)
            vec =
                splines.basis(bs.spline_basis, k_pivot, 1) +
                splines.basis(bs.spline_basis, k_pivot, 2) * delta +
                splines.basis(bs.spline_basis, k_pivot, 3) * delta * delta / T(2.0)
        elseif (ders == 2)
            vec =
                splines.basis(bs.spline_basis, k_pivot, 2) +
                splines.basis(bs.spline_basis, k_pivot, 3) * delta
        elseif (ders == 3)
            vec = splines.basis(bs.spline_basis, k_pivot, 3)
        else
            vec = k_pivot .* T(0)
        end
    else
        vec = basis(bs.spline_basis, x, ders)
    end
    if (!bs.intercept)
        vec = vec[2:length(vec)]
    end
    return vec
end

function basis(ns::NSplineBasis{T}, x::T, ders::Int=0) where {T<:Real}
    if (x < ns.b_spline_basis.boundary_knots[1])
        if (ders == 0)
            vec = ns.tl0 + (x - ns.b_spline_basis.boundary_knots[1]) * ns.tl1
        elseif (ders == 1)
            vec = ns.tl1
        else
            vec = ns.tl1 .* T(0)
        end
    elseif (x > ns.b_spline_basis.boundary_knots[2])
        if (ders == 0)
            vec = ns.tr0 + (x - ns.b_spline_basis.boundary_knots[2]) * ns.tr1
        elseif (ders == 1)
            vsc = ns.tr1
        else
            vec = ns.tr1 .* T(0)
        end
    else
        vec = ns.qmat * basis(ns.b_spline_basis, x, ders)
    end
    return vec
end

function basis(
    bs::AbstractSplineBasis{T}, x::AbstractArray{T,1}, ders::Int=0
) where {T<:Real}
    f(xi) = basis(bs, xi, ders)
    return copy(transpose(reduce(hcat, f.(x))))
end

# utility function for processing the spline arguments
function spline_args(
    x;
    boundary_knots::Union{Tuple{T,T},Nothing}=nothing,
    interior_knots::Union{Array{T,1},Nothing}=nothing,
    order::Int=4,
    intercept::Bool=false,
    df::Int=3 + Int(intercept),
    knots::Union{Array{T,1},Nothing}=nothing,
    knots_offset::Int=0,
) where {T<:Real}
    if (interior_knots !== nothing && boundary_knots !== nothing)
        # pass
    elseif (knots !== nothing) #TODO: Remove?
        boundary_knots = extrema(knots)
        interior_knots = length(knots) == 2 ? nothing : knots[2:(length(knots) - 1)]
    else
        if (boundary_knots === nothing)
            boundary_knots = extrema(x)
        end
        iKnots = df - order + knots_offset + 1 - Int(intercept)
        if (iKnots > 0)
            p = range(0; length=iKnots + 2, stop=1)[2:(iKnots + 1)]
            # index = (x .>= boundary_knots[1]) .* (x .<= boundary_knots[2]) RM!
            interior_knots = Statistics.quantile(x, p)
        end
    end
    return (boundary_knots, interior_knots)
end

"""
    bs_(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for B-splines and return a function with signature
`(x:: Array{T,1}; ders :: Int = 0)` for evaluation of `ders`
derivative for the splines at `x`.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int = order - 1 + Int(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines
"""
function bs_(
    x::Array{T,1};
    boundary_knots::Union{Tuple{T,T},Nothing}=nothing,
    interior_knots::Union{Array{T,1},Nothing}=nothing,
    order::Int=4,
    intercept::Bool=false,
    df::Int=order - 1 + Int(intercept),
    knots::Union{Array{T,1},Nothing}=nothing,
    centre::Union{T,Nothing}=nothing,
) where {T<:Real}
    (boundary_knots, interior_knots) = spline_args(
        x;
        boundary_knots=boundary_knots,
        interior_knots=interior_knots,
        order=order,
        intercept=intercept,
        df=df,
        knots=knots,
    )
    spline = BSplineBasis(boundary_knots, interior_knots, order, intercept)
    return spline
end

"""
    bs(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for B-splines.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int = order - 1 + Int(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines
- `ders :: Int = 0`: derivatives of the splines

"""
function bs(x::Array{T,1}; ders::Int=0, kwargs...) where {T<:Real}
    return bs_(x; kwargs...)(x; ders=ders)
end

"""
    ns_(x :: Array{T,1}; <keyword arguments>) where T<:Real

Calculate a basis for natural B-splines and return a function with signature
`(x:: Array{T,1}; ders :: Int = 0)` for evaluation of `ders`
derivative for the splines at `x`.

The keyword arguments include one of:
1. `df`, possibly in combination with `intercept`
2. `boundary_knots` and `interior_knots`
3. `knots`

# Arguments
- `boundary_knots :: Union{Tuple{T,T},Nothing} = nothing`: boundary knots
- `interior_knots :: Union{Array{T,1},Nothing} = nothing`: interior knots
- `order :: Int = 4`: order of the spline
- `intercept :: Bool = false`: bool for whether to include an intercept
- `df :: Int = order - 1 + Int(intercept)`: degrees of freedom
- `knots :: Union{Array{T,1}, Nothing} = nothing`: full set of knots (excluding repeats)
- `centre :: Union{T,Nothing} = nothing)`: value to centre the splines

"""
function ns_(
    x;
    boundary_knots::Union{Tuple{T,T},Nothing}=nothing,
    interior_knots::Union{Array{T,1},Nothing}=nothing,
    order::Int=4,
    intercept::Bool=false,
    df::Int=order - 3 + Int(intercept),
    knots::Union{Array{T,1},Nothing}=nothing,
    centre::Union{T,Nothing}=nothing,
) where {T<:Real}
    (boundary_knots, interior_knots) = spline_args(
        x;
        boundary_knots=boundary_knots,
        interior_knots=interior_knots,
        order=order,
        intercept=intercept,
        df=df,
        knots=knots,
        knots_offset=2,
    )
    spline = NSplineBasis(boundary_knots, interior_knots, order, intercept)
    return spline
end

function (spline::Union{BSplineBasis,NSplineBasis})(x; ders::Int=0)
    b = basis(spline, x, ders)
    #if (centre != nothing && ders==0)
    #    bc = basis(spline, centre, ders)
    #    for i=1:size(b,1)
    #        b[i,:] -= bc
    #    end
    #end
    return b
end

ext: jl


text:
module FindLocalIndex

using Gridap

import Base: size
import Base: getindex

export find_local_index
export get_local_item

function find_local_index(gids::CellNumber,lids_to_gids::CellArray)
  LocalIndices(Int,gids,lids_to_gids)
end

function get_local_item(a::IndexCellArray,litem::Int)
  LocalItems(a,litem)
end

struct LocalIndices{T,A,B} <: IndexCellNumber{T,1}
  gids::A
  lids_to_gids::B
end

function LocalIndices(
  T::Type,gids::IndexCellNumber,lids_to_gids::IndexCellArray)
  @assert T <: Integer
  A = typeof(gids)
  B = typeof(lids_to_gids)
  LocalIndices{T,A,B}(gids,lids_to_gids)
end

size(l::LocalIndices) = (length(l.gids),)

function getindex(l::LocalIndices{T},i::Integer) where T
  j = l.gids[i]
  lid_to_gid = l.lids_to_gids[j]
  k = zero(T)
  for g in lid_to_gid
    k += 1
    if g == i
      return k
    end
  end
  return zero(T)
end

struct LocalItems{T,A} <: IndexCellNumber{T,1}
  a::A
  lindex::Int
end

function LocalItems(a::IndexCellArray{T},lindex::Int) where T
  A = typeof(a)
  LocalItems{T,A}(a,lindex)
end

size(l::LocalItems) = (length(l.a),)

function getindex(l::LocalItems,i::Integer)
  v = l.a[i]
  v[l.lindex]
end

end # module


text:
using Test

println("Running tests:")

for x in (:CampbellCochrane, :Wachter, :BansalYaron, :GarleanuPanageas, :DiTella, :Haddad, :ArbitrageHoldingCosts)
	try
		include("../examples/AssetPricing/$(x).jl")
		@test distance <= 1e-5
		println("\t\033[1m\033[32mPASSED\033[0m: $(x)")
	catch e
		println("\t\033[1m\033[31mFAILED\033[0m: $(x)")
		showerror(stdout, e, backtrace())
		rethrow(e)
	end
end

# Test Campbell Cochrane with fixed time grid
m = CampbellCochraneModel()
stategrid = initialize_stategrid(m)
y0 = initialize_y(m, stategrid)
τs = range(1000, stop = 0, length = 50)
y, result, distance = pdesolve(m, stategrid, y0)
y2, result, distance = pdesolve(m, stategrid, y0, τs)
@test sum(abs2, y[:p] .- y2[:p][:, end]) <= 1e-10

for x in (:WangWangYang, :AchdouHanLasryLionsMoll_OneAsset, :AchdouHanLasryLionsMoll_TwoAssets)
	try
		include("../examples/ConsumptionProblem/$(x).jl")
		@test distance <= 1e-5
		println("\t\033[1m\033[32mPASSED\033[0m: $(x)")
	catch e
		println("\t\033[1m\033[31mFAILED\033[0m: $(x)")
		showerror(stdout, e, backtrace())
		rethrow(e)
	end
end

for x in (:BoltonChenWang, )
	try
		include("../examples/InvestmentProblem/$(x).jl")
		@test distance <= 1e-5
		println("\t\033[1m\033[32mPASSED\033[0m: $(x)")
	catch e
		println("\t\033[1m\033[31mFAILED\033[0m: $(x)")
		showerror(stdout, e, backtrace())
		rethrow(e)
	end
end


text:
# Various mathematical utils for cryptography processes.
# Copyright (C): Jaka Smrekar (vinctux) <vinctux@outlook.com>, 2016

using Primes
import Primes: isprime

# Get integer bitsize
function bitsize{T<:Integer}(N::T)
    b = BigInt(0)
    while N >> 1 > 0 b += 1 end
    return b
end

# Get reversed bit-order of x
# See: https://graphics.stanford.edu/~seander/bithacks.html
function rev{T<:Integer}(x::T)
    x = (((x & 0xAAAAAAAA) >> 1) | ((x & 0x55555555) << 1))
    x = (((x & 0xCCCCCCCC) >> 2) | ((x & 0x33333333) << 2))
    x = (((x & 0xF0F0F0F0) >> 4) | ((x & 0x0F0F0F0F) << 4))
    x = (((x & 0xFF00FF00) >> 8) | ((x & 0x00FF00FF) << 8))
    return ((x >> 16) | (x << 16))
end

# Get reversed bit-order for every N in X
function rev{T<:Integer}(X::Array{T, 1})
    A = Array{T, 1}()
    for N in X append!(A, rev(N)) end
    return A
end

# Cyclic bit-shift (rotation) to left = ROL
function rotl{T<:Integer}(u::T, c::Int)
    m = 8 * sizeof(typeof(u)) - 1
    c &= m
    return (u << c) | (u >> (m & (-c)))
end

# Cyclic bit-shift (rotation) to right = ROR
function rotr{T<:Integer}(u::T, c::Int)
    m = 8 * sizeof(typeof(u)) - 1
    c &= m
    return (u >> c) | (u << (m & (-c)))
end

# Flatten an nested array
# See http://codegolf.stackexchange.com/a/80140
flat(x, y = vcat(x...)) = x == y ? x : flat(y)

# Get first power of E, greater than N
function gpow{T<:Integer}(N::T, E::Integer = 2)
    if N <= 0 return 0 end
    x, I = N, 0
    while x <= N x *= E; I += 1 end
    return (x, I)
end

# Produce cryptographically-secure octet array of length b
function csrand{T<:Integer}(b::T)
    R, O = RandomDevice(), Array{UInt8, 1}()
    for i in 1:b append!(O, rand(R, 0o0:0o255)) end
    return O
end
csrand{T<:Integer}(t::UnitRange{T}) = rand(RandomDevice(), t)

function csrandbit{T<:Integer}(b::T)
    O = "0b"
    for i in 1:b O *= csrand(0:1) == 0 ? "0" : "1" end
    return parse(BigInt, O)
end
function randbit{T<:Integer}(b::T)
    O = "0b"
    for i in 1:b O *= rand(0:1) == 0 ? "0" : "1" end
    return parse(BigInt, O)
end

# Sieve of Atkin - fast prime generator
# Author of this piece of code: @VicDrastik
# For reference:
#   - https://en.wikipedia.org/wiki/Sieve_of_Atkin
#   - https://github.com/JuliaLang/julia/issues/11594#issuecomment-109753437
function e5(s::BitArray{1})
    n = length(s)
    s[:] = false
    s[2] = true
    s[3] = true
    s[5:6:n] = true
    s[7:6:n] = true
    for p = 5:2:isqrt(n) if s[p] s[p*p:2p:n] = false end end
end

# Generates a strong prime, according to 10.1.1.17.2713
function csprime{T<:Integer}(b::T, ERR::Int8 = Int8(75))
    sieve_base = falses(10^8)
    e5(sieve_base)
    X = csrand((div(1414213562373095 * 2^BigInt(b - 1), 10^15)):2^BigInt(b) - 1)
    y1, y2 = csrandbit(101), csrandbit(101)
    P = Dict(y1 => BigInt(0), y2 => BigInt(0))
    F = falses(5 * length(sieve_base))
    _S = [BigInt(i) for i in 1:length(sieve_base) if sieve_base[i]]
    for y in [y1, y2]
        # FIXME: Optimize with native zip - don't compute twice
        S = [(y % i, i) for i in 1:length(sieve_base) if sieve_base[i]]
        for i in 1:length(S) for j in (y + S[i][1] - S[i][2]):S[i][2]:length(F) F[j] = true end end
        for i in 1:length(F)
            if F[i] continue end
            t = y + i
            if isprime(t, ERR) P[y] = t; break end
        end
    end
    R = invmod(P[y2], P[y1]) * P[y2] - invmod(P[y1], P[y2]) * P[y1]
    Q = P[y1] * P[y2]
    Y0 = X + (R - (X % Q))
    while true
        pp = true
        for p in _S if Y0 % p == 0 pp = false; break end end
        if pp && isprime(Y0, ERR) break end
        Y0 += Q
    end
    return Y0
end
csprime(b::Integer) = csprime(b, Int8(75))

# Euler's totient function
# C++ reference:      https://gist.github.com/cslarsen/1635288
# Pythonic reference: https://stackoverflow.com/a/18114286/6489126
function phi{T<:Integer}(n::T, ERR::Int8 = Int8(75))
    if n < 0 n = -n end
    if n == 0 return 0
    elseif 0 < n <= 2 return 1 end
    if isprime(BigInt(n), ERR) return n - 1 end
    if n & 1 == 0 m = n >> 1; return ~(m & 1) == 0 ? phi(m) << 1 : phi(m) end
    O = 0
    for i in 1:n if gcd(n, i) == 1 O += 1 end end
    return O
end

# GDB function
# Pythonic reference: https://github.com/amintos/PyECC/blob/master/ecc/elliptic.py#L237-L244
function gdb{T<:Integer}(n::T)
    i = 1
    if n <= 0 return 0 end
    while n % i > 0 i <<= 1 end
    return i >> 2
end

# Transform n into a binary representation having signed bits.
# Pythonic reference: https://github.com/amintos/PyECC/blob/master/ecc/elliptic.py#L254-L271
function signbin{T<:Integer}(n::T)
    r = Array{Int8, 1}()
    while n > 1
        if n & 1 == 0w
            if gbd(n + 1) > gbd(n - 1) append!(r, -1); n += 1
            else append!(r, 1); n -= 1 end
        else append!(r, 0) end
        n >>= 1
    end
    r.append(n)
    return reverse(r)
end

# Find first prime greater or equal to N
function gprime(N::Integer)
    if iseven(N) N += 1 end
    while !isprime(N, 75) N += 2 end
    return N
end

# Find first prime greater or equal to N, while also N % M == R
function gprimemod(N::Integer, M::Integer, R::Integer = 1)
    while N % M != R && !isprime(N, 75) N = gprime(N + 2) end
    return N
end


text:
module TestSectorPerformance
using AlphaVantage
using Test

TEST_SLEEP_TIME =  parse(Float64, get(ENV, "TEST_SLEEP_TIME", "15"))
MAX_TESTS = parse(Int64, get(ENV, "MAX_TESTS", "1"))

@testset "Sector Performance" begin
    data = sector_performance()
    @test typeof(data) === Dict{String, Any}
    @test length(data) === 11
    sleep(TEST_SLEEP_TIME + 2*rand()) #as to not overload the API
end

end # module


text:
using Distributed
if nprocs() == 1
   # test with multiple processes, but use only 2
   # as travis testing has two processes
   addprocs(2)
end

# This assumes we are in the /src or /test folder.
@everywhere push!(LOAD_PATH, "../src")
import BoltzmannMachines
BMs = BoltzmannMachines

@everywhere push!(LOAD_PATH, "../test/")
import BMTest

#######################
# Deterministic tests #
#######################

BMTest.test_potentials()
BMTest.test_likelihoodbaserates()
BMTest.test_stackrbms_preparetrainlayers()
BMTest.test_dbmjoining()
BMTest.test_rbm_monitoring(BMs.GaussianBernoulliRBM)
BMTest.test_rbm_monitoring(BMs.GaussianBernoulliRBM2)
BMTest.test_mdbm_architectures()
BMTest.test_mdbm_monitoring()
BMTest.test_summing_out()
BMTest.test_unevenbatches()
BMTest.test_beam()
BMTest.test_softmaxvssigm()
BMTest.test_oneornone_encoding()
BMTest.test_intensities()
BMTest.test_monitored_fitrbm()
BMTest.test_monitored_stackrbms()
BMTest.test_monitored_fitdbm()
BMTest.test_batchsize_dbm_training()
BMTest.test_epochs_dbm_training()
BMTest.test_learningrates_dbm_training()
BMTest.test_top2latentdims()

println("Finished deterministic tests")

#########################################################################
# Tests which may fail sometimes due to stochastic nature of algorithms #
#########################################################################

nondeterministictestidx = 0

# Perform non-deterministic tests, allow failures up to n times
function softly(fun; nfail::Int = 7)

   treatasfail = [ErrorException("Not enough samples")]

   for i = 1:nfail
      pass = false
      try
         pass = fun()
      catch ex
         if ex in treatasfail
            pass = false
            @warn ex
         else
            rethrow(ex)
         end
      end

      if pass
         global nondeterministictestidx += 1
         # print some output to keep Travis going
         println("Non-deterministic test #$nondeterministictestidx passed" *
               " after $i trie(s).")
         return true
      end
   end
   false
end

using Test

@test softly(BMTest.check_rbm)
@test softly(BMTest.check_b2brbm)
@test softly(BMTest.check_likelihoodconsistency)
@test softly(() -> BMTest.check_gbrbm(BMs.GaussianBernoulliRBM))
@test softly(() -> BMTest.check_gbrbm(BMs.GaussianBernoulliRBM2))
@test softly(() ->
      BMTest.check_likelihoodconsistency_gaussian(BMs.GaussianBernoulliRBM))
@test softly(() ->
      BMTest.check_likelihoodconsistency_gaussian(BMs.GaussianBernoulliRBM2))
@test softly(BMTest.check_mdbm_rbm_b2brbm)
@test softly(BMTest.check_softmax0rbm)
@test softly(BMTest.check_softmax0dbm)
@test softly(BMTest.check_mdbm_gaussianvisibles)
@test softly(BMTest.check_exactsampling)
@test softly(BMTest.check_softmaxsampling)
@test softly(BMTest.check_sampling)

##########################################
# Run examples (without actual plotting) #
##########################################
@everywhere push!(LOAD_PATH, "../test/mock")
using BoltzmannMachinesPlots
include("../test/examples.jl")

@everywhere pop!(LOAD_PATH)
@everywhere pop!(LOAD_PATH)
@everywhere pop!(LOAD_PATH)


text:
# include("Load_LTSpice_Net.jl")
# include("SPICE2Matrix.jl")

"""
This function runs an AC analysis (freq. sweep) for an LTSPICE netlist file.

    if no inputs are given, it will prompt for the user to pick a file

    Inputs:
    FileName: Path to netlist (LTSPICE: View→SPICE netlist, and it will show a text file as well as populate in PWD)
    FreqList: List of frequencies to test (Hz), Default is 100 Hz to 100kHz in steps of 10
    inputs: vector of independent nodal current inputs (1: NumNodes) and independent V inputs (NumNodes+1:End)
            All voltage inputs default to 1
"""
function RunACAnalysis(FileName=nothing, FreqList = 100:10:100e3, inputs = nothing)
    if FileName===nothing
            FileName = open_dialog("Pick a file")
    end

    SPICE_DF,NodeList,InputList,NumVSources = SPICE2Matrix(FileName)
    println(InputList)
    if inputs===nothing
        println("No inputs given")
        inputs = zeros(length(InputList))
        inputs[end-(NumVSources-1):end] .= 1
    end

    ResultNodeNames = vcat("V(".*NodeList.*")", "I(".*(InputList[end-(NumVSources-1):end]).*")")

    Results = [ abs.(inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*FreqList[i],InputList))*inputs) for i in 1:length(FreqList)]
    Results = hcat(Results...)
    ResDict = Dict(ResultNodeNames[1] => Results[1,:])
    for i in 2:length(ResultNodeNames)

         merge!(ResDict,Dict(ResultNodeNames[i] => Results[i,:]))

    end

    return FreqList,Results, ResDict

end

"""
This function takes in:
    Filename to LTSPICE netlist (View→ SPICE netlist)
    DriveFreq, primary frequency of interest in Hz
    Θᵣ = 1, Θc = 1 are the lumped thermal resistance of inductors and resistors ( Θᵣ) and capacitors(Θc) in ᵒC per Watt
    TCᵣ = 0.004, TCc = 0.0003 are the temperature coefficients of the components-- TCᵣis the resistance tempco, and TCc is the capacitor tempco
        Units for TC is in UL per deg C

    InputScaling is the scaling of the input vector, for instance if you use the default input vec it will default to 1V at voltage sources, but this scalar can be used to modify that value

"""
function DetermineTempCo(FileName=nothing, DriveFreq = 25e3, ComponentName = "Ldrive";
    FreqList = 100:10:100e3,
    inputs = nothing,
    Θᵣ = 1, Θc = 1,
    TCᵣ = 0.004, TCc = 0.0003,
    InputScaling = 1)

    if FileName===nothing
        FileName = open_dialog("Pick a file")
    end

    SPICE_DF,NodeList,InputList,NumVSources = SPICE2Matrix(FileName)
    println(InputList)
    if inputs===nothing
        println("No inputs given")
        inputs = zeros(length(InputList))
        inputs[end-(NumVSources-1):end] .= (1*InputScaling)
    end

    ResultNodeNames = vcat("V(".*NodeList.*")", "I(".*(InputList[end-(NumVSources-1):end]).*")")

    Results = [ abs.(inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*FreqList[i],InputList))*inputs) for i in 1:length(FreqList)]
    Results = hcat(Results...)
    ResDict = Dict(ResultNodeNames[1] => Results[1,:])
    for i in 2:length(ResultNodeNames)

        merge!(ResDict,Dict(ResultNodeNames[i] => Results[i,:]))

    end
    CurVec = plotACElCurrent(SPICE_DF,FreqList,Results,ComponentName)
    plot(FreqList,abs.(CurVec))
    CurrResults =[ inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*DriveFreq,InputList))*inputs]
    CurrResults = hcat(CurrResults...)

    CurrentDF = getElementCurrents(SPICE_DF,CurrResults,DriveFreq)
    CurrentElIndex = findfirst(isequal(ComponentName),CurrentDF.Name)

    BaselineCurrent = CurrentDF.Current[CurrentElIndex]

    SPICE_DF.Drift = zeros(length(SPICE_DF.Name))
    SPICE_DF.TempRise = zeros(length(SPICE_DF.Name))
    SPICE_DF.Dissipation = zeros(length(SPICE_DF.Name))
    for i in 1:length(SPICE_DF.Name)

        if SPICE_DF.Type[i] =='R'
            Z = SPICE_DF.Value[i]
            TmpIndex = findfirst(isequal(SPICE_DF.Name[i]),CurrentDF.Name)

            SPICE_DF.Dissipation[i] = abs(CurrentDF.Current[TmpIndex])^2*Z #Power = I²R
            SPICE_DF.TempRise[i]    = Θᵣ*SPICE_DF.Dissipation[i]
            SPICE_DF.Drift[i]       = (1 + Θᵣ*SPICE_DF.Dissipation[i]*TCᵣ)
            SPICE_DF.Value[i]       = (1 + Θᵣ*SPICE_DF.Dissipation[i]*TCᵣ) * SPICE_DF.Value[i]

        elseif SPICE_DF.Type[i] =='L'
            Z = SPICE_DF.ESR[i]
            TmpIndex = findfirst(isequal(SPICE_DF.Name[i]),CurrentDF.Name)

            SPICE_DF.Dissipation[i] = abs(CurrentDF.Current[TmpIndex])^2*Z #Power = I²(ESR)
            SPICE_DF.TempRise[i]    = Θᵣ*SPICE_DF.Dissipation[i]
            SPICE_DF.Drift[i]       = (1 + Θᵣ*SPICE_DF.Dissipation[i]*TCᵣ)
            SPICE_DF.ESR[i]       = (1 + Θᵣ*SPICE_DF.Dissipation[i]*TCᵣ) * SPICE_DF.ESR[i]# for inductors, only the ESR will drift, not inductance
        elseif SPICE_DF.Type[i] =='C'
            Z = SPICE_DF.ESR[i]
            TmpIndex = findfirst(isequal(SPICE_DF.Name[i]),CurrentDF.Name)

            SPICE_DF.Dissipation[i] = abs(CurrentDF.Current[TmpIndex])^2*Z #Power = I²(ESR)
            SPICE_DF.TempRise[i]    = Θc*SPICE_DF.Dissipation[i]
            SPICE_DF.Drift[i]       = (1 + Θc*SPICE_DF.Dissipation[i]*TCc)
            SPICE_DF.Value[i]       = (1 + Θc*SPICE_DF.Dissipation[i]*TCc) * SPICE_DF.Value[i]

        end
    end

    NewCurrResults =[ inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*DriveFreq,InputList))*inputs]
    NewCurrResults = hcat(NewCurrResults...)

    NewCurrentDF = getElementCurrents(SPICE_DF,NewCurrResults,DriveFreq)
    PostHeatCurrent = NewCurrentDF.Current[CurrentElIndex]

    MagDriftPercent = abs((PostHeatCurrent - BaselineCurrent) / PostHeatCurrent)*100

    Results2 =[ abs.(inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*FreqList[i],InputList))*inputs) for i in 1:length(FreqList)]
    Results2 = hcat(Results2...)
    ResDict2 = Dict(ResultNodeNames[1] => Results2[1,:])
    for i in 2:length(ResultNodeNames)

        merge!(ResDict2,Dict(ResultNodeNames[i] => Results2[i,:]))

    end
    CurVec2 = plotACElCurrent(SPICE_DF,FreqList,Results2,ComponentName)
    plot(FreqList,abs.(CurVec2))

    return SPICE_DF,MagDriftPercent,NewCurrentDF
end

"""
This function takes in the SPICE dataframe, results vector(Nx1), and the frequency to test it at.
        This is needed because the Results vector is in terms of node voltages (and currents in voltage sources) but often the current in a single component is useful to know.

        This gets the voltage across two nodes and divides by the impednace of a component

"""
function getElementCurrents(SPICE_DF,Results,Freq)
    NumPasives = sum(SPICE_DF.Type .!= 'V')
    CurrentDF = DataFrame((Name = Any[], ΔV = Any[], Z = Any[],  Current = Any[]))
    for i in 1:length(SPICE_DF.Name)
        if SPICE_DF.Type[i] !='V'
            N1 = SPICE_DF.Node1[i]
            N2 = SPICE_DF.Node2[i]
            if (N1+N2)<999 #if both are non-GND
                ΔV  = Results[N1] - Results[N2]
            elseif N1==999
                ΔV  = Results[N2]
            else
                ΔV  = Results[N1]
            end

        end
        if SPICE_DF.Type[i] =='R'
            Z = SPICE_DF.Value[i]
            Cur = ΔV ./ Z
            push!(CurrentDF,[SPICE_DF.Name[i] ΔV Z Cur])

        elseif SPICE_DF.Type[i] =='L'
            Z = 2 .* π .* Freq.*im.*SPICE_DF.Value[i] .+ SPICE_DF.ESR[i]
            Cur = ΔV ./ Z
            push!(CurrentDF,[SPICE_DF.Name[i] ΔV Z Cur])

        elseif SPICE_DF.Type[i] =='C'
            Z = 1 ./( 2 .* π .* Freq.*im.*SPICE_DF.Value[i]) .+ SPICE_DF.ESR[i]
            Cur = ΔV ./ Z
            push!(CurrentDF,[SPICE_DF.Name[i] ΔV Z Cur])

        end

    end
    return CurrentDF
end

"""
Plot AC Element Current
similar to getElementCurrents, except the freq. is a vector, Results in Nx(FreqList Length), and only looks at a single component's current
"""
function plotACElCurrent(SPICE_DF,FreqList,Results,ComponentName)
    ElIndex = findfirst(isequal(ComponentName),SPICE_DF.Name)

            N1 = SPICE_DF.Node1[ElIndex]
            N2 = SPICE_DF.Node2[ElIndex]
            if (N1+N2)<999 #if both are non-GND
                ΔV  = Results[N1,:] - Results[N2,:]
            elseif N1==999
                ΔV  = Results[N2,:]
            else
                ΔV  = Results[N1,:]
            end

        if SPICE_DF.Type[ElIndex] =='R'
            Z = SPICE_DF.Value[ElIndex]
            CurrentVec = ΔV ./ Z

        elseif SPICE_DF.Type[ElIndex] =='L'
            Z = [2 .* π .* FreqList[i].*im.*SPICE_DF.Value[ElIndex] .+ SPICE_DF.ESR[ElIndex] for i in 1:length(FreqList)]
            CurrentVec = ΔV ./ Z

        elseif SPICE_DF.Type[ElIndex] =='C'
            Z = [1 ./( 2 .* π .* FreqList[i].*im.*SPICE_DF.Value[ElIndex]) .+ SPICE_DF.ESR[ElIndex] for i in 1:length(FreqList)]
            CurrentVec = ΔV ./ Z

        end

    return CurrentVec
end

"""
This function modifies the spice dataframe. Specifically, it updates a value of a component
    DF should be a SPICE Dataframe
    ComponentName should be a string of the name (e.g. ` "LDrive" `)
    NewVal should be a number, and the new value (e.g 0.005 for a 5mH inductor)

"""
function UpdateElementVal!(DF, ComponentName::String,NewVal)
    ElIndex = findfirst(isequal(ComponentName),DF.Name)
    DF.Value[ElIndex] = NewVal
end

"""
This function modifies the spice dataframe. Specifically, it updates a value of a component's ESR
    DF should be a SPICE Dataframe
    ComponentName should be a string of the name (e.g. ` "LDrive" `)
    NewVal should be a number, and the new value (e.g 0.01 for a 10mΩ ESR)

"""
function UpdateElementESR!(DF, ComponentName::String,NewVal)
    ElIndex = findfirst(isequal(ComponentName),DF.Name)
    DF.ESR[ElIndex] = NewVal
end

function GetElementVal(SPICE_DF,Name)
    SPICE_DF.Value[findfirst(isequal(Name),SPICE_DF.Name)]
end

function addPowerDiss_DF!(SPICE_DF,inputs,InputList,DriveFreq;TargetRMSCurrent=35,ComponentName="LDrive")

    CurrResults =[ inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*DriveFreq,InputList))*inputs]
    CurrResults = hcat(CurrResults...)

    CurrentDF = getElementCurrents(SPICE_DF,CurrResults,DriveFreq)
    CurrentElIndex = findfirst(isequal(ComponentName),CurrentDF.Name)

    BaselineCurrent = CurrentDF.Current[CurrentElIndex]

    UpdatedInputVoltage = TargetRMSCurrent/abs.(BaselineCurrent)

    CurrResults =[ inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*DriveFreq,InputList))*inputs*UpdatedInputVoltage]
    CurrResults = hcat(CurrResults...)

    CurrentDF = getElementCurrents(SPICE_DF,CurrResults,DriveFreq)

    SPICE_DF.Dissipation = zeros(length(SPICE_DF.Name))
    for i in 1:length(SPICE_DF.Name)

        if SPICE_DF.Type[i] =='R'
            Z = SPICE_DF.Value[i]
            TmpIndex = findfirst(isequal(SPICE_DF.Name[i]),CurrentDF.Name)

            SPICE_DF.Dissipation[i] = abs(CurrentDF.Current[TmpIndex])^2*Z #Power = I²R

        elseif SPICE_DF.Type[i] =='L'
            Z = SPICE_DF.ESR[i]
            TmpIndex = findfirst(isequal(SPICE_DF.Name[i]),CurrentDF.Name)

            SPICE_DF.Dissipation[i] = abs(CurrentDF.Current[TmpIndex])^2*Z #Power = I²(ESR)

        elseif SPICE_DF.Type[i] =='C'
            Z = SPICE_DF.ESR[i]
            TmpIndex = findfirst(isequal(SPICE_DF.Name[i]),CurrentDF.Name)

            SPICE_DF.Dissipation[i] = abs(CurrentDF.Current[TmpIndex])^2*Z #Power = I²(ESR)

        end

    end
    return CurrentDF
end

"""
This function takes in:

    DriveFreq, primary frequency of interest in Hz
    Θᵣ = 1, Θc = 1 are the lumped thermal resistance of inductors and resistors ( Θᵣ) and capacitors(Θc) in ᵒC per Watt
    TCᵣ = 0.004, TCc = 0.0003 are the temperature coefficients of the components-- TCᵣis the resistance tempco, and TCc is the capacitor tempco
        Units for TC is in UL per deg C

    InputScaling is the scaling of the input vector, for instance if you use the default input vec it will default to 1V at voltage sources, but this scalar can be used to modify that value

"""
function DetermineComponentsTempCoeffs(SPICE_DF,InputList,NumVSources,DriveFreq = 25e3, ComponentName = "LDrive";
    Θᵣ = 1, Θc = 1,
    TCᵣ = 0.004, TCc = 0.0003,
    InputScaling = 1,
    δppm = 900)

    δ = δppm*1e-6

    inputs = zeros(length(InputList))
    inputs[end-(NumVSources-1):end] .= (1*InputScaling)

    CurrResults =[ inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*DriveFreq,InputList))*inputs]
    CurrResults = hcat(CurrResults...)

    CurrentDF = getElementCurrents(SPICE_DF,CurrResults,DriveFreq)
    CurrentElIndex = findfirst(isequal(ComponentName),CurrentDF.Name)

    BaselineCurrent = CurrentDF.Current[CurrentElIndex]

    SPICE_DF.Drift = zeros(length(SPICE_DF.Name))
    SPICE_DF.TempRise = zeros(length(SPICE_DF.Name))
    SPICE_DF.Dissipation = zeros(length(SPICE_DF.Name))
    SPICE_DF.DriftCoeff = zeros(length(SPICE_DF.Name))
    for i in 1:length(SPICE_DF.Name)

        if SPICE_DF.Type[i] =='R'
            Z = SPICE_DF.Value[i]
            TmpIndex = findfirst(isequal(SPICE_DF.Name[i]),CurrentDF.Name)

            SPICE_DF.Dissipation[i] = abs(CurrentDF.Current[TmpIndex])^2*Z #Power = I²R
            SPICE_DF.TempRise[i]    = Θᵣ*SPICE_DF.Dissipation[i]
            SPICE_DF.Drift[i]       = (1 + Θᵣ*SPICE_DF.Dissipation[i]*TCᵣ)
            dVal = SPICE_DF.Value[i]*δ
            SPICE_DF.Value[i]       = dVal +  SPICE_DF.Value[i]

            NewCurrResults =[ inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*DriveFreq,InputList))*inputs]
            NewCurrResults = hcat(NewCurrResults...)

            NewCurrentDF = getElementCurrents(SPICE_DF,NewCurrResults,DriveFreq)
            PostHeatCurrent = NewCurrentDF.Current[CurrentElIndex]

            SPICE_DF.DriftCoeff[i] = abs((PostHeatCurrent - BaselineCurrent) / BaselineCurrent)/δ
            SPICE_DF.Value[i]       =  SPICE_DF.Value[i]-dVal

        elseif SPICE_DF.Type[i] =='L'
            Z = SPICE_DF.ESR[i]
            TmpIndex = findfirst(isequal(SPICE_DF.Name[i]),CurrentDF.Name)

            SPICE_DF.Dissipation[i] = abs(CurrentDF.Current[TmpIndex])^2*Z #Power = I²(ESR)
            SPICE_DF.TempRise[i]    = Θᵣ*SPICE_DF.Dissipation[i]
            SPICE_DF.Drift[i]       = (1 + Θᵣ*SPICE_DF.Dissipation[i]*TCᵣ)
            dVal = SPICE_DF.ESR[i]*δ
            SPICE_DF.ESR[i]       = dVal +  SPICE_DF.ESR[i]# for inductors, only the ESR will drift, not inductance
            NewCurrResults =[ inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*DriveFreq,InputList))*inputs]
            NewCurrResults = hcat(NewCurrResults...)

            NewCurrentDF = getElementCurrents(SPICE_DF,NewCurrResults,DriveFreq)
            PostHeatCurrent = NewCurrentDF.Current[CurrentElIndex]

            SPICE_DF.DriftCoeff[i] = abs((PostHeatCurrent - BaselineCurrent) / BaselineCurrent)/δ

            SPICE_DF.ESR[i]       =  SPICE_DF.ESR[i]-dVal

        elseif SPICE_DF.Type[i] =='C'
            # TempCo of celem PP caps is -300ppm/⁰
            lossTangent = 2e-4

            TmpIndex = findfirst(isequal(SPICE_DF.Name[i]),CurrentDF.Name)

            SPICE_DF.Dissipation[i] = abs(CurrentDF.Current[TmpIndex])*abs(CurrentDF.ΔV[TmpIndex])*lossTangent #Power = Total apparent power*tandδ
            SPICE_DF.TempRise[i]    = Θc*SPICE_DF.Dissipation[i]
            SPICE_DF.Drift[i]       = (1 + Θc*SPICE_DF.Dissipation[i]*TCc)
            dVal = SPICE_DF.Value[i]*δ
            SPICE_DF.Value[i]       = SPICE_DF.Value[i]-dVal #Subtract val because they lower capacitance as heat is applied
            NewCurrResults =[ inv(SPICE_DF2Matrix_ω(SPICE_DF,2*π*DriveFreq,InputList))*inputs]
            NewCurrResults = hcat(NewCurrResults...)

            NewCurrentDF = getElementCurrents(SPICE_DF,NewCurrResults,DriveFreq)
            PostHeatCurrent = NewCurrentDF.Current[CurrentElIndex]

            SPICE_DF.DriftCoeff[i] = abs((PostHeatCurrent - BaselineCurrent) / BaselineCurrent)/δ
            SPICE_DF.Value[i]       =  SPICE_DF.Value[i]+dVal

        end
    end

    return SPICE_DF
end


text:
using AbstractMCMC

"""
    AbstractProbabilisticProgram

Common base type for models expressed as probabilistic programs.
"""
abstract type AbstractProbabilisticProgram <: AbstractMCMC.AbstractModel end

"""
    logdensity(model, trace)

Evaluate the (possibly unnormalized) density of the model specified by the probabilistic program
in `model`, at specific values for the random variables given through `trace`.

`trace` can be of any supported internal trace type, or a fixed probability expression.

`logdensity` should interact with conditioning and deconditioning in the way required by probability
theory.
"""
function logdensity end

"""
    decondition(conditioned_model)

Remove the conditioning (i.e., observation data) from `conditioned_model`, turning it into a
generative model over prior and observed variables.

The invariant

```
m == condition(decondition(m), obs)
```

should hold for models `m` with conditioned variables `obs`.
"""
function decondition end

"""
    condition(model, observations)

Condition the generative model `model` on some observed data, creating a new model of the (possibly
unnormalized) posterior distribution over them.

`observations` can be of any supported internal trace type, or a fixed probability expression.

The invariant

```
m = decondition(condition(m, obs))
```

should hold for generative models `m` and arbitrary `obs`.
"""
function condition end


text:
"""
Simple server in Julia and client code in JS.

### Example client code (JS):
```http
<html>
<head>
    <meta charset="UTF-8">
    <title>Squaring numbers</title>
</head>
<body>
    <input id="number" placeholder="Input a number" type="number">
    <button id="submit">Square</button>
    <h4>Outputs</h4>
    <ul id="list"></ul>
</body>
<script>
    document.getElementById('submit').addEventListener('click', async function (event) {
        const list = document.getElementById('list');
        try {
            const r = await fetch('http://127.0.0.1:8080/api/square', {
                method: 'POST',
                body: document.getElementById('number').value
            });

            if (r.ok) {
                const body = await r.text()
                const newElement = document.createElement('li');
                newElement.textContent = body;
                list.insertBefore(newElement, list.firstChild);
            } else {
                console.error(r)
            };
        } catch (err) {
            console.error(err)
        }
    })
</script>
</html>
```

### Server code:
"""
using HTTP, JSON

const ROUTER = HTTP.Router()

function square(req::HTTP.Request)
    headers = [
        "Access-Control-Allow-Origin" => "*",
        "Access-Control-Allow-Methods" => "POST, OPTIONS"
    ]
    # handle CORS requests
    if HTTP.method(req) == "OPTIONS"
        return HTTP.Response(200, headers)
    end
    body = parse(Float64, String(HTTP.body(req)))
    square = body^2
    HTTP.Response(200, headers; body = string(square))
end

HTTP.@register(ROUTER, "POST", "/api/square", square)

HTTP.serve(ROUTER, "127.0.0.1", 8080)


text:
# Building OpenSpecFun from scratch

using BinDeps
using BinDeps: libdir, srcdir, includedir, depsdir, builddir
using Base.Math: libm
import BinaryProvider
using Libdl

BinDeps.@setup
const OSF_VERS = v"0.5.3"
openspecfun = library_dependency("libopenspecfun")

# If Julia is built with OpenLibm, we want to build OpenSpecFun with it as well.
# Unfortunately this requires a fair bit more work, as we need to link to the .so
# and to include the headers, which aren't readily available.
if libm == "libopenlibm"
    const OLM_VERS = v"0.5.4"
    const OLM_HASH = "9a8ae1d17825a4a6a4c013d36a7f4348b27c47eedb6549c521ecc9c79d021c13"
    const OLM_URL = "https://github.com/JuliaLang/openlibm/archive/v$OLM_VERS.tar.gz"
    use_openlibm = true

    if !isdir(libdir(openspecfun))
        mkpath(libdir(openspecfun))
    end

    openlibm_so = Libdl.dlpath(libm)

    # Copy over the OpenLibm .so
    for lib in readdir(dirname(openlibm_so))
        startswith(lib, "libopenlibm") || continue
        cp(joinpath(dirname(openlibm_so), lib), joinpath(libdir(openspecfun), lib),
           force=true, follow_symlinks=false)
    end

    # Grab and unpack the tarball so we can get the header filesa
    BinaryProvider.download_verify_unpack(OLM_URL, OLM_HASH, srcdir(openspecfun); verbose = true)

    # Copy over all of the OpenLibm headers
    cp(joinpath(srcdir(openspecfun), "openlibm-$OLM_VERS"), joinpath(srcdir(openspecfun), "openlibm"); force=true)
    openlibm_src = joinpath(srcdir(openspecfun), "openlibm")
    openlibm_include = joinpath(includedir(openspecfun), "openlibm")
    if !isdir(openlibm_include)
        mkpath(openlibm_include)
    end
    for f in readdir(joinpath(openlibm_src, "include"))
        cp(joinpath(openlibm_src, "include", f), joinpath(openlibm_include, f),
           force=true, follow_symlinks=true)
    end
    for f in readdir(joinpath(openlibm_src, "src"))
        if endswith(f, ".h")
            cp(joinpath(openlibm_src, "src", f), joinpath(openlibm_include, f),
               force=true, follow_symlinks=true)
        end
    end
else
    use_openlibm = false
end

fc = "gfortran"

# macOS has precompiled binaries, so it's just FreeBSD that should default to Clang
if Sys.KERNEL === :FreeBSD || Sys.KERNEL == :Darwin
    cc = "clang"
    use_clang = true
else
    cc = "gcc"
    use_clang = false
end

if Sys.ARCH in [:i386, :i387, :i486, :i586, :i686]
    cc *= " -m32"
    fc *= " -m32"
elseif Sys.ARCH === :x86_64
    cc *= " -m64"
    fc *= " -m64"
end

extra_ld = Sys.isapple() ? "" : "-Wl,-rpath,'\$\$ORIGIN' -Wl,-z,origin"

provides(Sources, URI("https://github.com/JuliaLang/openspecfun/archive/v$OSF_VERS.tar.gz"),
         openspecfun, unpacked_dir="openspecfun-$OSF_VERS")

provides(BuildProcess,
    (@build_steps begin
        GetSources(openspecfun)
        CreateDirectory(builddir(openspecfun))
        @build_steps begin
            ChangeDirectory(joinpath(srcdir(openspecfun), "openspecfun-$OSF_VERS"))
            FileRule(joinpath(libdir(openspecfun), "libopenspecfun." * Libdl.dlext),
                @build_steps begin
                    CreateDirectory(libdir(openspecfun))
                    ```
                    $MAKE_CMD install
                        ARCH="$(Sys.ARCH)"
                        CC="$cc"
                        FC="$fc"
                        USECLANG=$(Int(use_clang))
                        USEGCC=$(Int(!use_clang))
                        USE_OPENLIBM=$(Int(use_openlibm))
                        CFLAGS="-O3 -std=c99"
                        FFLAGS="-O2 -fPIC"
                        LDFLAGS="-L$(libdir(openspecfun)) $extra_ld"
                        DESTDIR=""
                        prefix=""
                        libdir="$(libdir(openspecfun))"
                        includedir="$(includedir(openspecfun))"
                        O=
                    ```
                end)
        end
    end), openspecfun)

BinDeps.@install Dict(:libopenspecfun => :openspecfun)

# BinDeps doesn't give us a `check_deps` function in the generated deps.jl file like
# BinaryProvider does. Instead, it uses a macro and checks immediately whether the
# library can be loaded. So let's just fake one.
depsjl = joinpath(depsdir(openspecfun), "deps.jl")
@assert isfile(depsjl)
if !any(line->occursin("check_deps", line), eachline(depsjl))
    open(depsjl, "a") do fh
        println(fh, """
            # NOTE: This function is a compatibility shim for BinaryProvider-based builds
            check_deps() = nothing
            """)
    end
end


text:
############################################################################################################################

# This file contains tests for Graph

############################################################################################################################

function graph_sanity_test(Nv::Int, Ne::Int, g::Graph)
   @test nv(g) == Nv
   @test ne(g) == Ne

   df1, df2 = g

   @test g == g
   @test copy(g) == g
   @test deepcopy(g) == g
end

function graph_label_test(Nv::Int, g::Graph, ls::Vector)
   @test nv(g) == Nv

   @test haslabel(g, rand(ls))

   u = rand(1 : Nv)
   v = rand(1 : Nv)

   @test decode(g, ls[v]) == v
   @test decode(g, ls) == collect(1:Nv)
   @test decode(g, Pair(ls[u], ls[v])) == EdgeID(u, v)

   @test encode(g, u) == ls[u]
   @test encode(g, 1 : 5) == ls[1 : 5]
   @test encode(g) == ls
   @test encode(g, EdgeID(u, v)) == Pair(ls[u], ls[v])
end

@testset "Graph Test" begin
   g = Graph(completeindxs(10))

   graph_sanity_test(10, 90, g)
   graph_label_test(10, g, collect(1:10))

   # Setlabel
   ls = map(string, 1:10)
   setlabel!(g, ls)
   graph_label_test(10, g, ls)

   # Relabel
   relabel!(g, 10, "TEN")
   @test decode(g, "TEN") == 10
   @test encode(g, 10) == "TEN"

   # Vector Relabel
   relabel!(g, [8,9], ["EIGHT", "NINE"])
   @test decode(g, ["EIGHT", "NINE", "TEN"]) == [8,9,10]
   @test encode(g, 8:10) == ["EIGHT", "NINE", "TEN"]

   # Disable labelling
   setlabel!(g)
   graph_label_test(10, g, collect(1:10))

   # Label on a property
   arr = [randstring() for i in 1:10]
   setvprop!(g, :, arr, :p1)
   setlabel!(g, :p1)
   graph_label_test(10, g, arr)
end


text:
function ffit(flist::Array{Function}, x::AbstractVector, y::AbstractVector)
    # TODO: figure out a fast and generic way of doing this
    A = eltype(x)[f(xi) for xi in x, f in flist]
    # TODO: check if this is faster than svd
    pinv(A)*y
end

# TODO: it would be cool to have a macro so that I could do
# @funclist [:1, :x, :x^2]
# and have it make the list of anon functions

x = [0.0, 1.0, 3.0, 5.0]
y = [1.0, 0.0, 2.0, 4.0]

ffit([(x)->1, (x)->x], x, y)
linreg(x, y)

ffit([(x)->1, (x)->x, (x)->x^2], x, y)

x = [-π, -π/2.0, 0, π/2, π]
y = [4.0, 0.0, 1.0, -1.0, -4.0]
ffit([(x)->sin(x/2), (x)->sin(x), (x)->sin(2*x)], x, y)

v = [5, 4, 3, 2, 1, 0, 1, 2, 3, 4, 5.0]
x = 1:length(v)
ffit([(x)->1,(x)->x,(x)->x^2], x, v)
ffit([(x)->1,(x)->x,(x)->x^2,(x)->x^3,(x)->x^4], x, v)

max_stars_repo_name: gabrielgellner/LeastSquares

max_issues_repo_name: gabrielgellner/LeastSquares

max_forks_repo_name: gabrielgellner/LeastSquares


text:
using Documenter, Hecke, Nemo, AbstractAlgebra, Markdown, Pkg

makedocs(
    modules = [Hecke, Nemo, AbstractAlgebra],
    clean   = true,
    format = :html,
    sitename = "Hecke",
    doctest = !false,
    pages = [
      "index.md",
      "number_fields/intro.md",
      "Orders" => [ "orders/introduction.md",
                    "orders/orders.md",
                    "orders/elements.md",
                    "orders/ideals.md",
                    "orders/frac_ideals.md"
                  ],
#      "Maximal Orders" => [ "MaximalOrders/Introduction.md",
#                            "MaximalOrders/Creation.md",
#                            "MaximalOrders/Elements.md",
#                            "MaximalOrders/Ideals.md"
#                          ],
      "abelian/introduction.md",
      "class_fields/intro.md",
      "sparse/intro.md",
      "FacElem.md"
      ]
)

# Hack around to get syntax highlighting working
#cd(joinpath(dirname(pathof(Hecke)), "..", "docs"))
#
#cp("application-f78e5cb881.palette.css", "build/application-f78e5cb881.palette.css", force = true)
#cp("application-e2807e330f.css", "build/application-e2807e330f.css", force = true)

deploydocs(
     deps = nothing,
#    deps = Deps.pip()#="pygments",
#                    "mkdocs==0.16.3",
#                    "python-markdown-math",
#                    "mkdocs-cinder")=#,
    repo = "github.com/thofma/Hecke.jl.git",
    target = "build",
    make = nothing,
    osname = "linux",
    julia = "1.0",
)

# Try out the following deps = Deps.pip("mkdocs==0.17.5", "mkdocs-material==2.9.4" ,"python-markdown-math", "pygments", "pymdown-extensions")


text:
# host array

export MtlArray

mutable struct MtlArray{T,N} <: AbstractGPUArray{T,N}
  buffer::MtlBuffer{T}
  dims::Dims{N}

  dev::MtlDevice
end

## constructors

# type and dimensionality specified, accepting dims as tuples of Ints
function MtlArray{T,N}(::UndefInitializer, dims::Dims{N}; storage=Shared) where {T,N}
    dev = device()
    # Check that requested size is not larger than maximum buffer size allowed
    sizeof(T) * prod(dims) > dev.maxBufferLength && error("Too large of Metal buffer requested of size $(Base.format_bytes(sizeof(T) * prod(dims))) (Max: $(Base.format_bytes(dev.maxBufferLength)))")
    buf = alloc(T, dev, prod(dims); storage=storage)

    obj = MtlArray{T,N}(buf, dims, dev)
    finalizer(obj) do arr
        free(arr.buffer)
    end
    return obj
end

# type and dimensionality specified, accepting dims as series of Ints
MtlArray{T,N}(::UndefInitializer, dims::Integer...) where {T,N} = MtlArray{T,N}(undef, Dims(dims))

# type but not dimensionality specified
MtlArray{T}(::UndefInitializer, dims::Dims{N}) where {T,N} = MtlArray{T,N}(undef, dims)
MtlArray{T}(::UndefInitializer, dims::Integer...) where {T} =
    MtlArray{T}(undef, convert(Tuple{Vararg{Int}}, dims))

# empty vector constructor
MtlArray{T,1}() where {T} = MtlArray{T,1}(undef, 0)

Base.similar(a::MtlArray{T,N}) where {T,N} = MtlArray{T,N}(undef, size(a))
Base.similar(a::MtlArray{T}, dims::Base.Dims{N}) where {T,N} = MtlArray{T,N}(undef, dims)
Base.similar(a::MtlArray, ::Type{T}, dims::Base.Dims{N}) where {T,N} = MtlArray{T,N}(undef, dims)

function Base.copy(a::MtlArray{T,N}) where {T,N}
  b = similar(a)
  @inbounds copyto!(b, a)
end

## array interface

Base.elsize(::Type{<:MtlArray{T}}) where {T} = sizeof(T)

Base.size(x::MtlArray) = x.dims
Base.sizeof(x::MtlArray) = Base.elsize(x) * length(x)

Base.pointer(x::MtlArray) = x.buffer
#Base.pointer(x::MtlArray, i::Integer) = x.ptr + (i-1) * Base.elsize(x)

## interop with other arrays

@inline function MtlArray{T,N}(xs::AbstractArray{T,N}) where {T,N}
  A = MtlArray{T,N}(undef, size(xs))
  copyto!(A, xs)
  return A
end

MtlArray{T,N}(xs::AbstractArray{S,N}) where {T,N,S} = MtlArray{T,N}(map(T, xs))

# underspecified constructors
MtlArray{T}(xs::AbstractArray{S,N}) where {T,N,S} = MtlArray{T,N}(xs)
(::Type{MtlArray{T,N} where T})(x::AbstractArray{S,N}) where {S,N} = MtlArray{S,N}(x)
MtlArray(A::AbstractArray{T,N}) where {T,N} = MtlArray{T,N}(A)

# idempotency
MtlArray{T,N}(xs::MtlArray{T,N}) where {T,N} = xs

## conversions

Base.convert(::Type{T}, x::T) where T <: MtlArray = x

## interop with C libraries

Base.unsafe_convert(::Type{Ptr{T}}, x::MtlArray{T}) where {T} = throw(ArgumentError("cannot take the host address of a $(typeof(x))"))
Base.unsafe_convert(::Type{Ptr{S}}, x::MtlArray{T}) where {S,T} = throw(ArgumentError("cannot take the host address of a $(typeof(x))"))

Base.unsafe_convert(t::Type{MTL.MTLBuffer}, x::MtlArray{T}) where {T}   = Base.unsafe_convert(t, pointer(x))

## interop with GPU arrays

# TODO Figure out global

function Base.convert(::Type{MtlDeviceArray{T,N,AS.Device}}, a::MtlArray{T,N}) where {T,N}
    MtlDeviceArray{T,N,AS.Device}(a.dims, reinterpret(Core.LLVMPtr{T, 1}, pointer(a).handle))
end

Adapt.adapt_storage(::Adaptor, xs::MtlArray{T,N}) where {T,N} =
  convert(MtlDeviceArray{T,N,AS.Device}, xs)

# Adapt.adapt_storage(::Adaptor, xs::MtlArray{T,N}) where {T,N} =
#   convert(Core.LLVMPtr{T,AS.Device}, xs)

function Base.convert(::Type{Core.LLVMPtr{T,AS.Device}}, a::MtlArray{T}) where {T}
    reinterpret(Core.LLVMPtr{T, 1}, a.buffer.handle)
end
## interop with CPU arrays

Base.unsafe_wrap(t::Type{<:Array}, arr::MtlArray, dims; own=false) = unsafe_wrap(t, arr.buffer, dims; own=own)

# We don't convert isbits types in `adapt`, since they are already
# considered GPU-compatible.

Adapt.adapt_storage(::Type{MtlArray}, xs::AbstractArray) =
  isbits(xs) ? xs : convert(MtlArray, xs)

# if an element type is specified, convert to it
Adapt.adapt_storage(::Type{<:MtlArray{T}}, xs::AbstractArray) where {T} =
  isbits(xs) ? xs : convert(MtlArray{T}, xs)

Adapt.adapt_storage(::Type{Array}, xs::MtlArray) = convert(Array, xs)

Base.collect(x::MtlArray{T,N}) where {T,N} = copyto!(Array{T,N}(undef, size(x)), x)

function Base.copyto!(dest::MtlArray{T}, doffs::Integer, src::Array{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(dest.dev, dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::MtlArray{T}, src::Array{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::Array{T}, doffs::Integer, src::MtlArray{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(src.dev, dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::Array{T}, src::MtlArray{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::MtlArray{T}, doffs::Integer, src::MtlArray{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  # TODO: which device to use here?
  if dest.dev == src.dev
    unsafe_copyto!(dest.dev, dest, doffs, src, soffs, n)
  else
    error("Copy between different devices not implemented")
  end
  return dest
end

Base.copyto!(dest::MtlArray{T}, src::MtlArray{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.unsafe_copyto!(dev::MtlDevice, dest::MtlArray{T}, doffs, src::Array{T}, soffs, n) where T
  GC.@preserve src dest begin
    unsafe_copyto!(dev, pointer(dest), doffs, pointer(src, soffs), n)
  end
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

function Base.unsafe_copyto!(dev::MtlDevice, dest::Array{T}, doff, src::MtlArray{T}, soff, n) where T
  GC.@preserve src dest begin
    unsafe_copyto!(dev, pointer(dest, doff), pointer(src), soff, n)
  end
  #GC.@preserve src dest unsafe_copyto!(dev, pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

function Base.unsafe_copyto!(dev::MtlDevice, dest::MtlArray{T}, doffs, src::MtlArray{T}, soffs, n) where T
  GC.@preserve src dest unsafe_copyto!(dev, pointer(dest), doffs, pointer(src), soffs, n)
#  GC.@preserve src dest unsafe_copyto!(dev, pointer(dest, doffs), pointer(src, soffs), n)
  if Base.isbitsunion(T)
    # copy selector bytes
    error("Not implemented")
  end
  return dest
end

## utilities

zeros(T::Type, dims...) = fill!(MtlArray{T}(undef, dims...), 0)
ones(T::Type, dims...) = fill!(MtlArray{T}(undef, dims...), 1)
zeros(dims...) = zeros(Float32, dims...)
ones(dims...) = Mtls(Float32, dims...)
fill(v, dims...) = fill!(MtlArray{typeof(v)}(undef, dims...), v)
fill(v, dims::Dims) = fill!(MtlArray{typeof(v)}(undef, dims...), v)

function Base.fill!(A::MtlArray{T}, val) where T
  B = convert(T, val)
  unsafe_fill!(A.dev, pointer(A), B, length(A))
  A
end


text:
module LargeVis
    using LightGraphs
    using LinearAlgebra
    using StatsBase

    include("grad.jl")

    export diff_likelihood
end


text:
@doc raw"""
```
set_seed_pixels!(labels::Array{Int,2}, img::AbstractArray; suppress_feedback::Bool=false)
```

A mutating function to allow for interactive selection of seed pixels to be used
with the segmentation algorithms. Mutates the labels variable to contain initial seeds.
Uses Makie to provide user interface.

# Details

The function allows for users to interactively select seed pixels and their regions.
This is achieved using mouse and keyboard interactions. This is described in more detail
below.

# Arguments

The function arguments are described in more detail below.

##  `labels`

An `Array{Int,2}` containing the labels to be assigned to each pixel. Should be
initialized to be the same size as image and passed containing only zeros or other
seeds.

##  `img`

An `AbstractArray` containing the image to be seeded.

##  `suppress_feedback`

A `Bool` keyword argument which controls if console feedback is to be returned.

# Controls

Left mouse button: select seeds.
Right mouse button: move image.
Scroll: zoom.
Left arrow: change region down.
Right arrow: change region up.

# Example

Select seeds in the "toucan" image.

```julia
using Images, PictureSegmentation
img = testimage("toucan")

#convert to RGB{Float64}
img=RGB{Float64}.(img)

#create array of seeds in itilized to 0
seeds = zeros(Int,axes(img))

#use set_seed_pixels to set the inital seeds
set_seed_pixels(clicks, img)

#you can now seed the image as you wish
```
"""
function set_seed_pixels!(labels::Array{Int,2}, img::AbstractArray; suppress_feedback::Bool=false)
    # setup initial scene
    scene = Scene()
    image!(img, show_axis = false)
    clicks = Node(Point2f0[])
    region = 1
    h, w = size(img)
    # place seed
    on(events(scene).mousedrag) do buttons
        if ispressed(scene, Mouse.left)
            pos = to_world(scene, Point2(scene.events.mouseposition[]))
            position = (round(Int, (pos[1] / w) * h), round(Int, (pos[2] / h) * w))
            if minimum(position) > 1 && position[1] < h-  1 && position[2] < w - 1
                push!(clicks, push!(clicks[], pos))
                labels[position[1], position[2]] = region
                if suppress_feedback == false
                    println("position, region: ", position, region)
                end
            end
        end
        return
    end
    # change region
    on(scene.events.keyboardbuttons) do buttons
        if ispressed(scene, Keyboard.right)
            region += 1
            if suppress_feedback == false
                println("region: ", region)
            end
        end
        if ispressed(scene, Keyboard.left)
            if region != 1
                region -= 1
                if suppress_feedback == false
                    println("region: ", region)
                end
            end
        end
        return
    end
    # display seeds
    scatter!(scene, clicks, color = :red, marker = '.', markersize = 3,show_axis = false)
end

@doc raw"""
```
set_seed_pixels!(labels::Array{Int,3}, img::AbstractArray, layer::Int = 1; suppress_feedback::Bool=false)
```

A mutating function to allow for interactive selection of seed pixels to be used
with segmentation algorithms. Mutates the labels variable to contain initial seeds.
Uses Makie to provide user interface.

# Details

The function allows for users to interactively select seed pixels and their regions.
This is achieved using mouse and keyboard interactions. This is described in more detail
below.

# Arguments

The function arguments are described in more detail below.

##  `labels`

An `Array{Int,3}` containing the labels to be assigned to each pixel. Should be
initialized to be the same size as the image and passed containing only zeros or other
seeds.

##  `img`

An `AbstractArray` containing the image to be seeded.

##  `layer`

An `Int` which controls the layer which is displayed initially. If not defined
then the first layer will be shown.

##  `suppress_feedback`

An `Bool` keyword argument which controls whether console feedback is to be returned.

# Controls

Left mouse button: select seeds.
Right mouse button: move image.
Scroll: zoom.
Left arrow: change region down.
Right arrow: change region up.
Up arrow: change layer up.
Down arrow: change layer down.

# Example

Select seeds in synthetic 3d image.

```julia
using Images, PictureSegmentation

# Create image of cone with multiple colours
img=Array{RGB{Float64},3}(undef,100,100,100)
for i in CartesianIndices(img)
    if (i[1]-50)^2 + (i[2]-50)^2 < (i[3]/4)^2
        img[i]=RGB(i[3]/100,0,1-(i[3]/100))
    else
        img[i]=RGB(0,0,0)
    end
end

# Create array of seeds in itilized to 0
seeds = zeros(Int,axes(img))

# Use set_seed_pixels to set the inital seeds
set_seed_pixels(clicks, img)

# You can now seed the image as you wish
```
"""
function set_seed_pixels!(labels::Array{Int,3}, img::AbstractArray, layer::Int = 1; suppress_feedback::Bool=false)
    # setup scene
    scene = Scene()
    clicks = Node(Point2f0[])
    region = 1
    h, w = size(img)
    image!(scene,img[:,:,layer], show_axis = false)
    # set seed
    on(events(scene).mousedrag) do buttons
        if ispressed(scene, Mouse.left)
            pos = to_world(scene, Point2(scene.events.mouseposition[]))
            position = (round(Int, (pos[1] / w) * h), round(Int, (pos[2] / h) * w))
            if minimum(position) > 1 && position[1] < h - 1 && position[2] < w - 1
                push!(clicks, push!(clicks[], pos))
                labels[position[1], position[2], layer] = region
                if suppress_feedback == false
                    println("position, layer, region: ", position, layer, region)
                end
            end
        end
        return
    end
    # change layer and region
    on(scene.events.keyboardbuttons) do buttons
        if ispressed(scene, Keyboard.right)
            region += 1
            if suppress_feedback == false
                println("region: ", region)
            end
        end
        if ispressed(scene, Keyboard.left)
            if region!=1
                region -= 1
                if suppress_feedback == false
                    println("region: ", region)
                end
            end
        end
        if ispressed(scene, Keyboard.up)
            if layer != length(img[1,1,:])
                layer += 1
                image!(scene,img[:,:,layer], show_axis = false)
                display(scene)
                if suppress_feedback == false
                    println("layer: ", layer)
                end
            end
        end
        if ispressed(scene, Keyboard.down)
            if layer != 1
                layer -= 1
                image!(scene,img[:,:,layer], show_axis = false)
                display(scene)
                if suppress_feedback == false
                    println("layer: ", layer)
                end
            end
        end
        return
    end
    scatter!(scene, clicks, color = :red, marker = '.', markersize = 3,show_axis = false)
end