text:
import Base: convert

export MPIScanner, MPIScannerGeneral, scannerBoreSize, scannerFacility,
       scannerManufacturer, scannerName, scannerTopology, scannerGradient,
       name, configDir, generalParams, getDevice, getDevices, getSequenceList,
       asyncMeasurement, SequenceMeasState, asyncProducer, prepareAsyncMeasurement,
       getProtocolList, setProtocol, getTransferFunctionList

"""
    $(SIGNATURES)

Recursively find all concrete types of the given type.
"""
function deepsubtypes(type::DataType)
  subtypes_ = subtypes(type)
  allSubtypes = subtypes_
  for subtype in subtypes_
    subsubtypes_ = deepsubtypes(subtype)
    allSubtypes = vcat(allSubtypes, subsubtypes_)
  end
  return allSubtypes
end

"""
    $(SIGNATURES)

Retrieve the concrete type of a given supertype corresponding to a given string.
"""
function getConcreteType(supertype_::DataType, type::String)
  knownTypes = deepsubtypes(supertype_)
  foundImplementation = nothing
  for Implementation in knownTypes
    if string(Implementation) == type
      foundImplementation = Implementation
    end
  end

  if !isnothing(foundImplementation)
    return foundImplementation
  else
    error("The type implied by the string `$type` could not be retrieved since its device struct was not found.")
  end
end

"""
    $(SIGNATURES)

Initiate devices from the given configuration dictionary.

The device types are referenced by strings matching their device struct name.
All device structs are supplied with the device ID and the corresponding
device configuration struct.
"""
function initiateDevices(devicesParams::Dict{String, Any})
  devices = Dict{String, Device}()

  # Get implementations for all devices in the specified order
  for deviceID in devicesParams["initializationOrder"]
    if haskey(devicesParams, deviceID)
      params = devicesParams[deviceID]
      deviceType = pop!(params, "deviceType")

      dependencies_ = Dict{String, Union{Device, Missing}}()
      if haskey(params, "dependencies")
        deviceDepencencies = pop!(params, "dependencies")
        for dependencyID in deviceDepencencies
          dependencies_[dependencyID] = missing
        end
      end

      DeviceImpl = getConcreteType(Device, deviceType)
      DeviceParamsImpl = getConcreteType(DeviceParams, deviceType*"Params") # Assumes the naming convention of ending with [...]Params!

      paramsInst = DeviceParamsImpl(params)
      devices[deviceID] = DeviceImpl(deviceID=deviceID, params=paramsInst, dependencies=dependencies_) # All other fields must have default values!
    else
      throw(ScannerConfigurationError("The device ID `$deviceID` was not found in the configuration. Please check your configuration."))
    end
  end

  # Set dependencies for all devices
  for device in values(devices)
    for dependencyID in keys(dependencies(device))
      device.dependencies[dependencyID] = devices[dependencyID]
    end

    if !checkDependencies(device)
      throw(ScannerConfigurationError("Unspecified dependency error in device with "
                                     *"ID `$(deviceID(device))`. The device depends "
                                     *"on the following device IDs: $(keys(device.dependencies))"))
    end
  end

  # Initiate all devices in the specified order
  for deviceID in devicesParams["initializationOrder"]
    init(devices[deviceID])
    if !isOptional(devices[deviceID]) && !isPresent(devices[deviceID])
      @error "The device with ID `$deviceID` should be present but isn't."
    end
  end

  return devices
end

"""
    $(SIGNATURES)

General description of the scanner.

Note: The fields correspond to the root section of an MDF file.
"""
Base.@kwdef struct MPIScannerGeneral
  "Bore size of the scanner."
  boreSize::typeof(1u"mm")
  "Facility where the scanner is located."
  facility::String
  "Manufacturer of the scanner."
  manufacturer::String
  "Name of the scanner"
  name::String
  "Topology of the scanner, e.g. FFL or FFP."
  topology::String
  "Gradient of the scanners selection field."
  gradient::typeof(1u"T/m")
  "Path of the dataset store."
  datasetStore::String
  "Default sequence of the scanner."
  defaultSequence::String = ""
  "Default protocol of the scanner."
  defaultProtocol::String = ""
  "Location of the scanner's transfer function."
  transferFunction::String = ""
  "Thread ID of the producer thread."
  producerThreadID::Int32 = 2
  "Thread ID of the consumer thread."
  consumerThreadID::Int32 = 3
  "Thread ID of the producer thread."
  protocolThreadID::Int32 = 4
end

abstract type AsyncBuffer end
abstract type AsyncMeasTyp end
struct FrameAveragedAsyncMeas <: AsyncMeasTyp end
struct RegularAsyncMeas <: AsyncMeasTyp end
# TODO Update
asyncMeasType(sequence::Sequence) = acqNumFrameAverages(sequence) > 1 ? FrameAveragedAsyncMeas() : RegularAsyncMeas()

mutable struct FrameAverageBuffer
  buffer::Array{Float32, 4}
  setIndex::Int
end
FrameAverageBuffer(samples, channels, periods, avgFrames) = FrameAverageBuffer(zeros(Float32, samples, channels, periods, avgFrames), 1)

mutable struct SequenceMeasState
  numFrames::Int
  nextFrame::Int
  channel::Union{Channel, Nothing}
  producer::Union{Task,Nothing}
  consumer::Union{Task, Nothing}
  asyncBuffer::AsyncBuffer
  buffer::Array{Float32,4}
  avgBuffer::Union{FrameAverageBuffer, Nothing}
  #temperatures::Matrix{Float64} temps are not implemented atm
  type::AsyncMeasTyp
end

"""
    $(SIGNATURES)

Basic description of a scanner.
"""
mutable struct MPIScanner
  "Name of the scanner"
  name::String
  "Path to the used configuration directory."
  configDir::String
  "General parameters of the scanner like its bore size or gradient."
  generalParams::MPIScannerGeneral
  "Device instances instantiated by the scanner from its configuration."
  devices::Dict{AbstractString, Device}
  seqMeasState::Union{SequenceMeasState, Nothing}

  """
    $(SIGNATURES)

  Initialize a scanner by its name.
  """
  function MPIScanner(name::AbstractString)
    # Search for scanner configurations of the given name in all known configuration directories
    # If you want to add a configuration directory, please use addConfigurationPath(path::String)
    filename = nothing
    configDir = nothing
    for path in scannerConfigurationPath
      configDir = joinpath(path, name)
      if isdir(configDir)
        filename = joinpath(configDir, "Scanner.toml")
        break
      end
    end

    if isnothing(filename)
      throw(ScannerConfigurationError("Could not find a valid configuration for scanner with name `$name`. Search path contains the following directories: $scannerConfigurationPath."))
    end

    @debug "Instantiating scanner `$name` from configuration file at `$filename`."

    params = TOML.parsefile(filename)
    generalParams = params_from_dict(MPIScannerGeneral, params["General"])
    @assert generalParams.name == name "The folder name and the scanner name in the configuration do not match."
    devices = initiateDevices(params["Devices"])

    scanner = new(name, configDir, generalParams, devices, nothing)

    return scanner
  end
end

"""
    $(SIGNATURES)

Close the devices when closing the scanner.
"""
function Base.close(scanner::MPIScanner)
  for device in getDevices(scanner, Device)
    close(device)
  end
end

"Name of the scanner"
name(scanner::MPIScanner) = scanner.name

"Path to the used configuration directory."
configDir(scanner::MPIScanner) = scanner.configDir

"General parameters of the scanner like its bore size or gradient."
generalParams(scanner::MPIScanner) = scanner.generalParams

"Location of the scanner's transfer function."
transferFunction(scanner::MPIScanner) = scanner.generalParams.transferFunction

"Check, whether the scanner has a transfer function defined."
hasTransferFunction(scanner::MPIScanner) = transferFunction(scanner) != ""

"""
    $(SIGNATURES)

Retrieve a device by its `deviceID`.
"""
getDevice(scanner::MPIScanner, deviceID::String) = scanner.devices[deviceID]

"""
    $(SIGNATURES)

Retrieve all devices of a specific `deviceType`. Returns an empty vector if none are found
"""
function getDevices(scanner::MPIScanner, deviceType::Type{<:Device})
  matchingDevices = Vector{Device}()
  for (deviceID, device) in scanner.devices
    if typeof(device) <: deviceType
      push!(matchingDevices, device)
    end
  end
  return matchingDevices
end
function getDevices(scanner::MPIScanner, deviceType::String)
  knownDeviceTypes = deepsubtypes(Device)
  deviceTypeSearched = knownDeviceTypes[findall(type->string(type)==deviceType, knownDeviceTypes)][1]
  return getDevices(scanner, deviceTypeSearched)
end

"""
$(SIGNATURES)

Retrieve all devices of a specific `deviceType` if it can be unambiguously retrieved. Returns nothing if no such device can be found and throws an error if multiple devices fit the type.
"""
function getDevice(scanner::MPIScanner, deviceType::Type{<:Device})
  devices = getDevices(scanner, deviceType)
  if length(devices) > 1
    error("The scanner has more than one $(string(deviceType)) device. Therefore, a single $(string(deviceType)) cannot be retrieved unambiguously.")
  elseif length(devices) == 0
    return nothing
  else
    return devices[1]
  end
end

"Bore size of the scanner."
scannerBoreSize(scanner::MPIScanner) = scanner.generalParams.boreSize

"Facility where the scanner is located."
scannerFacility(scanner::MPIScanner) = scanner.generalParams.facility

"Manufacturer of the scanner."
scannerManufacturer(scanner::MPIScanner) = scanner.generalParams.manufacturer

"Name of the scanner"
scannerName(scanner::MPIScanner) = scanner.generalParams.name

"Topology of the scanner, e.g. FFL or FFP."
scannerTopology(scanner::MPIScanner) = scanner.generalParams.topology

"Gradient of the scanners selection field."
scannerGradient(scanner::MPIScanner) = scanner.generalParams.gradient

"""
    $(SIGNATURES)

Retrieve a list of all sequences available for the scanner.
"""
function getSequenceList(scanner::MPIScanner)
  path = joinpath(configDir(scanner), "Sequences")
  if isdir(path)
    return String[ splitext(seq)[1] for seq in filter(a->contains(a,".toml"),readdir(path))]
  else
    return String[]
  end
end

"""
    $(SIGNATURES)

Constructor for a sequence of `name` from `configDir`.
"""
function MPIFiles.Sequence(configdir::AbstractString, name::AbstractString)
  path = joinpath(configdir, "Sequences", name*".toml")
  if !isfile(path)
    error("Sequence $(path) not available!")
  end
  return sequenceFromTOML(path)
end

"""
    $(SIGNATURES)

Constructor for a sequence of `name` from the configuration directory specified for the scanner.
"""
MPIFiles.Sequence(scanner::MPIScanner, name::AbstractString) = Sequence(configDir(scanner), name)

function getTransferFunctionList(scanner::MPIScanner)
  path = joinpath(configDir(scanner), "TransferFunctions")
  if isdir(path)
    return String[ splitext(seq)[1] for seq in filter(a->contains(a,".h5"),readdir(path))]
  else
    return String[]
  end
end

function MPIFiles.TransferFunction(configdir::AbstractString, name::AbstractString)
  path = joinpath(configdir, "TransferFunctions", name*".h5")
  if !isfile(path)
    error("TransferFunction $(path) not available!")
  end
  return TransferFunction(path)
end

MPIFiles.TransferFunction(scanner::MPIScanner) = TransferFunction(configDir(scanner),transferFunction(scanner))

#### Scanner Measurement Functions ####
####  Async version  ####
SequenceMeasState() = SequenceMeasState(0, 1, nothing, nothing, nothing, DummyAsyncBuffer(nothing), zeros(Float64,0,0,0,0), nothing, RegularAsyncMeas())

function asyncMeasurement(scanner::MPIScanner, sequence::Sequence)
  prepareAsyncMeasurement(scanner, sequence)
  scanner.seqMeasState.producer = @tspawnat scanner.generalParams.producerThreadID asyncProducer(scanner.seqMeasState.channel, scanner, sequence)
  bind(scanner.seqMeasState.channel, scanner.seqMeasState.producer)
  scanner.seqMeasState.consumer = @tspawnat scanner.generalParams.consumerThreadID asyncConsumer(scanner.seqMeasState.channel, scanner)
  return scanner.seqMeasState
end

function prepareAsyncMeasurement(scanner::MPIScanner, sequence::Sequence)
  daq = getDAQ(scanner)
  numFrames = acqNumFrames(sequence)
  rxNumSamplingPoints = rxNumSamplesPerPeriod(sequence)
  numPeriods = acqNumPeriodsPerFrame(sequence)
  frameAverage = acqNumFrameAverages(sequence)
  setup(daq, sequence)

  # Prepare buffering structures
  @info "Allocating buffer for $numFrames frames"
  # TODO implement properly with only RxMeasurementChannels
  buffer = zeros(Float32,rxNumSamplingPoints, length(rxChannels(sequence)),numPeriods,numFrames)
  #buffer = zeros(Float32,rxNumSamplingPoints,numRxChannelsMeasurement(daq),numPeriods,numFrames)
  avgBuffer = nothing
  if frameAverage > 1
    avgBuffer = FrameAverageBuffer(zeros(Float32, frameAverageBufferSize(daq, frameAverage)), 1)
  end
  channel = Channel{channelType(daq)}(32)

  # Prepare measState
  measState = SequenceMeasState(numFrames, 1, nothing, nothing, nothing, AsyncBuffer(daq), buffer, avgBuffer, asyncMeasType(sequence))
  measState.channel = channel

  scanner.seqMeasState = measState
end

function asyncProducer(channel::Channel, scanner::MPIScanner, sequence::Sequence; prepTx = true)
  su = getSurveillanceUnit(scanner) # Maybe also support multiple SU units?
  if !isnothing(su)
    enableACPower(su)
    # TODO Send expected enable time to SU
  end
  robots = getRobots(scanner)
  for robot in robots
    disable(robot)
  end

  try
    daq = getDAQ(scanner)
    asyncProducer(channel, daq, sequence, prepTx = prepTx)
  finally
    if !isnothing(su)
      disableACPower(su)
    end
    for robot in robots
      enable(robot)
    end
  end
end

# Default Consumer
function asyncConsumer(channel::Channel, scanner::MPIScanner)
  daq = getDAQ(scanner)
  measState = scanner.seqMeasState

  @info "Consumer start"
  while isopen(channel) || isready(channel)
    while isready(channel)
      chunk = take!(channel)
      updateAsyncBuffer!(measState.asyncBuffer, chunk)
      updateFrameBuffer!(measState, daq)
    end
    sleep(0.001)
  end
  @info "Consumer end"

  # TODO calibTemperatures is not filled in asyncVersion yet, would need own innerAsyncConsumer
  #if length(measState.temperatures) > 0
  #  params["calibTemperatures"] = measState.temperatures
  #end
end

#### Protocol ####
function getProtocolList(scanner::MPIScanner)
  path = joinpath(configDir(scanner), "Protocols/")
  if isdir(path)
    return String[ splitext(proto)[1] for proto in filter(a->contains(a,".toml"),readdir(path))]
  else
    return String[]
  end
end

function execute(scanner::MPIScanner, protocol::Protocol)
  return execute(protocol, scanner.generalParams.protocolThreadID)
end


text:
module TestSetters

include("preamble.jl")

@testset "nullsetter" begin
    @test set(1, nullsetter, 2) === 1
    @test string(nullsetter) === "nullsetter"
end

@testset "ToField" begin
    setter = (@lens _.x) ∘ ToField(@lens _.a)
    @test set((x = 1, y = 2), setter, (a = 10, b = 20)) === (x = 10, y = 2)
end

end  # module

ext: jl


text:
module bfga

    import Base.isless

    include("Types.jl")
    using .Types

    #using GeneticAlgorithms
    #include("GeneticAlgorithms.jl")
    #using .GeneticAlgorithms

    include("Bf.jl")
    import .BfInterpreter

    mutable struct bfMonster <: Entity
        dna::Array
        fitness
        bonus :: Float64
        age :: Int
        m_length :: Int
        program :: String

        # = new(create_entity(1, nbGenes).dna, nothing, 1, nbGenes)
        bfMonster(dna :: Array) = new(dna, 0, 0,  1, length(dna), join(genesToBfInstructions(dna), ""))
        bfMonster(dna :: Array, age :: Int) = new(dna, 0, 0, age, length(dna), join(genesToBfInstructions(dna), ""))
        #bfMonster(age :: Int) = create_entity(age,nbGenes)
    end

    #nbGenes = 150

    function isless(lhs::bfMonster, rhs::bfMonster)
        (lhs.fitness, rhs.bonus) < (rhs.fitness , lhs.bonus) # abs(lhs.fitness) > abs(rhs.fitness)
    end

    function create_entity(num, nbgenes)
        babybf = [rand() for i in 1:nbgenes]
        bfMonster(babybf, num)
    end

    #Constants
global instructions = ['<', '>', '+', '-', '.', ',', '[', ']'] #, '$', '!', '*']

    global instructions2 = ['<', '>', '+', '-', '.', ',', '[', ']','$', '!', '*' , '#' , '@', '/' , '%']

    function floatToBfInstruction( number ::Float64 )
        @assert (number >= 0.0 && number <= 1.0)
        memoire = number*length(instructions)
        ent = trunc(Int, memoire)
        if ent == memoire
            return instructions[ent]
        else
            return instructions[ent+1]
        end
    end

    function entityToBfInstructions!( ent )
        ent.program = join(genesToBfInstructions(ent.dna), "")
    end

    function genesToBfInstructions( genes :: Array{Float64,1} )
        instructions = []
        for num in genes
            push!(instructions,floatToBfInstruction(num))
        end
        instructions
    end

    function getBfCode(ent)
        ent.program
    end

    #=function group_entities!(pop)
        println("BEST: ", pop[1].fitness)
        @show pop[1].dna
        println("SECOND: ", pop[2].fitness)

        if pop[1].fitness == 0
            return
        end

        liste = []
        # simple naive groupings that pair the best entitiy with every other
        for i in 1:length(pop)
            #produce([1, i])
            prems = 1
            if rand() > 0.3
                prems = rand(2:10)
                #=if prems < 0
                    prems *= -1
                end
                prems = (prems % length(pop)) +1 =#
            end
            push!(liste , [prems,i])
        end
        liste
    end=#

    function getInstructionsSet()
        mem = BfInterpreter.getInstructionsDict()
        @show mem
        mem
    end

    function crossover(group)
        parent1 = group[1]
        parent2 = group[2]
        n1 = length(parent1.dna)
        n2 = length(parent2.dna)
        pos = trunc(Int, rand() * n1)
        minN = min(n1,n2)

        child1 = create_entity(parent1.age, minN)
        child2 = create_entity(parent1.age, minN)
        child1.age = parent1.age + 1
        child2.age = parent1.age + 1

        for i in 1:minN
            if i < pos
                child1.dna[i] = parent2.dna[i]
                child2.dna[i] = parent1.dna[i]
            else
                child2.dna[i] = parent2.dna[i]
                child1.dna[i] = parent1.dna[i]
            end
        end
        mutate(child1)
        mutate(child2)
        #clearCode!(child1)
        #clearCode!(child2)
        child1.m_length = minN #length(child1.dna)
        child2.m_length = minN #length(child2.dna)
        child1, child2
    end

    function mutate(ent)
        m_length = length(ent.dna)
        # Go through each bit.
        for pos in 1:m_length

            # Should this bit mutate?
            if (rand() < 0.06)
                # Select a mutation type.
                r = rand()
                if (r <= 0.25)
                    # Insertion mutation.
                    # Get shift index.
                    mutationIndex = pos;

                    # Make a copy of the current bit before we mutate it.
                    shiftBit = ent.dna[mutationIndex];

                    # Set random bit at mutation index.
                    ent.dna[mutationIndex] = rand();

                    # Bump bits up or down by 1.
                    up = rand() >= 0.5;
                    if up
                        # Bump bits up by 1.
                        for i in ((mutationIndex + 1 ):m_length)
                            nextShiftBit = ent.dna[i];

                            ent.dna[i] = shiftBit;

                            shiftBit = nextShiftBit;
                        end
                    else
                        # Bump bits down by 1.
                        for i in mutationIndex:-1:1
                            nextShiftBit = ent.dna[i];

                            ent.dna[i] = shiftBit;

                            shiftBit = nextShiftBit;
                        end
                    end
                elseif (r <= 0.5)
                    # Deletion mutation.
                    # Get deletion index.
                    mutationIndex = pos;

                    # Bump bits up or down by 1.
                    up = rand() >= 0.5;
                    if (up)
                        # Bump bits up by 1.
                        for i in mutationIndex:-1:2
                            ent.dna[i] = ent.dna[i - 1];
                        end

                        # Add a new mutation bit at front of genome to replace the deleted one.
                        ent.dna[1] = rand();
                    else
                        # Bump bits down by 1.
                        for i in mutationIndex:(m_length - 1)
                            ent.dna[i] = ent.dna[i + 1];
                        end

                        # Add a new mutation bit at end of genome to replace the deleted one.
                        ent.dna[m_length] = rand();
                    end
                elseif (r <= 0.75)
                    # Shift/rotation mutation.
                    # Bump bits up or down by 1.
                    up = rand() >= 0.5;
                    if (up)
                        # Bump bits up by 1. 1, 2, 3 => 3, 1, 2
                        shiftBit = ent.dna[1];

                        for i in 1:m_length
                            if (i > 1)
                                # Make a copy of the current bit.
                                temp = ent.dna[i];

                                # Set the current bit to the previous one.
                                ent.dna[i] = shiftBit;

                                # Select the next bit to be copied.
                                shiftBit = temp;
                            else
                                # Wrap last bit to front.
                                ent.dna[i] = ent.dna[m_length];
                            end
                        end
                    else
                        # Bump bits down by 1. 1, 2, 3 => 2, 3, 1
                        shiftBit = ent.dna[m_length ];

                        for i in m_length:-1:1
                            if (i < m_length)
                                # Make a copy of the current bit.
                                temp = ent.dna[i];

                                # Set the current bit to the previous one.
                                ent.dna[i] = shiftBit;

                                # Select the next bit to be copied.
                                shiftBit = temp;
                            else
                                # Wrap first bit to end.
                                ent.dna[i] = ent.dna[1];
                            end
                        end
                    end
                else
                    # Replacement mutation.
                    # Mutate bits.
                    mutation = rand();
                    ent.dna[pos] = mutation;
                end
            end
        end
    end

    function expand(ent , size :: Int)
        #println("EXPAND CAAALLLED")
        originalSize = length(ent.dna);
        difference = size - originalSize;

        # Resize the genome array.
        newGenes = []

        if (difference > 0)
            if (rand() < 0.5)
                # Extend at front.
                newGenes = [ [rand() for i in 1:difference]; deepcopy(ent.dna) ]
            else
                # Extend at back.
                newGenes = [ deepcopy(ent.dna); [rand() for i in 1:difference] ]
            end
            ent.dna = newGenes;
        else
            ent.dna = ent.dna[1:size]
            #Array.Resize(ref ent.dna, size);
        end
        ent.m_length = size;
    end

    function clearCode!(ent)
        n = ent.m_length
        bfCodeArray = [] #Array{Float64,1}
        #age = ent.age

        #=bfcode = join((genesToBfInstructions(bfgenes)), "")
        replace(bfcode, r"[^\+\-\<\>\.\,\[\]]" => s"") =#
        i= 1
        while i < n
            if floatToBfInstruction(ent.dna[i]) == '[' && floatToBfInstruction(ent.dna[i+1]) == ']'
                i += 1
            else
                append!(bfCodeArray, [ent.dna[i] ])
            end
            i+=1
        end
        # has to be changed, ugly convertion
        ent.dna = [ i::Float64 for i in bfCodeArray] :: Array{Float64,1}
        ent.m_length = length(bfCodeArray)
    end

end


text:
using GeoMakie
using CairoMakie
using GeoMakie.GeoInterface
using GeoJSON
using Test
using Downloads
using Proj4

@testset "GeoMakie" begin
    # select a coordinate projection, using a string that PROJ accepts
    # see e.g. https://proj.org/operations/projections/index.html
    source = "+proj=longlat +datum=WGS84"
    dest = "+proj=natearth2"
    trans = Proj4.Transformation(source, dest, always_xy=true)
    ptrans = Makie.PointTrans{2}(trans)

    fig = Figure()
    ax = Axis(fig[1,1], aspect = DataAspect())

    # all input data coordinates are projected using this function
    ax.scene.transformation.transform_func[] = ptrans

    # draw projected grid lines and set limits accordingly
    lats = -90:10.0:90
    lons = -180:10.0:180
    lons = collect(lons)
    lons[end] = prevfloat(lons[end])  # avoid PROJ wrapping 180 to -180
    field = [exp(cosd(l)) + 3(y/90) for l in lons, y in lats]
    points = map(CartesianIndices(size(field))) do xy
        x, y = Tuple(xy)
        Point2f0(lons[x], lats[y])
    end
    limits = FRect2D(Makie.apply_transform(ptrans, points))
    limits!(ax, limits)
    wireframe!(ax, lons, lats, field, color=(:gray, 0.2), transparency=true)

    # add black polygons for land area
    url = "https://raw.githubusercontent.com/nvkelso/natural-earth-vector/master/geojson/"
    land = Downloads.download(url * "ne_110m_land.geojson", IOBuffer())
    land_geo = GeoJSON.read(seekstart(land))
    poly!(ax, land_geo, color="black")

    # add grey dots for populated places
    pop = Downloads.download(url * "ne_10m_populated_places_simple.geojson", IOBuffer())
    pop_geo = GeoJSON.read(seekstart(pop))
    scatter!(ax, GeoMakie.geo2basic(pop_geo), color="lightgrey", markersize=1.2)

    hidedecorations!(ax)
    hidespines!(ax)
    display(fig)
    @test true # :)
end


text:
function round(x::FloatingPoint)
    isfinite(x) || return x+x
    t = trunc(x)
    if x > 0.0
        x-t >= half(x) ? t+one(t) : t
    else
        x-t <= -half(x) ? t-one(t) : t
    end
end

function floor(x::FloatingPoint)
    t = trunc(x)
    x < zero(x) && x!=t ? t-one(x) : t
end

function ceil(x::FloatingPoint)
    t = trunc(x)
    x > zero(x) && x!=t ? t+one(x) : t
end


text:
#EasyPlotQwt Display functionnality
#-------------------------------------------------------------------------------

#==Constants
===============================================================================#

const MIME2QWTFMT_MAP = Dict{String,String}(
	"image/png" => "png",
	"image/bmp" => "bmp",
	"image/tiff" => "tif",
	"image/jpeg" => "jpeg",
	"application/pdf" => "pdf",
)

#==Types
===============================================================================#
#TODO: support guimode=false

mutable struct PlotDisplay <: EasyPlot.EasyPlotDisplay #Don't export.  Qualify with Module
	guimode::Bool
	args::Tuple
	kwargs::Base.Iterators.Pairs
	PlotDisplay(args...; guimode=true, kwargs...) = new(guimode, args, kwargs)
end

#TODO: Generate from MIME2QWTFMT_MAP?
const SupportedMIME = Union{
	MIME"image/png",
	MIME"image/bmp",
	MIME"image/tiff",
	MIME"image/jpeg",
	MIME"application/pdf",
}

#==Helper functions
===============================================================================#
mimestr(::MIME{T}) where T = string(T)

#==Top-level rendering functions
===============================================================================#
#Do not overwrite Base.display... would circumvent display system.
function EasyPlot._display(fig::Figure)
	fig[:show]()

#=IMPORTANT:
	#PyCall.pygui_start(:qt_pyqt4) seems to handle QT events in background thread.
	#For blocking show(), the following can be used instead:
	app = GUICore.qapplication()
	app[:exec_]() #Process app events (modal)
=#
	nothing
end

function EasyPlot.render(d::PlotDisplay, eplot::EasyPlot.Plot)
	fig = Figure()
	return render(fig, eplot)
end

#Support _write function natively (use high-quality, by default):
function EasyPlot._write(filepath::String, mime::SupportedMIME, fig::Figure; draft::Bool=false)
	format = MIME2QWTFMT_MAP[mimestr(mime)]
	_save(fig, filepath, format, draft=draft)
end

#_save(fig::Figure, path::String, format::String, draft=false)
function Base.show(io::IO, mime::SupportedMIME, fig::Figure)
	tmpfile = "$(tempname())_export"
	EasyPlot._write(tmpfile, mime, fig, draft=true)
	data = readall(tmpfile)
	write(io, data)
	rm(tmpfile)
end

Base.showable(mime::MIME, eplot::EasyPlot.Plot, d::PlotDisplay) =
	hasmethod(show, (IO, typeof(mime), Figure))

#Last line


text:
################################################################################
#                           SUPPORT ITERATION METHODS
################################################################################
## Form a placeholder parameter reference given the object index
# IndependentParameterIndex
function _temp_parameter_ref(
    model::InfiniteOpt.InfiniteModel,
    idx::InfiniteOpt.IndependentParameterIndex
    )::InfiniteOpt.IndependentParameterRef
    return InfiniteOpt.IndependentParameterRef(model, idx)
end

# DependentParametersIndex
function _temp_parameter_ref(
    model::InfiniteOpt.InfiniteModel,
    idx::InfiniteOpt.DependentParametersIndex
    )::InfiniteOpt.DependentParameterRef
    idx = InfiniteOpt.DependentParameterIndex(idx, 1)
    return InfiniteOpt.DependentParameterRef(model, idx)
end

# Return the collected supports of an infinite parameter
function _collected_supports(
    pref::Union{InfiniteOpt.IndependentParameterRef, InfiniteOpt.DependentParameterRef}
    )::Vector
    supp_dict = InfiniteOpt._parameter_supports(pref)
    supp_list = collect(keys(supp_dict))
    # append a placeholder NaN support at the end to be used for efficient combinatorics
    return push!(supp_list, map(i -> NaN, first(supp_list)))
end

# Return the collected support labels of an infinite parameter
function _collected_support_labels(
    pref::Union{InfiniteOpt.IndependentParameterRef, InfiniteOpt.DependentParameterRef},
    supports::Vector
    )::Vector{Set{DataType}}
    supp_dict = InfiniteOpt._parameter_supports(pref)
    default = Set{DataType}()
    return map(k -> get(supp_dict, k, default), supports)
end

"""
    set_parameter_supports(trans_model::JuMP.Model,
                           inf_model::InfiniteOpt.InfiniteModel)::Nothing

Collect the infinite parameter supports stored in their respective dictionaries
form `inf_model` and process them into a tuple of vectors where each vector
contains the collected supports of a particular infinite parameter. These support
collections are ordered in accordance with the definition order of the
parameters (i.e., their object numbers). A support collection assocciated with
an independent will be a `Vector{Float64}` and a support collection associated
with a group of dependent parameters will be a `Vector{Vector{Float64}}`. Note
that each collection vector will include an extra final placeholder element
comprised of `NaN`s for convenience in generating support indices via
[`support_index_iterator`](@ref). This also gathers the associated support labels.

Before this is all done, `InfiniteOpt.add_generative_supports` is invoked as needed.
"""
function set_parameter_supports(
    trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    # gather the basic information
    param_indices = InfiniteOpt._param_object_indices(inf_model)
    prefs = map(idx -> _temp_parameter_ref(inf_model, idx), param_indices)
    data = transcription_data(trans_model)
    # check and add supports to prefs as needed
    for pref in prefs
        InfiniteOpt.add_generative_supports(pref)
        if InfiniteOpt.has_internal_supports(pref)
            data.has_internal_supports = true
        end
    end
    # build and add the support/label tuples
    supps = Tuple(_collected_supports(pref) for pref in prefs)
    labels = Tuple(_collected_support_labels(pref, supps[i])
                   for (i, pref) in enumerate(prefs))
    data.supports = supps
    data.support_labels = labels
    return
end

################################################################################
#                        VARIABLE INITIALIZATION METHODS
################################################################################
"""
    transcribe_finite_variables!(trans_model::JuMP.Model,
                               inf_model::InfiniteOpt.InfiniteModel)::Nothing

Create a transcription variable (i.e., a JuMP variable) for each `FiniteVariable`
stored in `inf_model` and add it to `trans_model`. The variable mapping is
also stored in `TranscriptionData.finvar_mappings` which enables
[`transcription_variable`](@ref) and [`lookup_by_support`](@ref).
"""
function transcribe_finite_variables!(
    trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (idx, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.FiniteVariable)
        hvref = InfiniteOpt._make_variable_ref(inf_model, idx)
        vref = JuMP.add_variable(trans_model,
                                 JuMP.ScalarVariable(object.variable.info),
                                 object.name)
        transcription_data(trans_model).finvar_mappings[hvref] = vref
    end
    return
end

# Return a proper scalar variable info object given on with a start value function
function _format_infinite_info(
    var::InfiniteOpt.InfiniteVariable,
    support::Vector{Float64}
    )::JuMP.VariableInfo
    # generate the start value
    info = var.info
    if var.is_vector_start
        start = info.start(support)
    else
        start = info.start(Tuple(support, var.parameter_refs)...)
    end
    # make the info and return
    return JuMP.VariableInfo(info.has_lb, info.lower_bound, info.has_ub,
                             info.upper_bound, info.has_fix, info.fixed_value,
                             info.has_start, start, info.binary, info.integer)
end

"""
    transcribe_infinite_variables!(trans_model::JuMP.Model,
                                   inf_model::InfiniteOpt.InfiniteModel)::Nothing

Create transcription variables (i.e., JuMP variables) for each `InfiniteVariable`
stored in `inf_model` and add them to `trans_model`. The variable mappings are
also stored in `TranscriptionData.infvar_mappings` in accordance with
`TranscriptionData.infvar_lookup` which enable [`transcription_variable`](@ref)
and [`lookup_by_support`](@ref). Note that the supports will not be generated
until `InfiniteOpt.variable_supports` is invoked via `InfiniteOpt.supports`.
Note that `TranscriptionData.infvar_support_labels` is also populated.
"""
function transcribe_infinite_variables!(
    trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (idx, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.InfiniteVariable)
        # get the basic variable information
        var = object.variable
        base_name = object.name
        param_nums = var.parameter_nums
        # prepare for iterating over its supports
        supp_indices = support_index_iterator(trans_model, var.object_nums)
        vrefs = Vector{JuMP.VariableRef}(undef, length(supp_indices))
        labels = Vector{Set{DataType}}(undef, length(supp_indices))
        lookup_dict = Dict{Vector{Float64}, Int}()
        # create a variable for each support
        for (counter, i) in enumerate(supp_indices)
            raw_supp = index_to_support(trans_model, i)
            supp = raw_supp[param_nums]
            info = _format_infinite_info(var, supp)
            var_name = string(base_name, "(support: ", counter, ")")
            @inbounds vrefs[counter] = JuMP.add_variable(trans_model,
                                                         JuMP.ScalarVariable(info),
                                                         var_name)
            lookup_dict[supp] = counter
            @inbounds labels[counter] = index_to_labels(trans_model, i)
        end
        # save the transcription information
        ivref = InfiniteOpt._make_variable_ref(inf_model, idx)
        data = transcription_data(trans_model)
        data.infvar_lookup[ivref] = lookup_dict
        data.infvar_mappings[ivref] = vrefs
        data.infvar_support_labels[ivref] = labels
    end
    return
end

# Return a proper scalar variable info object given on with a start value function
function _format_derivative_info(d::InfiniteOpt.Derivative,
    support::Vector{Float64}
    )::JuMP.VariableInfo
    # generate the start value
    info = d.info
    if d.is_vector_start
        start = info.start(support)
    else
        prefs = InfiniteOpt.raw_parameter_refs(d.variable_ref)
        start = info.start(Tuple(support, prefs)...)
    end
    # make the info and return
    return JuMP.VariableInfo(info.has_lb, info.lower_bound, info.has_ub,
                             info.upper_bound, info.has_fix, info.fixed_value,
                             info.has_start, start, info.binary, info.integer)
end

"""
    transcribe_derivative_variables!(trans_model::JuMP.Model,
                                     inf_model::InfiniteOpt.InfiniteModel)::Nothing

Create transcription variables (i.e., JuMP variables) for each `Derivative`
stored in `inf_model` and add them to `trans_model`. The variable mappings are
also stored in `TranscriptionData.infvar_mappings` in accordance with
`TranscriptionData.infvar_lookup` which enable [`transcription_variable`](@ref)
and [`lookup_by_support`](@ref). Note that the supports will not be generated
until `InfiniteOpt.variable_supports` is invoked via `InfiniteOpt.supports`. The
futher derivative evaluation constraints are added when
`transcribe_derivative_evaluations!` is invoked. Note that
`TranscriptionData.infvar_support_labels` is also populated.
"""
function transcribe_derivative_variables!(
    trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (idx, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.Derivative)
        # get the basic variable information
        dref = InfiniteOpt._make_variable_ref(inf_model, idx)
        d = object.variable
        base_name = InfiniteOpt.variable_string(JuMP.REPLMode, dispatch_variable_ref(dref))
        param_nums = InfiniteOpt._parameter_numbers(d.variable_ref)
        obj_nums = InfiniteOpt._object_numbers(d.variable_ref)
        # prepare for iterating over its supports
        supp_indices = support_index_iterator(trans_model, obj_nums)
        vrefs = Vector{JuMP.VariableRef}(undef, length(supp_indices))
        labels = Vector{Set{DataType}}(undef, length(supp_indices))
        lookup_dict = Dict{Vector{Float64}, Int}()
        # create a variable for each support
        for (counter, i) in enumerate(supp_indices)
            raw_supp = index_to_support(trans_model, i)
            supp = raw_supp[param_nums]
            info = _format_derivative_info(d, supp)
            deriv_name = string(base_name, "(support: ", counter, ")")
            @inbounds vrefs[counter] = JuMP.add_variable(trans_model,
                                                         JuMP.ScalarVariable(info),
                                                         deriv_name)
            lookup_dict[supp] = counter
            @inbounds labels[counter] = index_to_labels(trans_model, i)
        end
        # save the transcription information
        data = transcription_data(trans_model)
        data.infvar_lookup[dref] = lookup_dict
        data.infvar_mappings[dref] = vrefs
        data.infvar_support_labels[dref] = labels
    end
    return
end

# Setup the mapping for a given semi_infinite variable
function _set_semi_infinite_variable_mapping(
    trans_model::JuMP.Model,
    var::InfiniteOpt.SemiInfiniteVariable,
    rvref::InfiniteOpt.GeneralVariableRef,
    index_type
    )::Nothing
    param_nums = var.parameter_nums
    ivref = var.infinite_variable_ref
    ivref_param_nums = InfiniteOpt._parameter_numbers(ivref)
    eval_supps = var.eval_supports
    # prepare for iterating over its supports
    supp_indices = support_index_iterator(trans_model, var.object_nums)
    vrefs = Vector{JuMP.VariableRef}(undef, length(supp_indices))
    labels = Vector{Set{DataType}}(undef, length(supp_indices))
    lookup_dict = Dict{Vector{Float64}, Int}()
    counter = 1
    # map a variable for each support
    for i in supp_indices
        raw_supp = index_to_support(trans_model, i)
        # ensure this support is valid with the reduced restriction
        if any(!isnan(raw_supp[ivref_param_nums[k]]) && raw_supp[ivref_param_nums[k]] != v
               for (k, v) in eval_supps)
            continue
        end
        # map to the current transcription variable
        supp = raw_supp[param_nums]
        ivref_supp = [haskey(eval_supps, j) ? eval_supps[j] : raw_supp[k]
                      for (j, k) in enumerate(ivref_param_nums)]
        @inbounds vrefs[counter] = lookup_by_support(trans_model, ivref, ivref_supp)
        lookup_dict[supp] = counter
        @inbounds labels[counter] = index_to_labels(trans_model, i)
        counter += 1
    end
    # truncate vrefs if any supports were skipped because of dependent parameter supps
    deleteat!(vrefs, counter:length(vrefs))
    deleteat!(labels, counter:length(vrefs))
    # save the transcription information
    data = transcription_data(trans_model)
    data.infvar_lookup[rvref] = lookup_dict
    data.infvar_mappings[rvref] = vrefs
    data.infvar_support_labels[rvref] = labels
    return
end

# Empty mapping dispatch for infinite parameter functions
function _set_semi_infinite_variable_mapping(
    trans_model::JuMP.Model,
    var::InfiniteOpt.SemiInfiniteVariable,
    rvref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{InfiniteOpt.ParameterFunctionIndex}
    )::Nothing
    return
end

"""
    transcribe_semi_infinite_variables!(trans_model::JuMP.Model,
                                  inf_model::InfiniteOpt.InfiniteModel)::Nothing

Map each `SemiInfiniteVariable` in `inf_model` to transcription variables stored in
`trans_model`. The variable mappings are also stored in
`TranscriptionData.infvar_mappings` in accordance with
`TranscriptionData.infvar_lookup` which enable [`transcription_variable`](@ref)
and [`lookup_by_support`](@ref). Note that [`transcribe_infinite_variables!`](@ref)
must be called first. Note that the supports will not be generated
until `InfiniteOpt.variable_supports` is invoked via `InfiniteOpt.supports`.
Note that `TranscriptionData.infvar_support_labels` is also populated.
"""
function transcribe_semi_infinite_variables!(
    trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (idx, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.SemiInfiniteVariable)
        # get the basic variable information
        var = object.variable
        rvref = InfiniteOpt._make_variable_ref(inf_model, idx)
        # setup the mappings
        ivref = InfiniteOpt.infinite_variable_ref(rvref)
        _set_semi_infinite_variable_mapping(trans_model, var, rvref,
                                            InfiniteOpt._index_type(ivref))
    end
    return
end

# Override the info of the jump variable with the point variable's if it is different
function _update_point_info(
    gvref::InfiniteOpt.GeneralVariableRef,
    vref::JuMP.VariableRef
    )::Nothing
    pvref = InfiniteOpt.dispatch_variable_ref(gvref)
    if JuMP.has_lower_bound(pvref)
        if JuMP.is_fixed(vref)
            JuMP.unfix(vref)
        end
        JuMP.set_lower_bound(vref, JuMP.lower_bound(pvref))
    end
    if JuMP.has_upper_bound(pvref)
        if JuMP.is_fixed(vref)
            JuMP.unfix(vref)
        end
        JuMP.set_upper_bound(vref, JuMP.upper_bound(pvref))
    end
    if JuMP.is_fixed(pvref)
        JuMP.fix(vref, JuMP.fix_value(pvref), force = true)
    end
    if JuMP.is_binary(pvref) && !JuMP.is_binary(vref)
        if JuMP.is_integer(vref)
            JuMP.unset_integer(vref)
        end
        JuMP.set_binary(vref)
    elseif JuMP.is_integer(pvref) && !JuMP.is_integer(vref)
        if JuMP.is_binary(vref)
            JuMP.unset_binary(vref)
        end
        JuMP.set_integer(vref)
    end
    if JuMP.start_value(pvref) != JuMP.start_value(vref)
        JuMP.set_start_value(vref, JuMP.start_value(pvref))
    end
    return
end

"""
    transcribe_point_variables!(trans_model::JuMP.Model,
                               inf_model::InfiniteOpt.InfiniteModel)::Nothing

Map each `PointVariable` in `inf_model` to a transcription variable stored in
`trans_model`. The variable mapping is also stored in
`TranscriptionData.finvar_mappings` which enables
[`transcription_variable`](@ref) and [`lookup_by_support`](@ref). Note that
[`transcribe_infinite_variables!`](@ref) must be called first and that the
info constraints associated with the transcription variable will be updated in
accordance with the point variable.
"""
function transcribe_point_variables!(
    trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (idx, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.PointVariable)
        # get the basic variable information
        var = object.variable
        ivref = var.infinite_variable_ref
        supp = var.parameter_values
        # find the corresponding variable record the mapping
        vref = lookup_by_support(trans_model, ivref, supp)
        pvref = InfiniteOpt._make_variable_ref(inf_model, idx)
        transcription_data(trans_model).finvar_mappings[pvref] = vref
        # update the info constraints as needed
        _update_point_info(pvref, vref)
    end
    return
end

################################################################################
#                       TRANSCRIPTION EXPRESSION METHODS
################################################################################
"""
    transcription_expression(trans_model::JuMP.Model, expr, support::Vector{Float64})

Given the `expr` from an `InfiniteModel`, form its transcripted version in
accordance with the variable mappings available in `trans_model` defined at
`support`. This should only be used once all variables and measures have been
transcribed (e.g., via [`transcribe_finite_variables!`](@ref)).
"""
function transcription_expression(trans_model::JuMP.Model, expr, support)
    error("Unsupported expression type `$(typeof(expr))` for automated " *
          "transcription.")
end

# GeneralVariableRef
function transcription_expression(
    trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    support::Vector{Float64}
    )
    return transcription_expression(trans_model, vref,
                                    InfiniteOpt._index_type(vref), support)
end

# Infinite variables, infinite parameter functions, and measures
function transcription_expression(
    trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{V},
    support::Vector{Float64}
    ) where {V <: Union{InfVarIndex, InfiniteOpt.ParameterFunctionIndex, InfiniteOpt.MeasureIndex}}
    param_nums = InfiniteOpt._parameter_numbers(vref)
    return lookup_by_support(trans_model, vref, index_type, support[param_nums])
end

# Semi-Infinite variables
function transcription_expression(
    trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{InfiniteOpt.SemiInfiniteVariableIndex},
    support::Vector{Float64}
    )
    ivref = InfiniteOpt.infinite_variable_ref(vref)
    if InfiniteOpt._index_type(ivref) == InfiniteOpt.ParameterFunctionIndex
        prefs = InfiniteOpt.raw_parameter_refs(ivref)
        param_nums = InfiniteOpt._parameter_numbers(ivref)
        func = InfiniteOpt.raw_function(ivref)
        return func(Tuple(support[param_nums], prefs)...)
    else
        param_nums = InfiniteOpt._parameter_numbers(vref)
        return lookup_by_support(trans_model, vref, index_type, support[param_nums])
    end
end

# Point variables and finite variables
function transcription_expression(
    trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{V},
    support::Vector{Float64}
    ) where {V <: FinVarIndex}
    return lookup_by_support(trans_model, vref, index_type, support)
end

# Infinite parameters
function transcription_expression(
    trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{V},
    support::Vector{Float64}
    ) where {V <: InfiniteOpt.InfiniteParameterIndex}
    param_num = InfiniteOpt._parameter_number(vref)
    return support[param_num]
end

# Finite parameters
function transcription_expression(
    trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    index_type::Type{InfiniteOpt.FiniteParameterIndex},
    support::Vector{Float64}
    )
    return InfiniteOpt.parameter_value(vref)
end

# AffExpr and QuadExpr
function transcription_expression(
    trans_model::JuMP.Model,
    expr::Union{JuMP.GenericAffExpr, JuMP.GenericQuadExpr},
    support::Vector{Float64}
    )
    # TODO fix this temporary hack (need to handle NLP expressions better)
    try
        return InfiniteOpt.map_expression(
            v -> transcription_expression(trans_model, v, support),
            expr)
    catch
        return transcription_expression(trans_model,
            convert(InfiniteOpt.NLPExpr, expr),
            support)
    end
end

# NLPExpr
function transcription_expression(
    trans_model::JuMP.Model,
    nlp::InfiniteOpt.NLPExpr,
    support::Vector{Float64}
    )
    ast = InfiniteOpt.map_nlp_to_ast(
        v -> transcription_expression(trans_model, v, support),
        nlp)
    return JuMP.add_NL_expression(trans_model, ast)
end

# Real Number
function transcription_expression(
    trans_model::JuMP.Model,
    num::Real,
    support::Vector{Float64}
    )
    return zero(JuMP.AffExpr) + num
end

################################################################################
#                         MEASURE TRANSCRIPTION METHODS
################################################################################
"""
    transcribe_measures!(trans_model::JuMP.Model,
                         inf_model::InfiniteOpt.InfiniteModel)::Nothing

For each `Measure` in `inf_model` expand it via `InfiniteOpt.expand_measure` or
`analytic_expansion` as appropriate and transcribe the expanded expression via
[`transcription_expression`](@ref). Then store the measure to transcripted
expression mappings in `TranscriptionData.measure_mappings` and
`TranscriptionData.measure_lookup` to enable [`transcription_variable`](@ref)
and [`lookup_by_support`](@ref). Note that the supports will not be generated
until `InfiniteOpt.variable_supports` is invoked via `InfiniteOpt.supports`.
Note that `TranscriptionData.measure_support_labels` is also populated.
"""
function transcribe_measures!(
    trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )
    for (idx, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.Measure)
        # get the basic information
        meas = object.measure
        # expand the measure
        if meas.constant_func
            new_expr = InfiniteOpt.analytic_expansion(meas.func, meas.data, trans_model)
        else
            new_expr = InfiniteOpt.expand_measure(meas.func, meas.data, trans_model)
        end
        # prepare to transcribe over the supports
        supp_indices = support_index_iterator(trans_model, meas.object_nums)
        exprs = Vector{JuMP.AbstractJuMPScalar}(undef, length(supp_indices))
        labels = Vector{Set{DataType}}(undef, length(supp_indices))
        lookup_dict = Dict{Vector{Float64}, Int}()
        # map a variable for each support
        for (counter, i) in enumerate(supp_indices)
            raw_supp = index_to_support(trans_model, i)
            @inbounds exprs[counter] = transcription_expression(trans_model,
                                                                new_expr, raw_supp)
            supp = raw_supp[meas.parameter_nums]
            lookup_dict[supp] = counter
            @inbounds labels[counter] = index_to_labels(trans_model, i)
        end
        # save the transcription information
        mref = InfiniteOpt._make_variable_ref(inf_model, idx)
        data = transcription_data(trans_model)
        data.measure_lookup[mref] = lookup_dict
        data.measure_mappings[mref] = exprs
        data.measure_support_labels[mref] = labels
    end
    return
end

################################################################################
#                        OBJECTIVE TRANSCRIPTION METHODS
################################################################################
## Dispatch functions for setting the objective
# Normal Expr
function _set_objective(trans_model, sense, expr)
    return JuMP.set_objective(trans_model, sense, expr)
end

# NonlinearExpression
function _set_objective(trans_model, sense, expr::JuMP.NonlinearExpression)
    return JuMP.set_NL_objective(trans_model, sense, expr)
end

"""
    transcribe_objective!(trans_model::JuMP.Model,
                          inf_model::InfiniteOpt.InfiniteModel)::Nothing

Form the transcripted version of the objective stored in `inf_model` and add it
to `trans_model`. Note that all the variables and measures in `inf_model` must
by transcripted first (e.g., via [`transcribe_infinite_variables!`](@ref)).
"""
function transcribe_objective!(
    trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    expr = JuMP.objective_function(inf_model)
    sense = JuMP.objective_sense(inf_model)
    trans_expr = transcription_expression(trans_model, expr, Float64[])
    _set_objective(trans_model, sense, trans_expr)
    return
end

################################################################################
#                        CONSTRAINT TRANSCRIPTION METHODS
################################################################################
## Given a variable and its set from an info constraint, get the transcribed version
# Lower bound constraint
function _get_info_constr_from_var(
    trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    set::MOI.GreaterThan,
    support::Vector{Float64}
    )::Union{JuMP.ConstraintRef, Nothing}
    trans_vref = transcription_expression(trans_model, vref, support)
    return JuMP.has_lower_bound(trans_vref) ? JuMP.LowerBoundRef(trans_vref) : nothing
end

# Upper bound constraint
function _get_info_constr_from_var(
    trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    set::MOI.LessThan,
    support::Vector{Float64}
    )::Union{JuMP.ConstraintRef, Nothing}
    trans_vref = transcription_expression(trans_model, vref, support)
    return JuMP.has_upper_bound(trans_vref) ? JuMP.UpperBoundRef(trans_vref) : nothing
end

# Fix constraint
function _get_info_constr_from_var(
    trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    set::MOI.EqualTo,
    support::Vector{Float64}
    )::Union{JuMP.ConstraintRef, Nothing}
    trans_vref = transcription_expression(trans_model, vref, support)
    return JuMP.is_fixed(trans_vref) ? JuMP.FixRef(trans_vref) : nothing
end

# Binary constraint
function _get_info_constr_from_var(
    trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    set::MOI.ZeroOne,
    support::Vector{Float64}
    )::Union{JuMP.ConstraintRef, Nothing}
    trans_vref = transcription_expression(trans_model, vref, support)
    return JuMP.is_binary(trans_vref) ? JuMP.BinaryRef(trans_vref) : nothing
end

# Integer constraint
function _get_info_constr_from_var(
    trans_model::JuMP.Model,
    vref::InfiniteOpt.GeneralVariableRef,
    set::MOI.Integer,
    support::Vector{Float64}
    )::Union{JuMP.ConstraintRef, Nothing}
    trans_vref = transcription_expression(trans_model, vref, support)
    return JuMP.is_integer(trans_vref) ? JuMP.IntegerRef(trans_vref) : nothing
end

# Determine if a given raw support satisfies constraint domain restrictions
function _support_in_restrictions(
    support::Vector{Float64},
    indices::Vector{Int},
    domains::Vector{InfiniteOpt.IntervalDomain}
    )::Bool
    for i in eachindex(indices)
        s = support[indices[i]]
        if !isnan(s) && (s < JuMP.lower_bound(domains[i]) ||
            s > JuMP.upper_bound(domains[i]))
            return false
        end
    end
    return true
end

## Process constraint objects at a particular support value and make the new
## transcribed version
# JuMP.ScalarConstraint
function _process_constraint(
    trans_model::JuMP.Model,
    constr::JuMP.ScalarConstraint,
    func::JuMP.AbstractJuMPScalar,
    set::MOI.AbstractScalarSet,
    raw_supp::Vector{Float64},
    name::String
    )
    new_func = transcription_expression(trans_model, func, raw_supp)
    trans_constr = JuMP.build_constraint(error, new_func, set)
    return JuMP.add_constraint(trans_model, trans_constr, name)
end

# MOI.LessThan expr
function _make_constr_ast(ref, set::MOI.LessThan)
    return :($ref <= $(set.upper))
end

# MOI.LessGreat expr
function _make_constr_ast(ref, set::MOI.GreaterThan)
    return :($ref >= $(set.lower))
end

# MOI.EqualTo expr
function _make_constr_ast(ref, set::MOI.EqualTo)
    return :($ref == $(set.value))
end

# MOI.Interval expr
function _make_constr_ast(ref, set::MOI.Interval)
    return :($(set.lower) <= $ref <= $(set.upper))
end

# MOI.Set fallback
function _make_constr_ast(ref, set)
    error("TranscriptionOpt does not support constraint sets of type " *
          "`$(typeof(set))` for general nonlinear constraints because this " *
          "is not yet supported by JuMP.")
end

# JuMP.ScalarConstraint with NLPExpr
function _process_constraint(
    trans_model::JuMP.Model,
    constr::JuMP.ScalarConstraint,
    func::InfiniteOpt.NLPExpr,
    set::MOI.AbstractScalarSet,
    raw_supp::Vector{Float64},
    name::String
    )
    nlp_ref = transcription_expression(trans_model, func, raw_supp)
    return JuMP.add_NL_constraint(trans_model, _make_constr_ast(nlp_ref, set))
end

# JuMP.VectorConstraint
function _process_constraint(
    trans_model::JuMP.Model,
    constr::JuMP.VectorConstraint,
    func::Vector{<:JuMP.AbstractJuMPScalar},
    set::MOI.AbstractVectorSet,
    raw_supp::Vector{Float64},
    name::String
    )
    new_func = map(f -> transcription_expression(trans_model, f, raw_supp), func)
    if any(f -> f isa JuMP.NonlinearExpression, new_func)
        error("TranscriptionOpt does not support vector constraints of general " *
              "nonlinear expressions because this is not yet supported by JuMP.")
    end
    shape = JuMP.shape(constr)
    shaped_func = JuMP.reshape_vector(new_func, shape)
    shaped_set = JuMP.reshape_set(set, shape)
    trans_constr = JuMP.build_constraint(error, shaped_func, shaped_set)
    return JuMP.add_constraint(trans_model, trans_constr, name)
end

# Fallback
function _process_constraint(
    trans_model::JuMP.Model,
    constr,
    func,
    set,
    raw_supp,
    name
    )
    error("Transcribing constraints of type `$(typeof(constr))` is not ",
          "currently supported.")
end

"""
    transcribe_constraints!(trans_model::JuMP.Model,
                            inf_model::InfiniteOpt.InfiniteModel)::Nothing

For each constraint in `inf_model` form its transcripted version(s) and add them
to `trans_model`. The mappings are stored in `TranscriptionData.constr_mappings`
and the associated supports are stored in `TranscriptionData.constr_supports`
to enable [`transcription_constraint`](@ref) and `InfiniteOpt.constraint_supports`.
Note that variable info constraints are simply mapped to the existing info
constraints already generated along with the transcription variables. Note that
the variables and measures must all first be transcripted (e.g., via
[`transcribe_measures!`](@ref)). Note that
`TranscriptionData.constr_support_labels` is also populated.
"""
function transcribe_constraints!(
    trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    param_supps = parameter_supports(trans_model)
    for (idx, object) in inf_model.constraints
        # get the basic information
        constr = object.constraint
        func = JuMP.jump_function(constr)
        set = JuMP.moi_set(constr)
        obj_nums = object.object_nums
        cref = InfiniteOpt._make_constraint_ref(inf_model, idx)
        # prepare the iteration helpers
        supp_indices = support_index_iterator(trans_model, obj_nums)
        crefs = Vector{JuMP.ConstraintRef}(undef, length(supp_indices))
        supps = Vector{Tuple}(undef, length(supp_indices))
        labels = Vector{Set{DataType}}(undef, length(supp_indices))
        counter = 1
        # iterate over the support indices for the info constraints
        if object.is_info_constraint
            for i in supp_indices
                raw_supp = index_to_support(trans_model, i)
                info_ref = _get_info_constr_from_var(trans_model, func, set,
                                                     raw_supp)
                # not all supports may be defined if overwritten by a point variable
                if !isnothing(info_ref)
                    @inbounds crefs[counter] = info_ref
                    @inbounds supps[counter] = Tuple(param_supps[j][i[j]]
                                                     for j in obj_nums)
                    @inbounds labels[counter] = index_to_labels(trans_model, i)
                    counter += 1
                end
            end
        # iterate over the supports for regular constraints
        else
            # get basic setup information
            restrictions = InfiniteOpt.domain_restrictions(cref)
            prefs = collect(keys(restrictions))
            restrict_indices = map(p -> InfiniteOpt._parameter_number(p), prefs)
            restrict_domains = map(p -> restrictions[p], prefs)
            name = object.name
            for i in supp_indices
                raw_supp = index_to_support(trans_model, i)
                # ensure the support satisfies parameter bounds and then add it
                if _support_in_restrictions(raw_supp, restrict_indices,
                                            restrict_domains)
                    new_name = isempty(name) ? "" : string(name, "(support: ", counter, ")")
                    new_cref = _process_constraint(trans_model, constr, func,
                                                   set, raw_supp, new_name)
                    @inbounds crefs[counter] = new_cref
                    @inbounds supps[counter] = Tuple(param_supps[j][i[j]]
                                                     for j in obj_nums)
                    @inbounds labels[counter] = index_to_labels(trans_model, i)
                    counter += 1
                end
            end
        end
        # truncate the arrays in case not all the supports satisfied the bounds
        deleteat!(crefs, counter:length(crefs))
        deleteat!(supps, counter:length(supps))
        deleteat!(labels, counter:length(supps))
        # add the constraint mappings to the trans model
        data = transcription_data(trans_model)
        data.constr_mappings[cref] = crefs
        data.constr_supports[cref] = supps
        data.constr_support_labels[cref] = labels
    end
    return
end

################################################################################
#                      DERIVATIVE CONSTRAINT TRANSCRIPTION METHODS
################################################################################
"""
    transcribe_derivative_evaluations!(trans_model::JuMP.Model,
                                       inf_model::InfiniteOpt.InfiniteModel)::Nothing

Generate the auxiliary derivative evaluation equations and transcribe them
appropriately for all the derivatives in `inf_model`. These are in turn added to
`trans_model`. Note that no mapping information is recorded since the InfiniteModel
won't have any constraints that correspond to these equations. Also Note that
the variables and measures must all first be transcripted (e.g., via
[`transcribe_derivative_variables!`](@ref)).
"""
function transcribe_derivative_evaluations!(
    trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel
    )::Nothing
    for (idx, object) in InfiniteOpt._data_dictionary(inf_model, InfiniteOpt.Derivative)
        # get the basic variable information
        dref = InfiniteOpt._make_variable_ref(inf_model, idx)
        pref = dispatch_variable_ref(object.variable.parameter_ref)
        method = InfiniteOpt.derivative_method(pref)
        # generate the transcription constraints as needed
        if !InfiniteOpt.has_derivative_constraints(dref)
            # generate the evaluation expressions
            exprs = InfiniteOpt.evaluate_derivative(dref, method, trans_model)
            # prepare the iteration helpers
            param_obj_num = InfiniteOpt._object_number(pref)
            obj_nums = filter(!isequal(param_obj_num), InfiniteOpt._object_numbers(dref))
            supp_indices = support_index_iterator(trans_model, obj_nums)
            # transcribe the constraints
            set = MOI.EqualTo(0.0)
            for i in supp_indices
                raw_supp = index_to_support(trans_model, i)
                for expr in exprs
                    new_expr = transcription_expression(trans_model, expr, raw_supp)
                    trans_constr = JuMP.build_constraint(error, new_expr, set)
                    JuMP.add_constraint(trans_model, trans_constr) # TODO maybe add name?
                end
            end
        end
    end
    return
end

################################################################################
#                      INFINITEMODEL TRANSCRIPTION METHODS
################################################################################
"""
    build_transcription_model!(trans_model::JuMP.Model,
                               inf_model::InfiniteOpt.InfiniteModel;
                               [check_support_dims::Bool = true])::Nothing

Given an empty `trans_model` build it using the information stored in `inf_model`.
This is intended for a `TranscriptionModel` that serves as a internal optimizer model
of `inf_model`. This detail is important to correctly enable internally generated
semi-infinite variables during the transcription process such that `inf_model` is not
modified. Note that this will add supports to `inf_model` via
[`InfiniteOpt.fill_in_supports!`](@ref) for infinite parameters that contain
no supports. Also a warning is thrown when the transcription model contains
more than 15,000 support points to alert users when they may naively have
a few independent supports whose product quickly yields a very large grid.
For example having 3 independent parameters with 100 supports each would result
in 1,000,000 supports if all three are together in at least 1 constraint. This
behavior can be overcome using dependent parameters. The warning can be turned off
via `check_support_dims = false`.
"""
function build_transcription_model!(
    trans_model::JuMP.Model,
    inf_model::InfiniteOpt.InfiniteModel; check_support_dims::Bool = true
    )::Nothing
    # ensure there are supports to add and add them to the trans model
    InfiniteOpt.fill_in_supports!(inf_model, modify = false)
    set_parameter_supports(trans_model, inf_model)
    # check that there isn't a crazy amount of supports from taking the product
    supps = parameter_supports(trans_model)
    num_supps = isempty(supps) ? 0 : prod(length, supps)
    if check_support_dims && length(supps) > 1 && num_supps > 15000 # NOTE this is an arbitrary cutoff
        @warn("Due to necessarily considering the combinatorics of independent " *
              "parameter supports, the model will be transcripted over up to $(num_supps) " *
              "supports (if any constraint uses all the infinite parameters) "  *
              "and thus naive solution of the discretized problem may be slow. " *
              "This warning can be turned off via `check_support_dims = false`.")
    end
    # register functions as needed
    InfiniteOpt.add_registered_to_jump(trans_model, inf_model)
    # define the variables
    transcribe_finite_variables!(trans_model, inf_model)
    transcribe_infinite_variables!(trans_model, inf_model)
    transcribe_derivative_variables!(trans_model, inf_model)
    transcribe_semi_infinite_variables!(trans_model, inf_model)
    transcribe_point_variables!(trans_model, inf_model)
    transcribe_measures!(trans_model, inf_model)
    # define the objective
    transcribe_objective!(trans_model, inf_model)
    # define the constraints
    transcribe_constraints!(trans_model, inf_model)
    # define the derivative evaluation constraints
    transcribe_derivative_evaluations!(trans_model, inf_model)
    return
end

# TODO make TranscriptionModel(::InfiniteModel) --> have to deal with accessing semi_infinite vars


text:
using Yao

include("circuit.jl")
include("hamiltonians.jl")

export QNNL
export Layer

abstract type Layer{T<:Real} end

struct QNNL{T} <: Layer{T}
	encoder::YaoBlocks.ChainBlock
	transform::YaoBlocks.ChainBlock
	params::Array{T}
	Hami::Array{<:YaoBlocks.AbstractBlock, 1}
	bias::Array{T,1}
end

function QNNL{T}(he::HardwareEfficientAnsatz, nbit::Integer, in_dim::Integer,
	 			w::Vector{T}, Hami::Array{<:YaoBlocks.AbstractBlock, 1},
				b::Vector{T}) where {T}
	cir_e = YaoQDNN.encoder_circuit(he, nbit, in_dim)
	cir_t = YaoQDNN.transform_circuit(he, nbit, length(w))
	QNNL{T}(cir_e, cir_t, w, Hami, b)
end

function QNNL{T}(he::HardwareEfficientAnsatz, nbit::Integer, in_dim::Integer,
				npara::Integer, Hami::Array{<:YaoBlocks.AbstractBlock, 1};
				no_bias = false) where {T}
	w = T.((rand(npara) .- 0.5) * 2π)
	if no_bias
		b = T[]
	else
		b = T.(rand(length(Hami)) .- 0.5)
	end
	QNNL{T}(he, nbit, in_dim, w, Hami, b)
end

function QNNL{T}(he::HardwareEfficientAnsatz, nbit::Integer, in_dim::Integer,
				w::Vector{T}, Hs::String, b::Vector{T}) where {T}
	Hami = YaoBlocks.ChainBlock[]
	for c in Hs
		Hami = [Hami; genH(nbit, c)]
	end
	QNNL{T}(he, nbit, in_dim, w, Hami, b)
end

function QNNL{T}(he::HardwareEfficientAnsatz, nbit::Integer, in_dim::Integer,
				npara::Integer, Hs::String; no_bias = false) where {T}
	Hami = YaoBlocks.ChainBlock[]
	for c in Hs
		Hami = [Hami; genH(nbit, c)]
	end
	QNNL{T}(he, nbit, in_dim, npara, Hami; no_bias = no_bias)
end


text:
"""
## CircPlot

Plot circular data in a number of ways:

- `cplot_histogram`: Plot a polar histogram of angular data.
- `cplot_pts`: Plot angular data as points on the edge of a circle.
- `chistogram`: Get a set of `Plots.Shape`s which can be plotted as a polar
  histogram at any point on a map, for example.

By default, all input to routines is in radians, but all routines have a positional
argument to set `degrees` to `true` when not using radians.
"""
module CircPlot

import StatsBase

using Plots

export
    chistogram,
    cplot_histogram,
    cplot_pts

#=
Exported routines
=#
"""
    cplot_pts(a::Array, degrees=false; azimuth=false, hold=false, args...) -> ::Plots.Plot

Use Plots to plot a set of points on the circumference of the unit circle.

Use `azimuth=true` to plot angles clockwise from north.

Supply any additional arguments to `Plots.plot()` as keyword arguments; for example,
to set the size of points, do

    julia> cplot_pts(angles; markersize=20)
"""
function cplot_pts(a, degrees::Bool=false;
                   azimuth::Bool=false, axial::Bool=false, line::Bool=true,
                   padding=0.05, R=1, centre=(0,0), kwargs...)
    circum = 0:0.01:2π
    markerstyle = (0.1,:blue)
    degrees && (a = deg2rad.(a))
    azimuth && (a = π/2 .- a)
    R += padding
    x0, y0 = centre
    p = plot(aspect_ratio=:equal,
        xlim=(-R+x0,R+x0), ylim=(-R+y0,R+y0), legend=false)
    line && plot!(p, x0.+cos.(circum), y0.+sin.(circum), l=:black)
    scatter!(p, x0.+cos.(a), y0.+sin.(a), m=markerstyle; kwargs...)
    axial && scatter!(p, x0.+cos.(a+π), y0.+sin.(a+π), m=markerstyle; kwargs...)
    if azimuth
        xlabel!(p, "East")
        ylabel!(p, "North")
    else
        xlabel!(p, "x")
        ylabel!(p, "y")
    end
    p
end

"""
    cplot_histogram(a, binwidth, degrees=false; azimuth=false, axial=false, weights=ones(length(a)), circ=nothing; kwargs...) -> ::Plots.Plot

Plot a polar histogram of a set of angles `a` using bins `binwidth` wide.
`binwidth` should be given in the same convention (radians or °) as the data.
If π (or 180°) is not an integer multiple of `binwidth`, then the nearest value
for which this is true is selected.

Use `azimuth=true` to plot angles clockwise from north.

Use `axial=true` if data are π-ambiguous (i.e., are orientational rather than
directional data).

`weights` must have the same length as `a`, and if passed in will scale each
point `i` by `weights[i]`.

Set options for the surrounding circle using `circ=Dict(:opt=>val)`, where
the keys and values in the Dict will be passed to Plots for plotting the circle
shape around the plot.

Additional keyword arguments `kwargs...` are passed to `Plots.plot()` when
plotting the bars.
"""
function cplot_histogram(a, binwidth, degrees::Bool=false;
                         azimuth::Bool=false, axial::Bool=false,
                         weights=ones(eltype(a), length(a)), circ=nothing,
                         kwargs...)
    hist = fit_hist(a, binwidth, degrees, axial, weights=weights)
    n = length(hist.weights)
    edges = hist.edges[1]
    azimuth && (edges = π/2 .- collect(edges))
    R = maximum(hist.weights)
    p = plot(aspect_ratio=:equal, xlim=(-R,R), ylim=(-R,R), legend=false)
    circ != nothing && plot!(p, circle(R); circ...)
    plot!(p, sector.(edges[1:n], edges[2:n+1], hist.weights), fill=:black; kwargs...)
    p
end

"""
    chistogram(a, binwidth, degrees::Bool=true, x=0, y=0; axial=false, azimuth=false, scale=1, maxr=nothing) -> bars::Array{Shape}

Return a set of Shapes which show a polar histogram created from the set of
angles `a` with bins `binwidth` wide, centred at `x` and `y`.  These shapes can
be plotted with `plot(bars)`.

Set `axial` to `true` for data with π-ambiguity and `azimuth` to `true` if data are
clockwise from north.

`scale` determine what length to scale the radii by, or alternatively set `maxr` to
the radius (in plot units) of the maximum bin.  Set `minweight` to be the minimum
bar count required for the plot to fill the entire maximum radius from `maxr`.
"""
function chistogram(a, binwidth, degrees=false, x=0, y=0;
        weights=ones(a), axial=false, azimuth=false, scale=1, maxr=nothing,
        minweight=1, minobs=nothing)
    if minobs != nothing
        warning("`minobs` is deprecated: use `minweight` instead")
        minweight = minobs
    end
    hist = fit_hist(a, binwidth, degrees, axial, weights=weights)
    n = length(hist.weights)
    edges = hist.edges[1]
    R = maximum(hist.weights)
    maxr != nothing && (scale = maxr/max(minweight, R))
    azimuth && (edges = π/2 .- collect(edges))
    sector.(edges[1:n], edges[2:n+1], hist.weights*scale, x, y)
end

#=
Unexported routines
=#
"""    fit_hist(a, binwidth, degrees=false, axial=false; weights=ones(a))

Return a `StatsBase.Histogram` with the results of binning the angles `a`
into bins `binwidth` wide.  If `axial` is `true`, then directions have
π symmetry.  Optionally sypply a `weights` array containing the weighting
for each value in `a`."""
function fit_hist(a, binwidth, degrees::Bool, axial::Bool;
                  weights=ones(eltype(a), length(a)))
    binwidth isa Bool &&
        throw(ArgumentError("`binwidth` supplied as a `Bool`, but should be a real number"))
    length(weights) == length(a) ||
        throw(ArgumentError("`weights` must have same length as number of angles"))
    n = (degrees ? round(Int, 360/binwidth) : round(Int, 2pi/binwidth))
    n = max(1, n)
    binwidth = 2pi/n
    bins = range(0, stop=2pi, length=n+1)
    data = degrees ? deg2rad.(mod.(a, 360)) : mod.(a, 2pi)
    # FIXME: Many angles are given as whole numbers: this cludge helps
    # an issue where floating point values sometimes are binned one way,
    # sometimes another, giving polar histograms which depend on the range
    # in which angles are given.
    data .+= 10*eps(float(eltype(a)))
    axial && (data .= mod.(data, pi))
    h = StatsBase.fit(StatsBase.Histogram, data,
                      StatsBase.FrequencyWeights(weights), bins, closed=:left)
    axial && (h.weights[end÷2+1:end] .= h.weights[1:end÷2])
    h
end

"""Return a sensible number of segments between two angles in radians θ₁ and θ₂."""
number_of_segments(θ₁, θ₂) = max(2, round(Int, abs(θ₂ - θ₁)*90/π))

"""Return a circle shape with radius `r` centred at `x` and `y`."""
circle(r, x=0, y=0) = Shape(arc(0, 2π, r, x, y))

"""Return an arc (part of a circle's circumference) between `θ₁` and `θ₂` radians
anticlockwise from the +x axis, with radius `r`, centred at `x` and `y`."""
arc(θ₁, θ₂, r=1, x=0, y=0) = Tuple{Float64,Float64}[(x+r*cos(θ), y+r*sin(θ))
    for θ in range(θ₁, stop=θ₂, length=number_of_segments(θ₁, θ₂))]

"""Return a sector (pie wedge) between angles `θ₁` and `θ₂` radians anticlockwise
from the x-axis with radius `r`, centred at `x` and `y`."""
sector(θ₁, θ₂, r=1, x=0, y=0)::Shape =
    Shape(vcat(arc(θ₁, θ₂, r, x, y), Tuple{Float64,Float64}((x, y))))

"""Return a filled sector (arcshape) between two angles θ₁ and θ₂ radians
anticlockwise from the x-axis, and two radii, `r₁` and `r₂`, centred at `x` and `y`."""
arcshape(θ₁, θ₂, r₁, r₂, x=0, y=0) =
    Shape(vcat(arc(θ₁, θ₂, r₁, x, y), arc(θ₂, θ₁, r₂, x, y)))

end # module

max_stars_repo_name: anowacki/assorted-julia-modules

max_issues_repo_name: anowacki/assorted-julia-modules

max_forks_repo_name: anowacki/assorted-julia-modules


text:
using Test, StaticArrays
using CLIMA.VariableTemplates

struct TestModel{A,B,C}
  a::A
  b::B
  c::C
end

struct SubModelA
end
struct SubModelB
end
struct SubModelC{N}
end

function state(m::TestModel, T)
  @vars begin
    ρ::T
    ρu::SVector{3,T}
    ρe::T
    a::state(m.a,T)
    b::state(m.b,T)
    c::state(m.c,T)
    S::SHermitianCompact{3,T,6}
  end
end

state(m::SubModelA, T) = @vars()
state(m::SubModelB, T) = @vars(ρqt::T)
state(m::SubModelC{N}, T) where {N} = @vars(ρk::SVector{N,T})

model = TestModel(SubModelA(), SubModelB(), SubModelC{5}())

st = state(model, Float64)

@test varsize(st) == 17

v = Vars{st}(zeros(MVector{varsize(st),Float64}))
g = Grad{st}(zeros(MMatrix{3,varsize(st),Float64}))

@test v.ρ === 0.0
@test v.ρu === SVector(0.0, 0.0, 0.0)
v.ρu = SVector(1,2,3)
@test v.ρu === SVector(1.0, 2.0, 3.0)

@test v.b.ρqt === 0.0
v.b.ρqt = 12.0
@test v.b.ρqt === 12.0

@test v.S === zeros(SHermitianCompact{3,Float64,6})
v.S = SHermitianCompact{3,Float64,6}(1,2,3,2,3,4,3,4,5)
@test v.S[1,1] === 1.0
@test v.S[1,3] === 3.0
@test v.S[3,1] === 3.0
@test v.S[3,3] === 5.0

v.S = ones(SMatrix{3,3,Int64})
@test v.S[1,1] === 1.0
@test v.S[1,3] === 1.0
@test v.S[3,1] === 1.0
@test v.S[3,3] === 1.0

@test propertynames(v.a) == ()
@test propertynames(g.a) == ()

@test g.ρu == zeros(SMatrix{3,3,Float64})
g.ρu = SMatrix{3,3}(1:9)
@test g.ρu == SMatrix{3,3,Float64}(1:9)

@test size(v.c.ρk) == (5,)
@test size(g.c.ρk) == (3,5)

sv = similar(v)
@test typeof(sv) == typeof(v)
@test size(parent(sv)) == size(parent(v))

sg = similar(g)
@test typeof(sg) == typeof(g)
@test size(parent(sg)) == size(parent(g))

@test flattenednames(st) == ["ρ","ρu[1]","ρu[2]","ρu[3]","ρe",
                            "b.ρqt",
                            "c.ρk[1]","c.ρk[2]","c.ρk[3]","c.ρk[4]","c.ρk[5]",
                            "S[1,1]", "S[2,1]", "S[3,1]", "S[2,2]", "S[3,2]", "S[3,3]"]

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
__precompile__()

module LispREPL

import REPL: REPL, LineEdit
import LispSyntax, ParserCombinator

function valid_sexpr(s)
  try
    LispSyntax.read(String(take!(copy(LineEdit.buffer(s)))))
    true
  catch err
    isa(err, ParserCombinator.ParserException) || rethrow(err)
    false
  end
end

function initrepl(repl)

  text  = get(ENV, "LISP_PROMPT_TEXT", "lisp> ")
  color = Base.text_colors[Symbol(get(ENV, "LISP_PROMPT_COLOR", :magenta))]
  key   = ')'

  julia_mode = repl.interface.modes[1]

  prefix = repl.hascolor ? color : ""
  suffix = repl.hascolor ? (repl.envcolors ? Base.input_color : repl.input_color) : ""

  lisp_mode = LineEdit.Prompt(text;
    prompt_prefix    = prefix,
    prompt_suffix    = suffix,
    repl             = repl,
    on_enter         = valid_sexpr,
    complete         = REPL.REPLCompletionProvider(),
  )
  lisp_mode.on_done = REPL.respond(s -> :($(LispSyntax).@lisp($s)), repl, lisp_mode)

  push!(repl.interface.modes, lisp_mode)

  hp                     = julia_mode.hist
  hp.mode_mapping[:lisp] = lisp_mode
  lisp_mode.hist         = hp

  search_prompt, skeymap = LineEdit.setup_search_keymap(hp)

  mk = REPL.mode_keymap(julia_mode)

  lisp_keymap = Dict(
    key => (s, args...) ->
      if isempty(s) || position(LineEdit.buffer(s)) == 0
        buf = copy(LineEdit.buffer(s))
        LineEdit.transition(s, lisp_mode) do
          LineEdit.state(s, lisp_mode).input_buffer = buf
        end
      else
        LineEdit.edit_insert(s, key)
      end
  )

  lisp_mode.keymap_dict = LineEdit.keymap(Dict[
    skeymap,
    mk,
    LineEdit.history_keymap,
    LineEdit.default_keymap,
    LineEdit.escape_defaults,
  ])
  julia_mode.keymap_dict = LineEdit.keymap_merge(julia_mode.keymap_dict, lisp_keymap)

  nothing
end

function __init__()
  isdefined(Base, :active_repl) ? initrepl(Base.active_repl) : nothing
end

end # module


text:
using Test
using Config

function manager_tests()
    cfg = ConfigManager(joinpath(@__DIR__,"tdlambda.json"), @__DIR__)

    metasteps = [0.025, 0.05, 0.075, 0.1]
    lambdas = [0.0, 0.4, 0.8, 0.9]
    N = length(lambdas)
    cfgs = []
    for i=1:N
        for j=1:N
            push!(cfgs, (metasteps[j],lambdas[i]))
        end
    end

    for i=1:total_combinations(cfg)
        cfg = parse!(cfg, i)
        subconfig = get_subconfig(cfg, "agent")
        tru = cfgs[i]
        @test subconfig["metastep"] == tru[1]
        @test subconfig["lambda"] == tru[2]
    end
end

function manager_tests_toml()
    cfg = ConfigManager(joinpath(@__DIR__,"tdlambda.toml"), @__DIR__)

    metasteps = [0.025, 0.05, 0.075, 0.1]
    lambdas = [0.0, 0.4, 0.8, 0.9]
    N = length(lambdas)
    cfgs = []
    for i=1:N
        for j=1:N
            push!(cfgs, (metasteps[j],lambdas[i]))
        end
    end

    for i=1:total_combinations(cfg)
        cfg = parse!(cfg, i)
        subconfig = get_subconfig(cfg, "agent")
        tru = cfgs[i]
        @test subconfig["metastep"] == tru[1]
        @test subconfig["lambda"] == tru[2]
    end
end

@testset "Config.jl Tests:" begin
    manager_tests()
    manager_tests_toml()
end


text:
function to_rna(dna)
    function complement(n)
        n = uppercase(n)
        if n == 'G'
            return 'C'
        elseif n == 'C'
            return 'G'
        elseif n == 'T'
            return 'A'
        elseif n == 'A'
            return 'U'
        else
            throw(ErrorException("Invalid nucleotide $n in DNA strand"))
        end
    end
    return map(n -> complement(n), dna)
end

max_stars_repo_name: vladflore/exercism

max_issues_repo_name: vladflore/exercism

max_forks_repo_name: vladflore/exercism


text:
"""
    gaussiansmooth3d!(image, σ=[5,5,5]; mask=nothing, nbox=ifelse(isnothing(mask), 3, 6), weight=nothing, dims=1:min(ndims(image),3), boxsizes=getboxsizes.(σ, nbox))

Performs Gaussian smoothing on `image` with `σ` as standard deviation of the Gaussian.
By application of `nbox` times running average filters in each dimension.
The length of `σ` and the length of the `dims` that are smoothed have to match. (Default `3`)

Optional arguments:\n
`mask`: Smoothing can be performed using a mask to inter-/extrapolate missing values.\n
`nbox`: Number of box applications. Default is `3` for normal smoothing and `6` for masked smoothing.\n
`weight`: Apply weighted smoothing. Either weighted or masked smoothing can be porformed.\n
`dims`: Specify which dims should be smoothed. Corresponds to manually looping of the other dimensions.\n
`boxizes`: Manually specify the boxsizes, not using the provided σ. `length(boxsizes)==length(dims) && length(boxsizes[1])==nbox`
"""
function gaussiansmooth3d(image, σ=[5,5,5]; padding=false, kwargs...)
    if padding
        image = pad_image(image, σ)
    end
    smoothed = gaussiansmooth3d!(0f0 .+ copy(image), σ; kwargs...)
    if padding
        smoothed = remove_padding(smoothed, σ)
    end
    return smoothed
end

pad_image(image, σ) = PaddedView(0, image, Tuple(size(image) .+ 2σ), Tuple(σ .+ 1))
remove_padding(image, σ) = image[[σ[i]+1:size(image,i)-σ[i] for i in 1:ndims(image)]...]

function gaussiansmooth3d_phase(phase, σ=[5,5,5]; weight=1, kwargs...)
    return angle.(gaussiansmooth3d!(weight .* exp.(1im .* phase), σ; kwargs...))
end

function gaussiansmooth3d!(image, σ=[5,5,5]; mask=nothing, nbox=ifelse(isnothing(mask), 3, 4), weight=nothing, dims=1:min(ndims(image),3), boxsizes=getboxsizes.(σ, nbox))
    if length(σ) < length(dims) @error "Length of σ and dims does not match!" end
    if length(boxsizes) < length(dims) || length(boxsizes[1]) != nbox @error "boxsizes has wrong size!" end
    if typeof(mask) != Nothing
        image .*= ifelse.(mask .== 0, NaN, 1) # 0 in mask -> NaN in image
    end
    if typeof(weight) != Nothing
        weight = Float32.(weight)
        weight[weight .== 0] .= minimum(weight[weight .!= 0])
    end
    checkboxsizes!(boxsizes, size(image), dims)

    for ibox in 1:nbox, dim in dims
        bsize = boxsizes[dim][ibox]
        if size(image, dim) == 1 || bsize < 3
            continue
        end
        linefilter = getfilter(image, weight, mask, bsize, size(image, dim))
        K = ifelse(mask isa Nothing || isodd(ibox), :, size(image, dim):-1:1)

        # TODO parallel? -> Distributed arrays? -> use slices
        #loop = Iterators.product((size(image) |> sz -> (sz[1:(dim-1)], sz[(dim+1):end]) .|> CartesianIndices)...)
        #Threads.@threads for (I, J) in collect(loop)

        for J in CartesianIndices(size(image)[(dim+1):end])
            for I in CartesianIndices(size(image)[1:(dim-1)])
                w = if weight isa Nothing nothing else view(weight,I,:,J) end
                linefilter(view(image,I,K,J), w)
            end
        end
    end
    return image
end

## Calculate the filter sizes to achieve a given σ

function getboxsizes(σ, n)
    try
        wideal = √( (12σ^2 / n) + 1 )
        wl::Int = round(wideal - (wideal + 1) % 2) # next lower odd integer
        wu::Int = wl + 2

        mideal = (12σ^2 - n*wl.^2 - 4n*wl - 3n) / (-4wl - 4)
        m = round(mideal)

        [if i <= m wl else wu end for i in 1:n]
    catch
        zeros(n)
    end
end

function checkboxsizes!(boxsizes, sz, dims)
    for dim in dims
        bs = boxsizes[dim]
        for i in eachindex(bs)
            if iseven(bs[i])
                @warn "boxsize $i is even: $(bs[i]); it was changed to next bigger odd integer!"
                bs[i] += 1
            end
            if bs[i] > sz[dim] / 2
                val = sz[dim] ÷ 2
                if iseven(val) val += 1 end
                @warn "boxsize $i for dim $dim is limited to half the image; it was changed from $(bs[i]) to $(val)!"
                bs[i] = val
            end
        end
    end
end

## Function to initialize the filters

function getfilter(image, weight::Nothing, mask::Nothing, bsize, len)
    q = CircularBuffer{eltype(image)}(bsize)
    return (im, _) -> boxfilterline!(im, bsize, q)
end
function getfilter(image, weight, mask::Nothing, bsize, len)
    q = CircularBuffer{eltype(image)}(bsize)
    qw = CircularBuffer{eltype(weight)}(bsize)
    return (im, w) -> boxfilterline!(im, bsize, w, q, qw)
end
function getfilter(image, weight, mask, bsize, len)
    buffer = ones(eltype(image), len + bsize - 1) * NaN16
    return (im, _) -> nanboxfilterline!(im, bsize, buffer)
end

## Running Average Filters

function boxfilterline!(line::AbstractVector, boxsize::Int, q::CircularBuffer)
    r = div(boxsize, 2)
    initvals = view(line, 1:boxsize)
    lsum = sum(initvals)
    append!(q, initvals)

    @inbounds for i in (r+2):(length(line)-r)
        lsum += line[i+r] - popfirst!(q)
        push!(q, line[i+r])
        line[i] = lsum / boxsize
    end
end

function boxfilterline!(line::AbstractVector, boxsize::Int, weight::AbstractVector, lq::CircularBuffer, wq::CircularBuffer)
    r = div(boxsize, 2)

    wsmooth = wsum = sum = eps() # slightly bigger than 0 to avoid division by 0
    @inbounds for i in 1:boxsize
        sum += line[i] * weight[i]
        wsum += weight[i]
        wsmooth += weight[i]^2
        push!(lq, line[i])
        push!(wq, weight[i])
    end

    @inbounds for i in (r+2):(length(line)-r)
        w = weight[i+r]
        l = line[i+r]
        wold = popfirst!(wq)
        lold = popfirst!(lq)
        push!(wq, w)
        push!(lq, l)

        sum += l * w - lold * wold
        wsum += w - wold
        line[i] = sum / wsum
        wsmooth += w^2 - wold^2
        weight[i] = wsmooth / wsum
    end
end

function nanboxfilterline!(line::AbstractVector, boxsize::Int, orig::AbstractVector)
    n = length(line)
    r = div(boxsize, 2)
    maxfills = r

    orig[r+1:r+n] .= line
    orig[r+n+1:end] .= NaN

    lsum = sum(view(orig,r+1:2r))
    if isnan(lsum) lsum = 0. end
    nfills = 0
    nvalids = 0

    mode = :nan

    @inbounds for i in 1:length(line)
        if isnan(lsum) @warn "lsum nan"; break end

        # check for mode change
        if mode == :normal
            if isnan(orig[i+2r])
                mode = :fill
            end

        elseif mode == :nan
            if isnan(orig[i+2r])
                nvalids = 0
            else
                nvalids += 1
            end
            if nvalids == boxsize
                mode = :normal
                lsum = sum(view(orig,i:(i+2r)))
                line[i] = lsum / boxsize
                continue # skip to next loop iteration
            end

        elseif mode == :fill
            if isnan(orig[i+2r])
                nfills += 1
                if nfills > maxfills
                    mode = :nan
                    nfills = 0
                    lsum = 0
                    nvalids = 0
                end
            else
                mode = :normal
                nfills = 0
            end
        end

        # perform operation
        if mode == :normal
            lsum += orig[i+2r] - orig[i-1]
            line[i] = lsum / boxsize
        elseif mode == :fill
            lsum -= orig[i-1]
            line[i] = (lsum - orig[i]) / (boxsize - 2)

            orig[i+2r] = 2line[i] - line[i-r] # TODO maybe clamp the value
            if (i+r < n) line[i+r] = orig[i+2r] end
            lsum += orig[i+2r]
        end

    end
end


text:
function solve(days, counts)
    for day in 1:days
        first = popfirst!(counts)
        push!(counts, 0)
        counts[7] += first
        counts[9] += first
    end
    return sum(counts)
end

nums = (x -> parse(Int, x)).(split(readline(), ','))
counts = [sum(nums .== i) for i in 0:9]

println(solve(80, copy(counts)))
println(solve(256, copy(counts)))

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
include("benchmarks.jl")

include("plots.jl")

results = run(SUITE)
data = graphbench_data(results)
plot_all_subcats(data)

max_forks_repo_name: epatters/CompCat


text:
export CUDNNError

struct CUDNNError <: Exception
    code::cudnnStatus_t
end

Base.convert(::Type{cudnnStatus_t}, err::CUDNNError) = err.code

Base.showerror(io::IO, err::CUDNNError) =
    print(io, "CUDNNError: ", name(err), " (code $(reinterpret(Int32, err.code)))")

name(err::CUDNNError) = unsafe_string(cudnnGetErrorString(err))

## API call wrapper

# outlined functionality to avoid GC frame allocation
@noinline function throw_api_error(res)
    if res == CUDNN_STATUS_ALLOC_FAILED
        throw(OutOfGPUMemoryError())
    else
        throw(CUDNNError(res))
    end
end

macro check(ex, errs...)
    check = :(isequal(err, CUDNN_STATUS_ALLOC_FAILED))
    for err in errs
        check = :($check || isequal(err, $(esc(err))))
    end

    quote
        res = @retry_reclaim err->$check $(esc(ex))
        if res != CUDNN_STATUS_SUCCESS
            throw_api_error(res)
        end

        nothing
    end
end


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

module DateTests

#include("periods.jl")

for file in readlines(joinpath(@__DIR__, "testgroups"))
    include(file * ".jl")
end

end


text:
using RHEOS
using BenchmarkTools
using InverseLaplace
import MittagLeffler
using FunctionWrappers: FunctionWrapper

t=RheoFloat.(0.0:0.001:1.0);
t0=1.0

println("=========================")
println("Test on maths functions ")
println("=========================")

println()
println("Test on the MittalgLeffler function")
println()
println("Version from Julia repo")
@btime MittagLeffler.mittleff(0.5, t0)
println("Version from RHEOS (should be the same as above)")
@btime RHEOS.mittleff(0.5, t0)

println()
println()

println()
println("Test on the inverse Laplace transform")
println()
cₐ = 1.0; a=0.5; kᵦ = 2.0; kᵧ=3.0
println("Talbot algorithm on FractSLS_Zener creep function")
Ĵ(s) = (1/s)*(cₐ*s^a + kᵦ)/(cₐ*s^a*kᵦ + kᵧ*(cₐ*s^a + kᵦ))
print("  scalar value: ")
@btime InverseLaplace.talbot(s->(1/s)*(cₐ*s^a + kᵦ)/(cₐ*s^a*kᵦ + kᵧ*(cₐ*s^a + kᵦ)), t0)
print("  vector using talbot.: ")
@btime InverseLaplace.talbot.(s->(1/s)*(cₐ*s^a + kᵦ)/(cₐ*s^a*kᵦ + kᵧ*(cₐ*s^a + kᵦ)), t)
print("  vector using talbotarr: ")
@btime InverseLaplace.talbotarr(s -> Ĵ(s), t)
println()

println("Weeks (80) algorithm on FractSLS_Zener creep function")
print("  Initialisation: ")
@btime InverseLaplace.Weeks(Ĵ,80)
ft = InverseLaplace.Weeks(Ĵ,80)
print("  Scalar value: ")
@btime ft(t0)
print("  Vector values: ")
@btime ft.(t)
println()
println()

println("=========================")
println("Test on Moduli functions ")
println("=========================")

println()
println("Test on the Maxwell model")
println()
@btime relaxmod(Maxwell, t0, [2.0,3.0])
@btime relaxmod(Maxwell, t, [2.0,3.0]);

m = RheoModel(Maxwell,  (η=2.0, k=3.0))
@btime relaxmod(m, t0)
@btime relaxmod(m, t);

println()
println("Test on the FractSLS_Zener model relaxation modulus - with mittleff")
println()
cₐ = 1.0; a=0.5; kᵦ = 2.0; kᵧ=3.0
println("Direct function evaluation")
@btime kᵦ*MittagLeffler.mittleff(a, -(kᵦ/cₐ)*t0^a) + kᵧ
@btime kᵦ.*MittagLeffler.mittleff.(a, -(kᵦ./cₐ).*t.^a) .+ kᵧ;
println("From the RheoModelClass")
@btime relaxmod(FractSLS_Zener, t0, [1.0,0.5,2.0,3.0])
@btime relaxmod(FractSLS_Zener, t, [1.0,0.5,2.0,3.0])
m = RheoModel(FractSLS_Zener,  (cₐ = 1.0, a=0.5, kᵦ = 2.0, kᵧ=3.0))
println("From the RheoModel")
@btime relaxmod(m, t0)
@btime relaxmod(m, t)

println()
println("Test on the FractSLS_Zener model creep compliance - with inverse laplace")
println()

fwp = ((t,p)->(InverseLaplace.Weeks(s->(1/s)*(p[1]*s^p[2] + p[3])/(p[1]*s^p[2]*p[3] + p[4]*(p[1]*s^p[2] + p[3])),80))(t)) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat,Vector{RheoFloat}}}
fwpa = ((t,p)->(InverseLaplace.Weeks(s->(1/s)*(p[1]*s^p[2] + p[3])/(p[1]*s^p[2]*p[3] + p[4]*(p[1]*s^p[2] + p[3])),80)).(t)) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat},Vector{RheoFloat}}}
fw = InverseLaplace.Weeks(s->(1/s)*(1.0*s^0.5 + 2.0)/(1.0*s^0.5*2.0 + 3.0*(1.0*s^0.5 + 2.0)),80) |> FunctionWrapper{RheoFloat,Tuple{RheoFloat}}
fwa = t->fw.(t) |> FunctionWrapper{Vector{RheoFloat},Tuple{Vector{RheoFloat}}}
println()
println("function wrappers using Weeks")
@btime fwp(t0,[1.0,0.5,2.0,3.0])
@btime fwpa(t,[1.0,0.5,2.0,3.0])
@btime fw(t0)
@btime fwa(t)

println()
println("From the RheoModelClass")
@btime creepcomp(FractSLS_Zener, t0, [1.0,0.5,2.0,3.0])
@btime creepcomp(FractSLS_Zener, t, [1.0,0.5,2.0,3.0])
m = RheoModel(FractSLS_Zener,  (cₐ = 1.0, a=0.5, kᵦ = 2.0, kᵧ=3.0))
println("From the RheoModel")
@btime creepcomp(m, t0)
@btime creepcomp(m, t)
nothing

#=

Hardware:
Intel Core i5-7200U Processor (3MB Cache, up to 3.1Ghz)
8GB LPDDR3 1866

=========================
Test on maths functions
=========================

Test on the MittalgLeffler function

Version from Julia repo
  57.567 ns (1 allocation: 16 bytes)
Version from RHEOS (should be the same as above)
  58.195 ns (1 allocation: 16 bytes)

Test on the inverse Laplace transform

Talbot algorithm on FractSLS_Zener creep function
  scalar value:   21.206 μs (642 allocations: 18.75 KiB)
  vector using talbot.:   21.383 ms (642143 allocations: 18.33 MiB)
  vector using talbotarr:   3.418 ms (173413 allocations: 3.61 MiB)

Weeks (80) algorithm on FractSLS_Zener creep function
  Initialisation:   105.563 μs (2634 allocations: 110.75 KiB)
  Scalar value:   1.616 μs (1 allocation: 16 bytes)
  Vector values:   1.596 ms (2 allocations: 8.03 KiB)

=========================
Test on Moduli functions
=========================

Test on the Maxwell model

  67.185 ns (2 allocations: 112 bytes)
  7.559 μs (2 allocations: 8.09 KiB)
  33.592 ns (1 allocation: 16 bytes)
  8.153 μs (1 allocation: 8.00 KiB)

Test on the FractSLS_Zener model relaxation modulus - with mittleff

Direct function evaluation
  207.000 ns (7 allocations: 112 bytes)
  68.707 μs (15 allocations: 8.53 KiB)
From the RheoModelClass
  166.950 ns (4 allocations: 160 bytes)
  64.459 μs (2 allocations: 8.11 KiB)
From the RheoModel
  104.457 ns (1 allocation: 16 bytes)
  59.557 μs (1 allocation: 8.00 KiB)

Test on the FractSLS_Zener model creep compliance - with inverse laplace

function wrappers using Weeks
  62.121 μs (71 allocations: 30.80 KiB)
  1.660 ms (72 allocations: 38.81 KiB)
  1.636 μs (1 allocation: 16 bytes)
  1.610 ms (5 allocations: 8.12 KiB)

From the RheoModelClass
  10.674 μs (101 allocations: 2.64 KiB)
  10.721 ms (99102 allocations: 2.47 MiB)
From the RheoModel
  9.535 μs (1 allocation: 16 bytes)
  9.580 ms (1 allocation: 8.00 KiB)

=#


text:
@recipe function f(H::OpenQuantumBase.AbstractHamiltonian, s, lvl)
    y = []
    for x in s
        w, _ = eigen_decomp(H, x; lvl=lvl)
        push!(y, w)
    end
    y = hcat(y...)'
    lab = ["\$E_{$x}\$" for x in (0:lvl - 1)']
    label --> lab
    yguide --> "\$\\mathrm{GHz}\$"
    xguide --> "\$s\$"
    legend --> :best
    (s, y)
end


text:
# Project Files for Magnetic Rossby Wave project

#
# include project function files
include("pWaves.jl")
include("mgrid.jl")
include("modes_x.jl")
include("modes_y.jl")
include("transform.jl")
include("wave_equation.jl")
include("derivatives.jl")
include("ModelParameters.jl")
include("velocity.jl")
include("realization.jl")
include("divergence.jl")

max_stars_repo_name: bbuffett/MagRossby

max_issues_repo_name: bbuffett/MagRossby

max_forks_repo_name: bbuffett/MagRossby