text:
"""
    AbstractPattern{T} <: AbstractArray{T, 2}

`AbstractPatterns` are image-like array types which can be used to color
plottable objects. There are currently two subtypes: `LinePattern` and
`ImagePattern`. Any abstract pattern must implement the `to_image(pat)`
function, which must return a `Matrix{<: AbstractRGB}`.
"""
abstract type AbstractPattern{T} <: AbstractArray{T, 2} end

function Base.show(io::IO, ::MIME"text/plain", p::AbstractPattern)
    println(io, typeof(p))
    show(io, mime, image(p))
end
struct ImagePattern{T<:Colorant} <: AbstractPattern{T}
    img::Matrix{T}
end

"""
    Pattern(image)
    Pattern(mask[; color1, color2])

Creates an `ImagePattern` from an `image` (a matrix of colors) or a `mask`
(a matrix of real numbers). The pattern can be passed as a `color` to a plot to
texture it. If a `mask` is passed, one can specify to colors between which colors are
interpolated.
"""
Pattern(img::Array{<: Colorant, 2}) = ImagePattern(img)

function Pattern(mask::Matrix{<: Real}; color1=RGBA(0,0,0,1), color2=RGBA(1,1,1,0))
    img = map(x -> to_color(color1) * x + to_color(color2) * (1-x), mask)
    return ImagePattern(img)
end

to_image(p::ImagePattern) = p.img

struct LinePattern{T} <: AbstractPattern{T}
    dirs::Vector{Vec2}
    widths::Vector{AbstractFloat}
    shifts::Vector{Vec2}

    tilesize::NTuple{2, Integer}
    colors::NTuple{2, T}
end

"""
    LinePattern([; kwargs...])

Creates a `LinePattern` for the given keyword arguments:
- `direction`: The direction of the line.
- `width`: The width of the line
- `tilesize`: The size of the image on which the line is drawn. This should be
compatible with the direction.
- `shift`: Sets the starting point for the line.
- `linecolor`: The color with which the line is replaced.
- `background_color`:: The background color.

Multiple `direction`s, `width`s and `shift`s can also be given to create more
complex patterns, e.g. a cross-hatching pattern.
"""
function LinePattern(;
        direction = Vec2f0(1), width = 2f0, tilesize = (10,10),
        shift = map(w -> Vec2f0(0.5 - 0.5(w%2)), width),
        linecolor = RGBA(0,0,0,1), background_color = RGBA(1,1,1,0)
    )
    N = 1
    direction isa Vector{<:Vec2} && (N = length(direction))
    width isa Vector && (length(width) > N) && (N = length(width))
    shift isa Vector{<:Vec2} && (length(shift) > N) && (N = length(shift))

    dirs = direction isa Vector{<:Vec2} ? direction : [direction for _ in 1:N]
    widths = width isa Vector ? width : [width for _ in 1:N]
    shifts = shift isa Vector{<:Vec2} ? shift : [shift for _ in 1:N]
    colors = (to_color(linecolor), to_color(background_color))

    return LinePattern{eltype(colors)}(dirs, widths, shifts, tilesize, colors)
end

"""
    Pattern(style::String = "/"; kwargs...)
    Pattern(style::Char = '/'; kwargs...)

Creates a line pattern based on the given argument. Available patterns are
`'/'`, `'\\'`, `'-'`, `'|'`, `'x'`, and `'+'`. All keyword arguments correspond
to the keyword arguments for [`LinePattern`](@ref).
"""
Pattern(style::String; kwargs...) = Pattern(style[1]; kwargs...)
function Pattern(style::Char = '/'; kwargs...)
    if style == '/'
        LinePattern(direction=Vec2f0(1); kwargs...)
    elseif style == '\\'
        LinePattern(direction=Vec2f0(1, -1); kwargs...)
    elseif style == '-'
        LinePattern(direction=Vec2f0(1, 0); kwargs...)
    elseif style == '|'
        LinePattern(direction=Vec2f0(0, 1); kwargs...)
    elseif style == 'x'
        LinePattern(direction=[Vec2f0(1), Vec2f0(1, -1)]; kwargs...)
    elseif style == '+'
        LinePattern(direction=[Vec2f0(1, 0), Vec2f0(0, 1)]; kwargs...)
    else
        LinePattern(; kwargs...)
    end
end

function to_image(p::LinePattern)
    tilesize = p.tilesize
    full_mask = zeros(tilesize...)
    for (dir, width, shift) in zip(p.dirs, p.widths, p.shifts)
        mask = zeros(tilesize...)
        width -= 1 # take away center
        if abs(dir[1]) < abs(dir[2])
            # m = dx / dy; x = m * (y-1) + 1
            m = dir[1] / dir[2]
            for y in 1:tilesize[2]
                cx = m * (y-shift[2]) + shift[1]
                r = floor(Int64, cx-0.5width):ceil(Int64, cx+0.5width)
                for x in r[2:end-1]
                    mask[mod1(x, tilesize[1]), y] = 1.0
                end
                mask[mod1(r[1], tilesize[1]), y] = 1 - abs(cx-0.5width - r[1])
                mask[mod1(r[end], tilesize[1]), y] = 1 - abs(cx+0.5width - r[end])
            end
        else
            # m = dy / dx; y = m * (x-1) + 1
            m = dir[2] / dir[1]
            for x in 1:tilesize[1]
                cy = m * (x-shift[1]) + shift[2]
                r = floor(Int64, cy-0.5width):ceil(Int64, cy+0.5width)
                for y in r[2:end-1]
                    mask[x, mod1(y, tilesize[2])] = 1.0
                end
                mask[x, mod1(r[1], tilesize[2])] = 1 - abs(cy-0.5width - r[1])
                mask[x, mod1(r[end], tilesize[2])] = 1 - abs(cy+0.5width - r[end])
            end
        end
        full_mask .+= mask
    end

    return map(full_mask) do x
        p.colors[1] * clamp(x, 0, 1) + p.colors[2] * (1-clamp(x, 0, 1))
    end
end


text:
# === PRIVATE METHODS BELOW HERE ====================================================================================== #
function _basic_prediction(signalVector::Array{Int64,1}, strategyObject::VLAbstractGameStrategy)::Int64

    # we use binary 0,1 under the covers - so lets get rid of the -1 and replace w/0 -
    memory = length(signalVector)
    binarySignal = replace(signalVector, -1 => 0)
    iv = range(memory, stop=1, step=-1) |> collect

    # convert the binarySignal to an int -
    strategy_index = 0
    for index in iv
        value = binarySignal[index]
        strategy_index += value * 2^(memory - index)
    end

    # from the strategy object, get the actual impl of the strategy -
    strategyImpl = strategyObject.strategy

    # compute the predicted outcome -
    return strategyImpl[strategy_index + 1]
end

function _basic_minority(gameAgentArray::Array{VLBasicMinorityGameAgent,1}, signalVector::Array{Int64,1})::NamedTuple

    # iniialize -
    dim_output_array = Array{Int64,1}()
    numberOfAgents = length(gameAgentArray)
    actions = [-1,1]

    # ask the agents what their prediction will be, given this signal vector -
    agentPredictionArray = Array{Int64,1}(undef, numberOfAgents)
    for agent_index = 1:numberOfAgents

        # grab the agent -
        agentObject = gameAgentArray[agent_index]

        # get the strategy object -
        strategy_object = agentObject.bestAgentStrategy

        # what is predcited?
        predicted_action = _basic_prediction(signalVector, strategy_object)

        # grab -
        agentPredictionArray[agent_index] = predicted_action
    end

    # process each element of the alphabet -
    for value in actions

        # number of items -
        number_of_values = length(findall(x -> x == value, agentPredictionArray))

        # grab -
        push!(dim_output_array, number_of_values)
    end

    # find the argmin -
    arg_min_index = argmin(dim_output_array)
    minority_action_value = actions[arg_min_index]
    collective_action = sum(agentPredictionArray)

    # setup the return tuple -
    return_tuple = (winningAction = minority_action_value, sell = dim_output_array[1], buy = dim_output_array[2], collectiveAction = collective_action, agentActions = agentPredictionArray)

    # return -
    return return_tuple
end

function _gc_minority(gameAgentArray::Array{VLGCMinorityGameAgent,1}, signalVector::Array{Int64,1})::NamedTuple

     # iniialize -
    dim_output_array = Array{Int64,1}()
    numberOfAgents = length(gameAgentArray)
    actions = [-1,1]

    # ask the agents what their prediction will be, given this signal vector -
    agentPredictionArray = Array{Int64,1}(undef, numberOfAgents)
    for agent_index = 1:numberOfAgents

        # grab the agent -
        agentObject = gameAgentArray[agent_index]
        old_score = agentObject.score

         # compute a probability that the agen will skip this round, based on their score -
        P = 1 / (1 + exp(-0.01 * old_score)) # P of trading this round -
        #P = (1/√(2*π))*exp(-0.5*old_score^2)
        r = rand()
        if (P > r) # the agent trades -

            # get the strategy object -
            strategy_object = agentObject.bestAgentStrategy

            # what is predcited?
            predicted_action = _basic_prediction(signalVector, strategy_object)

            # grab -
            agentPredictionArray[agent_index] = predicted_action
        else

            # hold action - agent did not trade -
            agentPredictionArray[agent_index] = 0
        end
    end

    # process each element of the alphabet -
    for value in actions

        # number of items -
        number_of_values = length(findall(x -> x == value, agentPredictionArray))

        # grab -
        push!(dim_output_array, number_of_values)
    end

    # find the argmin -
    arg_min_index = argmin(dim_output_array)
    minority_action_value = actions[arg_min_index]
    collective_action = sum(agentPredictionArray)

    # setup the return tuple -
    buy_count = dim_output_array[2]
    sell_count = dim_output_array[1]
    hold_count = numberOfAgents - (buy_count + sell_count)
    return_tuple = (winningAction = minority_action_value, sell = sell_count, buy = buy_count, hold = hold_count,
        collectiveAction = collective_action, agentActions = agentPredictionArray)

    # return -
    return return_tuple
end

function _thermal_minority(gameAgentArray::Array{VLThermalMinorityGameAgent,1},
    signalVector::Array{Int64,1})::NamedTuple

    # iniialize -
    dim_output_array = Array{Int64,1}()
    numberOfAgents = length(gameAgentArray)
    actions = [-1,1]

    # ask the agents what their prediction will be, given this signal vector -
    agentPredictionArray = Array{Int64,1}(undef, numberOfAgents)
    for agent_index = 1:numberOfAgents

        # grab the agent -
        agentObject = gameAgentArray[agent_index]

        # Pick a strategy -
        agentStrategyCollection = agentObject.agentStrategyCollection
        strategy_probability_array = Array{Float64,1}()
        strategy_score_array = Array{Int64,1}()
        for strategy_object in agentStrategyCollection
            push!(strategy_probability_array, strategy_object.probability)
            push!(strategy_score_array, strategy_object.score)
        end

        # roll a rand number -
        r = rand()

        # how many strategies does this agent have?
        number_of_strategies = length(strategy_probability_array)

        # sort from the probabilities from smallest, to largest -
        idx_sort_probabilities = sortperm(strategy_probability_array)
        strategy_probability_array = strategy_probability_array[idx_sort_probabilities]

        # now, compute the partial sum -
        partial_sum_value = 0.0
        partial_sum_array = Array{Float64,1}()
        for strategy_index = 1:number_of_strategies

            partial_sum_value += strategy_probability_array[strategy_index]
            push!(partial_sum_array, partial_sum_value)
        end

        # compute the selected strategy index, and get the selected strategy -
        idx_selected_strategy = findfirst(x -> x >= r, partial_sum_array)
        if (isnothing(idx_selected_strategy) == true)
            @show r2, partial_sum_array, agentObject
        end

        actual_strategy_index = idx_sort_probabilities[idx_selected_strategy]
        strategy_object = agentObject.agentStrategyCollection[actual_strategy_index]

        # what is predcited?
        predicted_action = _basic_prediction(signalVector, strategy_object)

        # grab -
        agentPredictionArray[agent_index] = predicted_action
    end

    # process each element of the alphabet -
    for value in actions

        # number of items -
        number_of_values = length(findall(x -> x == value, agentPredictionArray))

        # grab -
        push!(dim_output_array, number_of_values)
    end

    # find the argmin -
    arg_min_index = argmin(dim_output_array)
    minority_action_value = actions[arg_min_index]
    collective_action = sum(agentPredictionArray)

    # setup the return tuple -
    return_tuple = (winningAction = minority_action_value, sell = dim_output_array[1], buy = dim_output_array[2], collectiveAction = collective_action, agentActions = agentPredictionArray)

    # return -
    return return_tuple
end

function _thermal_cooling_function(temperature::Float64)

    # default behavior -
    α = 0.9
    temperature =  α * temperature
    return temperature
end

# === PRIVATE METHODS ABOVE HERE ====================================================================================== #

# === PUBLIC METHODS BELOW HERE ======================================================================================= #
function execute(worldObject::VLThermalMinorityGameWorld, numberOfTimeSteps::Int64;
    voteManagerFunction::Function=_thermal_minority, predictionFuntion::Function=_basic_prediction,
    coolingFunction::Function=_thermal_cooling_function, temperatureMinimum::Float64=0.1,
    liquidity::Float64=10001.0, σ::Float64=0.0005)::VLMinorityGameSimulationResult

    try

        # get parameters from the gameworld -
        temperature = worldObject.temperature
        numberOfAgents = worldObject.numberOfAgents
        agentMemorySize = worldObject.agentMemorySize
        gameAgentArray = worldObject.gameAgentArray

        # actions -
        actions = [-1,1]    # we always use the binary model
        length_of_actions = length(actions)

        # initialize data structures -
        gameWorldMemoryBuffer = CircularBuffer{Int64}(agentMemorySize)
        game_state_table = DataFrame(temperature=Float64[], sell=Int[], buy=Int[], winner=Int[], A=Int[], price=Float64[])
        agent_state_table = Array{Int,2}(undef, numberOfTimeSteps, numberOfAgents)
        logAssetPriceArray = Array{Float64,1}(undef, (numberOfTimeSteps + 1))

        # initialze gameWorldMemoryBuffer -
        ra = rand(1:length_of_actions, (agentMemorySize + 1))
        for x = 1:(agentMemorySize + 1)
            push!(gameWorldMemoryBuffer, actions[ra[x]])
        end

        # initially asset price is 1.0 -
        logAssetPriceArray[1] = 1.0

        # add normal to perturbation -
        d = Normal(0, σ)

        # main run loop -
        # outer loop: we repeat this loop while T > T_min
        while (temperature > temperatureMinimum)

            # lets define the beta -
            β = (1 / temperature)

            # inner loop: we repeat this loop for 1 -> numberOfTimeSteps
            for time_step_index = 1:numberOfTimeSteps

                # grab the last agentMemorySize block -
                signalVector = gameWorldMemoryBuffer[(length(gameWorldMemoryBuffer) - (agentMemorySize - 1)):end]

                # ok, so let's ask all the agents what they voted to do, and get data on the minority position -
                results_tuple = voteManagerFunction(gameAgentArray, signalVector)
                winning_action = results_tuple.winningAction
                collective_action = results_tuple.collectiveAction
                sell = results_tuple.sell
                buy = results_tuple.buy
                agentPredictionArray = results_tuple.agentActions

                # update the price -
                old_price = logAssetPriceArray[time_step_index]
                r = ((1 / liquidity) * collective_action) + rand(d)
                logAssetPriceArray[time_step_index + 1] = old_price * (exp(r))

                # update the agents data with the winning outcome -
                for agent_index = 1:numberOfAgents

                    # grab the agent -
                    agentObject = gameAgentArray[agent_index]

                    # what did this agent predict?
                    agentPrediction = agentPredictionArray[agent_index]

                    # capture the current score, then update -
                    # if the agent predicted correctly, then increase personal score -
                    agent_state_table[time_step_index, agent_index] = agentObject.score
                    agentObject.score += -1 * sign(agentPrediction * collective_action)

                    # so now we need to compute the probabaility for each strategy for this
                    # agent -
                    agentStrategyCollection = agentObject.agentStrategyCollection
                    numberOfStrategiesPerAgent = length(agentStrategyCollection)
                    factor_array = Array{Float64,1}(undef, numberOfStrategiesPerAgent)
                    for sindex = 1:numberOfStrategiesPerAgent

                        strategy_object = agentStrategyCollection[sindex]

                        # update the score for this strategy -
                        strategy_prediction = predictionFuntion(signalVector, strategy_object)

                        # update the score -
                        strategy_object.score += -1 * sign(strategy_prediction * collective_action)
                        if (strategy_object.score > 2 * agentMemorySize)
                            strategy_object.score = 2 * agentMemorySize
                        elseif (strategy_object.score < -2 * agentMemorySize)
                            strategy_object.score = -2 * agentMemorySize
                        end

                        # compute the factors -
                        factor_array[sindex] = exp(β * strategy_object.score)
                    end

                    # ok, so know that we that updated scores, lets compute the facrors -
                    denom_value = sum(factor_array)
                    for agent_strategy_index = 1:numberOfStrategiesPerAgent

                        strategy_object = agentStrategyCollection[agent_strategy_index]

                        # compute the probabaility -
                        probability = factor_array[agent_strategy_index] / denom_value
                        if (probability == 0.0)
                            @show factor_array, agent_strategy_index, strategy_object.score
                        end

                        # for this agent, get the strategy -
                        strategy_object.probability = probability
                    end
                end # end agent update for loop

                # add the winning outcome to the system memory -
                push!(gameWorldMemoryBuffer, winning_action)

                # populate state table -
                data_row = [temperature, sell, buy, winning_action, collective_action, logAssetPriceArray[time_step_index]]
                push!(game_state_table, data_row)

            end # end time for loop

            # after each time range, update the temperature using the cooling function
            temperature = coolingFunction(temperature)
        end # end Temp while loop

        # package up and return -
        return VLMinorityGameSimulationResult(game_state_table)

    catch error
        rethrow(error)
    end
end

function execute(worldObject::VLBasicMinorityGameWorld, numberOfTimeSteps::Int64;
    voteManagerFunction::Function=_basic_minority, predictionFuntion::Function=_basic_prediction,
    liquidity::Float64=10001.0, σ::Float64=0.0005)::NamedTuple

    try

        # setup -
        numberOfAgents = worldObject.numberOfAgents
        agentMemorySize = worldObject.agentMemorySize
        gameAgentArray = worldObject.gameAgentArray
        actions = [-1,1]    # we always use the binary model
        length_of_actions = length(actions)

        # initialize data structures -
        gameWorldMemoryBuffer = CircularBuffer{Int64}(agentMemorySize)
        game_state_table = DataFrame(sell=Int[], buy=Int[], winner=Int[], A=Int[], price=Float64[])
        agent_state_table = Array{Int,2}(undef, numberOfTimeSteps, numberOfAgents)
        logAssetPriceArray = Array{Float64,1}(undef, (numberOfTimeSteps + 1))

        # add normal to perturbation -
        d = Normal(0, σ)

        # initialize -
        for _ = 1:(agentMemorySize + 1)

            # generate a random number {0,1}
            r = rand(1:length_of_actions)

            # capture -
            push!(gameWorldMemoryBuffer, actions[r])
        end

        # initially asset price is 1.0 -
        logAssetPriceArray[1] = 1.0

        # main loop -
        for time_step_index = 1:numberOfTimeSteps

            # grab the last agentMemorySize block -
            signalVector = gameWorldMemoryBuffer[(length(gameWorldMemoryBuffer) - (agentMemorySize - 1)):end]

            # ok, so let's ask all the agents what they voted to do, and get data on the minority position -
            results_tuple = voteManagerFunction(gameAgentArray, signalVector)
            winning_action = results_tuple.winningAction
            collective_action = results_tuple.collectiveAction
            sell = results_tuple.sell
            buy = results_tuple.buy
            agentPredictionArray = results_tuple.agentActions

            # update the price -
            old_price = logAssetPriceArray[time_step_index]
            r = ((1 / liquidity) * collective_action) + rand(d)
            logAssetPriceArray[time_step_index + 1] = old_price * (exp(r))

            # update the agents data with the winning outcome -
            for agent_index = 1:numberOfAgents

                # grab the agent -
                agentObject = gameAgentArray[agent_index]

                # what did this agent predict?
                agentPrediction = agentPredictionArray[agent_index]

                # old score -
                old_score = agentObject.score
                agent_state_table[time_step_index, agent_index] = old_score

                # if the agent predicted correctly, then increase personal score -
                agentObject.score += -1 * sign(agentPrediction * collective_action)

                # let's re-rank all the strategies for this agent, and then put the new scores in an array that we can sort them
                agentStrategyCollection = agentObject.agentStrategyCollection
                tmp_score_array = Array{Int64,1}()
                for strategy_object in agentStrategyCollection

                    # compute the prediction for this strategy -
                    strategy_prediction = predictionFuntion(signalVector, strategy_object)

                    # update the score -
                    old_score = strategy_object.score
                    ΔS = -1 * sign(strategy_prediction * collective_action)
                    new_score = old_score + ΔS
                    strategy_object.score = new_score

                    # cache new score -
                    push!(tmp_score_array, new_score)
                end

                # sort the scores -
                idx_sort_score = sortperm(tmp_score_array)

                # ok, so grab the best strategy, and update the best strategy pointer
                agentObject.bestAgentStrategy = agentStrategyCollection[last(idx_sort_score)]
            end

            # add the winning outcome to the system memory -
            push!(gameWorldMemoryBuffer, winning_action)

            # populate state table -
            data_row = [sell, buy, winning_action, collective_action, logAssetPriceArray[time_step_index]]
            push!(game_state_table, data_row)
        end

        # return tuple -
        return_tuple = (market_table = game_state_table, agent_score_table = agent_state_table)

        # return -
        return return_tuple
    catch error
        # if we get here ... something bad has happend. do nothing for now
        rethrow(error)
    end
end

function execute(worldObject::VLGCMinorityGameWorld; kwargs...)::VLMinorityGameSimulationResult

    try

        # ok: so we need a bunch of stuff. if this stuff is not in the kwargs dictionary, then we need to provide default values -
        kwargs_dictionary = Dict(kwargs)
        numberOfTimeSteps = get(kwargs_dictionary,:numberOfTimeSteps, 100) # check: number of steps -
        voteManagerFunction = get(kwargs_dictionary, :voteManagerFunction, _gc_minority) # check: voteManagerFunction -
        predictionFuntion = get(kwargs_dictionary, :predictionFuntion, _basic_prediction) # check: predictionFuntion -
        liquidity = get(kwargs_dictionary, :liquidity, 10001.0)  # check: liquidity -
        σ = get(kwargs_dictionary, :σ, 0.0005) # check: sigma -
        agent_score_refresh_rate = get(kwargs_dictionary, :agentScoreRefreshRate, 50) # check: sigma -

        # initialize -
        numberOfAgents = worldObject.numberOfAgents
        agentMemorySize = worldObject.agentMemorySize
        gameAgentArray = worldObject.gameAgentArray
        actions = [-1,1]    # we always use the binary model
        length_of_actions = length(actions)

        # initialize data structures -
        gameWorldMemoryBuffer = CircularBuffer{Int64}(agentMemorySize)
        game_state_table = DataFrame(sell=Int[], buy=Int[], hold=[], winner=Int[], A=Int[], price=Float64[])
        agent_state_table = Array{Int,2}(undef, numberOfTimeSteps, numberOfAgents)
        logAssetPriceArray = Array{Float64,1}(undef, (numberOfTimeSteps + 1))

        # add normal to perturbation -
        d = Normal(0, σ)

        # initialize -
        for _ = 1:(agentMemorySize + 1)

            # generate a random number {0,1}
            r = rand(1:length_of_actions)

            # capture -
            push!(gameWorldMemoryBuffer, actions[r])
        end

        # initially asset price is 1.0 -
        logAssetPriceArray[1] = 1.0

        # main loop -
        for time_step_index = 1:numberOfTimeSteps

            # grab the last agentMemorySize block -
            signalVector = gameWorldMemoryBuffer[(length(gameWorldMemoryBuffer) - (agentMemorySize - 1)):end]

            # ok, so let's ask all the agents what they voted to do, and get data on the minority position -
            results_tuple = voteManagerFunction(gameAgentArray, signalVector)
            winning_action = results_tuple.winningAction
            collective_action = results_tuple.collectiveAction
            sell = results_tuple.sell
            buy = results_tuple.buy
            hold = results_tuple.hold
            agentPredictionArray = results_tuple.agentActions

            # update the price -
            old_price = logAssetPriceArray[time_step_index]
            r = ((1 / liquidity) * collective_action) + rand(d)
            logAssetPriceArray[time_step_index + 1] = old_price * (exp(r))

            # update the agents data with the winning outcome -
            for agent_index = 1:numberOfAgents

                # grab the agent -
                agentObject = gameAgentArray[agent_index]

                # get the agents score -
                old_score = agentObject.score
                agent_state_table[time_step_index, agent_index] = old_score

                # what did this agent predict?
                agentPrediction = agentPredictionArray[agent_index]

                # if the agent predicted correctly, then increase personal score -
                if (mod(time_step_index,agent_score_refresh_rate) == 0)
                    agentObject.score = 0
                else
                    agentObject.score += -1 * sign(agentPrediction * collective_action)
                end

                # let's re-rank all the strategies for this agent, and then put the new scores in an array that we can sort them
                agentStrategyCollection = agentObject.agentStrategyCollection
                tmp_score_array = Array{Int64,1}()
                for strategy_object in agentStrategyCollection

                    # compute the prediction for this strategy -
                    strategy_prediction = predictionFuntion(signalVector, strategy_object)

                    # update the score -
                    old_score = strategy_object.score
                    ΔS = -1 * sign(strategy_prediction * collective_action)
                    new_score = old_score + ΔS
                    strategy_object.score = new_score

                    # cache new score -
                    push!(tmp_score_array, new_score)
                end

                # sort the scores -
                idx_sort_score = sortperm(tmp_score_array)

                # ok, so grab the best strategy, and update the best strategy pointer
                agentObject.bestAgentStrategy = agentStrategyCollection[last(idx_sort_score)]
            end

            # add the winning outcome to the system memory -
            push!(gameWorldMemoryBuffer, winning_action)

            # populate state table -
            data_row = [sell, buy, hold, winning_action, collective_action, logAssetPriceArray[time_step_index]]
            push!(game_state_table, data_row)
        end

        # package up and return -
        return VLMinorityGameSimulationResult(game_state_table)
    catch error
        rethrow(error)
    end
end
# === PUBLIC METHODS ABOVE HERE ======================================================================================= #

ext: jl


text:
export crelu

doc"""
    crelu(x::Var)

Concatenated Rectified Linear Unit.
The output is twice the size of the input.

```math
f(x) = (\max(0,x), \max(0,-x))
```

# References
* Shang et al., ["Understanding and Improving Convolutional Neural Networks via Concatenated Rectified Linear Units"](https://arxiv.org/abs/1603.05201), arXiv 2016.
"""
crelu(x::Var) = Var(crelu(x.data), (crelu,x))
crelu(x::Node) = Node(crelu, x)

function crelu(x::Array{T}) where T
    y = Array{T}(2size(x,1), Base.tail(size(x))...)
    @inbounds for i = 1:length(x)
        k = (i-1)*2 + 1
        y[k] = max(x[i], zero(T))
        y[k+1] = max(-x[i], zero(T))
    end
    y
end

function addgrad!(y::Var, ::typeof(crelu), x::Var)
    isvoid(x.grad) && return
    ∇crelu!(y.grad, x.data, x.grad)
end

function ∇crelu!(gy::Array{T}, x::Array{T}, gx::Array{T}) where T
    @inbounds for i = 1:length(x)
        k = (i-1)*2 + 1
        gx[i] += x[i] > zero(T) ? gy[k] : -gy[k+1]
    end
end


text:
abstract type AbstractPointFunction end

struct PointFunction{F} <: AbstractPointFunction
    f::F
end

struct FilterPointFunction{F,T} <: AbstractPointFunction
    f::F
    default::T
end

@inline function (vf::AbstractPointFunction)(
    m::AbstractMetricParams{T},
    gp::GeodesicBase.AbstractGeodesicPoint{T},
    max_time;
    kwargs...,
)::T where {T}
    convert(T, vf.f(m, gp, max_time; kwargs...))
end

@inline function apply(
    vf::AbstractPointFunction,
    rc::SolutionRenderCache{M,T,G};
    kwargs...,
) where {M,T,G}
    map(sol -> vf.f(rc.m, get_endpoint(m, sol), rc.max_time; kwargs...), rc.geodesics)
end

@inline function apply(
    vf::AbstractPointFunction,
    rc::EndpointRenderCache{M,T,P};
    kwargs...,
) where {M,T,P}
    map(gp -> vf.f(rc.m, gp, rc.max_time; kwargs...), rc.points)
end

@inline function Base.:∘(vf1::AbstractPointFunction, vf2::AbstractPointFunction)
    PointFunction((m, gp, max_time; kwargs...) -> vf1.f(vf2.f(m, gp, max_time; kwargs...)))
end

@inline function Base.:∘(vf1::AbstractPointFunction, vf2::FilterPointFunction{F}) where {F}
    PointFunction(
        (m, gp, max_time; kwargs...) -> begin
            pass_on = vf2.f(m, gp, max_time; kwargs...)
            if pass_on
                vf1.f(m, gp, max_time; kwargs...)
            else
                vf2.default
            end
        end,
    )
end

module ConstPointFunctions
import ..GeodesicRendering: PointFunction, FilterPointFunction

const filter_early_term =
    FilterPointFunction((m, gp, max_time; kwargs...) -> gp.t < max_time, NaN)

const filter_intersected =
    FilterPointFunction((m, gp, max_time; kwargs...) -> gp.retcode == :Intersected, NaN)

const affine_time = PointFunction((m, gp, max_time; kwargs...) -> gp.t)

const shadow = affine_time ∘ filter_early_term
end # module

export PointFunction, FilterPointFunction, ConstPointFunctions


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg
using Base.BinaryPlatforms

const YGGDRASIL_DIR = "../.."
include(joinpath(YGGDRASIL_DIR, "fancy_toys.jl"))
include(joinpath(YGGDRASIL_DIR, "platforms", "llvm.jl"))

name = "faust"
version = v"2.40.10"

# Collection of sources required to build faust
sources = [
    GitSource("https://github.com/grame-cncm/faust.git",
              "04bc38cd56ba5f5b10af1c86d0534364f0d5cf62"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/faust

CMAKE_FLAGS=()
CMAKE_FLAGS+=(-DINCLUDE_HTTP=OFF)
CMAKE_FLAGS+=(-DINCLUDE_ITP=OFF)
CMAKE_FLAGS+=(-DINCLUDE_OSC=OFF)
CMAKE_FLAGS+=(-DINCLUDE_STATIC=OFF)
CMAKE_FLAGS+=(-DITPDYNAMIC=OFF)

CMAKE_TARGET=${target}

if [[ "${target}" == *musl* ]]; then
    export CXXFLAGS="-DALPINE"
fi

if [[ "${target}" == *apple* ]]; then
    export MACOSX_DEPLOYMENT_TARGET=10.9
    export LDFLAGS="${LDFLAGS} -mmacosx-version-min=10.9"
    # If we're building for Apple, CMake gets confused with `aarch64-apple-darwin` and instead prefers
    # `arm64-apple-darwin`.  If this issue persists, we may have to change our triplet printing.
    if [[ "${target}" == aarch64* ]]; then
        CMAKE_TARGET=arm64-${target#*-}
    fi
fi

if [[ "${target}" == *mingw* ]]; then
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/ws2.patch
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/ws2_libraries.patch
fi

if [[ "${target}" == *freebsd* ]]; then
    export LDFLAGS="${LDFLAGS} -lexecinfo"
fi

atomic_patch -p1 ${WORKSPACE}/srcdir/patches/unset_llvm_libs.patch

if [[ "${bb_full_target}" == x86_64-linux-musl-* ]]; then
    llvm_lib=$(basename /workspace/destdir/lib/libLLVM-??jl.so)
    LLVM_LIB="-l${llvm_lib:3:9}"
elif [[ "${bb_full_target}" == *mingw*+13 ]]; then
    LLVM_LIB="-lLLVM-13jl.dll"
elif [[ "${target}" == *mingw* ]]; then
    LLVM_LIB="-lLLVM.dll"
else
    LLVM_LIB="-lLLVM"
fi

LLVM_LIBS="${LLVM_LIB} -lstdc++"

CMAKE_FLAGS+=(-DLLVM_LIBS=\'${LLVM_LIBS}\')

CMAKE_FLAGS+=(-DCMAKE_C_COMPILER_TARGET=${CMAKE_TARGET})
CMAKE_FLAGS+=(-DCMAKE_CXX_COMPILER_TARGET=${CMAKE_TARGET})

CMAKE_FLAGS+=(-DCMAKE_INSTALL_PREFIX=${prefix})
CMAKE_FLAGS+=(-DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN})
CMAKE_FLAGS+=(-DCMAKE_BUILD_TYPE=Release)

CMAKE_FLAGS+=(-DUSE_LLVM_CONFIG=OFF)
CMAKE_FLAGS+=(-DLLVM_DIR=${prefix}/lib/cmake/llvm)

export CMAKEOPT="${CMAKE_FLAGS[@]}"
export PREFIX=$prefix

make -j${nproc} all
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms(; experimental=true))

# The products that we will ensure are always built
products = [
    ExecutableProduct("faust", :faust),
    LibraryProduct("libfaust", :libfaust),
    LibraryProduct("libHTTPDFaust", :libHTTPDFaust),
    LibraryProduct("libOSCFaust", :libOSCFaust),

    # Generated list of Faust scripts using:
    # find $prefix/bin -maxdepth 1 -name 'faust2*' -printf '    FileProduct(joinpath("bin", "%f"), :%f),\n' | sort
    FileProduct(joinpath("bin", "faust2alqt"), :faust2alqt),
    FileProduct(joinpath("bin", "faust2alsaconsole"), :faust2alsaconsole),
    FileProduct(joinpath("bin", "faust2alsa"), :faust2alsa),
    FileProduct(joinpath("bin", "faust2android"), :faust2android),
    FileProduct(joinpath("bin", "faust2androidunity"), :faust2androidunity),
    FileProduct(joinpath("bin", "faust2api"), :faust2api),
    FileProduct(joinpath("bin", "faust2atomsnippets"), :faust2atomsnippets),
    FileProduct(joinpath("bin", "faust2au"), :faust2au),
    FileProduct(joinpath("bin", "faust2bela"), :faust2bela),
    FileProduct(joinpath("bin", "faust2caqt"), :faust2caqt),
    FileProduct(joinpath("bin", "faust2caqtios"), :faust2caqtios),
    FileProduct(joinpath("bin", "faust2csound"), :faust2csound),
    FileProduct(joinpath("bin", "faust2csvplot"), :faust2csvplot),
    FileProduct(joinpath("bin", "faust2dplug"), :faust2dplug),
    FileProduct(joinpath("bin", "faust2dssi"), :faust2dssi),
    FileProduct(joinpath("bin", "faust2dummy"), :faust2dummy),
    FileProduct(joinpath("bin", "faust2dummymem"), :faust2dummymem),
    FileProduct(joinpath("bin", "faust2eps"), :faust2eps),
    FileProduct(joinpath("bin", "faust2esp32"), :faust2esp32),
    FileProduct(joinpath("bin", "faust2faustvst"), :faust2faustvst),
    FileProduct(joinpath("bin", "faust2firefox"), :faust2firefox),
    FileProduct(joinpath("bin", "faust2gen"), :faust2gen),
    FileProduct(joinpath("bin", "faust2graph"), :faust2graph),
    FileProduct(joinpath("bin", "faust2graphviewer"), :faust2graphviewer),
    FileProduct(joinpath("bin", "faust2ios"), :faust2ios),
    FileProduct(joinpath("bin", "faust2jackconsole"), :faust2jackconsole),
    FileProduct(joinpath("bin", "faust2jack"), :faust2jack),
    FileProduct(joinpath("bin", "faust2jackrust"), :faust2jackrust),
    FileProduct(joinpath("bin", "faust2jackserver"), :faust2jackserver),
    FileProduct(joinpath("bin", "faust2jaqtchain"), :faust2jaqtchain),
    FileProduct(joinpath("bin", "faust2jaqt"), :faust2jaqt),
    FileProduct(joinpath("bin", "faust2juce"), :faust2juce),
    FileProduct(joinpath("bin", "faust2ladspa"), :faust2ladspa),
    FileProduct(joinpath("bin", "faust2linuxunity"), :faust2linuxunity),
    FileProduct(joinpath("bin", "faust2lv2"), :faust2lv2),
    FileProduct(joinpath("bin", "faust2mathdoc"), :faust2mathdoc),
    FileProduct(joinpath("bin", "faust2mathviewer"), :faust2mathviewer),
    FileProduct(joinpath("bin", "faust2max6"), :faust2max6),
    FileProduct(joinpath("bin", "faust2md"), :faust2md),
    FileProduct(joinpath("bin", "faust2msp"), :faust2msp),
    FileProduct(joinpath("bin", "faust2netjackconsole"), :faust2netjackconsole),
    FileProduct(joinpath("bin", "faust2netjackqt"), :faust2netjackqt),
    FileProduct(joinpath("bin", "faust2nodejs"), :faust2nodejs),
    FileProduct(joinpath("bin", "faust2octave"), :faust2octave),
    FileProduct(joinpath("bin", "faust2osxiosunity"), :faust2osxiosunity),
    FileProduct(joinpath("bin", "faust2owl"), :faust2owl),
    FileProduct(joinpath("bin", "faust2paqt"), :faust2paqt),
    FileProduct(joinpath("bin", "faust2pdf"), :faust2pdf),
    FileProduct(joinpath("bin", "faust2plot"), :faust2plot),
    FileProduct(joinpath("bin", "faust2png"), :faust2png),
    FileProduct(joinpath("bin", "faust2portaudiorust"), :faust2portaudiorust),
    FileProduct(joinpath("bin", "faust2puredata"), :faust2puredata),
    FileProduct(joinpath("bin", "faust2pure"), :faust2pure),
    FileProduct(joinpath("bin", "faust2raqt"), :faust2raqt),
    FileProduct(joinpath("bin", "faust2ros"), :faust2ros),
    FileProduct(joinpath("bin", "faust2rosgtk"), :faust2rosgtk),
    FileProduct(joinpath("bin", "faust2rpialsaconsole"), :faust2rpialsaconsole),
    FileProduct(joinpath("bin", "faust2rpinetjackconsole"), :faust2rpinetjackconsole),
    FileProduct(joinpath("bin", "faust2sam"), :faust2sam),
    FileProduct(joinpath("bin", "faust2sc"), :faust2sc),
    FileProduct(joinpath("bin", "faust2sig"), :faust2sig),
    FileProduct(joinpath("bin", "faust2sigviewer"), :faust2sigviewer),
    FileProduct(joinpath("bin", "faust2smartkeyb"), :faust2smartkeyb),
    FileProduct(joinpath("bin", "faust2sndfile"), :faust2sndfile),
    FileProduct(joinpath("bin", "faust2soul"), :faust2soul),
    FileProduct(joinpath("bin", "faust2supercollider"), :faust2supercollider),
    FileProduct(joinpath("bin", "faust2svg"), :faust2svg),
    FileProduct(joinpath("bin", "faust2svgviewer"), :faust2svgviewer),
    FileProduct(joinpath("bin", "faust2teensy"), :faust2teensy),
    FileProduct(joinpath("bin", "faust2unity"), :faust2unity),
    FileProduct(joinpath("bin", "faust2vcvrack"), :faust2vcvrack),
    FileProduct(joinpath("bin", "faust2vst"), :faust2vst),
    FileProduct(joinpath("bin", "faust2vsti"), :faust2vsti),
    FileProduct(joinpath("bin", "faust2w32max6"), :faust2w32max6),
    FileProduct(joinpath("bin", "faust2w32msp"), :faust2w32msp),
    FileProduct(joinpath("bin", "faust2w32puredata"), :faust2w32puredata),
    FileProduct(joinpath("bin", "faust2w32vst"), :faust2w32vst),
    FileProduct(joinpath("bin", "faust2w64max6"), :faust2w64max6),
    FileProduct(joinpath("bin", "faust2w64vst"), :faust2w64vst),
    FileProduct(joinpath("bin", "faust2wasm"), :faust2wasm),
    FileProduct(joinpath("bin", "faust2webaudiowasm"), :faust2webaudiowasm),
    FileProduct(joinpath("bin", "faust2webaudiowast"), :faust2webaudiowast),
    FileProduct(joinpath("bin", "faust2winunity"), :faust2winunity),

    FileProduct(joinpath("bin","faustoptflags"), :faustoptflags),
    FileProduct(joinpath("bin","usage.sh"), :usage_sh),

    FileProduct(joinpath("share", "faust", "stdfaust.lib"), :stdfaust_lib),
]

llvm_versions = [v"11.0.1", v"12.0.1", v"13.0.1"]

sources = Dict(
    v"11.0.1" => [sources; ],
    v"12.0.1" => [sources; ],
    v"13.0.1" => [sources; ],
)

base_dependencies = [
    Dependency("libmicrohttpd_jll"),
    Dependency("libsndfile_jll"),
    BuildDependency("Ncurses_jll"),
    BuildDependency("XML2_jll"),
    BuildDependency("Zlib_jll"),
]

augment_platform_block = """
    using Base.BinaryPlatforms
    $(LLVM.augment)
    function augment_platform!(platform::Platform)
        augment_llvm!(platform)
    end
"""

# determine exactly which tarballs we should build
builds = []
for llvm_version in llvm_versions
    # Dependencies that must be installed before this package can be built
    llvm_name = "LLVM_full_jll"
    dependencies = [
	base_dependencies;
        BuildDependency(PackageSpec(name=llvm_name, version=llvm_version))
    ]

    for platform in platforms
        augmented_platform = deepcopy(platform)
        augmented_platform[LLVM.platform_name] = LLVM.platform(llvm_version, false)

        should_build_platform(triplet(augmented_platform)) || continue
        push!(builds, (;
            dependencies,
            sources=sources[llvm_version],
            platforms=[augmented_platform],
        ))
    end
end

# don't allow `build_tarballs` to override platform selection based on ARGS.
# we handle that ourselves by calling `should_build_platform`
non_platform_ARGS = filter(arg -> startswith(arg, "--"), ARGS)

# `--register` should only be passed to the latest `build_tarballs` invocation
non_reg_ARGS = filter(arg -> arg != "--register", non_platform_ARGS)

for (i,build) in enumerate(builds)
    build_tarballs(i == lastindex(builds) ? non_platform_ARGS : non_reg_ARGS,
                   name, version, build.sources, script,
                   build.platforms, products, build.dependencies;
                   preferred_gcc_version=v"8", julia_compat="1.6",
                   augment_platform_block, lazy_artifacts=true)
end


text:
@testset "Atom Type" begin
    atom = Atom("He")

    @test mass(atom) ≈ 4.002602 atol=1e-6
    @test charge(atom) == 0
    @test name(atom) == "He"
    @test atom_type(atom) == "He"

    set_mass!(atom, 678)
    @test mass(atom) == 678
    set_charge!(atom, -1.5)
    @test charge(atom) == -1.5
    set_name!(atom, "Zn")
    @test name(atom) == "Zn"

    @test fullname(atom) == "Helium"
    set_atom_type!(atom, "Zn")
    @test atom_type(atom) == "Zn"
    @test fullname(atom) == "Zinc"

    @test vdw_radius(atom) ≈ 2.1 atol=1e-1
    @test covalent_radius(atom) ≈ 1.31 atol=1e-2
    @test atomic_number(atom) == 30

    copy = deepcopy(atom)
    @test name(copy) == "Zn"

    set_name!(copy, "I")
    @test name(copy) == "I"
    @test name(atom) == "Zn"
end


text:
using Flux
import Flux: trainable
using Zygote: hook
using Statistics: mean, std

"""
    Hook(layer, [forward, backward, state = Dict()])

Wraps a Flux layer, and holds a Dict with state (`Hook.state`)

After the forward pass of the wrapped layer calls `Hook.forward(state, layeroutput)`
After the gradient calculation of the wrapped layer calls `Hook.backward(state, layergradient)`

For adding hooks to existing models, see [`addhook`](@ref)
"""
struct Hook
    layer
    forward
    backward
    state
end

Hook(layer; forward = const_, backward = const_, state = Dict()) = Hook(
        layer, forward, backward, state)
const_(_, x) = x

(hook::Hook)(x) = hook.forward(hook.state, Zygote.hook(x̄ -> hook.backward(hook.state, x̄), hook.layer(x)))
trainable(hook::Hook) = (hook.layer,)

"""
    addhook(model; [forward, backward, filter_fn])

Add a [`Hook`](@ref) with callbacks `forward` and `backward` to all layers
of `model` where `filter_fn(layer) == true`.

See [`Hook`](@ref)
"""
function addhook(model; forward = const_, backward = const_, filter_fn = _ -> true)
    mapmodel(model) do layer
        filter_fn(layer) ? Hook(layer, forward = forward, backward = backward) : layer
    end
end

function removehooks(model)
    mapmodel(l -> l isa Hook ? l.layer : l, model)
end

# hook functions

function forwardstatshook(state::Dict, output)
    state["activation_mean"] = mean(output)
    state["activation_std"] = std(output)
    return output
end

function hooksaveactivations(state::Dict, output)
    state["activation"] = output
    return output
end

function hooksavegrads(state::Dict, gradients)
    state["gradient"] = gradients
    return gradients
end

function backwardstatshook(state::Dict, gradients)
    state["gradient_mean"] = mean(gradients)
    state["gradient_std"] = std(gradients)
    return gradients
end


text:
#
# OptimPack.jl --
#
# Optimization for Julia.
#
# ----------------------------------------------------------------------------
#
# This file is part of OptimPack.jl which is licensed under the MIT
# "Expat" License:
#
# Copyright (C) 2014-2016, Éric Thiébaut.
#
# ----------------------------------------------------------------------------

module OptimPack

export nlcg, vmlm, spg2

export fzero, fmin, fmin_global

# Functions must be imported to be extended with new methods.
import Base.size
import Base.length
import Base.eltype
import Base.ndims
import Base.copy
import Base.dot

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("OptimPack not properly installed. Please run Pkg.build(\"OptimPack\")")
end

"""
`Float` is any floating point type supported by the library.
"""
typealias Float Union{Cfloat,Cdouble}

cint(i::Integer) = convert(Cint, i)
cuint(i::Integer) = convert(Cuint, i)

const SUCCESS = cint( 0)
const FAILURE = cint(-1)

#------------------------------------------------------------------------------
# ERROR MANAGEMENT
#
# We must provide an appropriate error handler to OptimPack in order to throw
# an error exception and avoid aborting the program in case of misuse of the
# library.

__error__(ptr::Ptr{UInt8}) = (ErrorException(bytestring(ptr)); nothing)

const __cerror__ = cfunction(__error__, Void, (Ptr{UInt8},))

function __init__()
    ccall((:opk_set_error_handler, opklib), Ptr{Void}, (Ptr{Void},),
          __cerror__)
    nothing
end

__init__()

#------------------------------------------------------------------------------
# OBJECT MANAGEMENT

"""
OptimPack Objects
=================

All concrete types derived from the abstract `Object` type have a `handle`
member which stores the address of the OptimPack object.  To avoid conflicts
with Julia `Vector` type, an OptimPack vector (*i.e.* `opk_vector_t`)
corresponds to the type `Variable` in Julia.
"""
abstract Object

"""
Reference Counting
==================

OptimPack use reference counting for managing the memory.  The number of
references of an object is given by `references(obj)`.

Two *private* functions are used for managing the reference count:
`__hold_object__(ptr)` set a reference on an OptimPack object while
`__drop_object__(ptr)` discards a reference on an OptimPack object.  The
argument `ptr` is the address of the object.  This functions shall not be
directly called by a user of the code.
"""
function references(obj::Object)
    ccall((:opk_get_object_references, opklib), Cptrdiff_t, (Ptr{Void},),
          obj.handle)
end

function __hold_object__(ptr::Ptr{Void})
    ccall((:opk_hold_object, opklib), Ptr{Void}, (Ptr{Void},), ptr)
end

function __drop_object__(ptr::Ptr{Void})
    ccall((:opk_drop_object, opklib), Void, (Ptr{Void},), ptr)
end

@doc (@doc references) __hold_object__
@doc (@doc references) __drop_object__

#------------------------------------------------------------------------------
# VARIABLE SPACE

abstract VariableSpace <: Object
"""
Variable Space
==============
Abstract type `VariableSpace` corresponds to a *vector space* (type
`opk_vspace_t`) in OptimPack.
"""

type DenseVariableSpace{T,N} <: VariableSpace
    handle::Ptr{Void}
    eltype::Type{T}
    size::NTuple{N,Int}
    length::Int
end

# Extend basic functions for arrays.
length(vsp::DenseVariableSpace) = vsp.length
eltype(vsp::DenseVariableSpace) = vsp.eltype
size(vsp::DenseVariableSpace) = vsp.size
size(vsp::DenseVariableSpace, n::Integer) = vsp.size[n]
ndims(vsp::DenseVariableSpace) = length(vsp.size)

DenseVariableSpace(T::Union{Type{Cfloat},Type{Cdouble}},
                   dims::Int...) = DenseVariableSpace(T, dims)

function checkdims{N}(dims::NTuple{N,Int})
    number::Int = 1
    for dim in dims
        if dim < 1
            error("invalid dimension")
        end
        number *= dim
    end
    return number
end

for (T, f) in ((Cfloat, :opk_new_simple_float_vector_space),
               (Cdouble, :opk_new_simple_double_vector_space))
    @eval begin
        function DenseVariableSpace{N}(::Type{$T}, dims::NTuple{N,Int})
            length::Int = checkdims(dims)
            ptr = ccall(($(string(f)), opklib), Ptr{Void}, (Cptrdiff_t,), length)
            systemerror("failed to create vector space", ptr == C_NULL)
            obj = DenseVariableSpace{$T,N}(ptr, $T, dims, length)
            finalizer(obj, obj -> __drop_object__(obj.handle))
            return obj
        end
    end
end

#------------------------------------------------------------------------------
# VARIABLES

abstract Variable <: Object
"""
Variables
=========
Abstract type `Variable` correspond to *vectors* (type `opk_vector_t`) in OptimPack.
"""

# Note: There are no needs to register a reference for the owner of a
# vector (it already owns one internally).
type DenseVariable{T<:Float,N} <: Variable
    handle::Ptr{Void}
    owner::DenseVariableSpace{T,N}
    array::Union{Array,Void}
end

length(v::DenseVariable) = length(v.owner)
eltype(v::DenseVariable) = eltype(v.owner)
size(v::DenseVariable) = size(v.owner)
size(v::DenseVariable, n::Integer) = size(v.owner, n)
ndims(v::DenseVariable) = ndims(v.owner)

# FIXME: add means to wrap a Julia array around this or (better?, simpler?)
#        just use allocate a Julia array and wrap a vector around it?
"""
`v = create(s)` creates a new variable of the variable space `s`.
"""
function create{T<:Float,N<:Integer}(space::DenseVariableSpace{T,N})
    ptr = ccall((:opk_vcreate, opklib), Ptr{Void}, (Ptr{Void},), space.handle)
    systemerror("failed to create vector", ptr == C_NULL)
    obj = DenseVariable{T,N}(ptr, space, nothing)
    finalizer(obj, obj -> __drop_object__(obj.handle))
    return obj
end

for (T, wrap, rewrap) in ((Cfloat, :opk_wrap_simple_float_vector,
                           :opk_rewrap_simple_float_vector),
                          (Cdouble, :opk_wrap_simple_double_vector,
                           :opk_rewrap_simple_double_vector))
    @eval begin
        function wrap{N}(s::DenseVariableSpace{$T,N}, a::DenseArray{$T,N})
            assert(size(a) == size(s))
            ptr = ccall(($(string(wrap)), opklib), Ptr{Void},
                        (Ptr{Void}, Ptr{$T}, Ptr{Void}, Ptr{Void}),
                        s.handle, a, C_NULL, C_NULL)
            systemerror("failed to wrap vector", ptr == C_NULL)
            obj = DenseVariable{$T,N}(ptr, s, a)
            finalizer(obj, obj -> __drop_object__(obj.handle))
            return obj
        end

        function wrap!{N}(v::DenseVariable{$T,N}, a::DenseArray{$T,N})
            assert(size(a) == size(v))
            assert(v.array != nothing)
            status = ccall(($(string(rewrap)), opklib), Cint,
                           (Ptr{Void}, Ptr{$T}, Ptr{Void}, Ptr{Void}),
                           v.handle, a, C_NULL, C_NULL)
            systemerror("failed to re-wrap vector", status != SUCCESS)
            v.array = a
            return v
        end
    end
end

"""
`v = wrap(s, a)` wraps the Julia array `a` into a variable of the space `s` and
returns the resulting variable `v`.  Array `a` must have the correct dimensions
and element type.
""" wrap

"""
`wrap!(v, a)` rewraps the Julia array `a` into the variable `v` and returns
`v`.  Array `a` must have the correct dimensions and element type.
""" wrap

#------------------------------------------------------------------------------
# OPERATIONS ON VARIABLES

"""
`norm1(v)` returns the L1 norm (sum of absolute values) ov *variables* `v`.
"""
function norm1(v::Variable)
    ccall((:opk_vnorm1, opklib), Cdouble, (Ptr{Void},), v.handle)
end

"""
`norm2(v)` returns the Euclidean (L2) norm (square root of the sum of squared
values) of *variables* `v`.
"""
function norm2(v::Variable)
    ccall((:opk_vnorm2, opklib), Cdouble, (Ptr{Void},), v.handle)
end

"""
`norminf(v)` returns the infinite norm (maximum absolute value) of *variables*
`v`.
"""
function norminf(v::Variable)
    ccall((:opk_vnorminf, opklib), Cdouble, (Ptr{Void},), v.handle)
end

"""
`zero!(v)` fills *variables* `v` with zeros.
"""
function zero!(v::Variable)
    ccall((:opk_vzero, opklib), Void, (Ptr{Void},), v.handle)
end

"""
`fill!(v, alpha)` fills *variables* `v` with value `alpha`.
"""
function fill!(v::Variable, alpha::Real)
    ccall((:opk_vfill, opklib), Void, (Ptr{Void},Cdouble),
          v.handle, alpha)
end

"""
`copy!(dst, src)` copies source *variables* `src` into the destination
*variables* `dst`.
"""
function copy!(dst::Variable, src::Variable)
    ccall((:opk_vcopy, opklib), Void, (Ptr{Void},Ptr{Void}),
          dst.handle, src.handle)
end

"""
`scale!(dst, alpha, src)` stores `alpha` times the source *variables* `src`
into the destination *variables* `dst`.
"""
function scale!(dst::Variable, alpha::Real, src::Variable)
    ccall((:opk_vscale, opklib), Void, (Ptr{Void},Cdouble,Ptr{Void}),
          dst.handle, alpha, src.handle)
end

"""
`swap!(x, y)` exchanges the contents of *variables* `x` and `y`.
"""
function swap!(x::Variable, y::Variable)
    ccall((:opk_vswap, opklib), Void, (Ptr{Void},Ptr{Void}),
          x.handle, y.handle)
end

"""
`dot(x, y)` returns the inner product of *variables* `x` and `y`.
"""
function dot(x::Variable, y::Variable)
    ccall((:opk_vdot, opklib), Cdouble, (Ptr{Void},Ptr{Void}),
          x.handle, y.handle)
end

"""
`combine!(dst, alpha, x, beta, y)` stores into the destination `dst` the linear
combination `alpha*x + beta*y`.

`combine!(dst, alpha, x, beta, y, gamma, z)` stores into the destination `dst`
the linear combination `alpha*x + beta*y + gamma*z`.
"""
function combine!(dst::Variable,
                  alpha::Real, x::Variable,
                  beta::Real,  y::Variable)
    ccall((:opk_vaxpby, opklib), Void,
          (Ptr{Void},Cdouble,Ptr{Void},Cdouble,Ptr{Void}),
          dst.handle, alpha, x.handle, beta, y.handle)
end

function combine!(dst::Variable,
                  alpha::Real, x::Variable,
                  beta::Real,  y::Variable,
                  gamma::Real, z::Variable)
    ccall((:opk_vaxpbypcz, opklib), Void,
          (Ptr{Void},Cdouble,Ptr{Void},Cdouble,Ptr{Void},Cdouble,Ptr{Void}),
          dst.handle, alpha, x.handle, beta, y.handle, gamma, y.handle)
end

#------------------------------------------------------------------------------
# OPERATORS

abstract Operator <: Object

for s in (:apply_direct, :apply_adoint, :apply_inverse)
    @eval begin
        function $(Symbol(string(s)*"!"))(op::Operator,
                                          dst::Variable,
                                          src::Variable)
            status = ccall(($("opk_"*string(s)), opklib), Cint,
                           (Ptr{Void},Ptr{Void},Ptr{Void}),
                           op.handle, dst.handle, src.handle)
            if status != SUCCESS
                error("something wrong happens")
            end
            nothing
        end
    end
end

#------------------------------------------------------------------------------
# LINE SEARCH METHODS

abstract LineSearch <: Object

type ArmijoLineSearch <: LineSearch
    handle::Ptr{Void}
    ftol::Cdouble
    function ArmijoLineSearch(;ftol::Real=1e-4)
        assert(0.0 <= ftol < 1.0)
        ptr = ccall((:opk_lnsrch_new_backtrack, opklib), Ptr{Void},
                    (Cdouble,), ftol)
        systemerror("failed to create linesearch", ptr == C_NULL)
        obj = new(ptr, ftol)
        finalizer(obj, obj -> __drop_object__(obj.handle))
        return obj
    end
end

type MoreThuenteLineSearch <: LineSearch
    handle::Ptr{Void}
    ftol::Cdouble
    gtol::Cdouble
    xtol::Cdouble
    function MoreThuenteLineSearch(;ftol::Real=1e-4, gtol::Real=0.9,
                                   xtol::Real=eps(Cdouble))
        assert(0.0 <= ftol < gtol < 1.0)
        assert(0.0 <= xtol < 1.0)
        ptr = ccall((:opk_lnsrch_new_csrch, opklib), Ptr{Void},
                (Cdouble, Cdouble, Cdouble), ftol, gtol, xtol)
        systemerror("failed to create linesearch", ptr == C_NULL)
        obj = new(ptr, ftol, gtol, xtol)
        finalizer(obj, obj -> __drop_object__(obj.handle))
        return obj
    end
end

type NonmonotoneLineSearch <: LineSearch
    handle::Ptr{Void}
    mem::Int
    ftol::Cdouble
    amin::Cdouble
    amax::Cdouble
    function NonmonotoneLineSearch(;mem::Integer=10, ftol::Real=1e-4,
                                   amin::Real=0.1, amax::Real=0.9)
        assert(mem >= 1)
        assert(0.0 <= ftol < 1.0)
        assert(0.0 < amin < amax < 1.0)
        ptr = ccall((:opk_lnsrch_new_nonmonotone, opklib), Ptr{Void},
                (Cptrdiff_t, Cdouble, Cdouble, Cdouble), mem, ftol, amin, amax)
        systemerror("failed to create nonmonotone linesearch", ptr == C_NULL)
        obj = new(ptr, mem, ftol, amin, amax)
        finalizer(obj, obj -> __drop_object__(obj.handle))
        return obj
    end
end

const LNSRCH_ERROR_ILLEGAL_ADDRESS                    = cint(-12)
const LNSRCH_ERROR_CORRUPTED_WORKSPACE                = cint(-11)
const LNSRCH_ERROR_BAD_WORKSPACE                      = cint(-10)
const LNSRCH_ERROR_STP_CHANGED                        = cint( -9)
const LNSRCH_ERROR_STP_OUTSIDE_BRACKET                = cint( -8)
const LNSRCH_ERROR_NOT_A_DESCENT                      = cint( -7)
const LNSRCH_ERROR_STPMIN_GT_STPMAX                   = cint( -6)
const LNSRCH_ERROR_STPMIN_LT_ZERO                     = cint( -5)
const LNSRCH_ERROR_STP_LT_STPMIN                      = cint( -4)
const LNSRCH_ERROR_STP_GT_STPMAX                      = cint( -3)
const LNSRCH_ERROR_INITIAL_DERIVATIVE_GE_ZERO         = cint( -2)
const LNSRCH_ERROR_NOT_STARTED                        = cint( -1)
const LNSRCH_SEARCH                                   = cint(  0)
const LNSRCH_CONVERGENCE                              = cint(  1)
const LNSRCH_WARNING_ROUNDING_ERRORS_PREVENT_PROGRESS = cint(  2)
const LNSRCH_WARNING_XTOL_TEST_SATISFIED              = cint(  3)
const LNSRCH_WARNING_STP_EQ_STPMAX                    = cint(  4)
const LNSRCH_WARNING_STP_EQ_STPMIN                    = cint(  5)

function start!(ls::LineSearch, f0::Real, df0::Real,
                stp1::Real, stpmin::Real, stpmax::Real)
    ccall((:opk_lnsrch_start, opklib), Cint,
          (Ptr{Void}, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble),
          ls, f0, df0, stp1, stpmin, stpmax)
end

function iterate!(ls::LineSearch, stp::Real, f::Real, df::Real)
    _stp = Cdouble[stp]
    task = ccall((:opk_lnsrch_iterate, opklib), Cint,
                 (Ptr{Void}, Ptr{Cdouble}, Cdouble, Cdouble),
                 ls, _stp, f, df)
    return (task, _stp[1])
end

get_step(ls::LineSearch) = ccall((:opk_lnsrch_get_step, opklib),
                                 Cdouble, (Ptr{Void}, ), ls)
get_status(ls::LineSearch) = ccall((:opk_lnsrch_get_status, opklib),
                                   Cint, (Ptr{Void}, ), ls)
has_errors(ls::LineSearch) = (ccall((:opk_lnsrch_has_errors, opklib),
                                    Cint, (Ptr{Void}, ), ls) != 0)
has_warnings(ls::LineSearch) = (ccall((:opk_lnsrch_has_warnings, opklib),
                                      Cint, (Ptr{Void}, ), ls) != 0)
converged(ls::LineSearch) = (ccall((:opk_lnsrch_converged, opklib),
                                   Cint, (Ptr{Void}, ), ls) != 0)
finished(ls::LineSearch) = (ccall((:opk_lnsrch_finished, opklib),
                                  Cint, (Ptr{Void}, ), ls) != 0)
get_ftol(ls::LineSearch) = ls.ftol
get_gtol(ls::MoreThuenteLineSearch) = ls.gtol
get_xtol(ls::MoreThuenteLineSearch) = ls.xtol

#------------------------------------------------------------------------------
# NON LINEAR OPTIMIZERS

# Codes returned by the reverse communication version of optimzation
# algorithms.

const TASK_ERROR       = cint(-1) # An error has ocurred.
const TASK_PROJECT_X   = cint( 0) # Caller must project variables x.
const TASK_COMPUTE_FG  = cint( 1) # Caller must compute f(x) and g(x).
const TASK_PROJECT_D   = cint( 2) # Caller must project the direction d.
const TASK_FREE_VARS   = cint( 3) # Caller must update the subspace of free variables.
const TASK_NEW_X       = cint( 4) # A new iterate is available.
const TASK_FINAL_X     = cint( 5) # Algorithm has converged, solution is available.
const TASK_WARNING     = cint( 6) # Algorithm terminated with a warning.

abstract Optimizer <: Object
abstract Options

const NLCG_FLETCHER_REEVES        = cuint(1)
const NLCG_HESTENES_STIEFEL       = cuint(2)
const NLCG_POLAK_RIBIERE_POLYAK   = cuint(3)
const NLCG_FLETCHER               = cuint(4)
const NLCG_LIU_STOREY             = cuint(5)
const NLCG_DAI_YUAN               = cuint(6)
const NLCG_PERRY_SHANNO           = cuint(7)
const NLCG_HAGER_ZHANG            = cuint(8)
const NLCG_POWELL                 = cuint(1<<8) # force beta >= 0
const NLCG_SHANNO_PHUA            = cuint(1<<9) # compute scale from previous iterations

# For instance: (NLCG_POLAK_RIBIERE_POLYAK | NLCG_POWELL) merely
# corresponds to PRP+ (Polak, Ribiere, Polyak) while (NLCG_PERRY_SHANNO |
# NLCG_SHANNO_PHUA) merely corresponds to the conjugate gradient method
# implemented in CONMIN.

type NLCGoptions <: Options
    gatol::Cdouble
    grtol::Cdouble
    stpmin::Cdouble
    stpmax::Cdouble
    method::UInt
    function NLCGoptions(;gatol::Real=0.0, grtol::Real=1E-6,
                         stpmin::Real=1E-20, stpmax::Real=1E+20,
                         method::Integer=(NLCG_HAGER_ZHANG|NLCG_SHANNO_PHUA))
        new(convert(Cdouble, gatol),
            convert(Cdouble, grtol),
            convert(Cdouble, stpmin),
            convert(Cdouble, stpmax),
            convert(Cptrdiff_t, method))
    end
    NLCGoptions(method::Integer) = new(convert(UInt, method))
end

"""
Default settings for non linear conjugate gradient (should correspond to the
method which is, in general, the most successful).
"""
const NLCG_DEFAULT = NLCGoptions()

type NLCG <: Optimizer
    handle::Ptr{Void}
    space::VariableSpace
    method::Cuint
    lnsrch::LineSearch
    function NLCG(options::NLCGoptions,
                  space::VariableSpace,
                  lnsrch::LineSearch)
        ptr = ccall((:opk_new_nlcg_optimizer_with_line_search, opklib),
                    Ptr{Void}, (Ptr{Void}, Cuint, Ptr{Void}),
                    space.handle, options.method, lnsrch.handle)
        systemerror("failed to create optimizer", ptr == C_NULL)
        obj = new(ptr, space, options.method, lnsrch)
        finalizer(obj, obj -> __drop_object__(obj.handle))
        set_gatol!(obj, options.gatol)
        set_grtol!(obj, options.grtol)
        set_stpmin_and_stpmax!(obj, options.stpmin, options.stpmax)
        return obj
    end
end

const SCALING_NONE             = cint(0)
const SCALING_OREN_SPEDICATO   = cint(1) # gamma = <s,y>/<y,y>
const SCALING_BARZILAI_BORWEIN = cint(2) # gamma = <s,s>/<s,y>

type VMLMoptions <: Options
    gatol::Cdouble
    grtol::Cdouble
    stpmin::Cdouble
    stpmax::Cdouble
    mem::Cptrdiff_t
    scaling::Cint
    function VMLMoptions(;gatol::Real=0.0, grtol::Real=1E-6,
                         stpmin::Real=1E-20, stpmax::Real=1E+20,
                         mem::Integer=3,
                         scaling::Integer=SCALING_OREN_SPEDICATO)
        new(convert(Cdouble, gatol),
            convert(Cdouble, grtol),
            convert(Cdouble, stpmin),
            convert(Cdouble, stpmax),
            convert(Cptrdiff_t, mem),
            convert(Cint, scaling))
    end
end

"""
Default settings for non linear conjugate gradient (should correspond to the
method which is, in general, the most successful).
"""
const VMLM_DEFAULT = VMLMoptions()

type VMLM <: Optimizer
    handle::Ptr{Void}
    space::VariableSpace
    mem::Cptrdiff_t
    lnsrch::LineSearch
    function VMLM(options::VMLMoptions,
                  space::VariableSpace,
                  lnsrch::LineSearch)
        mem = options.mem
        mem ≥ 1 || error("illegal number of memorized steps")
        mem = min(mem, length(space))
        ptr = ccall((:opk_new_vmlm_optimizer_with_line_search, opklib),
                    Ptr{Void}, (Ptr{Void}, Cptrdiff_t, Ptr{Void},
                                Cdouble, Cdouble, Cdouble),
                    space.handle, mem, lnsrch.handle, 0.0, 0.0, 0.0)
        systemerror("failed to create optimizer", ptr == C_NULL)
        obj = new(ptr, space, mem, lnsrch)
        finalizer(obj, obj -> __drop_object__(obj.handle))
        set_gatol!(obj, options.gatol)
        set_grtol!(obj, options.grtol)
        set_stpmin_and_stpmax!(obj, options.stpmin, options.stpmax)
        set_scaling!(obj, options.scaling)
        return obj
    end
end

for (T, start, iterate,
     get_task,
     get_iterations,
     get_evaluations,
     get_restarts,
     get_gatol, get_grtol,
     set_gatol, set_grtol,
     get_stpmin, get_stpmax,
     set_stpmin_and_stpmax) in ((NLCG, :opk_start_nlcg, :opk_iterate_nlcg,
                                 :opk_get_nlcg_task, :opk_get_nlcg_iterations,
                                 :opk_get_nlcg_evaluations,
                                 :opk_get_nlcg_restarts,
                                 :opk_get_nlcg_gatol, :opk_get_nlcg_grtol,
                                 :opk_set_nlcg_gatol, :opk_set_nlcg_grtol,
                                 :opk_get_nlcg_stpmin, :opk_get_nlcg_stpmax,
                                 :opk_set_nlcg_stpmin_and_stpmax),
                                (VMLM, :opk_start_vmlm, :opk_iterate_vmlm,
                                 :opk_get_vmlm_task, :opk_get_vmlm_iterations,
                                 :opk_get_vmlm_evaluations,
                                 :opk_get_vmlm_restarts,
                                 :opk_get_vmlm_gatol, :opk_get_vmlm_grtol,
                                 :opk_set_vmlm_gatol, :opk_set_vmlm_grtol,
                                 :opk_get_vmlm_stpmin, :opk_get_vmlm_stpmax,
                                 :opk_set_vmlm_stpmin_and_stpmax))
    @eval begin

        start!(opt::$T) = ccall(($(string(start)), opklib), Cint,
                                  (Ptr{Void},), opt.handle)

        function iterate!(opt::$T, x::Variable, f::Real, g::Variable)
            ccall(($(string(iterate)), opklib), Cint,
                  (Ptr{Void}, Ptr{Void}, Cdouble, Ptr{Void}),
                  opt.handle, x.handle, f, g.handle)
        end

        get_task(opt::$T) = ccall(($(string(get_task)), opklib),
                                    Cint, (Ptr{Void},), opt.handle)

        iterations(opt::$T) = ccall(($(string(get_iterations)), opklib),
                                    Cptrdiff_t, (Ptr{Void},), opt.handle)

        evaluations(opt::$T) = ccall(($(string(get_evaluations)), opklib),
                                     Cptrdiff_t, (Ptr{Void},), opt.handle)

        restarts(opt::$T) = ccall(($(string(get_restarts)), opklib),
                                  Cptrdiff_t, (Ptr{Void},), opt.handle)

        get_gatol(opt::$T) = ccall(($(string(get_gatol)), opklib),
                                   Cdouble, (Ptr{Void},), opt.handle)

        get_grtol(opt::$T) = ccall(($(string(get_grtol)), opklib),
                                   Cdouble, (Ptr{Void},), opt.handle)

        function set_gatol!(opt::$T, gatol::Real)
            if ccall(($(string(set_gatol)), opklib),
                     Cint, (Ptr{Void},Cdouble), opt.handle, gatol) != SUCCESS
                e = errno()
                if e == Base.EINVAL
                    error("invalid value for parameter gatol")
                else
                    error("unexpected error while setting parameter gatol")
                end
            end
        end

        function set_grtol!(opt::$T, grtol::Real)
            if ccall(($(string(set_grtol)), opklib),
                     Cint, (Ptr{Void},Cdouble), opt.handle, grtol) != SUCCESS
                e = errno()
                if e == Base.EINVAL
                    error("invalid value for parameter grtol")
                else
                    error("unexpected error while setting parameter grtol")
                end
            end
        end

        get_stpmin(opt::$T) = ccall(($(string(get_stpmin)), opklib),
                                    Cdouble, (Ptr{Void},), opt.handle)

        get_stpmax(opt::$T) = ccall(($(string(get_stpmax)), opklib),
                                    Cdouble, (Ptr{Void},), opt.handle)

        function set_stpmin_and_stpmax!(opt::$T, stpmin::Real, stpmax::Real)
            if ccall(($(string(set_stpmin_and_stpmax)), opklib),
                     Cint, (Ptr{Void},Cdouble,Cdouble),
                     opt.handle, stpmin, stpmax) != SUCCESS
                e = errno()
                if e == Base.EINVAL
                    error("invalid values for parameters stpmin and stpmax")
                else
                    error("unexpected error while setting parameters stpmin and stpmax")
                end
            end
        end
    end
end

"""
`task = start!(opt)` starts optimization with the nonlinear optimizer
`opt` and returns the next pending task.
""" start!

"""
`task = iterate!(opt, x, f, g)` performs one optimization step with
the nonlinear optimizer `opt` for variables `x`, function value `f`
and gradient `g`.  The method returns the next pending task.
""" iterate!

"""
`get_task(opt)` returns the current pending task for the nonlinear
optimizer `opt`.
""" get_task

"""
`iterations(opt)` returns the number of iterations performed by the
nonlinear optimizer `opt`.
""" iterations

"""
`evaluations(opt)` returns the number of function (and gradient)
evaluations requested by the nonlinear optimizer `opt`.
""" evaluations

"""
`restarts(opt)` returns the number of restarts performed by the
nonlinear optimizer `opt`.
""" restarts

"""
`get_gatol(opt)` returns the absolute gradient threshold used by the
nonlinear optimizer `opt` to check for the convergence.
""" get_gatol

"""
`get_grtol(opt)` returns the relative gradient threshold used by the
nonlinear optimizer `opt` to check for the convergence.
""" get_grtol

"""
`set_gatol!(opt, gatol)` set the absolute gradient threshold used by the
nonlinear optimizer `opt` to check for the convergence.
""" set_gatol!

"""
`set_grtol!(opt, grtol)` set the relative gradient threshold used by the
nonlinear optimizer `opt` to check for the convergence.
""" set_grtol!

"""
`get_stpmin(opt)` returns the minimum relative step length used by the
nonlinear optimizer `opt` during line searches.
""" get_stpmin

"""
`get_stpmax(opt)` returns the maximum relative step length used by the
nonlinear optimizer `opt` during line searches.
""" get_stpmax

"""
`set_stpmin_and_stpmax!(opt, stpmin, stpmax)` set the minimum and the
maximum relative step length used by the nonlinear optimizer `opt` during
line searches.
""" set_stpmin_and_stpmax!

get_scaling(opt::VMLM) = ccall((:opk_get_vmlm_scaling, opklib),
                                      Cint, (Ptr{Void},), opt.handle)

function set_scaling!(opt::VMLM, scaling::Integer)
    if ccall((:opk_set_vmlm_scaling, opklib),
             Cint, (Ptr{Void},Cint), opt.handle, scaling) != SUCCESS
        error("unexpected error while setting scaling scaling")
    end
end

#------------------------------------------------------------------------------
# DRIVERS FOR NON-LINEAR OPTIMIZATION

default_nlcg_line_search() = MoreThuenteLineSearch(ftol=1E-4, gtol=0.1)
default_vmlm_line_search() = MoreThuenteLineSearch(ftol=1E-4, gtol=0.9)

"""
Nonlinear Conjugate Gradient
============================

Minimizing the smooth mulivariate function `f(x)` by a nonlinear conjugate
gradient methods is done by:
```
x = nlcg(fg!, x0, method)
```
where `fg!` implements the objective function (and its gradient), `x0` gives
the initial value of the variables (as well as the data type and dimensions of
the solution) and optional argument `method` may be used to choose a specific
conjugate gradient method.

See `vmlm` for more details.
"""
function nlcg{T<:Float,N}(fg!::Function, x0::DenseArray{T,N},
                          method::Integer=NLCG_DEFAULT.method;
                          lnsrch::LineSearch=default_nlcg_line_search(),
                          gatol::Real=NLCG_DEFAULT.gatol,
                          grtol::Real=NLCG_DEFAULT.grtol,
                          stpmin::Real=NLCG_DEFAULT.stpmin,
                          stpmax::Real=NLCG_DEFAULT.stpmax,
                          maxeval::Integer=-1, maxiter::Integer=-1,
                          verb::Bool=false, debug::Bool=false)
    # Create an optimizer and solve the problem.
    space = DenseVariableSpace(T, size(x0))
    options = NLCGoptions(gatol=gatol, grtol=grtol, stpmin=stpmin,
                          stpmax=stpmax, method=method)
    opt = NLCG(options, space, lnsrch)
    solve(opt, fg!, x0, maxeval=maxeval, maxiter=maxiter,
          verb=verb, debug=debug)
end

"""
Limited Memory Variable Metric
==============================

Minimizing the smooth mulivariate function `f(x)` by a limited-memory version
of the LBFGS variable metric method is done by:
```
x = vmlm(fg!, x0, mem)
```
where `fg!` implements the objective function (see below), `x0` gives
the initial value of the variables (as well as the data type and dimensions of
the solution) and optional argument `mem` is the number of previous steps to
memorize (by default `mem = 3`).

The objective function is implemented by `fg!` which is called as:
```
f = fg!(x, g)
```
with `x` the current variables and `g` a Julia array (of same type and
simensions as `x`) to store the gradient of the function.  The value returned
by `fg!` is `f(x)`.
"""
function vmlm{T<:Float,N}(fg!::Function, x0::DenseArray{T,N},
                          mem::Integer=VMLM_DEFAULT.mem;
                          lnsrch::LineSearch=default_vmlm_line_search(),
                          scaling::Integer=VMLM_DEFAULT.scaling,
                          gatol::Real=VMLM_DEFAULT.gatol,
                          grtol::Real=VMLM_DEFAULT.grtol,
                          stpmin::Real=VMLM_DEFAULT.stpmin,
                          stpmax::Real=VMLM_DEFAULT.stpmax,
                          maxeval::Integer=-1, maxiter::Integer=-1,
                          verb::Bool=false, debug::Bool=false)
    # Create an optimizer and solve the problem.
    options = VMLMoptions(gatol=gatol, grtol=grtol, stpmin=stpmin,
                          stpmax=stpmax, mem=mem, scaling=scaling)
    space = DenseVariableSpace(T, size(x0))
    opt = VMLM(options, space, lnsrch)
    solve(opt, fg!, x0, maxeval=maxeval, maxiter=maxiter,
          verb=verb, debug=debug)
end

function solve(opt::Optimizer, fg!::Function, x0::DenseArray;
               maxeval::Integer=-1, maxiter::Integer=-1,
               verb::Bool=false, debug::Bool=false)
    if debug
        @printf("gatol=%E; grtol=%E; stpmin=%E; stpmax=%E\n",
                get_gatol(opt), get_grtol(opt),
                get_stpmin(opt), get_stpmax(opt))
    end
    dims = size(x0)
    space = opt.space
    x = copy(x0)
    g = similar(x)
    wx = wrap(space, x)
    wg = wrap(space, g)
    task = start!(opt)
    while true
        if task == TASK_COMPUTE_FG
            f = fg!(x, g)
        elseif task == TASK_NEW_X || task == TASK_FINAL_X
            iter = iterations(opt)
            eval = evaluations(opt)
            if verb
                if iter == 0
                    @printf("%s\n%s\n",
                            " ITER   EVAL  RESTARTS         F(X)             ||G(X)||",
                            "--------------------------------------------------------")
                end
                @printf("%5d  %5d  %5d  %24.16E %10.3E\n",
                        iter, eval, restarts(opt), f, norm2(wg))
            end
            if task == TASK_FINAL_X
                return x
            end
            if maxiter >= 0 && iter >= maxiter
                warn("exceeding maximum number of iterations ($maxiter)")
                return x
            end
            if maxeval >= 0 && eval >= maxeval
                warn("exceeding maximum number of evaluations ($eval >= $maxeval)")
                return x
            end
        elseif task == TASK_WARNING
            @printf("some warnings...\n")
            return x
        elseif task == TASK_ERROR
            @printf("some errors...\n")
            return nothing
        else
            @printf("unexpected task...\n")
            return nothing
        end
        task = iterate!(opt, wx, f, wg)
    end
end

# Load other components.
include("Brent.jl")
include("Powell.jl")
include("spg2.jl")

end # module


text:
import Observables
import AbstractPlotting

scene = AbstractPlotting.Scene()
len = 1 + rand()
width = 3
color = :gold
transparency = false
tail, head = ℝ³(rand(3)), ℝ³(rand(3))
arrow = Arrow(tail,
              head,
              scene,
              width = width,
              color = color)
tail1 = Observables.to_value(arrow.tail)[1]
head1 = Observables.to_value(arrow.head)[1]
tail, head = ℝ³(rand(3)), ℝ³(rand(3))
update(arrow, tail, head)
tail2 = Observables.to_value(arrow.tail)[1]
head2 = Observables.to_value(arrow.head)[1]

@test isapprox(tail1, tail2) == false
@test isapprox(head1, head2) == false


text:
# export ckepler, ckepler!, pkepler, pkepler!
export vdp, vdp!, basic!, basic

m = 1;
k = 1.0168951928e3;
kh = k/m;

function ckepler!(dx, x, p, t)
    dx[1] = x[2]
    dx[2] = -kh/(norm([x[1], x[3]])^3) * x[1]
    dx[3] = x[4]
    dx[4] = -kh/(norm([x[1], x[3]])^3) * x[3]
end
function ckepler(x, p, t)
    dx = similar(x)
    ckepler!(dx, x, p, t)
    return dx
end

function pkepler!(dx, x, p, t)
    dx[1] = x[2]
    dx[2] = x[1]*x[4]^2 - kh/(x[1]^2)
    dx[3] = x[4]
    dx[4] = -2*x[2]*x[4]/x[1]
end
function pkepler(x, p, t)
    dx = similar(x)
    pkepler!(dx, x, p, t)
    return dx
end

function vdp!(dx, x, p, t)
    dx[1] = x[2]
    dx[2] = (1-x[1]^2)*x[2] - x[1]
end
function vdp(x, p, t)
    dx = similar(x)
    vdp!(dx, x, p, t)
    return dx
end

function basic!(dx, x, p, t)
    A = [0. -0.95; 0.95 0.]
    dx[:] = A * x
end

function basic(x, p, t)
    dx = similar(x)
    basic!(dx, x, p, t)
    return dx
end


text:
"""
    struct Ellipsoid{T,CO,TR,TP,TT} <: AbstractVolumePrimitive{T, CO}

Volume primitive describing an [Ellipsoid](@ref).

## Parametric types
* `T`: Precision type.
* `CO`: Describes whether the surface belongs to the primitive.
    It can be `ClosedPrimitive`, i.e. the surface points belong to the primitive,
    or `OpenPrimitive`, i.e. the surface points do not belong to the primitive.
* `TR`: Type of the radius `r`.
    * `TR == T`: Sphere (constant radius `r` along all axes).
* `TP`: Type of the azimuthial angle `φ`.
    * `TP == Nothing`: Full 2π in `φ`.
* `TT`: Type of the polar angle `θ`.
    * `TT == Nothing`: Full 2π in `θ`.

## Fields
* `r::TR`: Definition of the radius of the `Ellipsoid` (in m).
* `φ::TP`: Range in azimuthial angle `φ` of the `Ellipsoid`.
* `θ::TT`: Range in polar angle `θ` of the `Ellipsoid`.
* `origin::CartesianPoint{T}`: The position of the center of the `Ellipsoid`.
* `rotation::SMatrix{3,3,T,9}`: Matrix that describes a rotation of the `Ellipsoid` around its `origin`.

## Definition in Configuration File

So far, the only `Ellipsoid` implemented so far is a `FullSphere`.
A `FullSphere` is defined in the configuration file as part of the `geometry` field
of an object through the field `sphere`.

Example definitions of a `FullSphere` looks like this:
```yaml
sphere:
  r: 2
```
This is a full sphere with radius 2.

To define a sphere with inner cut-out, use [`CSGDifference`](@ref):
```yaml
difference:
  - sphere:
      r: 2
  - sphere:
      r: 1
```
This is a sphere with inner radius 1 and outer radius 2.

See also [Constructive Solid Geometry (CSG)](@ref).
"""
struct Ellipsoid{T,CO,TR,TP,TT} <: AbstractVolumePrimitive{T, CO}
    r::TR
    φ::TP
    θ::TT

    origin::CartesianPoint{T}
    rotation::SMatrix{3,3,T,9}
end

#Type conversion happens here
function Ellipsoid{T,CO}(r, origin, rotation) where {T,CO}
    _r = _csg_convert_args(T, r)
    Ellipsoid{T,CO,typeof(_r),Nothing,Nothing}(_r, nothing, nothing, origin, rotation)
end

#Type promotion happens here
function Ellipsoid(CO, r::TR, origin::PT, rotation::ROT) where {TR, PT, ROT}
    eltypes = _csg_get_promoted_eltype.((TR, PT, ROT))
    T = float(promote_type(eltypes...))
    Ellipsoid{T,CO}(r, origin, rotation)
end

function Ellipsoid(::Type{CO}=ClosedPrimitive;
    r = 1,
    origin = zero(CartesianPoint{Int64}),
    rotation = one(SMatrix{3, 3, Int64, 9})
) where {CO}
    Ellipsoid(CO, r, origin, rotation)
end

function Ellipsoid{T}(::Type{CO}=ClosedPrimitive;
    r = 1.0,
    origin = zero(CartesianPoint{Float64}),
    rotation = one(SMatrix{3, 3, Float64, 9})
) where {T, CO}
    Ellipsoid{T,CO}( r, origin, rotation)
end

Ellipsoid{T,CO,TR,TP,TT}( e::Ellipsoid{T,CO,TR,TP,TT}; COT = CO,
            origin::CartesianPoint{T} = e.origin,
            rotation::SMatrix{3,3,T,9} = e.rotation) where {T,CO<:Union{ClosedPrimitive, OpenPrimitive},TR,TP,TT} =
    Ellipsoid{T,COT,TR,TP,TT}(e.r,e.φ,e.θ,origin,rotation)

const Sphere{T,CO,TP,TT} = Ellipsoid{T,CO,T,TP,TT}
const FullSphere{T,CO} = Ellipsoid{T,CO,T,Nothing,Nothing}
const FullEllipsoid{T,CO} = Ellipsoid{T,CO,NTuple{3,T},Nothing,Nothing}

function Geometry(::Type{T}, ::Type{Ellipsoid}, dict::AbstractDict, input_units::NamedTuple, transformations::Transformations{T}) where {T}
    length_unit = input_units.length
    angle_unit = input_units.angle
    origin = get_origin(T, dict, length_unit)
    rotation = get_rotation(T, dict, angle_unit)

    r = parse_r_of_primitive(T, dict, input_units.length)
    φ = parse_φ_of_primitive(T, dict, angle_unit)
    if !(φ === nothing)
        error("Partial Ellipsoid (`φ = φ`) is not yet supported.")
    end
    θ = parse_θ_of_primitive(T, dict, angle_unit)
    if !(θ === nothing)
        error("Partial Ellipsoid (`θ = θ`) is not yet supported.")
    end
    e = if r isa Tuple{T,T}
        Ellipsoid{T}(ClosedPrimitive,
            r = r[2],
            origin = origin,
            rotation = rotation,
        ) - Ellipsoid{T}(OpenPrimitive,
            r = r[1],
            origin = origin,
            rotation = rotation,
        )
    else
        Ellipsoid{T}(ClosedPrimitive,
            r = r,
            origin = origin,
            rotation = rotation,
        )
    end
    transform(e, transformations)
end

function Dictionary(e::Ellipsoid{T})::OrderedDict{String, Any} where {T}
    dict = OrderedDict{String, Any}()
    dict["r"] = e.r # always a Real
    if !isnothing(e.φ) error("Partial Ellipsoid (`φ = φ`) is not yet supported.") end
    if !isnothing(e.θ) error("Partial Ellipsoid (`θ = θ`) is not yet supported.") end
    if e.origin != zero(CartesianVector{T}) dict["origin"] = e.origin end
    if e.rotation != one(SMatrix{3,3,T,9}) dict["rotation"] = Dictionary(e.rotation) end
    OrderedDict{String, Any}("sphere" => dict)
end

_in(pt::CartesianPoint{T}, s::FullSphere{<:Any, ClosedPrimitive}; csgtol::T = csg_default_tol(T)) where {T} =
    hypot(pt.x, pt.y, pt.z) <= s.r + csgtol

_in(pt::CartesianPoint{T}, s::FullSphere{<:Any, OpenPrimitive}; csgtol::T = csg_default_tol(T)) where {T} =
    hypot(pt.x, pt.y, pt.z) < s.r - csgtol

function surfaces(e::Ellipsoid{T,ClosedPrimitive}) where {T}
    em = EllipsoidMantle{T,typeof(e.r),typeof(e.φ),typeof(e.θ),:inwards}(e.r, e.φ, e.θ, e.origin, e.rotation)
    (em,)
end
function surfaces(e::Ellipsoid{T,OpenPrimitive}) where {T}
    em = EllipsoidMantle{T,typeof(e.r),typeof(e.φ),typeof(e.θ),:outwards}(e.r, e.φ, e.θ, e.origin, e.rotation)
    (em,)
end

extremum(e::Ellipsoid{T}) where {T} = max(e.r...)


text:
"""
    traintensor([T = N0f8], [indices]; [dir]) -> Array{T}

Returns the MNIST **training** images corresponding to the given
`indices` as a multi-dimensional array of eltype `T`.

The image(s) is/are returned in the horizontal-major
memory layout as a single numeric array. If `T <: Integer`, then
all values will be within `0` and `255`, otherwise the values are
scaled to be between `0` and `1`.

If the parameter `indices` is omitted or an `AbstractVector`, the
images are returned as a 3D array (i.e. a `Array{T,3}`), in
WHN format (width, height, #images).
For integer `indices` instead, a 2D array in WH format is returned.

```julia-repl
julia> MNIST.traintensor() # load all training images
28×28×60000 Array{N0f8,3}:
[...]

julia> MNIST.traintensor(Float32, 1:3) # first three images as Float32
28×28×3 Array{Float32,3}:
[...]
```

If `indices` is an `Integer`, the single image is returned as
`Matrix{T}`.

```julia-repl
julia> MNIST.traintensor(1) # load first training image
28×28 Array{N0f8,2}:
[...]
```

You can use the utility function
[`convert2image`](@ref) to convert an MNIST array into a
vertical-major Julia image with the corrected color values.

```julia-repl
julia> MNIST.convert2image(MNIST.traintensor(1)) # convert to column-major colorant array
28×28 Array{Gray{N0f8},2}:
[...]
```

$(download_docstring("MNIST", DEPNAME))
"""
function traintensor(::Type{T}, args...; dir = nothing) where T
    path = datafile(DEPNAME, TRAINIMAGES, dir)
    images = Reader.readimages(path, args...)
    bytes_to_type(T, images)
end

function traintensor(args...; dir = nothing)
    traintensor(N0f8, args...; dir = dir)
end

"""
    testtensor([T = N0f8], [indices]; [dir]) -> Array{T}

Returns the MNIST **test** images corresponding to the given
`indices` as a multi-dimensional array of eltype `T`.

The image(s) is/are returned in the horizontal-major
memory layout as a single numeric array. If `T <: Integer`, then
all values will be within `0` and `255`, otherwise the values are
scaled to be between `0` and `1`.

If the parameter `indices` is omitted or an `AbstractVector`, the
images are returned as a 3D array (i.e. a `Array{T,3}`), in
WHN format (width, height, #images).
For integer `indices` instead, a 2D array in WH format is returned.

```julia-repl
julia> MNIST.testtensor() # load all test images
28×28×10000 Array{N0f8,3}:
[...]

julia> MNIST.testtensor(Float32, 1:3) # first three images as Float32
28×28×3 Array{Float32,3}:
[...]
```

If `indices` is an `Integer`, the single image is returned as
`Matrix{T}`.

```julia-repl
julia> MNIST.testtensor(1) # load first test image
28×28 Array{N0f8,2}:
[...]
```

You can use the utility function
[`convert2image`](@ref) to convert an MNIST array into a
vertical-major Julia image with the corrected color values.

```julia-repl
julia> MNIST.convert2image(MNIST.testtensor(1)) # convert to column-major colorant array
28×28 Array{Gray{N0f8},2}:
[...]
```

$(download_docstring("MNIST", DEPNAME))
"""
function testtensor(::Type{T}, args...; dir = nothing) where T
    path = datafile(DEPNAME, TESTIMAGES, dir)
    images = Reader.readimages(path, args...)
    bytes_to_type(T, images)
end

function testtensor(args...; dir = nothing)
    testtensor(N0f8, args...; dir = dir)
end

"""
    trainlabels([indices]; [dir])

Returns the MNIST **trainset** labels corresponding to the given
`indices` as an `Int` or `Vector{Int}`. The values of the labels
denote the digit that they represent. If `indices` is omitted,
all labels are returned.

```julia-repl
julia> MNIST.trainlabels() # full training set
60000-element Array{Int64,1}:
 5
 0
 ⋮
 6
 8

julia> MNIST.trainlabels(1:3) # first three labels
3-element Array{Int64,1}:
 5
 0
 4

julia> MNIST.trainlabels(1) # first label
5
```

$(download_docstring("MNIST", DEPNAME))
"""
function trainlabels(args...; dir = nothing)
    path = datafile(DEPNAME, TRAINLABELS, dir)
    Vector{Int}(Reader.readlabels(path, args...))
end

function trainlabels(index::Integer; dir = nothing)
    path = datafile(DEPNAME, TRAINLABELS, dir)
    Int(Reader.readlabels(path, index))
end

"""
    testlabels([indices]; [dir])

Returns the MNIST **testset** labels corresponding to the given
`indices` as an `Int` or `Vector{Int}`. The values of the labels
denote the digit that they represent. If `indices` is omitted,
all labels are returned.

```julia-repl
julia> MNIST.testlabels() # full test set
10000-element Array{Int64,1}:
 7
 2
 ⋮
 5
 6

julia> MNIST.testlabels(1:3) # first three labels
3-element Array{Int64,1}:
 7
 2
 1

julia> MNIST.testlabels(1) # first label
7
```

$(download_docstring("MNIST", DEPNAME))
"""
function testlabels(args...; dir = nothing)
    path = datafile(DEPNAME, TESTLABELS, dir)
    Vector{Int}(Reader.readlabels(path, args...))
end

function testlabels(index::Integer; dir = nothing)
    path = datafile(DEPNAME, TESTLABELS, dir)
    Int(Reader.readlabels(path, index))
end

"""
    traindata([T = N0f8], [indices]; [dir]) -> Tuple

Returns the MNIST **trainingset** corresponding to the given
`indices` as a two-element tuple. If `indices` is omitted the
full trainingset is returned. The first element of three return
values will be the images as a multi-dimensional array, and the
second element the corresponding labels as integers.

The image(s) is/are returned in the horizontal-major
memory layout as a single numeric array of eltype `T`. If `T <:
Integer`, then all values will be within `0` and `255`, otherwise
the values are scaled to be between `0` and `1`. The integer
values of the labels correspond 1-to-1 the digit that they
represent.

```julia
train_x, train_y = MNIST.traindata() # full datatset
train_x, train_y = MNIST.traindata(2) # only second observation
train_x, train_y = MNIST.traindata(dir="./MNIST") # custom folder
```

$(download_docstring("MNIST", DEPNAME))

Take a look at [`MNIST.traintensor`](@ref) and
[`MNIST.trainlabels`](@ref) for more information.
"""
function traindata(::Type{T}, args...; dir = nothing) where T
    (traintensor(T, args...; dir = dir),
     trainlabels(args...; dir = dir))
end

traindata(args...; dir = nothing) = traindata(N0f8, args...; dir = dir)

"""
    testdata([T = N0f8], [indices]; [dir]) -> Tuple

Returns the MNIST **testset** corresponding to the given
`indices` as a two-element tuple. If `indices` is omitted the
full testset is returned. The first element of three return
values will be the images as a multi-dimensional array, and the
second element the corresponding labels as integers.

The image(s) is/are returned in the horizontal-major
memory layout as a single numeric array of eltype `T`. If `T <:
Integer`, then all values will be within `0` and `255`, otherwise
the values are scaled to be between `0` and `1`. The integer
values of the labels correspond 1-to-1 the digit that they
represent.

```julia
test_x, test_y = MNIST.testdata() # full datatset
test_x, test_y = MNIST.testdata(2) # only second observation
test_x, test_y = MNIST.testdata(dir="./MNIST") # custom folder
```

$(download_docstring("MNIST", DEPNAME))

Take a look at [`MNIST.testtensor`](@ref) and
[`MNIST.testlabels`](@ref) for more information.
"""
function testdata(::Type{T}, args...; dir = nothing) where T
    (testtensor(T, args...; dir = dir),
     testlabels(args...; dir = dir))
end

testdata(args...; dir = nothing) = testdata(N0f8, args...; dir = dir)


text:
import Base: -, +, *, filter
import LinearAlgebra: norm

abstract type AbstractProperty end

@inline +(φ1::T, φ2::T) where {T <: AbstractProperty} = T( φ1.val + φ2.val )
@inline -(φ1::T, φ2::T) where {T <: AbstractProperty} = T( φ1.val - φ2.val )
@inline -(φ::T) where {T <: AbstractProperty} = T( -φ.val)
@inline *(a::Union{Number, AbstractMatrix}, φ::T) where {T <: AbstractProperty} =
      T(a * φ.val)
@inline *(φ::T, a::Union{Number, AbstractMatrix}) where {T <: AbstractProperty} =
      T(φ.val * a)
@inline norm(φ::T) where {T <: AbstractProperty} = norm(φ.val)
@inline Base.length(φ::AbstractProperty) = length(φ.val)
@inline Base.size(φ::AbstractProperty) = size(φ.val)
@inline Base.zero(φ::T) where {T <: AbstractProperty} = T(zero(φ.val))
@inline Base.zero(::Type{T}) where {T <: AbstractProperty} = zero(T())

Base.isapprox(φ1::T, φ2::T) where {T <: AbstractProperty} =
      isapprox(φ1.val, φ2.val)

"""
`struct Invariant{D}` : specifies that the output of an ACE is
an invariant scalar.
"""
struct Invariant{T} <: AbstractProperty
   val::T
end

Invariant{T}() where {T <: Number} = Invariant{T}(zero(T))

Invariant(T::DataType = Float64) = Invariant{T}()

filter(φ::Invariant, b::Array) = ( length(b) <= 1 ? true :
     iseven(sum(bi.l for bi in b)) && iszero(sum(bi.m for bi in b))  )

rot3Dcoeffs(::Invariant, T=Float64) = Rot3DCoeffs(T)

# TODO: this is a naive implementation of differentiation.
#       cf https://github.com/ACEsuit/ACE.jl/issues/27
#       for further discussion
*(φ::Invariant, dAA::SVector) = φ.val * dAA

@doc raw"""
`struct EuclideanVector{D, T}` : specifies that the output $\varphi$ of an
ACE is an equivariant $\varphi \in \mathbb{R}^{3}$, i.e., it transforms under
$O(3)$ as
```math
      \varphi \circ Q = Q \cdot \varphi,
```
where $\cdot$ denotes the standard matrix-vector product.
"""
struct EuclideanVector{T} <: AbstractProperty
   val::SVector{3, T}
end

EuclideanVector{T}() where {T <: Number} = EuclideanVector{T}(zero(SVector{3, T}))

EuclideanVector(T::DataType=Float64) = EuclideanVector{T}()

filter(φ::EuclideanVector, b::Array) = ( length(b) <= 1 ? true :
             isodd( sum(bi.l for bi in b)) &&
            (abs(sum(bi.m for bi in b)) <= 1) )

rot3Dcoeffs(::EuclideanVector,T=Float64) = Rot3DCoeffsEquiv{T,1}(Dict[], ClebschGordan(T))

# differentiation - cf #27
*(φ::EuclideanVector, dAA::SVector) = φ.val * dAA'

struct SphericalVector{L, LEN, T} <: AbstractProperty
   val::SVector{LEN, T}
   _valL::Val{L}
end

# differentiation - cf #27
*(φ::SphericalVector, dAA::SVector) = φ.val * dAA'

getL(φ::SphericalVector{L}) where {L} = L

# L = 0 -> (0,0)
# L = 1 -> (0,0), (1,-1), (1,0), (1,1)  -> 4
# L = 3 ->  ... + 5 -> 9
# 1 + 3 + 5 + ... + 2*L+1
# = L + 2 * (1 + ... + L) = L+1 + 2 * L * (L+1) / 2 = (L+1)^2
function SphericalVector(L::Integer; T = Float64)
   LEN = 2L+1   # length of SH basis up to L
   return SphericalVector( zero(SVector{LEN, T}), Val{L}() )
end

function SphericalVector{L, LEN, T}(x::AbstractArray) where {L, LEN, T}
   @assert length(x) == LEN
   SphericalVector{L, LEN, T}( SVector{LEN, T}(x...), Val(L) )
end

SphericalVector{L, LEN, T}()  where {L, LEN, T} =
      SphericalVector( zero(SVector{LEN, T}), Val{L}() )

filter(φ::SphericalVector, b::Array) = ( length(b) <= 1 ? true :
        ( ( iseven(sum(bi.l for bi in b)) == iseven(getL(φ)) ) &&
         ( abs(sum(bi.m for bi in b)) <= getL(φ) )  ) )

rot3Dcoeffs(::SphericalVector, T::DataType=Float64) = Rot3DCoeffs(T)

struct SphericalMatrix{L1, L2, LEN1, LEN2, T} <: AbstractProperty
   val::SMatrix{LEN1, LEN2, T}
   _valL1::Val{L1}
   _valL2::Val{L2}
end

# differentiation - cf #27
# actually this here appears to be the generic form how to do the
# differentiation for arbtirary order tensors.
*(φ::SphericalMatrix{L1, L2, LEN1, LEN2}, dAA::SVector{N}
      ) where {L1, L2, LEN1, LEN2, N} =
      reshape(φ.val[:] * dAA', Size(LEN1, LEN2, N))

getL(φ::SphericalMatrix{L1,L2}) where {L1,L2} = L1, L2

# L = 0 -> (0,0)
# L = 1 -> (0,0), (1,-1), (1,0), (1,1)  -> 4
# L = 3 ->  ... + 5 -> 9
# 1 + 3 + 5 + ... + 2*L+1
# = L + 2 * (1 + ... + L) = L+1 + 2 * L * (L+1) / 2 = (L+1)^2
function SphericalMatrix(L1::Integer, L2::Integer; T = Float64)
   LEN1 = 2L1+1   # length of SH basis up to L
   LEN2 = 2L2+1
   return SphericalMatrix( zero(SMatrix{LEN1, LEN2, T}), Val{L1}(), Val{L2}() )
end

function SphericalMatrix{L1, L2, LEN1, LEN2, T}(x::AbstractMatrix) where {L1, L2, LEN1, LEN2, T}
   @assert size(x) == (LEN1, LEN2)
   SphericalMatrix{L1, L2, LEN1, LEN2, T}( SMatrix{LEN1, LEN2, T}(x...), Val(L1), Val(L2) )
end

SphericalMatrix{L1, L2, LEN1, LEN2, T}()  where {L1, L2, LEN1, LEN2, T} =
      SphericalMatrix( zero(SMatrix{LEN1, LEN2, T}), Val{L1}(), Val{L2}() )

filter(φ::SphericalMatrix, b::Array) = ( length(b) < 1 ? true :
        ( ( iseven(sum(bi.l for bi in b)) == iseven(sum(getL(φ))) ) &&
         ( abs(sum(bi.m for bi in b)) <= sum(getL(φ)) )  ) )

rot3Dcoeffs(::SphericalMatrix, T::DataType=Float64) = Rot3DCoeffs(T)


text:
function __solve(prob::OptimizationProblem, opt::MultistartOptimization.TikTak;
                 local_method, local_maxiters = nothing,
                 progress = false, kwargs...)
    local x, _loss

    if !(isnothing(local_maxiters)) && local_maxiters <= 0.0
        error("The number of local_maxiters has to be a non-negative and non-zero number.")
    else !(isnothing(local_maxiters))
        local_maxiters = convert(Int, local_maxiters)
    end

    _loss = function(θ)
        x = prob.f(θ, prob.p)
        return first(x)
    end

    t0 = time()

    P = MultistartOptimization.MinimizationProblem(_loss, prob.lb, prob.ub)
    multistart_method = opt
    if !(isnothing(local_maxiters))
        local_method = MultistartOptimization.NLoptLocalMethod(local_method, maxeval = local_maxiters)
    else
        local_method = MultistartOptimization.NLoptLocalMethod(local_method)
    end
    p = MultistartOptimization.multistart_minimization(multistart_method, local_method, P)

    t1 = time()

    SciMLBase.build_solution(prob, opt, p.location, p.value; original=p)
end


text:
include("MultiCellArrays.jl")
@reexport using Gridap.MultiCellArrays

include("MultiCellMaps.jl")
@reexport using Gridap.MultiCellMaps

include("MultiFESpaces.jl")
@reexport using Gridap.MultiFESpaces

include("MultiAssemblers.jl")
@reexport using Gridap.MultiAssemblers

include("MultiFEFunctions.jl")
@reexport using Gridap.MultiFEFunctions

include("MultiFEBases.jl")
@reexport using Gridap.MultiFEBases

include("MultiFEOperators.jl")
@reexport using Gridap.MultiFEOperators