text:
#
# Advent of Code
# Day 11: Radioisotope Thermoelectric Generators
#

println("Advent of Code")
println("Day 11: Radioisotope Thermoelectric Generators")
println("")

DEBUG = false

BIG = typemax(Int64)

# Test inputs
# The first floor contains a hydrogen-compatible microchip and a lithium-compatible microchip.
# The second floor contains a hydrogen generator.
# The third floor contains a lithium generator.
# The fourth floor contains nothing relevant.
#

# Test setup

MH = 1
MLi = 2
GH = -1
GLi = -2

floors = Vector{Vector{Int}}(4)
floors[1] = [MH, MLi]
floors[2] = [GH]
floors[3] = [GLi,]
floors[4] = []

elevator = [0,0]

ALLCHIPS = IntSet([MH, MLi])

# Part 1 inputs
# The first floor contains a promethium generator and a promethium-compatible microchip.
# The second floor contains a cobalt generator, a curium generator, a ruthenium generator, and a plutonium generator.
# The third floor contains a cobalt-compatible microchip, a curium-compatible microchip, a ruthenium-compatible microchip, and a plutonium-compatible microchip.
# The fourth floor contains nothing relevant.

# Part  1 setup

# MPr = 1
# GPr = -1
# MCo = 2
# GCo = -2
# MCu = 3
# GCu = -3
# MRu = 4
# GRu = -4
# MPl = 5
# GPl = -5
#
# floors = Vector{Vector{Int}}(4)
# floors[1] = [GPr, MPr]
# floors[2] = [GCo, GCu, GRu, GPl]
# floors[3] = [MCo, MCu, MRu, MPl]
# floors[4] = []
#
# elevator = [0,0]

elevatorFloor = 1
elevatorTrips = 0

# The goal is to get everything to the 4th floor with the fewest elevator moves
# But:
#  The elevator stops at each floor
#  The elevator can only hold two items at once
#  The elevator can only move if it has a generator or a chip
#  Chips get fried in the presence of a incompatible generator
#

# Brute force this.
# Somehow.

function changeConfiguration(unloadCounter1, unloadCounter2, loadCounter1, loadCounter2, startingFloors, startingElevator, elevatorFloor)

  newFloors = deepcopy(startingFloors)
  newElevator = deepcopy(startingElevator)
  startingElevatorLoad = length(startingElevator)
  startingFloorLoad = length(startingFloors[elevatorFloor])

  # println("Changing configuration")
  # println("Unload: $unloadCounter1, $unloadCounter2")
  # println("Load: $loadCounter1, $loadCounter2")
  # println("Floor: $elevatorFloor")
  # println("Starting Floor: ", startingFloors[elevatorFloor])
  # println("Starting Elevator: ", startingElevator)
  # println("")

  abort = false

  # wait a minute, we can't remove the same thing from the elevator
  # or load the same thing on to the elevator
  #
  # But doing nothing is an option.
  #

  if (unloadCounter1 == unloadCounter2) && unloadCounter1 != 0
    if DEBUG println("Abort: unloading the same thing twice ($unloadCounter1, $unloadCounter2)") end
    abort = true
  elseif (loadCounter1 == loadCounter2) && loadCounter1 != 0
    if DEBUG println("Abort: loading same thing twice ($loadCounter1, $loadCounter2)") end
    abort = true
  end

  if !abort
    # unload 1
    if unloadCounter1 > 0
      if newElevator[unloadCounter1] != 0
        push!(newFloors[elevatorFloor], newElevator[unloadCounter1])
        newElevator[unloadCounter1] = 0
      else
        if DEBUG println("Abort: Nothing to unload (1)") end
        abort = true
      end
    end

    # unload 2
    if unloadCounter2 > 0 && !abort
      if newElevator[unloadCounter2] != 0
        push!(newFloors[elevatorFloor], newElevator[unloadCounter2])
        newElevator[unloadCounter2] = 0
      else
        if DEBUG println("Abort: Nothing to unload (2)") end
        abort = true
      end
    end
  end

  if !abort
    # only load if there is space
    availableSpaces = IntSet([1,2])

    if newElevator[1] != 0
      delete!(availableSpaces, 1)
    end

    if newElevator[2] != 0
      delete!(availableSpaces, 2)
    end

    # load 1
    if loadCounter1 != 0
      if length(availableSpaces) > 0 && length(newFloors[elevatorFloor]) > 0
        elevatorSpot = pop!(availableSpaces)
        thing = startingFloors[elevatorFloor][loadCounter1]
        #println("Putting ", thing, " in the elevator spot $elevatorSpot")
        newElevator[elevatorSpot] = thing
        thingIndex = findin(newFloors[elevatorFloor], thing)[1]
        #println("Trying to remove ", thing, " from ", newFloors[elevatorFloor], " at index ", thingIndex)
        deleteat!(newFloors[elevatorFloor], thingIndex)
      else
        if DEBUG println("Abort: No room to load elevator (1)") end
        abort = true
      end
    end

    # load 2
    if !abort
      if loadCounter2 != 0
        if length(availableSpaces) > 0 && length(newFloors[elevatorFloor]) > 0
          elevatorSpot = pop!(availableSpaces)
          thing = startingFloors[elevatorFloor][loadCounter2]
          #println("Putting ", thing, " in the elevator spot $elevatorSpot")
          newElevator[elevatorSpot] = thing
          thingIndex = findin(newFloors[elevatorFloor], thing)[1]
          #println("Trying to remove ", thing, " from ", newFloors[elevatorFloor], " at index ", thingIndex)
          deleteat!(newFloors[elevatorFloor], thingIndex)
        else
          if DEBUG println("Abort: No room to load elevator (2)") end
          abort = true
        end
      end
    end

    if !abort
      # ok, we've loaded and unloaded the elevator, let's check to see if we have a valid configuration

      # first, the elevator needs something to make it go
      if newElevator[1] == 0 && newElevator[2] == 0
        if DEBUG println("Abort: Nothing in the elevator") end
        abort = true
      end

      # things in the elevator must be compatible
      if newElevator[1] != 0 && newElevator[2] != 0
        if newElevator[1] != -newElevator[2]
          if DEBUG println("Abort: Incompatible items in elevator (", newElevator[1],",",newElevator[2],")") end
          abort = true
        end
      end

      # things left on the floor must be compatible
      # which means no chips on this floor with an incompatible generator

      for item1 in newFloors[elevatorFloor]
        if item1 > 0 # it's a chip find all generators
          for item2 in newFloors[elevatorFloor]
            if item2 < 0 && item1 != -item2
              if DEBUG println("Abort: Incompatible items on floor (", item1,",",item2,")") end
              abort = true
              break
            end
          end
          if abort
            break
          end
        end
      end
    end
  end

  if abort
    isValidConfiguration = false
  else
    isValidConfiguration = true

    # println("Valid change")
    # println("Starting Floor: ", startingFloors[elevatorFloor])
    # println("Starting Elevator: ", startingElevator)
    # println("Unload: $unloadCounter1, $unloadCounter2")
    # println("Load: $loadCounter1, $loadCounter2")
    # println("Ending Floor: ", newFloors[elevatorFloor])
    # println("Ending Elevator: ", newElevator)
  end

  return (isValidConfiguration, newFloors, newElevator)
end

function printFloors(floors::Vector{Vector{Int}}, elevatorFloor::Int)
  for counter in 4:-1:1
    if counter == elevatorFloor
      line = " E "
    else
      line = "   "
    end

    for item in floors[counter]
      line *= string(item)
      line *= " "
    end
    println("$counter: $line")
  end
end

function findMinimumTrips(recursionDepth::Int, elevatorFloor::Int, startingFloors::Vector{Vector{Int}}, startingElevator::Vector{Int})

  minimumTrips = BIG

  recursionDepth += 1

  # println("Recursion depth: $recursionDepth")

  if recursionDepth > 15
    # too many steps, try something else
    recursionDepth -= 1
    return minimumTrips
  end

  # Current options include:
  #  * unload zero, one, or two things from the elevator
  #  * load zero, one, or two compatible things in the elevator
  # Then
  #  * move the elevator up or down

  newElevator = deepcopy(startingElevator)
  newFloors = deepcopy(floors)
  startingElevatorLoad = length(startingElevator)
  currentFloorLoad = length(startingFloors[elevatorFloor])

  found = false
  abort = false

  for unloadCounter1 in 0:2
    for unloadCounter2 in 0:2
      for loadCounter1 in 0:currentFloorLoad
        for loadCounter2 in 0:currentFloorLoad
          (isValidConfiguration, newFloors, newElevator) = changeConfiguration(unloadCounter1, unloadCounter2, loadCounter1, loadCounter2, startingFloors, startingElevator, elevatorFloor)

          if isValidConfiguration
            #
            # Before moving the elevator, check to see if we have all the chips on the top floor
            #

            # println("Valid change")
            # println("Starting Floor: ", startingFloors[elevatorFloor])
            # println("Starting Elevator: ", startingElevator)
            # println("Unload: $unloadCounter1, $unloadCounter2")
            # println("Load: $loadCounter1, $loadCounter2")
            # println("Ending Floor: ", newFloors[elevatorFloor])
            # println("Ending Elevator: ", newElevator)
            #
            # println("")
            #
            # printFloors(newFloors, elevatorFloor)
            #
            # println("")

            # found = true

            # for item in ALLCHIPS
            #   if !(item in newFloors[4])
            #     found = false
            #   end
            # end

            if length(newFloors[4]) == 2
              if MH in newFloors[4]
                if MLi in newFloors[4]
                  found = true
                end
              end
            end

            if found
              # VICTORY
              println("Found one!")
              exit()
            else
              # No victory yet, move again
              minTripsUp = BIG
              minTripsDown = BIG

              println("\u1b[1H")
              printFloors(newFloors, elevatorFloor)

              if elevatorFloor < 4
                minTripsUp = findMinimumTrips(recursionDepth, elevatorFloor + 1, newFloors, newElevator)
              end

              if elevatorFloor > 1
                minTripsDown = findMinimumTrips(recursionDepth, elevatorFloor - 1, newFloors, newElevator)
              end

              currentTrips = min(minTripsUp, minTripsDown) + 1

              if currentTrips > minimumTrips
                abort = true
              else
                minimumTrips = currentTrips
              end
            end
          end
          if found || abort break end
        end
        if found || abort break end
      end
      if found || abort break end
    end
    if found || abort break end
  end

  recursionDepth -= 1
  return minimumTrips
end

recursionDepth = 0

overallMinimumTrips = findMinimumTrips(recursionDepth, elevatorFloor, floors, elevator)

println("The minimum number of trips is: $overallMinimumTrips")


text:
module TestWithBase

using ..Utils: scan_functions_with_kwargs
using Test

function test_with_base(scan!, op, xs)
    @test scan!(op, copy(xs)) == accumulate(op, xs)
end

function test_with_base(scan!)
    @testset "length(xs) = $n" for n in [2^8, 2^8 + 1]
        xs = rand(1:2^10, n)
        @testset "$op" for op in [+, &, xor]
            test_with_base(scan!, op, xs)
        end
    end
end

function test_with_base()
    @testset "$scan!" for scan! in scan_functions_with_kwargs()
        test_with_base(scan!)
    end
end

end  # module

ext: jl


text:
using Router, ChirpsController

route("/") do
  # Router.serve_static_file("welcome.html")
  respond_with_html("home", "home_page.md")
end

route("/chirps", ChirpsController.index)
route("/chirps/new", ChirpsController.new)
route("/chirps", ChirpsController.create, method = POST)

max_stars_repo_name: essenciary/chirper

max_issues_repo_name: essenciary/chirper

max_forks_repo_name: essenciary/chirper


text:
export WeightedExplorer

using Random
using StatsBase: sample, Weights

"""
    WeightedExplorer(;is_normalized::Bool)

`is_normalized` is used to indicate if the feeded action values
are alrady normalized to have a sum of `1.0`.
"""
struct WeightedExplorer{T,R<:AbstractRNG} <: AbstractExplorer
    rng::R
end

function WeightedExplorer(; is_normalized::Bool = false, seed = nothing)
    rng = MersenneTwister(seed)
    WeightedExplorer{is_normalized,typeof(rng)}(rng)
end

(s::WeightedExplorer{true})(values::AbstractVector{T}) where {T} =
    sample(s.rng, Weights(values, one(T)))
(s::WeightedExplorer{false})(values::AbstractVector) =
    sample(s.rng, Weights(values, sum(values)))


text:
function getFileName(a1::CXFile,)
    c = CXString()
    ccall( (:wci_getFileName, libwci), Void, (CXFile,Ptr{Void},), a1, c.data,)
    return get_string(c)
end
function getIncludedFile(cu::InclusionDirective)
    getFileName(ccall( (:wci_getIncludedFile, libwci),
            Ptr{Void},
            (Ptr{UInt8},),
            pointer(cu.data)))
end
function getNullLocation()
    c = CXSourceLocation()
    ccall( (:wci_getNullLocation, libwci), Void, (Ptr{Void},),    c.data,)
    return c
end
function equalLocations(a1::CXSourceLocation,a2::CXSourceLocation,)
    ccall( (:wci_equalLocations, libwci), UInt32, (Ptr{Void},Ptr{Void},), a1.data,a2.data, )
end
function getNullRange()
    c = CXSourceRange()
    ccall( (:wci_getNullRange, libwci), Void, (Ptr{Void},),    c.data,)
    return c
end
function getRange(a1::CXSourceLocation,a2::CXSourceLocation,)
    c = CXSourceRange()
    ccall( (:wci_getRange, libwci), Void, (Ptr{Void},Ptr{Void},Ptr{Void},), a1.data,a2.data, c.data,)
    return c
end
function equalRanges(a1::CXSourceRange,a2::CXSourceRange,)
    ccall( (:wci_equalRanges, libwci), UInt32, (Ptr{Void},Ptr{Void},), a1.data,a2.data, )
end
function Range_isNull(a1::CXSourceRange,)
    ccall( (:wci_Range_isNull, libwci), Int32, (Ptr{Void},), a1.data, )
end
function getNullCursor()
    c = TmpCursor()
    ccall( (:wci_getNullCursor, libwci), Void, (Ptr{Void},),    c.data,)
    return CXCursor(c)
end
function getTranslationUnitCursor(a1::CXTranslationUnit,)
    c = TmpCursor()
    ccall( (:wci_getTranslationUnitCursor, libwci), Void, (CXTranslationUnit,Ptr{Void},), a1, c.data,)
    return CXCursor(c)
end
function equalCursors(a1::CLCursor,a2::CLCursor,)
    ccall( (:wci_equalCursors, libwci), UInt32, (Ptr{Void},Ptr{Void},), a1.data,a2.data, )
end
function Cursor_isNull(a1::CLCursor,)
    ccall( (:wci_Cursor_isNull, libwci), Int32, (Ptr{Void},), a1.data, )
end
function hashCursor(a1::CLCursor,)
    ccall( (:wci_hashCursor, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function getCursorKind(a1::CLCursor,)
    ccall( (:wci_getCursorKind, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function isDeclaration(a1::CXCursorKind,)
    ccall( (:wci_isDeclaration, libwci), UInt32, (CXCursorKind,), a1, )
end
function isReference(a1::CXCursorKind,)
    ccall( (:wci_isReference, libwci), UInt32, (CXCursorKind,), a1, )
end
function getCursorLinkage(a1::CLCursor,)
    ccall( (:wci_getCursorLinkage, libwci), CXLinkageKind, (Ptr{Void},), a1.data, )
end
function getCursorAvailability(a1::CLCursor,)
    ccall( (:wci_getCursorAvailability, libwci), CXAvailabilityKind, (Ptr{Void},), a1.data, )
end
function getCursorLanguage(a1::CLCursor,)
    ccall( (:wci_getCursorLanguage, libwci), CXLanguageKind, (Ptr{Void},), a1.data, )
end
function getCursorSemanticParent(a1::CLCursor,)
    c = TmpCursor()
    ccall( (:wci_getCursorSemanticParent, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXCursor(c)
end
function getCursorLexicalParent(a1::CLCursor,)
    c = TmpCursor()
    ccall( (:wci_getCursorLexicalParent, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXCursor(c)
end
function getCursorType(a1::CLCursor,)
    c = TmpType()
    ccall( (:wci_getCursorType, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXType(c)
end
function getTypedefDeclUnderlyingType(a1::CLCursor,)
    c = TmpType()
    ccall( (:wci_getTypedefDeclUnderlyingType, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXType(c)
end
function getEnumDeclIntegerType(a1::CLCursor,)
    c = TmpType()
    ccall( (:wci_getEnumDeclIntegerType, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXType(c)
end
function getEnumConstantDeclValue(a1::CLCursor,)
    ccall( (:wci_getEnumConstantDeclValue, libwci), Int64, (Ptr{Void},), a1.data, )
end
function getEnumConstantDeclUnsignedValue(a1::CLCursor,)
    ccall( (:wci_getEnumConstantDeclUnsignedValue, libwci), UInt64, (Ptr{Void},), a1.data, )
end
function Cursor_getNumArguments(a1::CLCursor,)
    ccall( (:wci_Cursor_getNumArguments, libwci), Int32, (Ptr{Void},), a1.data, )
end
function Cursor_getArgument(a1::CLCursor,a2::Int32,)
    c = TmpCursor()
    ccall( (:wci_Cursor_getArgument, libwci), Void, (Ptr{Void},Int32,Ptr{Void},), a1.data,a2, c.data,)
    return CXCursor(c)
end
function equalTypes(a1::CLType,a2::CLType,)
    ccall( (:wci_equalTypes, libwci), UInt32, (Ptr{Void},Ptr{Void},), a1.data,a2.data, )
end
function getCanonicalType(a1::CLType,)
    c = TmpType()
    ccall( (:wci_getCanonicalType, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXType(c)
end
function isConstQualifiedType(a1::CLType,)
    ccall( (:wci_isConstQualifiedType, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function isVolatileQualifiedType(a1::CLType,)
    ccall( (:wci_isVolatileQualifiedType, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function isRestrictQualifiedType(a1::CLType,)
    ccall( (:wci_isRestrictQualifiedType, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function getPointeeType(a1::CLType,)
    c = TmpType()
    ccall( (:wci_getPointeeType, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXType(c)
end
function getTypeDeclaration(a1::CLType,)
    c = TmpCursor()
    ccall( (:wci_getTypeDeclaration, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXCursor(c)
end
function getDeclObjCTypeEncoding(a1::CLCursor,)
    c = CXString()
    ccall( (:wci_getDeclObjCTypeEncoding, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return get_string(c)
end
function getTypeKindSpelling(a1::CXTypeKind,)
    c = CXString()
    ccall( (:wci_getTypeKindSpelling, libwci), Void, (CXTypeKind,Ptr{Void},), a1, c.data,)
    return get_string(c)
end
function getFunctionTypeCallingConv(a1::CLType,)
    ccall( (:wci_getFunctionTypeCallingConv, libwci), Int32, (Ptr{Void},), a1.data, )
end
function getResultType(a1::CLType,)
    c = TmpType()
    ccall( (:wci_getResultType, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXType(c)
end
function getNumArgTypes(a1::CLType,)
    ccall( (:wci_getNumArgTypes, libwci), Int32, (Ptr{Void},), a1.data, )
end
function getArgType(a1::CLType,a2::UInt32,)
    c = TmpType()
    ccall( (:wci_getArgType, libwci), Void, (Ptr{Void},UInt32,Ptr{Void},), a1.data,a2, c.data,)
    return CXType(c)
end
function isFunctionTypeVariadic(a1::CLType,)
    ccall( (:wci_isFunctionTypeVariadic, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function getCursorResultType(a1::CLCursor,)
    c = TmpType()
    ccall( (:wci_getCursorResultType, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXType(c)
end
function isPODType(a1::CLType,)
    ccall( (:wci_isPODType, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function getElementType(a1::CLType,)
    c = TmpType()
    ccall( (:wci_getElementType, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXType(c)
end
function getNumElements(a1::CLType,)
    ccall( (:wci_getNumElements, libwci), Int64, (Ptr{Void},), a1.data, )
end
function getArrayElementType(a1::CLType,)
    c = TmpType()
    ccall( (:wci_getArrayElementType, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXType(c)
end
function getArraySize(a1::CLType,)
    ccall( (:wci_getArraySize, libwci), Int64, (Ptr{Void},), a1.data, )
end
function isVirtualBase(a1::CLCursor,)
    ccall( (:wci_isVirtualBase, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function getCXXAccessSpecifier(a1::CLCursor,)
    ccall( (:wci_getCXXAccessSpecifier, libwci), Int64, (Ptr{Void},), a1.data, )
end
function getNumOverloadedDecls(a1::CLCursor,)
    ccall( (:wci_getNumOverloadedDecls, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function getOverloadedDecl(a1::CLCursor,a2::UInt32,)
    c = TmpCursor()
    ccall( (:wci_getOverloadedDecl, libwci), Void, (Ptr{Void},UInt32,Ptr{Void},), a1.data,a2, c.data,)
    return CXCursor(c)
end
function getCursorUSR(a1::CLCursor,)
    c = CXString()
    ccall( (:wci_getCursorUSR, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return get_string(c)
end
function getCursorSpelling(a1::CLCursor,)
    c = CXString()
    ccall( (:wci_getCursorSpelling, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return get_string(c)
end
function getCursorDisplayName(a1::CLCursor,)
    c = CXString()
    ccall( (:wci_getCursorDisplayName, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return get_string(c)
end
function getCursorReferenced(a1::CLCursor,)
    c = TmpCursor()
    ccall( (:wci_getCursorReferenced, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXCursor(c)
end
function getCursorDefinition(a1::CLCursor,)
    c = TmpCursor()
    ccall( (:wci_getCursorDefinition, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXCursor(c)
end
function isCursorDefinition(a1::CLCursor,)
    ccall( (:wci_isCursorDefinition, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function getCanonicalCursor(a1::CLCursor,)
    c = TmpCursor()
    ccall( (:wci_getCanonicalCursor, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXCursor(c)
end
function CXXMethod_isStatic(a1::CLCursor,)
    ccall( (:wci_CXXMethod_isStatic, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function CXXMethod_isVirtual(a1::CLCursor,)
    ccall( (:wci_CXXMethod_isVirtual, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function CXXMethod_isProtected(a1::CLCursor,)
    ccall( (:wci_CXXMethod_isProtected, libwci), UInt32, (Ptr{Void},), a1.data, )
end
function getTemplateCursorKind(a1::CLCursor,)
    ccall( (:wci_getTemplateCursorKind, libwci), Int32, (Ptr{Void},), a1.data, )
end
function getSpecializedCursorTemplate(a1::CLCursor,)
    c = TmpCursor()
    ccall( (:wci_getSpecializedCursorTemplate, libwci), Void, (Ptr{Void},Ptr{Void},), a1.data, c.data,)
    return CXCursor(c)
end
function getTokenKind(a1::CXToken,)
    ccall( (:wci_getTokenKind, libwci), Int32, (Ptr{Void},), a1.data, )
end
function getTokenSpelling(a1::CXTranslationUnit,a2::CXToken,)
    c = CXString()
    ccall( (:wci_getTokenSpelling, libwci), Void, (CXTranslationUnit,Ptr{Void},Ptr{Void},), a1,a2.data, c.data,)
    return get_string(c)
end
function getCursorKindSpelling(a1::CXCursorKind,)
    c = CXString()
    ccall( (:wci_getCursorKindSpelling, libwci), Void, (CXCursorKind,Ptr{Void},), a1, c.data,)
    return get_string(c)
end
function getClangVersion()
    c = CXString()
    ccall( (:wci_getClangVersion, libwci), Void, (Ptr{Void},),    c.data,)
    return get_string(c)
end
function getCursorExtent(a1::CLCursor)
    c = CXSourceRange()
    ccall( (:wci_getCursorExtent, libwci), Void, (Ptr{Void},Ptr{Void}), a1.data, c.data)
    return c
end
function disposeTokens(tu::CXTranslationUnit, tokens::Ptr{_CXToken}, numtokens::Cuint)
    ccall( (:clang_disposeTokens, :libclang), Void, (Ptr{Void}, Ptr{_CXToken}, Cuint),
           tu, tokens, numtokens)
end
function tokenize(tu::CXTranslationUnit, sr::CXSourceRange)
    tokens = Ptr{_CXToken}[C_NULL]
    numtokens = ccall( (:wci_tokenize, libwci), Cuint, (Ptr{Void}, Ptr{Void}, Ptr{_CXToken}),
            tu, sr.data, pointer(tokens))
    return TokenList(tokens[1], convert(Cuint, numtokens - 1), tu)
end
function Cursor_getTranslationUnit(cu::CLCursor)
    ccall( (:wci_Cursor_getTranslationUnit, libwci), CXTranslationUnit, (Ptr{Void},), cu.data)
end


text:
#Quantities: A Julian way to deal with quantities of arbitrary units
#-------------------------------------------------------------------------------

module Quantities

include("base.jl")
include("show.jl")

#==Exported symbols
===============================================================================#
export Length, Dist, Position

#Useful typealiases (sadly, names conflict with constructors):
export _Length, _Dist, _Position

#Accessor functions:
export value #High-collision WARNING: other modules probably want to export "value"

end

#Last line


text:
"""Abstract type of all possible patterns.
"""

abstract type Pattern end

struct CreateInfo <: Pattern
    sdef::SDefinition
end

struct Create <: Pattern
    f::FDefinition
    type

    function Create(decl)
        type = remove_parts(decl.name, discarded_parts=[1, 2])
        new(decl, type)
    end
end

struct Command <: Pattern
    f::FDefinition
end

struct Finalizer <: Pattern
    f::FDefinition
    x
    function Finalizer(decl)
        x = remove_parts(decl.name, discarded_parts=[1, 2], convert_to=SnakeCaseLower)
        new(decl, x)
    end
end

const patterns_dict = OrderedDict(
    (decl -> startswith(decl.name, "create_")) => Create,
    (decl -> occursin("CreateInfo", decl.name)) => CreateInfo,
    (decl -> startswith(decl.name, "vkCmd")) => Command,
    (decl -> startswith(decl.name, "vkDestroy")) => Finalizer,
)

function patterns(decl::Declaration)
    patterns_found = []
    for (ispattern, pat) ∈ patterns_dict
        ispattern(decl) ? push!(patterns_found, pat(decl)) : nothing
    end
    patterns_found
end


text:
# draw histgram of sensor data

module DrawHistogram
    using Plots, DataFrames, CSV
    pyplot()

    function main(is_test=false)
        data_path = joinpath(split(@__FILE__, "src")[1], "data/sensor_data_200.txt")
        df_200_mm = CSV.read(data_path, DataFrame,
                             header=["date", "time", "ir", "lidar"],
                             delim=' ')

        bin_min_max = maximum(df_200_mm.lidar) - minimum(df_200_mm.lidar)
        histogram(df_200_mm.lidar, bins=bin_min_max, label="histogram")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/freq_dist/histogram/histogram_200_mm.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
function sample(range::UnitRange{Int64}, n; nc::Bool = false)
    N = range.stop - range.start + 2
    result = Array{Int64}(undef, n)
    m = 1
    ncm = 1
    t = 1
    if nc
        nc_result = Array{Int64}(undef, N-n-1)
        while m <= n
            if (N-t) * rand() >= (n - m+1)
                nc_result[ncm] = t
                ncm+=1
            else
                result[m] = t
                m+=1
            end
            t+=1
        end
        nc_result[ncm:end] = collect(t:range.stop)
        return result, nc_result
    else
        while m <= n
            if (N-t) * rand() < (n - m+1)
                result[m] = t
                m+=1
            end
            t+=1
        end

        return result
    end
end


text:
@testset "solve_diffusion.jl" begin
    domain = ClosedIntervalDomain(0.0, 1.0)
    a, b = domain.left, domain.right

    T = 1.0

    # Exact solutions (from Borgaard 2006)
    u(ξ, t) = t + sin(2π * ξ) + sin(8π * ξ)
    u_int(ξ, t) = t * ξ - 1 / 2π * cos(2π * ξ) - 1 / 8π * cos(8π * ξ)
    f(ξ, t) = 1 + 4π^2 * sin(2π * ξ) + 64π^2 * sin(8π * ξ)
    g_a(t) = t
    g_b(t) = t

    tols = (; abstol = 1e-6, reltol = 1e-4)

    # Number of mesh points
    NN = floor.(Int, 10 .^ LinRange(2.5, 4.0, 4))

    err_ref = 0
    err_bar_ref = 0
    err_allbar_ref = 0
    n_ref = 0
    for (i, N) ∈ enumerate(NN)
        M = N

        # Discretization
        x = discretize(domain, M)
        ξ = discretize(domain, N)
        Δx = (domain.right - domain.left) / N

        # Filter
        h(x) = Δx / 2

        # Exact filtered solution
        function ū(x, t)
            α = max(a, x - h(x))
            β = min(b, x + h(x))
            1 / (β - α) * (u_int(β, t) - u_int(α, t))
        end

        # Equations
        equation = DiffusionEquation(domain, IdentityFilter(), f, g_a, g_b)
        equation_filtered = DiffusionEquation(domain, TopHatFilter(h), f, g_a, g_b)

        # Solve discretized problem
        sol = solve(
            equation,
            ξ -> u(ξ, 0.0),
            (0.0, T),
            M,
            N;
            method = "discretizefirst",
            boundary_conditions = "derivative",
            tols...,
        )

        # Solve filtered-then-discretized problem
        @test_broken sol_bar = solve(
            equation_filtered,
            ξ -> u(ξ, 0.0),
            (0.0, T),
            M,
            N;
            method = "filterfirst",
            boundary_conditions = "derivative",
            tols...,
        )

        # Solve discretized-then-filtered problem
        sol_allbar = solve(
            equation_filtered,
            ξ -> u(ξ, 0.0),
            (0.0, T),
            M,
            N;
            method = "discretizefirst",
            boundary_conditions = "derivative",
            tols...,
        )

        ## Relative error
        u_exact = u.(ξ, T)
        ū_exact = ū.(x, T)
        err = norm(sol(T) - u_exact) / norm(u_exact)
        # err_bar = norm(sol_bar(T) - ū_exact) / norm(ū_exact)
        err_allbar = norm(sol_allbar(T) - ū_exact) / norm(ū_exact)

        if i == 1
            err_ref = err
            # err_bar_ref = err_bar
            err_allbar_ref = err_allbar
            n_ref = N
        end

        # Test for quadratic convergence
        @test err < 1.6 * err_ref * (n_ref / N)^2
        # @test err_bar < 1.6 * err_bar_ref * (n_ref / N)^2
        @test err_allbar < 1.6 * err_allbar_ref * (n_ref / N)^2
    end

    # Test exact boundary conditions
    N = 100
    M = N
    equation = DiffusionEquation(domain, IdentityFilter(), f, g_a, g_b)

    # Solve discretized problem with DAE formulation
    sol_exact = solve(
        equation,
        ξ -> u(ξ, 0.0),
        (0.0, T),
        M,
        N;
        method = "discretizefirst",
        boundary_conditions = "exact",
        tols...,
    )

    # Solve discretized problem with derivative BC
    sol_derivative = solve(
        equation,
        ξ -> u(ξ, 0.0),
        (0.0, T),
        M,
        N;
        method = "discretizefirst",
        boundary_conditions = "derivative",
        tols...,
    )

    @test sol_exact(T) ≈ sol_derivative(T)
end


text:
module AdvancedVI

using Random: AbstractRNG

using Distributions, DistributionsAD, Bijectors
using DocStringExtensions

using ProgressMeter, LinearAlgebra

using ForwardDiff
using Tracker

const PROGRESS = Ref(true)
function turnprogress(switch::Bool)
    @info("[AdvancedVI]: global PROGRESS is set as $switch")
    PROGRESS[] = switch
end

const DEBUG = Bool(parse(Int, get(ENV, "DEBUG_ADVANCEDVI", "0")))

include("ad.jl")

using Requires
function __init__()
    @require Flux="587475ba-b771-5e3f-ad9e-33799f191a9c" begin
        apply!(o, x, Δ) = Flux.Optimise.apply!(o, x, Δ)
        Flux.Optimise.apply!(o::TruncatedADAGrad, x, Δ) = apply!(o, x, Δ)
        Flux.Optimise.apply!(o::DecayedADAGrad, x, Δ) = apply!(o, x, Δ)
    end
    @require Zygote = "e88e6eb3-aa80-5325-afca-941959d7151f" begin
        include("compat/zygote.jl")
        export ZygoteAD

        function AdvancedVI.grad!(
            vo,
            alg::VariationalInference{<:AdvancedVI.ZygoteAD},
            q,
            model,
            θ::AbstractVector{<:Real},
            out::DiffResults.MutableDiffResult,
            args...
        )
            f(θ) = if (q isa Distribution)
                - vo(alg, update(q, θ), model, args...)
            else
                - vo(alg, q(θ), model, args...)
            end
            y, back = Zygote.pullback(f, θ)
            dy = first(back(1.0))
            DiffResults.value!(out, y)
            DiffResults.gradient!(out, dy)
            return out
        end
    end
    @require ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" begin
        include("compat/reversediff.jl")
        export ReverseDiffAD

        function AdvancedVI.grad!(
            vo,
            alg::VariationalInference{<:AdvancedVI.ReverseDiffAD{false}},
            q,
            model,
            θ::AbstractVector{<:Real},
            out::DiffResults.MutableDiffResult,
            args...
        )
            f(θ) = if (q isa Distribution)
                - vo(alg, update(q, θ), model, args...)
            else
                - vo(alg, q(θ), model, args...)
            end
            tp = AdvancedVI.tape(f, θ)
            ReverseDiff.gradient!(out, tp, θ)
            return out
        end
    end
end

export
    vi,
    ADVI,
    ELBO,
    elbo,
    TruncatedADAGrad,
    DecayedADAGrad,
    VariationalInference

abstract type VariationalInference{AD} end

getchunksize(::Type{<:VariationalInference{AD}}) where AD = getchunksize(AD)
getADtype(::VariationalInference{AD}) where AD = AD

abstract type VariationalObjective end

const VariationalPosterior = Distribution{Multivariate, Continuous}

"""
    grad!(vo, alg::VariationalInference, q, model::Model, θ, out, args...)

Computes the gradients used in `optimize!`. Default implementation is provided for
`VariationalInference{AD}` where `AD` is either `ForwardDiffAD` or `TrackerAD`.
This implicitly also gives a default implementation of `optimize!`.

Variance reduction techniques, e.g. control variates, should be implemented in this function.
"""
function grad! end

"""
    vi(model, alg::VariationalInference)
    vi(model, alg::VariationalInference, q::VariationalPosterior)
    vi(model, alg::VariationalInference, getq::Function, θ::AbstractArray)

Constructs the variational posterior from the `model` and performs the optimization
following the configuration of the given `VariationalInference` instance.

# Arguments
- `model`: `Turing.Model` or `Function` z ↦ log p(x, z) where `x` denotes the observations
- `alg`: the VI algorithm used
- `q`: a `VariationalPosterior` for which it is assumed a specialized implementation of the variational objective used exists.
- `getq`: function taking parameters `θ` as input and returns a `VariationalPosterior`
- `θ`: only required if `getq` is used, in which case it is the initial parameters for the variational posterior
"""
function vi end

function update end

# default implementations
function grad!(
    vo,
    alg::VariationalInference{<:ForwardDiffAD},
    q,
    model,
    θ::AbstractVector{<:Real},
    out::DiffResults.MutableDiffResult,
    args...
)
    f(θ_) = if (q isa Distribution)
        - vo(alg, update(q, θ_), model, args...)
    else
        - vo(alg, q(θ_), model, args...)
    end

    chunk_size = getchunksize(typeof(alg))
    # Set chunk size and do ForwardMode.
    chunk = ForwardDiff.Chunk(min(length(θ), chunk_size))
    config = ForwardDiff.GradientConfig(f, θ, chunk)
    ForwardDiff.gradient!(out, f, θ, config)
end

function grad!(
    vo,
    alg::VariationalInference{<:TrackerAD},
    q,
    model,
    θ::AbstractVector{<:Real},
    out::DiffResults.MutableDiffResult,
    args...
)
    θ_tracked = Tracker.param(θ)
    y = if (q isa Distribution)
        - vo(alg, update(q, θ_tracked), model, args...)
    else
        - vo(alg, q(θ_tracked), model, args...)
    end
    Tracker.back!(y, 1.0)

    DiffResults.value!(out, Tracker.data(y))
    DiffResults.gradient!(out, Tracker.grad(θ_tracked))
end

"""
    optimize!(vo, alg::VariationalInference{AD}, q::VariationalPosterior, model::Model, θ; optimizer = TruncatedADAGrad())

Iteratively updates parameters by calling `grad!` and using the given `optimizer` to compute
the steps.
"""
function optimize!(
    vo,
    alg::VariationalInference,
    q,
    model,
    θ::AbstractVector{<:Real};
    optimizer = TruncatedADAGrad()
)
    # TODO: should we always assume `samples_per_step` and `max_iters` for all algos?
    alg_name = alg_str(alg)
    samples_per_step = alg.samples_per_step
    max_iters = alg.max_iters

    num_params = length(θ)

    # TODO: really need a better way to warn the user about potentially
    # not using the correct accumulator
    if (optimizer isa TruncatedADAGrad) && (θ ∉ keys(optimizer.acc))
        # this message should only occurr once in the optimization process
        @info "[$alg_name] Should only be seen once: optimizer created for θ" objectid(θ)
    end

    diff_result = DiffResults.GradientResult(θ)

    i = 0
    prog = if PROGRESS[]
        ProgressMeter.Progress(max_iters, 1, "[$alg_name] Optimizing...", 0)
    else
        0
    end

    # add criterion? A running mean maybe?
    time_elapsed = @elapsed while (i < max_iters) # & converged
        grad!(vo, alg, q, model, θ, diff_result, samples_per_step)

        # apply update rule
        Δ = DiffResults.gradient(diff_result)
        Δ = apply!(optimizer, θ, Δ)
        @. θ = θ - Δ

        AdvancedVI.DEBUG && @debug "Step $i" Δ DiffResults.value(diff_result)
        PROGRESS[] && (ProgressMeter.next!(prog))

        i += 1
    end

    return θ
end

# objectives
include("objectives.jl")

# optimisers
include("optimisers.jl")

# VI algorithms
include("advi.jl")

end # module


text:
struct AddCostSpec
    variable_type::Type
    component_type::Type
    has_status_variable::Bool
    has_status_parameter::Bool
    sos_status::SOS_STATUS_VARIABLE
    multiplier::Float64
    variable_cost::Union{Nothing, Function}
    start_up_cost::Union{Nothing, Function}
    shut_down_cost::Union{Nothing, Function}
    fixed_cost::Union{Nothing, Function}
    addtional_linear_terms::Dict{String, Symbol}
end

function AddCostSpec(;
    variable_type,
    component_type,
    has_status_variable = false,
    has_status_parameter = false,
    sos_status = NO_VARIABLE,
    multiplier = OBJECTIVE_FUNCTION_POSITIVE,
    variable_cost = nothing,
    start_up_cost = nothing,
    shut_down_cost = nothing,
    fixed_cost = nothing,
    addtional_linear_terms = Dict{String, Symbol}(),
)
    return AddCostSpec(
        variable_type,
        component_type,
        has_status_variable,
        has_status_parameter,
        sos_status,
        multiplier,
        variable_cost,
        start_up_cost,
        shut_down_cost,
        fixed_cost,
        addtional_linear_terms,
    )
end

function AddCostSpec(
    ::Type{<:T},
    ::Type{<:U},
    ::PSIContainer,
) where {T <: PSY.Component, U <: AbstractDeviceFormulation}
    error("AddCostSpec is not implemented for $T / $U")
end

set_addtional_linear_terms!(spec::AddCostSpec, key, value) =
    spec.addtional_linear_terms[key] = value

function add_service_variables!(spec::AddCostSpec, services)
    for service in services
        name = PSY.get_name(service)
        set_addtional_linear_terms!(spec, name, make_variable_name(name, typeof(service)))
    end
    return
end

"""
Add variables to the PSIContainer for a service.
"""
function cost_function!(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{T},
    ::DeviceModel{T, U},
    ::Type{<:PM.AbstractPowerModel},
    feedforward::Union{Nothing, AbstractAffectFeedForward} = nothing,
) where {T <: PSY.Component, U <: AbstractDeviceFormulation}
    for d in devices
        spec = AddCostSpec(T, U, psi_container)
        @debug T, spec
        services = PSY.get_services(d)
        add_service_variables!(spec, services)
        add_to_cost!(psi_container, spec, PSY.get_operation_cost(d), d)
    end
    return
end

function add_to_cost_expression!(
    psi_container::PSIContainer,
    cost_expression::JuMP.AbstractJuMPScalar,
)
    T_ce = typeof(cost_expression)
    T_cf = typeof(psi_container.cost_function)
    if T_cf <: JuMP.GenericAffExpr && T_ce <: JuMP.GenericQuadExpr
        psi_container.cost_function += cost_expression
    else
        JuMP.add_to_expression!(psi_container.cost_function, cost_expression)
    end
    return
end

function has_on_variable(
    psi_container::PSIContainer,
    ::Type{T};
    variable_type = OnVariable,
) where {T <: PSY.Component}
    #get_variable can't be used because the default behavior is to error if variables is not present
    return !isnothing(get(
        psi_container.variables,
        make_variable_name(variable_type, T),
        nothing,
    ))
end

function has_on_parameter(psi_container::PSIContainer, ::Type{T}) where {T <: PSY.Component}
    if !model_has_parameters(psi_container)
        return false
    end
    return !isnothing(get(
        psi_container.parameters,
        encode_symbol(OnVariable, string(T)),
        nothing,
    ))
end

function _get_pwl_vars_container(psi_container::PSIContainer)
    if !haskey(psi_container.variables, :PWL_cost_vars)
        time_steps = model_time_steps(psi_container)
        contents = Dict{Tuple{String, Int, Int}, Any}()
        container = JuMP.Containers.SparseAxisArray(contents)
        assign_variable!(psi_container, :PWL_cost_vars, container)
    else
        container = get_variable(psi_container, :PWL_cost_vars)
    end
    return container
end

function slope_convexity_check(slopes::Vector{Float64})
    flag = true
    for ix in 1:(length(slopes) - 1)
        if slopes[ix] > slopes[ix + 1]
            @debug slopes
            return flag = false
        end
    end
    return flag
end

@doc raw"""
Returns True/False depending on compatibility of the cost data with the linear implementation method

Returns ```flag```

# Arguments

* cost_ : container for quadratic and linear factors
"""
function pwlparamcheck(cost_)
    slopes = PSY.get_slopes(cost_)
    # First element of the return is the average cost at P_min.
    # Shouldn't be passed for convexity check
    return slope_convexity_check(slopes[2:end])
end

function linear_gen_cost!(
    psi_container::PSIContainer,
    var_name::Symbol,
    component_name::String,
    linear_term::Float64,
    time_period::Int,
)
    resolution = model_resolution(psi_container)
    dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR
    variable = get_variable(psi_container, var_name)[component_name, time_period]
    gen_cost = sum(variable) * linear_term
    add_to_cost_expression!(psi_container, gen_cost * dt)
    return
end

@doc raw"""
Returns piecewise cost expression using SOS Type-2 implementation for psi_container model.

# Equations

``` variable = sum(sos_var[i]*cost_data[2][i])```

``` gen_cost = sum(sos_var[i]*cost_data[1][i]) ```

# LaTeX

`` variable = (sum_{i\in I} c_{2, i} sos_i) ``

`` gen_cost = (sum_{i\in I} c_{1, i} sos_i) ``

Returns ```gen_cost```

# Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* variable::JuMP.Containers.DenseAxisArray{JV} : variable array
* cost_data::PSY.VariableCost{NTuple{2, Float64}} : container for quadratic and linear factors
"""
function pwl_gencost_sos!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    component_name::String,
    cost_data::Vector{NTuple{2, Float64}},
    time_period::Int,
)
    var_name = make_variable_name(spec.variable_type, spec.component_type)
    variable = get_variable(psi_container, var_name)[component_name, time_period]
    settings_ext = get_ext(get_settings(psi_container))
    export_pwl_vars = get_export_pwl_vars(psi_container.settings)
    @debug export_pwl_vars
    total_gen_cost = JuMP.AffExpr(0.0)

    if spec.sos_status == NO_VARIABLE
        bin = 1.0
        @debug("Using Piecewise Linear cost function but no variable/parameter ref for ON status is passed. Default status will be set to online (1.0)")
    elseif spec.sos_status == PARAMETER
        param_key = encode_symbol(OnVariable, string(spec.component_type))
        bin =
            get_parameter_container(psi_container, param_key).parameter_array[component_name]
        @debug("Using Piecewise Linear cost function with parameter $(param_key)")
    elseif spec.sos_status == VARIABLE
        var_key = make_variable_name(OnVariable, spec.component_type)
        bin = get_variable(psi_container, var_key)[component_name, time_period]
        @debug("Using Piecewise Linear cost function with variable $(var_key)")
    else
        @assert false
    end

    gen_cost = JuMP.AffExpr(0.0)
    pwlvars = Array{JuMP.VariableRef}(undef, length(cost_data))
    for i in 1:length(cost_data)
        pwlvars[i] = JuMP.@variable(
            psi_container.JuMPmodel,
            base_name = "{$(variable)}_{sos}",
            start = 0.0,
            lower_bound = 0.0,
            upper_bound = 1.0
        )
        if export_pwl_vars
            container = _get_pwl_vars_container(psi_container)
            container[(component_name, time_period, i)] = pwlvars[i]
        end
        JuMP.add_to_expression!(gen_cost, cost_data[i][1] * pwlvars[i])
    end
    JuMP.@constraint(psi_container.JuMPmodel, sum(pwlvars) == bin)
    JuMP.@constraint(
        psi_container.JuMPmodel,
        pwlvars in MOI.SOS2(collect(1:length(pwlvars)))
    )
    JuMP.@constraint(
        psi_container.JuMPmodel,
        variable == sum([var_ * cost_data[ix][2] for (ix, var_) in enumerate(pwlvars)])
    )
    JuMP.add_to_expression!(total_gen_cost, gen_cost)

    return total_gen_cost
end

@doc raw"""
Returns piecewise cost expression using linear implementation for psi_container model.

# Equations

``` 0 <= pwl_var[i] <= (cost_data[2][i] - cost_data[2][i-1])```

``` variable = sum(pwl_var[i])```

``` gen_cost = sum(pwl_var[i]*cost_data[1][i]/cost_data[2][i]) ```

# LaTeX
`` 0 <= pwl_i <= (c_{2, i} - c_{2, i-1})``

`` variable = (sum_{i\in I} pwl_i) ``

`` gen_cost = (sum_{i\in I}  pwl_i) c_{1, i}/c_{2, i} ``

Returns ```gen_cost```

# Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* variable::JuMP.Containers.DenseAxisArray{JV} : variable array
* cost_data::Vector{NTuple{2, Float64}} : container for quadratic and linear factors
"""
function pwl_gencost_linear!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    component_name::String,
    cost_data::Vector{NTuple{2, Float64}},
    time_period::Int,
)
    var_name = make_variable_name(spec.variable_type, spec.component_type)
    variable = get_variable(psi_container, var_name)[component_name, time_period]
    settings_ext = get_ext(get_settings(psi_container))
    export_pwl_vars = get_export_pwl_vars(psi_container.settings)
    @debug export_pwl_vars
    total_gen_cost = JuMP.AffExpr(0.0)

    gen_cost = JuMP.AffExpr(0.0)
    total_gen = JuMP.AffExpr(0.0)
    for i in 1:length(cost_data)
        pwlvar = JuMP.@variable(
            psi_container.JuMPmodel,
            base_name = "{$(variable)}_{pwl}",
            start = 0.0,
            lower_bound = 0.0,
            upper_bound = PSY.get_breakpoint_upperbounds(cost_data)[i]
        )
        if export_pwl_vars
            container = _get_pwl_vars_container(psi_container)
            container[(component_name, time_period, i)] = pwlvar
        end
        JuMP.add_to_expression!(gen_cost, PSY.get_slopes(cost_data)[i] * pwlvar)
        JuMP.add_to_expression!(total_gen, pwlvar)
    end
    JuMP.@constraint(psi_container.JuMPmodel, variable == total_gen)
    JuMP.add_to_expression!(total_gen_cost, gen_cost)
    return total_gen_cost
end

"""
Adds to the models costs represented by PowerSystems TwoPart costs
"""
function add_to_cost!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    cost_data::PSY.TwoPartCost,
    component::PSY.Component,
)
    component_name = PSY.get_name(component)
    time_steps = model_time_steps(psi_container)
    @debug "TwoPartCost" component_name
    if !(spec.variable_cost === nothing)
        variable_cost = spec.variable_cost(cost_data)
        for t in time_steps
            variable_cost!(psi_container, spec, component_name, variable_cost, t)
        end
    else
        @warn "No variable cost defined for $component_name"
    end

    if !isnothing(spec.fixed_cost) && spec.has_status_variable
        @debug "Fixed cost" component_name
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                make_variable_name(OnVariable, spec.component_type),
                component_name,
                spec.fixed_cost,
                t,
            )
        end
    end
    return
end

"""
Adds to the models costs represented by PowerSystems ThreePart costs
"""
function add_to_cost!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    cost_data::PSY.ThreePartCost,
    component::PSY.Component,
)
    component_name = PSY.get_name(component)
    @debug "ThreePartCost" component_name
    resolution = model_resolution(psi_container)
    dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR
    variable_cost = PSY.get_variable(cost_data)
    time_steps = model_time_steps(psi_container)
    for t in time_steps
        variable_cost!(psi_container, spec, component_name, variable_cost, t)
    end

    if !isnothing(spec.start_up_cost)
        @debug "Start up cost" component_name
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                make_variable_name(StartVariable, spec.component_type),
                component_name,
                spec.start_up_cost(cost_data) * spec.multiplier,
                t,
            )
        end
    end

    if !isnothing(spec.shut_down_cost)
        @debug "Shut down cost" component_name
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                make_variable_name(StopVariable, spec.component_type),
                component_name,
                spec.shut_down_cost(cost_data) * spec.multiplier,
                t,
            )
        end
    end

    if !isnothing(spec.fixed_cost) && spec.has_status_variable
        @debug "Fixed cost" component_name
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                make_variable_name(OnVariable, spec.component_type),
                component_name,
                spec.fixed_cost(cost_data) * spec.multiplier,
                t,
            )
        end
    end

    return
end

"""
Adds to the models costs represented by PowerSystems Multi-Start costs
"""
function add_to_cost!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    cost_data::PSY.MultiStartCost,
    component::PSY.Component,
)
    component_name = PSY.get_name(component)
    resolution = model_resolution(psi_container)
    dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR
    time_steps = model_time_steps(psi_container)

    if !isnothing(spec.fixed_cost) && spec.has_status_variable
        @debug "Fixed cost" component_name
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                make_variable_name(OnVariable, spec.component_type),
                component_name,
                spec.fixed_cost(cost_data) * spec.multiplier,
                t,
            )
        end
    end

    if !isnothing(spec.shut_down_cost)
        @debug "Shut down cost" component_name
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                make_variable_name(StopVariable, spec.component_type),
                component_name,
                spec.shut_down_cost(cost_data) * spec.multiplier,
                t,
            )
        end
    end

    # Original implementation had SOS by default
    variable_cost_data = PSY.get_cost(PSY.get_variable(cost_data))
    if !all(iszero.(last.(variable_cost_data)))
        for t in time_steps
            gen_cost =
                pwl_gencost_sos!(psi_container, spec, component_name, variable_cost_data, t)
            add_to_cost_expression!(psi_container, spec.multiplier * gen_cost * dt)
        end
    else
        @debug "No Variable Cost associated with $(component_name)"
    end

    # Start-up costs
    start_cost_data = PSY.get_start_up(cost_data)
    for (st, var_type) in
        enumerate((HotStartVariable, WarmStartVariable, ColdStartVariable))
        var_name = make_variable_name(var_type, spec.component_type)
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                var_name,
                component_name,
                start_cost_data[st] * spec.multiplier,
                t,
            )
        end
    end

    return
end

"""
Adds to the models costs represented by PowerSystems Market-Bid costs
"""
function add_to_cost!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    cost_data::PSY.MarketBidCost,
    component::PSY.Component,
)
    component_name = PSY.get_name(component)
    @debug "Market Bid" component_name
    resolution = model_resolution(psi_container)
    dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR
    time_steps = model_time_steps(psi_container)
    initial_time = model_initial_time(psi_container)
    # TODO: Use mthods from PowerSystems to eliminate this step
    variable_cost_forecast = get_time_series(psi_container, component, "variable_cost")
    variable_cost_forecast = map(PSY.VariableCost, variable_cost_forecast)
    for t in time_steps
        variable_cost!(psi_container, spec, component_name, variable_cost_forecast[t], t)
    end

    if !isnothing(spec.start_up_cost)
        @debug "Start up cost" component_name
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                make_variable_name(StartVariable, spec.component_type),
                component_name,
                spec.start_up_cost(cost_data) * spec.multiplier,
                t,
            )
        end
    end

    if spec.has_status_variable
        @debug "no_load cost" component_name
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                make_variable_name(OnVariable, spec.component_type),
                component_name,
                PSY.get_no_load(cost_data) * spec.multiplier,
                t,
            )
        end
    end

    if !isnothing(spec.shut_down_cost)
        @debug "Shut down cost" component_name
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                make_variable_name(StopVariable, spec.component_type),
                component_name,
                spec.shut_down_cost(cost_data) * spec.multiplier,
                t,
            )
        end
    end

    #Service Cost Bid
    ancillary_services = PSY.get_ancillary_services(cost_data)
    for service in ancillary_services
        add_service_bid_cost!(psi_container, spec, component, service)
    end
    return
end

function add_service_bid_cost!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    component::PSY.Component,
    service::PSY.Service,
)
    return
end

function add_service_bid_cost!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    component::PSY.Component,
    service::PSY.Reserve{T},
) where {T <: PSY.ReserveDirection}
    # TODO: Use mthods from PowerSystems to eliminate this step
    forecast_data = get_time_series(psi_container, component, PSY.get_name(service))
    forecast_data = map(PSY.VariableCost, forecast_data)
    time_steps = model_time_steps(psi_container)
    if eltype(forecast_data) == PSY.VariableCost{Float64}
        for t in time_steps
            linear_gen_cost!(
                psi_container,
                spec.addtional_linear_terms[PSY.get_name(service)],
                PSY.get_name(component),
                forecast_data,
                t,
            )
        end
    else
        error("Current version only supports linear cost bid for services, please change the forecast data for $(PSY.get_name(service))")
    end
    return
end

function add_service_bid_cost!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    component::PSY.Component,
    service::PSY.ReserveDemandCurve{T},
) where {T <: PSY.ReserveDirection}
    error("Current version doesn't supports cost bid for  ReserveDemandCurve services, please change the forecast data for $(PSY.get_name(service))")
    return
end

@doc raw"""
Adds to the cost function cost terms for sum of variables with common factor to be used for cost expression for psi_container model.

    # Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* var_name::Symbol: The variable name
* component_name::String: The component_name of the variable container
* cost_component::PSY.VariableCost{Float64} : container for cost to be associated with variable
"""
function variable_cost!(
    ::PSIContainer,
    ::AddCostSpec,
    component_name::String,
    ::Nothing,
    ::Int,
)
    @debug "Empty Variable Cost" component_name
    return
end

@doc raw"""
Adds to the cost function cost terms for sum of variables with common factor to be used for cost expression for psi_container model.

    # Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* var_name::Symbol: The variable name
* component_name::String: The component_name of the variable container
* cost_component::PSY.VariableCost{Float64} : container for cost to be associated with variable
"""
function variable_cost!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    component_name::String,
    cost_component::PSY.VariableCost{Float64},
    time_period::Int,
)
    @debug "Linear Variable Cost" component_name
    var_name = make_variable_name(spec.variable_type, spec.component_type)
    cost_data = PSY.get_cost(cost_component)
    linear_gen_cost!(
        psi_container,
        var_name,
        component_name,
        cost_data * spec.multiplier,
        time_period,
    )
    return
end

@doc raw"""
Adds to the cost function cost terms for sum of variables with common factor to be used for cost expression for psi_container model.

# Equation

``` gen_cost = dt*sign*(sum(variable.^2)*cost_data[1] + sum(variable)*cost_data[2]) ```

# LaTeX

`` cost = dt\times sign (sum_{i\in I} c_1 v_i^2 + sum_{i\in I} c_2 v_i ) ``

for quadratic factor large enough. If the first term of the quadratic objective is 0.0, adds a
linear cost term `sum(variable)*cost_data[2]`

# Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* var_name::Symbol: The variable name
* component_name::String: The component_name of the variable container
* cost_component::PSY.VariableCost{NTuple{2, Float64}} : container for quadratic and linear factors
"""
function variable_cost!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    component_name::String,
    cost_component::PSY.VariableCost{NTuple{2, Float64}},
    time_period::Int,
)
    var_name = make_variable_name(spec.variable_type, spec.component_type)
    cost_data = PSY.get_cost(cost_component)
    if cost_data[1] >= eps()
        @debug "Quadratic Variable Cost" component_name
        resolution = model_resolution(psi_container)
        dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR
        variable = get_variable(psi_container, var_name)[component_name, time_period]
        gen_cost = sum(variable .^ 2) * cost_data[1] + sum(variable) * cost_data[2]
        add_to_cost_expression!(psi_container, spec.multiplier * gen_cost * dt)
    else
        @debug "Quadratic Variable Cost with only linear term" component_name
        linear_gen_cost!(
            psi_container,
            var_name,
            component_name,
            cost_data[2] * spec.multiplier,
            time_period,
        )
    end
    return
end

@doc raw"""
Creates piecewise linear cost function using a sum of variables and expression with sign and time step included.

# Expression

```JuMP.add_to_expression!(gen_cost, c)```

Returns sign*gen_cost*dt

# LaTeX

``cost = sign\times dt \sum_{v\in V} c_v``

where ``c_v`` is given by

`` c_v = \sum_{i\in Ix} \frac{y_i - y_{i-1}}{x_i - x_{i-1}} v^{p.w.}_i ``

# Arguments

* psi_container::PSIContainer : the psi_container model built in PowerSimulations
* var_name::Symbol: The variable name
* component_name::String: The component_name of the variable container
* cost_component::PSY.VariableCost{Vector{NTuple{2, Float64}}}
"""
function variable_cost!(
    psi_container::PSIContainer,
    spec::AddCostSpec,
    component_name::String,
    cost_component::PSY.VariableCost{Vector{NTuple{2, Float64}}},
    time_period::Int,
)
    @debug "PWL Variable Cost" component_name
    resolution = model_resolution(psi_container)
    dt = Dates.value(Dates.Second(resolution)) / SECONDS_IN_HOUR
    # If array is full of tuples with zeros return 0.0
    cost_data = PSY.get_cost(cost_component)
    if all(iszero.(last.(cost_data)))
        @debug "All cost terms for component $(component_name) are 0.0"
        return JuMP.AffExpr(0.0)
    end

    var_name = make_variable_name(spec.variable_type, spec.component_type)
    if !pwlparamcheck(cost_component)
        @warn("The cost function provided for $(var_name) device is not compatible with a linear PWL cost function.
        An SOS-2 formulation will be added to the model.
        This will result in additional binary variables added to the model.")
        gen_cost =
            pwl_gencost_sos!(psi_container, spec, component_name, cost_data, time_period)
    else
        gen_cost =
            pwl_gencost_linear!(psi_container, spec, component_name, cost_data, time_period)
    end
    add_to_cost_expression!(psi_container, spec.multiplier * gen_cost * dt)
    return
end


text:
module Sesam

using CSV
using DataFrames
using JuMP
using Gurobi

# utility functions
include("util.jl")

# power plants
include("powerplants.jl")

# renewables
include("renewables.jl")

# nodes
include("nodes.jl")

# lines
include("lines.jl")

### models
# dcopf
include("dcopf.jl")

# scopf
include("scopf.jl")

# economic dispatch
include("economic_dispatch.jl")

# congestion management
include("congestion_management.jl")

# congestion management with PtG
include("congestion_management_ptg.jl")

export
    CSV,
    DataFrames,
    JuMP,
    Gurobi,
    ProgressMeter,
    ElectronDisplay,
    Nodes,
    Lines,
    PowerPlants,
    Renewables,
    DCOPF,
    SCOPF,
    EconomicDispatch,
    CongestionManagement,
    CongestionManagementPtG
end


text:
package java.io;

public class ByteArrayOutputStream extends OutputStream {
    private byte[] bytes;

    public ByteArrayOutputStream() {
	bytes = new byte[1];
    }

    public synchronized byte[] toByteArray() {
	return bytes;
    }

    public synchronized void write(int i) {
	bytes[0] = (byte)i;
    }
}

max_stars_repo_name: HarvardPL/cryptoerase

max_issues_repo_name: HarvardPL/cryptoerase

max_forks_repo_name: HarvardPL/cryptoerase


text:
using Test
using TestSetExtensions

using SIAN
using ModelingToolkit

@info "Testing started"

@testset "All the tests" begin
    @includetests ARGS
end

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
##
##  m i s c . j l
##

#-- Algebraic-geometric mean
function agm(a::Number, b::Number)
    local a0 = a, b0 = b
    local c  = max(abs(a), abs(b))
    if typeof(a) <: Real && a < 0; a0 = complex(a, 0); end
    if typeof(b) <: Real && b < 0; b0 = complex(b, 0); end

    while abs(b0-a0) >= eps(c)
        a0, b0 = (a0 + b0)/2.0, sqrt(a0 * b0)
    end
    return (a0 + b0) / 2.0
end

#-- Arc length of a curve defined through f:R --> R
function arc_length(f::Function, a::Real, b::Real;
                        nmax::Int=25, tol::Real=2.5e-15)
    local k::Int
    h = b - a;
    fa = f(a); fb = f(b)
    A = zeros(nmax, nmax)
    A[1, 1] = sqrt(sum((fb - fa)^2))
    for k = 1:(nmax-1)
        h /= 2
        x = linspace(a, b, 2^k + 1)
        y = map(f, x)
        X = [x y]
        dX = diff(X)
        A[k+1, 1] = sum(sqrt(sum(dX.^2, 2)))
        for l in 1:k
            A[k+1, l+1] = (4^l * A[k+1, l] - A[k, l]) / (4^l - 1)
        end
        if abs(A[k+1, k+1] - A[k, k]) < tol && k > 3
            break
        end
    end
    err = abs(A[k+1, k+1] - A[k, k])
    return A[k+1, k+1], k+1, err
end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libtestu01"], :libtestu01),
    LibraryProduct(prefix, String["libprobdist"], :libprobdist),
    LibraryProduct(prefix, String["libtestu01extractors"], :libtestu01extractors),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/andreasnoack/TestU01Builder/releases/download/v1.2.3-0"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/TestU01.aarch64-linux-gnu.tar.gz", "205e8212a2da9e4d2d9cc9a400ffa7766062ec2dc245a076c01e2c9e51184228"),
    Linux(:aarch64, :musl) => ("$bin_prefix/TestU01.aarch64-linux-musl.tar.gz", "842ca0cc5f5b3edc40da6d84ff12689c0554012251a99edd96892ab953cd191d"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/TestU01.arm-linux-gnueabihf.tar.gz", "069dbc418562b385125dd05289c807443869c6dc5fa731bbb3dd7f57d302b2ca"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/TestU01.arm-linux-musleabihf.tar.gz", "2f0ff74623bcd5fbff1f2fa4833376f42a235a407da2ac1192b71a0f1eeff4f9"),
    Linux(:i686, :glibc) => ("$bin_prefix/TestU01.i686-linux-gnu.tar.gz", "616971bc49f3f7066fa77b7890f94750ba1666c10130c3e5b93b2dea9c6cf3b3"),
    Linux(:i686, :musl) => ("$bin_prefix/TestU01.i686-linux-musl.tar.gz", "db381b8e8193f90c3485e514dc5d20682b3941ddebbb0fdfb8aa6fc246866e19"),
    Windows(:i686) => ("$bin_prefix/TestU01.i686-w64-mingw32.tar.gz", "1500b9b1115cdbd07c1c28b744bf5ec3cfc815caff84dacfc816d99b93432c64"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/TestU01.powerpc64le-linux-gnu.tar.gz", "64ba48e9fad7c26c267ed338bb47edf9a1020c854cf80f7e31a238d51689e80a"),
    MacOS(:x86_64) => ("$bin_prefix/TestU01.x86_64-apple-darwin14.tar.gz", "95b73e9dd271db69f14401474363e45ec401a34a608207276c4e03d2e4c5320b"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/TestU01.x86_64-linux-gnu.tar.gz", "187f70dbc7293d47b5a9bc586eec0ebd99475205155a94045a15ded1b48ce40f"),
    Linux(:x86_64, :musl) => ("$bin_prefix/TestU01.x86_64-linux-musl.tar.gz", "dc1e796ad5dce302d61507fe896f911aa3e17148739f955e195b84fe6709271a"),
    FreeBSD(:x86_64) => ("$bin_prefix/TestU01.x86_64-unknown-freebsd11.1.tar.gz", "65d0ff7a7e5e30781aa5bcc4fd16e8b51eddec10fa2555ce0eb2b759b1333827"),
    Windows(:x86_64) => ("$bin_prefix/TestU01.x86_64-w64-mingw32.tar.gz", "8609cb84cafbb4866cb9c80c5c654a0188e986b750376c761bf153dc8fdd63d4"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Function to run the MGEO Real.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Lopes, I. M. L (2013). Uma abordagem multi-objetivo para a otimização de
#       trajetórias de uma vela solar (INPE-15384-TDI/1403). PhD Thesis,
#       Instituto Nacional de Pesquisas Espaciais, São José dos Campos, SP,
#       Brazil (in Portuguese).
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export mgeo_run

"""
    mgeo_run(mgeod::MGEO_Structure{Nv, Nf, Design_Variable_MGEO_Real}, f_obj::Function, show_debug::Bool = false) where {Nv, Nf}

Run the MGEO Real configured in `mgeod` using the objective functions `f_obj`.
The algorithm was based in [1].

# Keywords

* `show_debug`: Print debug information (**Default** = `true`).

# Returns

The Pareto frontier, which is an array of `Pareto_Point`.

# Remarks

The objective function must have the following signature

        function f_obj(vars)
            1) Compute f[1], f[2], ..., f[nf] using vars[1], vars[2], ... var[n]
            2) return (valid_point, f)
        end

where `nf` is the number of objective functions, `n` is the number of design
variables, and `valid_point` is a boolean value that indicates if vars yield to
a valid point.

"""
function mgeo_run(mgeod::MGEO_Structure{Nv, Nf, Design_Variable_MGEO_Real},
                  f_obj::Function,
                  show_debug::Bool = false) where {Nv, Nf}
    # Pareto frontier.
    pareto = Pareto_Point{Nv,Nf}[]

    # Maximum number of generations per run.
    ngen_max_per_run = floor(mgeod.ngen_max/mgeod.run_max)

    # Vector to store the population.
    vars = zeros(MVector{Nv,Float64})

    # Initialization of Pareto Frontier
    # =================================

    # Initialize the string with one valid point in the Pareto Frontier.
    @inbounds for run = 1:mgeod.run_max
        # Sample a new population using a uniform distribution between the
        # minimum and maximum allowed range for each design variable.
        vars .= map(x->x.min + rand()*(x.max-x.min), mgeod.design_vars)

        # Call the objective functions.
        (valid, f) = f_obj(vars)

        # Check if the number of objective functions are valid.
        (length(f) != Nf) && @error("The number of objective functions returned by f_obj is not equal to the one specified in mgeod.")

        if valid
            # Add the point to the Pareto Frontier.
            push!(pareto, Pareto_Point{Nv, Nf}(SVector(vars),
                                               SVector{Nf, Float64}(f)))
            break
        end
    end

    # Check if a valid point was found.
    if (length(pareto) == 0)
        return pareto
    end

    # Loop - Independent Runs
    # =======================

    # List of all points created after the modification of population.
    candidate_points = Vector{Pareto_Point{Nv,Nf}}(undef,Nv)

    # Array to sort the candidate points.
    f_rank = Vector{sRank_Real}(undef, Nv)

    @inbounds for run = 1:mgeod.run_max
        if show_debug
            println("--------------------------------------------------------------------------------")
            println("RUN = $run")
            println("")
            println("Generations ($ngen_max_per_run):")
            print("    0%")
        end

        # Sample a new string if it is not the first run.
        if run > 1
            vars .= map(x->x.min + rand()*(x.max-x.min), mgeod.design_vars)
        end

        # Loop - MGEO Generations
        # =======================

        genDebug = 0

        @inbounds for ngen_per_run = 1:ngen_max_per_run
            if (show_debug)
                genDebug += 1

                if (genDebug >= 0.1*ngen_max_per_run)
                    genDebug = 0
                    print("...$(ngen_per_run/ngen_max_per_run*100)%")
                end
            end

            # Choose which objective function will be used to compute the
            # adaptability and to assemble the rank.
            chosen_func = rand(1:Nf)

            # Evaluate the objective functions using parallel computing.
            Threads.@threads for i=1:Nv
                # Modify the i-th design variable.
                @inbounds if !mgeod.design_vars[i].perc_perturbation
                    new_var_i = vars[i] + mgeod.design_vars[i].σ*randn()
                else
                    new_var_i = vars[i] + (mgeod.design_vars[i].σ*vars[i])*randn()
                end

                # Make sure that the new value lies within the defined interval.
                @inbounds if new_var_i < mgeod.design_vars[i].min
                    new_var_i = mgeod.design_vars[i].min
                elseif new_var_i > mgeod.design_vars[i].max
                    new_var_i = mgeod.design_vars[i].max
                end

                vars_i              = copy(vars)
                @inbounds vars_i[i] = new_var_i

                # Evaluate the objective functions.
                (valid, f) = f_obj(vars_i)

                # Check if the solution is valid.
                if valid
                    # Create the candidate point.
                    candidate_point =
                        Pareto_Point{Nv, Nf}(SVector(vars_i),
                                             SVector{Nf, Float64}(f))

                    # Add the result to the rank.
                    @inbounds f_rank[i] = sRank_Real(true,
                                                     i,
                                                     vars_i[i],
                                                     f[chosen_func])
                    @inbounds candidate_points[i] = candidate_point
                else
                    @inbounds f_rank[i] = sRank_Real(false, i, vars_i[i], 0.0)
                end
            end

            # Add the points to the Pareto frontier.
            num_valid_points = 0

            @inbounds for i=1:length(candidate_points)
                if f_rank[i].valid
                    num_valid_points += 1
                    check_dominance!(pareto, candidate_points[i], mgeod)
                end
            end

            # If there are no valid points, an independent run must be called.
            (num_valid_points == 0) && break

            # Ranking
            # ========
            #
            # The points will be ranked according to the selected objective
            # function (`chosen_func`).
            #
            # Notice that the invalid points will be placed after the valid
            # ones.

            sort!(f_rank, lt=(a,b)->begin
                      # Check if data is valid.
                      if (a.valid && b.valid)
                          return a.f < b.f
                      elseif (a.valid)
                          return true
                      else
                          return false
                      end
                  end
                 )

            # Choose a variable to be changed for the new generation.
            var_accepted = false

            while !var_accepted
                # Sample a valid point.
                v_sample = rand(1:num_valid_points)

                # Accept the change with probability r^(-τ), where r is the
                # rank of the sampled variable.
                Pk = v_sample^(-mgeod.τ)

                if rand() <= Pk
                    # The bit is accepted, then exit the loop.
                    @inbounds vars[f_rank[v_sample].index] =
                        f_rank[v_sample].var_value
                    var_accepted = true
                end
            end
        end

        if show_debug
            println("")
            println("")
            println("Pareto frontier information:")
            println("    Number of points: $(length(pareto))")
            println("")
        end
    end

    pareto
end


text:
using Jecco
using LinearAlgebra
using SparseArrays
using Plots

source(x,y) = exp(-x^2 - y^2) * (-4 + 3 * (x^2 + y^2) + 4 * x^2 * y^2)

# returns axx Dxx + ayy Dyy + axy Dxy + bx Dx + by Dy + cc. note that this
# function overwrites the input matrices to save memory
function build_operator(Dxx::SparseMatrixCSC, Dyy::SparseMatrixCSC, Dxy::SparseMatrixCSC,
                        Dx::SparseMatrixCSC, Dy::SparseMatrixCSC,
                        axx::Diagonal, ayy::Diagonal, axy::Diagonal,
                        bx::Diagonal, by::Diagonal, cc::Diagonal)
    Dxx_ = similar(Dxx)
    Dyy_ = similar(Dyy)
    Dxy_ = similar(Dxy)
    Dx_  = similar(Dx)
    Dy_  = similar(Dy)

    mul!(Dxx_, axx, Dxx)
    mul!(Dyy_, ayy, Dyy)
    mul!(Dxy_, axy, Dxy)
    mul!(Dx_, bx, Dx)
    mul!(Dy_, by, Dy)

    Dxx_ + Dyy_ + Dxy_ + Dx_ + Dy_ + cc
end

#=
use the Kronecker product (kron) to build the 2-dimensional derivation matrices
from the 1-dimensional ones. see for instance:

  https://en.wikipedia.org/wiki/Kronecker_product

  https://arxiv.org/pdf/1801.01483.pdf (section 5)
=#
function deriv_operators(hx, hy, Nx::Int, Ny::Int, ord::Int)
    Dx_op  = CenteredDiff{1}(1, ord, hx, Nx)
    Dxx_op = CenteredDiff{1}(2, ord, hx, Nx)

    Dy_op  = CenteredDiff{2}(1, ord, hy, Ny)
    Dyy_op = CenteredDiff{2}(2, ord, hy, Ny)

    Dx  = kron(I(Ny), SparseMatrixCSC(Dx_op))
    Dxx = kron(I(Ny), SparseMatrixCSC(Dxx_op))
    Dy  = kron(SparseMatrixCSC(Dy_op), I(Nx))
    Dyy = kron(SparseMatrixCSC(Dyy_op), I(Nx))

    Dx, Dy, Dxx, Dyy, Dx * Dy
end

x_min    = -5.0
x_max    =  5.0
x_nodes  =  128
y_min    = -5.0
y_max    =  5.0
y_nodes  =  64

ord = 4

xcoord  = Cartesian{1}("x", x_min, x_max, x_nodes, endpoint=false)
ycoord  = Cartesian{2}("y", y_min, y_max, y_nodes, endpoint=false)

hx = Jecco.delta(xcoord)
hy = Jecco.delta(ycoord)

Nx = xcoord.nodes
Ny = ycoord.nodes

f_exact = [exp(-xcoord[i]^2 - ycoord[j]^2) for i in 1:Nx, j in 1:Ny]

Dx, Dy, Dxx, Dyy, Dxy = deriv_operators(hx, hy, Nx, Ny, ord)

f0    = zeros(Nx,Ny)
ind2D = LinearIndices(f0)

M = Nx * Ny

b_vec = zeros(M)

axx     = ones(M)
ayy     = ones(M)
axy     = zeros(M)
bx      = zeros(M)
by      = zeros(M)
cc      = zeros(M)

for j in 1:Ny, i in 1:Nx
    idx = ind2D[i,j]

    xi  = xcoord[i]
    yi  = ycoord[j]

    axy[idx] = xi * yi
    bx[idx]  = xi
    by[idx]  = yi
    cc[idx]  = xi^2 + yi^2

    b_vec[idx] = source(xi, yi)
end

# build operator A = Dxx + Dyy + x y Dxy + x Dx + y Dy + (x^2 + y^2)
A_mat = build_operator(Dxx, Dyy, Dxy, Dx, Dy, Diagonal(axx), Diagonal(ayy),
                       Diagonal(axy), Diagonal(bx), Diagonal(by), Diagonal(cc))

A_fact = factorize(A_mat)
sol    = A_fact \ b_vec

@inbounds for idx in eachindex(f0)
    f0[idx] = sol[idx]
end

j_slice = div(Ny,2) + 1
x = xcoord[:]

plot(x, f_exact[:,j_slice])
scatter!(x, f0[:,j_slice])


text:
using Base.Test

@testset "RAD" begin
end


text:
function get_derivatives(u::Vector{Float64},p::Vector{Float64})
    # derivatives: dF/d[bifurcation_param]
    dFdp::Vector{Float64} = zeros(V.NUM)

    Wee1p = 1 - u[V.Wee1] - u[V.Wee1pp]
    Cdc25p = 1 - u[V.Cdc25] - u[V.Cdc25pp]

    dFdp[V.Subp] = - (1/p[C.Kd])*(p[C.kcBc1Sub]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*(p[C.SubT]-u[V.Subp])
    dFdp[V.PP1] = (1/p[C.Kd])*(p[C.kipp1] + p[C.kipp1C]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*u[V.PP1]
    dFdp[V.Gwlp] = - (1/p[C.Kd])*(p[C.kcBc1G]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*(p[C.Gwtot] - u[V.Gwlp])
    dFdp[V.Wee1] = (1/p[C.Kd])*(p[C.kcBc1W1]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*u[V.Wee1]
    dFdp[V.Wee1pp] = - (1/p[C.Kd])*(p[C.kcBc1W1]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*Wee1p
    dFdp[V.Cdc25] = (1/p[C.Kd])*(p[C.kcBc125]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*u[V.Cdc25]
    dFdp[V.Cdc25pp] = - (1/p[C.Kd])*(p[C.kcBc125]*u[V.CycBCdk1]/(1 + (p[C.InhCDK]/p[C.Kd]))^2)*(Cdc25p)

    return dFdp
end

max_stars_repo_name: himoto/Bifurcation

max_issues_repo_name: himoto/Bifurcation

max_forks_repo_name: himoto/Bifurcation


text:
function BayesianBlocks(M,ndim)

    results = zeros(size(M))

    for i in 1:ndim

        if length(unique(M[i,:]))<3
            res_temp = zeros(length(M[i,:]))

            for u in 1:length(unique(M[i,:]))
                res_temp[M[i,:] .== unique(M[i,:])[u]] .= u
            end

            results[i,:] = res_temp
        else

            disc = LinearDiscretizer(binedges(DiscretizeBayesianBlocks(), M[i,:]))
            results[i,:] = [encode(disc,one) for one in M[i,:]]

        end

    end

    return(DataFrame(results))

end

function BayesianBlocks_cutpoints(M,ndim)

    results = DataFrame()

    for i in 1:ndim

        if length(unique(M[i,:]))>3
            cuts = binedges(DiscretizeBayesianBlocks(), M[i,:])
            disc = LinearDiscretizer(cuts)
            results[!,i] = [encode(disc,one) for one in M[i,:]]

            results[!,i] = [encode(disc,one) for one in M[i,:]]
            results[!,i] = float(results[!,i])

            for j in 1:length(cuts)
                results[!,i][results[!,i].==j] .= cuts[j]
            end

        end

    end

    return(results)

end


text:
type AssignmentModel <: RedeployModel
    model::Gurobi.Model
    lambda::Float64

    hosp2stn::Matrix{Float64}
    stn2stn::Matrix{Float64}

    assignment::Vector{Int} # which location the ambulance is assigned to
    ambulances::Vector{Vector{Int}} # list of ambulances assigned to each location
    status::Vector{Symbol} # the current status of the ambulance
        # possible statuses: :available, :responding, :atscene, :conveying, :returning
    fromtime::Vector{Int} # the time it started the new status
    hospital::Vector{Int} # the hospital the ambulance is at (0 otherwise)

    soln::Vector{Float64} # buffer for storing dynamic assignment solutions
end

function AssignmentModel(
        p::DeploymentProblem,
        available::Vector{Int},
        # utilization::Vector{Float64},
        hospitals::DataFrame,
        stations::DataFrame;
        lambda::Float64 = 100.0
    )

    nambulances = sum(available)
    nlocations = length(available)
    assignment = zeros(Int, nambulances)

    hosp2stn = convert(Matrix{Float64},
        hcat([hospitals[Symbol("stn$(i)_min")] for i in 1:nlocations]...)
    )
    stn2stn =  convert(Matrix{Float64},
        hcat([stations[Symbol("stn$(i)_min")] for i in 1:nlocations]...)
    )

    k = 1
    ambulances = [Int[] for i in 1:nlocations]
    for i in eachindex(available), j in 1:available[i]
        assignment[k] = i
        push!(ambulances[i], k)
        k += 1
    end
    @assert k == nambulances + 1

    @assert sum(length(a) for a in ambulances) == nambulances
    status = fill(:available, nambulances)
    fromtime = zeros(Int, nambulances)
    hospital = zeros(Int, nambulances)

    m = Gurobi.Model(Gurobi.Env(), "redeploy", :minimize)
    Gurobi.setparam!(m, "OutputFlag", 0)
    for a in 1:nambulances, i in 1:nlocations # w variables
        Gurobi.add_bvar!(m, 0.)
    end
    for i in 1:nlocations # eta1
        Gurobi.add_cvar!(m, 0., -Inf, Inf)
    end
    for a in 1:nambulances, i in 1:nlocations # eta2 variables
        Gurobi.add_cvar!(m, lambda, 0., Inf)
    end
    for i in 1:nlocations # eta3 := eta1^2
        Gurobi.add_cvar!(m, 1., 0., Inf)
    end

    # η₁[i] >= available[i] - sum(w[a,i] for a in 1:nambulances)
    #     reformulated to
    # sum(w[a,i] for a in 1:nambulances) + η₁[i] >= available[i]
    for i in 1:nlocations
        Gurobi.add_constr!(m,
            [((1:nambulances)-1)*nlocations+i; nambulances*nlocations + i], # inds
            ones(nambulances + 1), # coeffs
            '>', Float64(available[i]))
    end

    # sum(w[a,i] for i in 1:nlocations) == 1       [a=1:nambulances]
    for a in 1:nambulances
        Gurobi.add_constr!(m,
            collect((a-1)*nlocations + (1:nlocations)), # inds
            ones(nlocations), # coeffs
            '=', 1.)
    end

    # eta2[a,i] >= |w[a,i] - (assignment[a] == i)|   [a=1:nambulances, i=1:nlocations]
    #     reformulated to
    # eta2[a,i] >= w[a,i] - (assignment[a] == i)   i.e. eta2[a,i] - w[a,i] >= -(assignment[a] == i)
    # eta2[a,i] >= - w[a,i] + (assignment[a] == i) i.e. eta2[a,i] + w[a,i] >=  (assignment[a] == i)
    for a in 1:nambulances, i in 1:nlocations
        offset = (a-1)*nlocations + i
        inds = [(nambulances+1)*nlocations + offset, offset]
        Gurobi.add_constr!(m, inds, [1., -1.], '>', - Float64(assignment[a] == i))
        Gurobi.add_constr!(m, inds, [1., 1.], '>', Float64(assignment[a] == i))
    end

    # eta3[i] = eta1[i]^2
    #     reformulated into
    # eta3[i] >= 0
    # eta3[i] >= 1 + 2*eta1[i]
    # eta3[i] >= 4 + 4*eta1[i]
    # eta3[i] >= 9 + 6*eta1[i]
    # ...
    for i in 1:nlocations
        for k in 0.1:0.1:0.9
            Gurobi.add_constr!(m,
                [2*nambulances*nlocations + nlocations + i, nambulances*nlocations + i], # inds
                [1., -2*k], # coeffs
                '>', Float64(k^2))
        end
        for k in 1:3
            Gurobi.add_constr!(m,
                [2*nambulances*nlocations + nlocations + i, nambulances*nlocations + i], # inds
                [1., -Float64(2*k)], # coeffs
                '>', Float64(k^2))
        end
    end

    AssignmentModel(m, lambda, hosp2stn, stn2stn, assignment, ambulances,
                    status, fromtime, hospital, zeros(nambulances*nlocations))
end

function reassign_ambulances!(
        ems,
        problem::DispatchProblem,
        redeploy::AssignmentModel,
        t::Int
    )
    nlocations = length(redeploy.ambulances)
    nambulances = length(redeploy.assignment)
    # DEBUG
    # for a in 1:nambulances
    #     @assert t >= redeploy.fromtime[a]
    # end
    cost(a,i) = if redeploy.status[a] == :available
             0                                 + redeploy.stn2stn[redeploy.assignment[a],i]
        elseif redeploy.status[a] == :responding
            55 - (t - redeploy.fromtime[a])/60 + redeploy.stn2stn[redeploy.assignment[a],i] # + redeploy.hosp2stn[redeploy.hospital[a],i]
        elseif redeploy.status[a] == :atscene
            45 - (t - redeploy.fromtime[a])/60 + redeploy.stn2stn[redeploy.assignment[a],i] # + redeploy.hosp2stn[redeploy.hospital[a],i]
        elseif redeploy.status[a] == :conveying
            30 - (t - redeploy.fromtime[a])/60 + redeploy.stn2stn[redeploy.assignment[a],i] # + redeploy.hosp2stn[redeploy.hospital[a],i]
        elseif redeploy.status[a] == :returning
            15 - (t - redeploy.fromtime[a])/60 + redeploy.stn2stn[redeploy.assignment[a],i]
        elseif redeploy.status[a] == :redeploying
            10 - (t - redeploy.fromtime[a])/60 + redeploy.stn2stn[redeploy.assignment[a],i]
        end
    # (1) Optimize Dynamic Assignment Problem
    # (1i) change RHS to account for wait_queue
    nwait = [length(problem.wait_queue[nbhd]) for nbhd in 1:size(problem.coverage,1)]
    Gurobi.set_dblattrarray!(redeploy.model, "RHS", 1, nlocations, [
        Float64(problem.deployment[i] + sum(
            nwait[nbhd] for nbhd in 1:size(problem.coverage,1)
            if problem.coverage[nbhd,i]
        ))
        for i in 1:nlocations
    ])
    # (1ii) change coeffs to account for traveling time
    let con = Cint[0], ind = Cint[0], val = Float64[0.0]
        for i in 1:nlocations
            con[1] = i
            for a in 1:nambulances
                ind[1] = (a-1)*nlocations + i
                val[1] = 1-max(min(60,cost(a,i)),0)/45
                Gurobi.chg_coeffs!(redeploy.model, con, ind, val)
            end
        end
    end
    Gurobi.update_model!(redeploy.model)
    Gurobi.optimize(redeploy.model)
    # @show Gurobi.get_objval(redeploy.model)

    # (2) Reassign ambulances based on optimal solution
    Gurobi.get_dblattrarray!(redeploy.soln, redeploy.model, "X", 1)
    for a in 1:nambulances
        let stn = redeploy.assignment[a]
            for i in 1:nlocations
                if stn != i && redeploy.soln[(a-1)*nlocations + i] > 0.5
                    # redeploy an existing ambulance
                    push!(problem.redeploy_events, (a,redeploy.assignment[a],i,t))
                    if redeploy.status[a] == :available
                        @assert problem.available[redeploy.assignment[a]] > 0
                        problem.available[stn] -= 1
                        t_end = t + ceil(Int, 0*60*redeploy.stn2stn[stn,i])
                        redeploying_to!(redeploy, a, i, t)
                        enqueue!(ems.eventqueue, (:done, 0, t_end, a), t_end)
                    else
                        println("redeploying amb $a from $(redeploy.assignment[a]) to $i")
                        ambulances = redeploy.ambulances[redeploy.assignment[a]]
                        @assert !in(a, redeploy.ambulances[i])
                        @assert !in(a, ambulances)
                        in(a, ambulances) && deleteat!(ambulances, findfirst(ambulances, a))
                        @assert !in(a, ambulances)
                        redeploy.assignment[a] = i
                    end
                end
            end
        end
    end
end


text:
@doc raw"""
    Sphere

An unconstrained implementation of the Sphere function defined as:

```math
f(\mathbf{x}) = \sum_{i=1}^{d} x_i^2
```

where ``d`` is the dimension of the input vector ``\mathbf{x}``.
"""
struct Sphere <: Unconstrained end

_sphere(x) = sum(x.^2)

@doc raw"""
    Easom

An unconstrained implementation of the 2-dimensional
Easom function defined as:

```math
f(\mathbf{x}) = -\cos{(x_1)} \cos{(x_2)} \exp{[-(x_1 - \pi)^2 - (x_2 - \pi)^2]}
```

where ``x_1`` and ``x_2`` refer to the first and second element of the
input vector ``\mathbf{x}``.
"""
struct Easom <: Unconstrained end

@inline function _easom(x)
    @assert length(x) == 2 "This is a 2D function"

    term_1 = -cos(x[1]) * cos(x[2])
    term_2 = exp(-(x[1] - π)^2 - (x[2] - π)^2)

    return term_1 * term_2
end

@doc raw"""
    Ackley

An unconstrained implementation of the d-dimensional
Ackley function defined as:

```math
f(\mathbf{x}) = -20 e^{ -0.02 \sqrt{\frac{1}{d}\sum_{i=1}^{d}{x_i^2}}} - e^{\frac{1}{d}\sum_{i=1}^{d}{\cos{(2 \pi x_i)}}} + 20 + e
```

where ``d`` is the dimension of the input vector ``\mathbf{x}``.
"""
struct Ackley <: Unconstrained end

@inline function _ackley(x)
    dimension = length(x)
    @assert dimension > 0 "Must have positive dimension"

    term_1 = exp(-0.02 * sqrt(_sphere(x) / dimension))
    term_2 = exp(sum(cos.(2.0 * π * x)) / dimension)

    return -20.0 * term_1 - term_2 + 20.0 + exp(1.0)
end

@doc raw"""
    Rosenbrock

An unconstrained implementation of the d-dimensional
Rosenbrock function defined as:

```math
f(\mathbf{x}) = \sum_{i=1}^{N-1} \left[100(x_{i-1}-x_i^2)^2 +(1-x_i)^2 \right]
```

where ``N`` is the dimension of the input vector ``\mathbf{x}``.
"""
struct Rosenbrock <: Unconstrained end

@inline function _rosenbrock(x)
    dimension = length(x)
    @assert dimension >= 2 "Must be at least 2D"

    total = 0.0
    for i in 1:dimension - 1
        @inbounds total += 100.0 * (x[i + 1] - x[i]^2)^2 + (1.0 - x[i])^2
    end

    return total
end  # function _rosenbrock

@doc raw"""
    Goldstein-Price

An unconstrained implementation of the d-dimensional
Goldstein-Price function defined as:

```math
f(x,y)=[1 + (x + y + 1)^2(19 − 14x+3x^2− 14y + 6xy + 3y^2)] \times \\
[30 + (2x − 3y)^2(18 − 32x + 12x^2 + 4y − 36xy + 27y^2)]
```
"""
struct GoldsteinPrice <: Unconstrained end

@inline function _goldprice(x)
    @assert length(x) == 2 "Exactly 2D"

    term_1 = (x[1] + x[2] + 1.0)^2
    term_1 *= 19.0 - 14.0 * x[1] + 3.0 * x[1]^2 - 14.0 * x[2] + 6.0 * x[1] * x[2] + 3.0 * x[2]^2
    term_1 += 1.0
    term_2 = (2.0 * x[1] - 3.0 * x[2])^2
    term_2 *= 18.0 - 32.0 * x[1] + 12.0 * x[1]^2 + 48.0 * x[2] - 36.0 * x[1] * x[2] + 27.0 * x[2]^2
    term_2 += 30.0

    return term_1 * term_2
end  # function _goldprice

@doc raw"""
    Beale

An unconstrained implementation of the d-dimensional
Beale function defined as:

```math
f(x, y) = (1.5-x+xy)^2+(2.25-x+xy^2)^2+(2.625-x+xy^3)^2
```
"""
struct Beale <: Unconstrained end

@inline function _beale(x)
    @assert length(x) == 2 "Exactly 2D"

    term_1 = (1.5 - x[1] + x[1] * x[2])^2
    term_2 = (2.25 - x[1] + x[1] * x[2]^2)^2
    term_3 = (2.625 - x[1] + x[1] * x[2]^3)^2

    return term_1 + term_2 + term_3
end  # function _beale

@doc raw"""
    Lévy

An unconstrained implementation of the d-dimensional
Lévy function defined as:

```math
f(\mathbf{x}) = \sin^{2}{\pi w_1} + \sum_{i=1}^{d-1} (w_i-1)^2 [1+10\sin^{2}{\pi w_1 + 1}]
+ (w_d-1)^2 [1+\sin^{2}{2\pi w_d}]
```

where

``w_i = 1 + \frac{x_i-1}{4}`` and ``d`` is the dimension of the vector.
"""
struct Levy <: Unconstrained end

@inline function _levy(x)
    ω = @. 1.0 + (x - 1.0) / 4.0

    term_1 = (sin(π * ω[1]))^2

    term_2 = sum(@. (ω - 1.0)^2 * (1.0 + 10.0 * (sin(π * ω + 1.0))^2))
    term_3 = (ω[end] - 1.0)^2 * (1.0 + (sin(2.0 * π * ω[end]))^2)

    return term_1 + term_2 + term_3
end  # function _levy

# Build a dictionary of test functions and their implementations
test_functions = Dict([
    :Sphere => :_sphere,
    :Easom => :_easom,
    :Ackley => :_ackley,
    :Rosenbrock => :_rosenbrock,
    :GoldsteinPrice => :_goldprice,
    :Beale => :_beale,
    :Levy => :_levy
])

for (k, v) in test_functions
    # Create the methods for the given test functions
    @eval $k(x::T) where T = $v(x)
    # Create a special method for the `TestFunction` type
    @eval evaluate(b::$k, x::T) where T = $v(x)
end


text:
#=
main:
- Julia version: 1.6.3
- Author: lqxc
- Date: 2021-10-19
=#

# D means the whole dataset, which should be Vector{Vector{T}}
# ! count(X) = filter(D, X ⊆ _).count()
# support(X => Y) = count(X ∪ Y) / |D|
# confidence(X => Y) = count(X ∪ Y) / count(X)

function rules(
    dicts::Vector{Dict{Set{T},Float64}},
    support::Dict{Set{T}, Float64},
    ς::Float64
   )::Dict{Pair{Set{T},Set{T}}, Tuple{Float64, Float64}} where {T}
    result = Dict{Pair{Set{T}, Set{T}},Tuple{Float64, Float64}}()
    for i ∈ 2:(dicts |> length)
        for x ∈ 1:i-1
            for set ∈ dicts[x] |> keys
                for n_set ∈ dicts[i] |> keys
                    if set ⊆ n_set
                        diff = setdiff(n_set, set)
                        # set ∪ diff = n_set
                        # confidence(set => diff) = count(n_set) / count(set)
                        # (n_set, set => diff, support[n_set] / support[set]) |> println
                        confidence = support[n_set] / support[set]
                        if confidence >= ς
                            result[set => diff] = confidence, confidence / support[diff]
                        end
                    end
                end
            end
        end
    end
    result
end

#=
# scan the data with each item set in item set set,
# and get each support
=#

function scan(
    data::Vector{Vector{T}},
    C::Set{Set{T}},
    ε::Float64,
)::Dict{Set{T},Float64} where {T}
    dict = Dict{Set{T},Float64}()
    data .|> row -> C .|> set -> if set ⊆ row # check each set is subset of each row in data which in C
        v = get!(dict, set, 0) + 1 # Int, yes, set default if is missing, and get value
        dict[set] = v
    end
    len = data |> length # Int, get data set size
    for (key, value) ∈ dict # Vector{T} => Float64, cal each sets support%
        support = value / len
        dict[key] = support
    end
    dict |> keys .|> key -> if dict[key] < ε # remove set which support less than ε
        delete!(dict, key)
    end
    dict
end

function apriori(
    data::Vector{Vector{T}};
    ε::Float64 = 0.5,
    ς::Float64 = ε,
)::Tuple{
    Vector{Dict{Set{T}, Float64}},
    Dict{Pair{Set{T},Set{T}}, Tuple{Float64, Float64}}
} where {T}
    n = 2
    C1::Set{Set{T}} = (∪(data...) .|> v -> Set([v])) |> Set
    support = Dict{Set{T}, Float64}()
    dict = scan(data, C1, ε)
    merge!(+, support, dict)
    L::Vector{Set{Set{T}}} = [[(dict |> keys)...] |> Set]
    D::Vector{Dict{Set{T},Float64}} = [dict]
    while L[n-1] |> length != 0
        C = Set{Set{T}}()
        # len = L[n-1] |> length
        for p ∈ L[n-1], q ∈ L[n-1]
            if p != q
                push!(C, p ∪ q)
            end
        end
        dict = scan(data, C, ε)
        merge!(+, support, dict)
        push!(D, dict)
        push!(L, [keys(dict)...] |> Set)
        n += 1
    end
    D, rules(D, support, ς)
end

# data = readlines("files/groceries.csv") .|> s -> split(s, ",")

# apriori(data; ε = 0.03, ς = 0.2) .|>
# [
#    dicts -> dicts .|> dict -> dict |> keys .|> key -> println(key, " -> ", dict[key]),
#     dict -> dict |> keys .|> key -> println(key, " -> ", dict[key])
# ]
apriori([
    [:I1, :I2, :I5],
    [:I2, :I4],
    [:I2, :I3],
    [:I1, :I2, :I4],
    [:I1, :I3],
    [:I2, :I3],
    [:I1, :I3],
    [:I1, :I2, :I3, :I5],
    [:I1, :I2, :I3]
]; ε = 0.22,ς = 0.7) .|>
[
    dict -> dict .|> dict -> dict |> keys .|> key -> println(key, " -> ", dict[key]),
    dict -> dict |> keys .|> key -> println(key, " -> ", dict[key])
]

max_stars_repo_name: lovebaihezi/lab

max_issues_repo_name: lovebaihezi/lab

max_forks_repo_name: lovebaihezi/lab


text:
compiler_env, = ARGS

if VERSION < v"0.7-"
    error("Unsupported Julia version $VERSION")
end

using Pkg

function cat_build_log(pkg)
    modpath = Base.locate_package(pkg)
    if modpath !== nothing
        logfile = joinpath(dirname(modpath), "..", "deps", "build.log")
        if isfile(logfile)
            print(stderr, read(logfile, String))
            return
        end
    end
    @error "build.log for $pkg not found"
end

Pkg.activate(compiler_env)
@info "Installing PackageCompiler..."

Pkg.add([
    PackageSpec(
        name = "PackageCompiler",
        uuid = "9b87118b-4619-50d2-8e1e-99f35a4d4d9d",
        version = "1",
    )
])
cat_build_log(Base.PkgId(
    Base.UUID("9b87118b-4619-50d2-8e1e-99f35a4d4d9d"),
    "PackageCompiler"))

@info "Loading PackageCompiler..."
using PackageCompiler

@info "PackageCompiler is successfully installed at $compiler_env"

max_stars_repo_name: grahamgill/pyjulia

max_issues_repo_name: grahamgill/pyjulia

max_forks_repo_name: grahamgill/pyjulia