text:
"""

## Envelope

Holds a single MLine (often an *upper envelope*), and a vector of `Point`s marking the intersections between MLines.

### Fields

* `env    `: Cleaned upper envelope
* `isects `: Vector of intersecting `Point`s (new points in clean env)
* `removed`: Vector of indices of Points removed from `dirty` MLine during Envelope assembly
* `vbound` : Value on lower bound of asset domain

"""
mutable struct Envelope{T<:Number}
    env       :: MLine{T}
    isects    :: Vector{Point{T}}
    removed   :: Vector{Int}
    vbound    :: T
    # Envelope(1) builds a test object
    # function Envelope(x::T) where {T<:Number}
    #     this = new{T}()
    #     this.dirty = MLine([typemin(T)],[typemin(T)])
    #     this.env = MLine([typemin(T)],[typemin(T)])
    #     this.isects = Point{T}[]
    #     this.removed = Int[ ]
    #     this.vbound = zero(T)
    #     return this
    # end
    # Envelope(MLine) builds an object with the MLine as envelope
    function Envelope(d::MLine{T}) where {T<:Number}
        this = new{T}()
        this.env = deepcopy(d)
        this.isects = Point{T}[]
        this.removed = Int[ ]
        this.vbound = zero(T)
        return this
    end
end
function show(io::IO, ::MIME"text/plain", en::Envelope{T}) where {T<:Number}
    print(io,"$T Envelope\n")
    print(io,"num of intersections: $(length(en.isects))\n")
    print(io,"num of pts removed: $(length(en.removed))\n")
end
function show(io::IO, en::Envelope{T}) where {T<:Number}
    print(io,length(en.env),"-point $T Envelope")
end

size(e::Envelope) = size(e.L)
eltype(e::Envelope) = eltype(e.L)
bound(e::Envelope) = e.vbound
gets(en::Envelope) = en.isects
getr(en::Envelope) = en.removed
Base.getindex(en::Envelope,id::Int) = en.L[id]

"""
    removed!(e::Envelope,L::MLine)

Find which points from initial `MLine` did not end up in the
`env` and write them to `removed`.
"""
function removed!(e::Envelope,d::MLine)
    # find all points in dirty that did not make it into clean
    rmvd = findall( d .∉ Ref(e.env.v))
    if length(rmvd) > 0
        push!(e.removed,rmvd...)   # splice array into this call
    end
    return nothing
end

"""
    splitLine(m::MLine)

splits a `MLine` object at wrong EGM solution points. Wrong solutions appear in kinked regions.
Returns an array of `MLine`.
"""
function splitLine(o::MLine{T}) where T<:Number

    # 1) find all jump-backs in x-grid
    # println(o.x)
    xvec = getx(o)
    ii = xvec[2:end].>xvec[1:end-1]

    # 2) if no backjumps at all, exit
    if all(ii)
        # return as is
        return [o]   # always return a vector or MLine
    else
    # 3) else, identify subsets
        i = 1
        sections = MLine{T}[]  # an array of Lines
        new_sections = MLine{T}[]  # an array of Lines
        while true
            j = findfirst(ii .!= ii[1])  # identifies all indices within kinked region from left to right until the first kink
            # if no more kinks
            if isnothing(j)
                if i > 1
                    # add remaining MLine
                    push!(sections,o)
                end
                # then break
                break
            end
            newm,o = splitat(o,j)  # split old MLine at j
            push!(sections,newm)
            ii = ii[j:end] # chop off from jump index
            i += 1
        end

        for s in sections
            # fedor keeps those
            # if (!issorted(s.v)) && (length(s.v)==2)
                # disregard those
            # else
                if !issorted(s)
                    sortx!(s)
                end
                push!(new_sections,s)
            # end
        end
        return new_sections
    end
end

"""
    upper_env!(e::Envelope{T}) where T<:Number

### Outline

This function computes the *upper envelope* over a set of split lines.
In particular:

1. generates a common support `xx` by concatenating the `x` coords of each `MLine` in `e.L`
2. interpolates all `MLine` in `e.L` on that support `xx`
3. Remembers which indices of each line have been extrapolated, so they cannot be the optimum
4. finds the index of the `MLine` with maximal `y`-value for each `xx` (among non extrapolated ones)
5. assembles the *upper envelope* by just combining all values from 3. into a new `MLine`. This means to collect intersections between lines (i.e. new points) along the way as well.
"""
function upper_env(L::Vector{MLine{T}}) where T<:Number
    # compute upper envelope of all lines
        # - get all x's from all s and sort into a vector xx
        # - interpolate(extrapolate) all s on xx
        # - how to deal with points at which some MLine is infeasible? => set to -Inf

    if length(L)<2
        error("an upper envelope requires by definition at least 2 lines.")
    end

    # - get all x's from all Lines and sort into a vector xx
    xx = sort(unique(reduce(vcat,[getx(l) for l in L])))

    # - interpolate(extrapolate) all Ls on xx
    # this returns an array (length(L)) of interpolated MLines
    yy = interp(L,xx,extrap=false)  # set a line to -Inf if xx not in its domain

    n = length(xx)   # consider all points

    r_idx = linemax(yy)  # get row indices only: the row index tells us which MLine was optimal at that point.
    # envelope over all lines: just pick max points for each xx
    top = [yy[r_idx[i]][i] for i in 1:length(r_idx)]

    # build up upper envelope one by one point in xx
    # ==============================================

    env = Point{T}[]
    push!(env,top[1])  # first point in envelope point vector
    isec = Point[]     # empty vector of point to collect intersections

    k0 = r_idx[1]  # start index
    for j in 2:n
        k1 = r_idx[j]
        if k0 != k1
            # @debug "switching from $k0 to $k1 at index $j"
            # is there an intersection between both lines?
            ln1 = k0; ln2 = k1  # candidate line indices
            L1 = L[ln1]; L2 = L[ln2]  # candidate lines
            pt1 = xx[j-1]; pt2 = xx[j]  # intersection ∈ [pt1,pt2]

            y1 = interp(L[ln1],[pt1,pt2],extrap=true)
            y2 = interp(L[ln2],[pt1,pt2],extrap=true)

            # check neither y1 nor y2 have any extrapolated points
            neither = (length(y1.iextrap) == 0) && ( length(y2.iextrap) == 0 )
            # neither = (y1.iextrap != [1,2]) && (y2.iextrap != [1,2])
            # and that we have different signs on both ends of search interval
            f_closure(z) = interp(L1,[z])[1].y - interp(L2,[z])[1].y
            diffsign = f_closure(pt1) * f_closure(pt2) < 0

            if neither && diffsign

                # find intersection point
                while true
                    # @debug "checking lines ln1=$ln1 ln2=$ln2 in [$(round(pt1,digits=3)),$(round(pt1,digits=3))]"
                    pt3 = fzero(x -> interp(L[ln1],[x])[1].y - interp(L[ln2],[x])[1].y,pt1,pt2)
                    y3  = interp(L[ln1],[pt3])[1].y  # get function value of ln1 at intersection
                    # @debug "found intersection" pt3=pt3 y3=y3

                    # are there other lines above this intersection point?
                    # interpolate *all* lines in new point pt3
                    yy2 = interp(L,[pt3],extrap=false)
                    y3, ln3 = findmax([yy2[i].v[1].y for i in 1:length(L)])

                    if ln3 == ln2 || ln3 == ln1
                        # @debug "no higher line at" pt3=pt3
                        # no new lines above at pt3!
                        # add intersection point to the end of the env points
                        push!(env, Point(pt3,y3))
                        # add intersection to intersections container
                        push!(isec, Point(pt3,y3))
                        # additional intersections before next point?
                        if ln2 == k1
                            # no other intersection.
                            # if ln2 is new optimal line and there is no ln3 above,
                            # there cannot be an additional intersection to the right of pt3.
                            break

                        else
                            #
                            # @debug "additional iscecs. updating" ln1=ln2 pt1=pt3 ln2=k1 pt2=xx[j]
                            ln1=ln2
                            pt1=pt3
                            ln2=k1
                            pt2=xx[j]
                        end
                    else
                        # there is indeed a third line over point pt3
                        # what we have is not the upper envelope yet.
                        # redo search.
                        # there must be an intersection to the left of pt3
                        # @debug "additional isec to left. updating" ln2=ln3 pt2=pt3
                        ln2 = ln3 # new candidate line's index
                        pt2 = pt3 # new right boundary of search interval

                    end
                end
            else
                # @debug "no valid intersection" notboth=notboth diffsign=diffsign
            end
        end # if k0 != k1
        # regular addition of points
        # if current xx[j] is a member of line k1, add the point
        # if j is the last index, add the point (no intersections to the right)
        if (xx[j] ∈ L[k1].v) || j==n
            # @debug "adding regular point" x=xx[j] Lx=getx(L[k1].v)
            # look how beautiful that is: x ∈ v where v is an array of x-y points
            push!(env, top[j])  # push that point onto the envelope
        end
        k0 = k1 # update current index
    end # end j in 2:n
    e = Envelope(MLine(env))
    e.isects = isec
    return e
end

"""
    secondary_envelope(L::MLine)

* computes secondary envelope of a `MLine`
* splits line at backward bends and computes upper envelope over those
* returns an `Envelope` object with fields `removed` giving the indices of removed points for each constituting line.
"""
function secondary_envelope(L::MLine)

    # identify loop-backs and split line
    A = splitLine(L)

    if length(A) > 1
        # compute upper envelope over array of lines
        e = upper_env(A)
    else
        e = Envelope(A[1])
    end

    # record which points in original L did not make it into final e.env
    removed!(e,L)

    return e

end


text:
root_dir = "../."
using Pkg; Pkg.activate(root_dir)
using Strathreads, Temporal, Indicators, Dates
using DotEnv

cfg = DotEnv.config(path="$root_dir/.env")
pathCoins = cfg["INPUT_DIR"]
#assets = ["EOS-USDT"]
#A
d= Dict()
_assets = readlines("../data/211223.trade_list")
for asset in _assets
    assets=[asset]
    universe = Universe(assets)
    function datasource(asset::String)::TS
        savedata_path = joinpath(pathCoins, "$asset.csv")
        return Temporal.tsread(savedata_path, indextype=UInt64, format="yyyy-mm-dd HH:MM:SS")
    end
    universe = Universe(assets)
    universe = gather(assets, source=datasource)

    # define indicators and parameter space
    function fun(x::TS; args...)::TS
        close_prices = x[:close]
        macd_ma = macd(close_prices; args...)
        output = [close_prices macd_ma]
        #output.fields = [:close, :MACD] -> shit renames fields
        return output
    end
    arg_names = [:nfast, :nslow]
    arg_defaults = [12, 26]
    arg_ranges = [5:1:12, 15:1:35]
    paramset = ParameterSet(arg_names, arg_defaults, arg_ranges)
    indicator = Indicator(fun, paramset)

    # define signals that will trigger trading decisions
    siglong = @signal  MACD ↑ Signal
    sigexit = @signal MACD == Signal
    sigshort = @signal MACD ↓ Signal

    # define the trading rules
    longrule = @rule siglong → long 100
    shortrule = @rule sigshort → short 100
    exitrule = @rule sigexit → liquidate 1.0
    rules = (longrule, shortrule, exitrule)

    # run strategy
    strat = Strategy(universe, indicator, rules)
    bt = backtest(strat, px_trade=:open, px_close=:close, verbose=true)
    opt = optimize(strat, px_trade=:open, px_close=:close)

    ##visualize
    #using Plots
    #gr()
    #(x, y, z) = (opt[:,i] for i in 1:3)
    #surface(x, y, z)

    #eval best
    opt_idx = findmax(opt[:,3])[2]
    sets = Dict(zip(arg_names, opt[opt_idx, 1:2]))
    println(opt[opt_idx, :])
    d[asset] = opt[opt_idx, :]
end

ext: jl

max_stars_repo_name: lithiferous/Strathreads

max_issues_repo_name: lithiferous/Strathreads

max_forks_repo_name: lithiferous/Strathreads


text:
using ForSyDeIO
using Documenter

DocMeta.setdocmeta!(ForSyDeIO, :DocTestSetup, :(using ForSyDeIO); recursive=true)

makedocs(;
    modules=[ForSyDeIO],
    authors="Rodolfo Jordão <jordao@kth.se> and contributors",
    repo="https://github.com/Rojods/ForSyDeIO.jl/blob/{commit}{path}#{line}",
    sitename="ForSyDeIO.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Rojods.github.io/ForSyDeIO.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Rojods/ForSyDeIO.jl",
)

max_stars_repo_name: forsyde/forsyde-io

max_issues_repo_name: forsyde/forsyde-io

max_forks_repo_name: forsyde/forsyde-io


text:
using NNlib, Test

include("activation.jl")
include("conv.jl")
include("pooling.jl")


text:
function g(x, t1, t2, t3, t4)
    return t1 * x^3 + t2 * x^2 + t3 * x + t4
end

h(x) = g(x, 4, -3, 2, 10)

h(1.2)


text:
"""
*Extensions for the Julia docsystem.*

# Introduction

This package provides a collection of useful extensions for Julia's built-in docsystem.
These are features that are still regarded as "experimental" and not yet mature enough to be
considered for inclusion in `Base`, or that have sufficiently niche use cases that including
them with the default Julia installation is not seen as valuable enough at this time.

Currently `DocStringExtensions.jl` exports a collection of so-called "abbreviations", which
can be used to add useful automatically generated information to docstrings. These include
information such as:

  * simplified method signatures;
  * documentation for the individual fields of composite types;
  * import and export lists for modules;
  * and source-linked lists of methods applicable to a particular docstring.

Users are most welcome to suggest additional abbreviation ideas, or implement and submit
them themselves. Julia's strong support for program introspection makes this a reasonably
straight forward process.

Details of the currently available abbreviations can be viewed in their individual
docstrings listed below in the "Exports" section.

# Examples

In simple terms an abbreviation can be used by simply interpolating it into a suitable
docstring. For example:

```julia
using DocStringExtensions

\"""
A short summary of `func`...

\$(SIGNATURES)

where `x` and `y` should both be positive.

# Details

Some details about `func`...
\"""
func(x, y) = x + y
```

`\$(SIGNATURES)` will be replaced in the above docstring with

````markdown
# Signatures

```julia
func(x, y)
```
````

The resulting generated content can be viewed via Julia's `?` mode or, if `Documenter.jl` is
set up, the generated external documentation.

The advantage of using [`SIGNATURES`](@ref) (and other abbreviations) is that docstrings are
less likely to become out-of-sync with the surrounding code. Note though that references to
the argument names `x` and `y` that have been manually embedded within the docstring are, of
course, not updated automatically.

# Exports
$(EXPORTS)

# Imports
$(IMPORTS)

"""
module DocStringExtensions

# Imports.

import LibGit2

# Exports.

export @template, FIELDS, TYPEDFIELDS, EXPORTS, METHODLIST, IMPORTS
export SIGNATURES, TYPEDSIGNATURES, TYPEDEF, DOCSTRING, FUNCTIONNAME
export README, LICENSE

# Includes.

include("utilities.jl")
include("abbreviations.jl")
include("templates.jl")

#
# Bootstrap abbreviations.
#
# Within the package itself we would like to be able to use the abbreviations that have been
# implemented. To do this we need to delay evaluation of the interpolated abbreviations
# until they have all been defined. We use `Symbol`s in place of the actual constants, such
# as `METHODLIST` which is written as `:METHODLIST` instead.
#
# The docstring for the module itself, defined at the start of the file, does not need to
# use `Symbol`s since with the way `@doc` works the module docstring gets inserted at the
# end of the module definition and so has all the definitions already defined.
#
let λ = s -> isa(s, Symbol) ? getfield(DocStringExtensions, s) : s
    for (binding, multidoc) in Docs.meta(DocStringExtensions)
        for (typesig, docstr) in multidoc.docs
            docstr.text = Core.svec(map(λ, docstr.text)...)
        end
    end
end

__init__() = (hook!(template_hook); nothing)

end # module


text:
module ChainRules

using Cassette
using LinearAlgebra
using Base.Broadcast: materialize, materialize!, broadcasted, Broadcasted, broadcastable

import NaNMath, SpecialFunctions, LinearAlgebra, LinearAlgebra.BLAS

export AbstractRule, Rule, frule, rrule

include("differentials.jl")
include("rules.jl")
include("rules/base.jl")
include("rules/array.jl")
include("rules/broadcast.jl")
include("rules/linalg/dense.jl")
include("rules/linalg/diagonal.jl")
include("rules/linalg/symmetric.jl")
include("rules/linalg/factorization.jl")
include("rules/blas.jl")
include("rules/nanmath.jl")
include("rules/specialfunctions.jl")

end # module


text:
import Base.+, Base.-, Base./, Base.*
import SchumakerSpline.evaluate
function evaluate(sf::Sum_Of_Functions, point::Real)
    total = 0.0
    for func in sf.functions_
        total = total + evaluate(func,point)
    end
    return total
end
function (s::Sum_Of_Functions)(x::Union{Real,Date,DateTime,DatePeriod})
    return evaluate(s, x)
end

function derivative(f::Sum_Of_Functions)
    return Sum_Of_Functions( derivative.(f.functions_))
end

function indefinite_integral(f::Sum_Of_Functions)
    return Sum_Of_Functions(indefinite_integral.(f.functions_))
end

function +(f::Sum_Of_Functions,number::Real)
    constant_function = PE_Function(number, 0.0,0.0,0)
    return Sum_Of_Functions(vcat(f.functions_, [constant_function]))
end
function -(f::Sum_Of_Functions, number::Real)
    return +(f, -number)
end
function *(f::Sum_Of_Functions, number::Real)
    funcs = deepcopy(f.functions_)
    for i in 1:length(funcs)
        funcs[i] = funcs[i] * number
    end
    return Sum_Of_Functions(funcs)
end
function /(f::Sum_Of_Functions, number::Real)
    return *(f, 1/number )
end

function +(f1::Sum_Of_Functions, f2::Sum_Of_Functions)
    return Sum_Of_Functions([f1,f2])
end
function +(f1::Sum_Of_Functions, f2::Piecewise_Function)
    added_functions = (f1 .+ f2.functions_)::Array
    return Piecewise_Function(f2.starts_,added_functions)
end

function +(f1::Piecewise_Function, f2::Sum_Of_Functions)
    return +(f2,f1)
end

function -(f1::Sum_Of_Functions, f2::Sum_Of_Functions)
    return Sum_Of_Functions([f1,-1*f2])
end
function -(f1::Sum_Of_Functions, f2::Piecewise_Function)
    return +(f1, -1 * f2)
end

function -(f1::Piecewise_Function, f2::Sum_Of_Functions)
    return +(f1, -1*f2)
end

function *(f1::Sum_Of_Functions,f2::Sum_Of_Functions)
    results = Array{Sum_Of_Functions}(undef, length(f1.functions_))
    for i in 1:length(f1.functions_)
        new_funcs = f1.functions_[i] * f2
        results[i] = new_funcs
    end
    return Sum_Of_Functions(results)
end
function *(f1::Sum_Of_Functions, f2::Piecewise_Function)
    return Piecewise_Function(f2.starts_, f1 .* f2.functions_)
end

function *(f1::Piecewise_Function, f2::Sum_Of_Functions)
    return *(f2,f1)
end


text:
#OVERVIEW
#----------------
#Any JuMP model that describes an optimization problem must have four parts:

#Model Object,
#Variables,
#Objective,
#Constraints.

#PACKAGE IMPORT
#----------------
println("Initializing...")

using JuMP  # Need to say it whenever we use JuMP

using GLPK

using CSV # Import data

using DataFrames # For results export to CSV

using Plots

#=
electrolyzer_tea.jl is a function used to evaluate the TEA of the electrochemical
system, including separations, stack cost, and other operational parameters.
This function takes a series of inputs specific to an electrochemical pathway,
as shown below:

function electrolyzer_TEA(
    flow, product flow in [kg/year]
    product_phase, phase of product ["Gas" or "Liquid"]
    product_mw, molecular weight of product [g/mol]
    product_density, product density [kg/m^3]
    electrons_product, [# e-]
    electrons_co2, [#e-/co2]
    system_potential, whole cell potential [V]
    selectivity, selectivity, or faradaic efficiency of CO2R system to product [%]
    hours_year, hours in a given year [hours/year]
    current_density, [A/cm^2]
    single_pass_conversion, single pass conversion of CO2 to product [%]
    electrolyzer_cost, electrolyzer capital cost [$/m^2]
    electricity_price, [$/kWh]
    water_price, process water cost [$/kg]
    water_rate, rate of water consusumption [kg water/kg product]
    crf capital recovery factor [1/years]
    )
=#
include("./electrolyzer_tea.jl")

#=
The electrolyzer tea function operates identically to the electrolyzer tea function
discussed above, however it returns an array of results for costs rather
than a single value. This allows for a more detailed analysis of costs
included in the TEA. Input values are identical, output values are as follows:

    return:
    electrolyzer_capex + electrolyzer_BOP, electrolyzer costs [$/year]
    PSA_capex + PSA_opex + distil_capex + distil_opex, separations costs [$/year]
    battery_capex, battery energy storage costs [$/year]
    maintenance_cost, system maintenance costs [$/year]
    water_opex annual process water costs [$/year]
=#
include("./electrolyzer_tea_array.jl")

__precompile__()

# Define function to perform optimization routine - decreases runtime
function co2_optimizer()

    #DATA IMPORT
    #--------------------
    input_filename = "input_data_full.csv"
    input_filepath = joinpath(@__DIR__, input_filename)

    # Import point source and point sink data
    data = DataFrame!(CSV.File(input_filepath))

    #=
    Option to exclude specific point types from either sink or
    source data in optimizer. To use, fill the array with the proper name
    of the point source to be excluded, e.g. ["COAL", "GAS"]
    =#
    #--------------------
    source_exclusions =  ["COAL","GAS"]#["None"]#["None"]#["None"]
    sink_exclusions = ["None"]#["Import","Export","Acetic Acid Plant", "Petroleum Terminal"]#["Ethane Crackers"]#["Import","Export"]#["None"]

    for i = 1:length(source_exclusions)

        data = data[data[!,8].!= source_exclusions[i] ,:]

    end

    for i = 1:length(sink_exclusions)

        data = data[data[!,8].!= sink_exclusions[i] ,:]

    end

    println("Source and sink points to be considered:")
    println(size(data,1))

    #VARIABLES and generic functions
    #---------

    # A variable is modeled using @defVar(name of the model object, variable name and bound, variable type)
    # Bound can be lower bound, upper bound or both. If no variable type is defined, then it is treated as
    #real. For binary variable write Bin and for integer use Int.

    # Sink/source point datasets
    #_____________________________________________________________________________

    #= Isolate consumption point data by including only points where
    demand for the molecule is not equal to zero =#
    sinks_C2H4 = data[data[:,4].!= 0  ,:]
    sinks_HCOO = data[data[:,5].!= 0  ,:]
    sinks_CO = data[data[:,6].!= 0  ,:]

    #= Isolate CO2 point source data by including only points where
    feedstock supply is not equal to zero =#
    source_pts = data[data[:,7].!= 0  ,:]
    sink_pts = data[data[:,7].== 0  ,:]

    # Number of chemical products considered
    num_species = 3

    # Find number of CO2 point sources
    num_source_pts = size(data[data[:,7].!= 0  ,:],1)

    # Find number of sink points
    num_sink_pts = size(data[data[:,7].== 0  ,:],1)

    # Find total number of points
    num_pts = size(data,1)

    # Faradays constant
    F = 96485 # [C/mol]

    # Sink/source distance evaluation datasets
    # Goal: build an array showing distance between sink and source points
    #_____________________________________________________________________________

    # Array of lat/long for sources
    lat_source = source_pts[:, 2]
    long_source = source_pts[:, 3]

    # Array of lat/long for sinks
    lat_sink = sink_pts[:, 2]
    long_sink = sink_pts[:, 3]

    # Molecular weights
    #_____________________________________________________________________________
    C_mw = 12.011 # [g/mol]
    CO2_mw = 44 # [g/mol]
    C2H4_mw = 28.05 # [g/mol]
    HCOO_mw = 45.017 # [g/mol]
    CO_mw = 28.01 # [g/mol]
    H2O_mw = 18 # [g/mol]

    product_mw = [C2H4_mw , HCOO_mw, CO_mw]

    # Molecule carbons
    #_____________________________________________________________________________
    C2H4_carbons = 2 # [carbons/molecule]
    HCOO_carbons = 1 # [carbons/molecule]
    CO_carbons = 1 # [carbons/molecule]

    carbons = [C2H4_carbons , HCOO_carbons, CO_carbons]

    # Electrons per mol product
    #_____________________________________________________________________________
    C2H4_electrons_product = 12  # [mol e- / mol product]
    HCOO_electrons_product = 2 # [mol e- / mol product]
    CO_electrons_product = 2 # [mol e- / mol product]
    H2O_electrons = 4 # [mol e- / mol product]

    electrons_product = [C2H4_electrons_product , HCOO_electrons_product, CO_electrons_product]

    # Product state
    #_____________________________________________________________________________
    C2H4_phase = "Gas"
    HCOO_phase = "Liquid"
    CO_phase = "Gas"

    product_phase = [C2H4_phase, HCOO_phase, CO_phase]

    # Product densities
    #_____________________________________________________________________________
    C2H4_density = 1.18 # [kg/m^3]
    HCOO_density = 1221 # [kg/m^3]
    CO_density = 1.14 # [kg/m^3]

    product_density = [C2H4_density , HCOO_density, CO_density]

    # Electrons per mol CO2
    #_____________________________________________________________________________
    C2H4_electrons_co2 = 6  # [mol e- / mol co2]
    HCOO_electrons_co2 = 2 # [mol e- / mol co2]
    CO_electrons_co2 = 2 # [mol e- / mol co2]

    electrons_co2 = [C2H4_electrons_co2 , HCOO_electrons_co2, CO_electrons_co2]

    # Current densities
    #_____________________________________________________________________________
    C2H4_cd = 1.55 # [A/cm^2] 10.1126/science.aay4217
    HCOO_cd = 0.5 # [A/cm^2] 10.1021/acsenergylett.0c00860
    CO_cd = 0.35 # [A/cm^2] http://xlink.rsc.org/?DOI=C9EE02410G

    current_density = [C2H4_cd , HCOO_cd, CO_cd]

    # single pass conversion
    #_____________________________________________________________________________
    C2H4_spc = 0.5 # [%] assumed
    HCOO_spc = 0.5 # [%] assumed
    CO_spc = 0.5 # [%] assumed

    single_pass_conversion = [C2H4_spc , HCOO_spc, CO_spc]

    # Selectivities
    #___________________________________________________________________________
    C2H4_selectivity = 0.6 # 10.1126/science.aay4217
    HCOO_selectivity = 0.93 #  10.1021/acsenergylett.0c00860
    CO_selectivity = 0.95  #  http://xlink.rsc.org/?DOI=C9EE02410G

    selectivity = [C2H4_selectivity , HCOO_selectivity, CO_selectivity]

    # Potentials
    #_____________________________________________________________________________
    C2H4_system_potential = 3.9 #[V] 10.1126/science.aay4217 taken from fig s33, durability study
    HCOO_system_potential = 5.8 #[V] 10.1021/acsenergylett.0c00860 table 1
    CO_system_potential = 3.0 #[V] 10.1021/acs.iecr.7b03514

    system_potential = [C2H4_system_potential , HCOO_system_potential, CO_system_potential]

    # Market Values
    #_____________________________________________________________________________
    C2H4_value = 0.555 # [$/kg] C2H4 value from PEP yearbook and: https://www.spglobal.com/platts/en/market-insights/podcasts/focus/040720-north-america-lng-exports-price-plunge
    HCOO_value = 0.735 # [$/kg] CO value https://pubs.acs.org/doi/10.1021/acs.iecr.7b03514
    CO_value = 0.6 # [$/kg] CO value https://pubs.acs.org/doi/10.1021/acs.iecr.7b03514

    value = [C2H4_value , HCOO_value, CO_value]

    # Market Demand
    #_____________________________________________________________________________
    C2H4_demand = sink_pts[:,4]  # consumption at each point in [kg/year]
    HCOO_demand = sink_pts[: ,5] # consumption at each point in [kg/year]
    CO_demand = sink_pts[: ,6] # consumption at each point in [kg/year]

    demand = hcat(C2H4_demand, HCOO_demand, CO_demand)

    # Water rate of consumption
    #_____________________________________________________________________________
    C2H4_water = 1.28 #[kg H2O/kg product]
    HCOO_water = 0.391 #[kg H2O/kg product]
    CO_water = 0.6426 #[kg H2O/kg product]

    water_rate = [C2H4_water, HCOO_water, CO_water]

    #_____________________________________________________________________________
    capacity = data[data[:,7].!= 0  ,7] * 1e3 * (C_mw / CO2_mw) # Feedstock carbon at each point in [kg C/year]

    transport_A = 1e-5 # Best fit line for CO2 transportation costs from IPCC (see Excel sheet)
    transport_B = -0.0002 # Best fit line for CO2 transportation costs from IPCC (see Excel sheet)

    #CO2 flue gas capture cost for each location [$/kg C]
    CO2_cap_cost = data[data[:,7].!= 0  ,12] / 1e3 * (C_mw / CO2_mw) #[$/kg C]

    # Spatially distinct electricity prices in [$/kWh] (sink pts)
    electricity_price = data[data[:,7].!= 0  ,15]

    # Economic assumptions
    #_____________________________________________________________________________

    hours_year = 8000 # [hours of operation annually]
    electrolyzer_cost = 10000 # [$/m^2]
    interest_rate = 0.08 # [%]
    facility_lifetime = 15 #[years]
    water_price = 8.929e-4 # [$/kg] https://www.energy.gov/sites/prod/files/2017/10/f38/water_wastewater_escalation_rate_study.pdf
    crf = (interest_rate * (1 + interest_rate)^facility_lifetime) / ((1 + interest_rate)^facility_lifetime - 1)

    # Define an array with rows = number of sources, columns = number of sinks, and sheets = number of products
    D = zeros(Float16, num_source_pts, num_sink_pts, num_species)

    # Distance constraint evaluation.
    #_____________________________________________________________________________
    #=
    Define an array with rows = source points, columns = sink points,
    sheets = number of species. Default values are zeros. The loop below changes
    the value to 1 if distance between sink/source is <1,000 km.
    This is then incorporated in the flow variable definition as an upper
    bound on the variable.
    =#
    up_bound = zeros(num_source_pts,num_sink_pts,num_species)
    for i in 1:num_source_pts
        for j in 1:num_sink_pts
            for k in 1:num_species
                pt1 = [long_source[i]; lat_source[i]]
                pt2 = [long_sink[j]; lat_sink[j]]

                D[i, j, k] = (abs(pt1[1] - pt2[1]) + abs(pt1[2] - pt2[2])) / 1e3

                if (abs(source_pts[i,3] - sink_pts[j,3]) + abs(source_pts[i,2] - sink_pts[j,2])) / 1e3 < 1000
                    up_bound[i, j, k] = 1.0

                end
            end
        end
    end

    #MODEL CONSTRUCTION
    #--------------------

    # Preparing an optimization model
    eco2RR = Model(GLPK.Optimizer)

    #OPTIMIZATION VARIABLES
    #--------------------
    #Initialize flow array [kg/year]
    @variable(eco2RR, flow[i = 1:num_source_pts , j = 1:num_sink_pts , k = 1:num_species] >= 0, upper_bound = up_bound[i, j, k].*1e10) # distance constraint
    #@variable(eco2RR, flow[i = 1:num_source_pts , j = 1:num_sink_pts , k = 1:num_species] >= 0)

    #OPTIMIZATION CONSTRAINTS
    #--------------------

    # Ship no more than plant capacity, in units of [kg carbon]
    @constraint(eco2RR, capacity_con[i in 1:num_source_pts],
    + sum(flow[i, : , 1]) * carbons[1] * C_mw / product_mw[1] / selectivity[1]
    + sum(flow[i, : , 2]) * carbons[2] * C_mw / product_mw[2] / selectivity[2]
    + sum(flow[i, : , 3]) * carbons[3] * C_mw / product_mw[3] / selectivity[3]
    <= capacity[i])

    #=
    Constrain total electrochemical system power consumption at a source point
    to less than 300MW. Provides a realistic constraint on system size

    Modifying this equality only applies when the system is feedstock constrained
    i.e. total demand is greater than carbon feedstock available

    Two key options are available:

    == equality - model required to export 300MW worth at each CO2 point.

    <= equality - model can provide anything up to 300MW system at a CO2 point source.

    =#

    @constraint(eco2RR, elec_con[i in 1:num_source_pts],
    sum(
    + sum(flow[i, :, 1]) / hours_year / 3600 * 1e3 / product_mw[1] * electrons_product[1] * F * system_potential[1] / selectivity[1] / 1e6
    + sum(flow[i, :, 2]) / hours_year / 3600 * 1e3 / product_mw[2] * electrons_product[2] * F * system_potential[2] / selectivity[2] / 1e6
    + sum(flow[i, :, 3]) / hours_year / 3600 * 1e3 / product_mw[3] * electrons_product[3] * F * system_potential[3] / selectivity[3] / 1e6
    )
    <= 500)

    #=
    Choose how the model attempts to meet demand for different molecules.
    This constraint can cause the model to be unable to find a SOLUTION
    if feedstock inputs are not sufficient to meet demand.

    Modifying this equality only applies when the system is demand constrained
    i.e. total demand is less than carbon feedstock available

    Two key options are available:

    == equality - model required to meet demand at each point. In this scenario,
        if feedstock inputs are not sufficient then optimizer will not find a
        solution

    <= equality - model can provide anything up to market demand at a point.
        Using this equality allows the optimizer to only make flows it deems
        economically sustainable (positive NPV).

    =#
    @constraint(eco2RR, demand_con[j in 1:num_sink_pts, k in 1:num_species],
    sum(flow[:, j , k]) == demand[j,k])

    println("Constraints Defined")

    #OPTIMIZATION OBJECTIVE FUNCTION
    #--------------------
    println("Beginning Optimization...")

    @objective(
        eco2RR, # Designate model name
        Max, # Tell JuMP to minimize/optimize
        sum(
        (

        - (D[i, j, k] * transport_A + transport_B) * flow[i, j, k] # Transportation cost
        - electrolyzer_TEA(
            flow[i, j, k],
            product_phase[k],
            product_mw[k],
            product_density[k],
            electrons_product[k],
            electrons_co2[k],
            system_potential[k],
            selectivity[k],
            hours_year,
            current_density[k],
            single_pass_conversion[k],
            electrolyzer_cost,
            electricity_price[i],
            water_price,
            water_rate[k],
            crf) # System TEA
        - flow[i, j, k] * carbons[k] * C_mw / product_mw[k] * CO2_cap_cost[i] # Carbon capture cost
        - flow[i, j, k] / hours_year / 3600 * 1e3 / product_mw[k] * electrons_product[k] * F * system_potential[k] / selectivity[k] / 1e3 * hours_year * electricity_price[i] # Electrochemical system electricity cost
        + flow[i, j, k] * value[k] # Product value
        )

        for
            i in 1:num_source_pts, # Iterate over source points
            j in 1: num_sink_pts, # Iterate over sink points
            k in 1: num_species#, # Iterate over number of species
    ))

    # Solving the optimization problem
    optimize!(eco2RR)

    # Set objective value
    obj_value = JuMP.objective_value(eco2RR)

    #OPTIMIZATION SOLUTION WRANGLING
    #--------------------
    println("Optimizer termination status:")
    println(termination_status(eco2RR))
    println(primal_status(eco2RR))
    println("Objective function value:")
    println(obj_value)

    #Output filepath and variable name
    #--------------------

    C2H4_result_filepath = joinpath(@__DIR__, "results/flows/C2H4_flows.csv")
    HCOO_result_filepath = joinpath(@__DIR__, "results/flows/HCOO_flows.csv")
    CO_result_filepath = joinpath(@__DIR__, "results/flows/CO_flows.csv")

    C2H4_NPV_filepath = joinpath(@__DIR__, "results/npv/C2H4_NPV.csv")
    HCOO_NPV_filepath = joinpath(@__DIR__, "results/npv/HCOO_NPV.csv")
    CO_NPV_filepath = joinpath(@__DIR__, "results/npv/CO_NPV.csv")

    C2H4_trans_filepath = joinpath(@__DIR__, "results/econ/trans/C2H4_trans.csv")
    HCOO_trans_filepath = joinpath(@__DIR__, "results/econ/trans/HCOO_trans.csv")
    CO_trans_filepath = joinpath(@__DIR__, "results/econ/trans/CO_trans.csv")

    C2H4_tea_filepath = joinpath(@__DIR__, "results/econ/tea/C2H4_tea.csv")
    HCOO_tea_filepath = joinpath(@__DIR__, "results/econ/tea/HCOO_tea.csv")
    CO_tea_filepath = joinpath(@__DIR__, "results/econ/tea/CO_tea.csv")

    C2H4_ccu_filepath = joinpath(@__DIR__, "results/econ/ccu/C2H4_ccu.csv")
    HCOO_ccu_filepath = joinpath(@__DIR__, "results/econ/ccu/HCOO_ccu.csv")
    CO_ccu_filepath = joinpath(@__DIR__, "results/econ/ccu/CO_ccu.csv")

    C2H4_elec_filepath = joinpath(@__DIR__, "results/econ/electricity/C2H4_elec.csv")
    HCOO_elec_filepath = joinpath(@__DIR__, "results/econ/electricity/HCOO_elec.csv")
    CO_elec_filepath = joinpath(@__DIR__, "results/econ/electricity/CO_elec.csv")

    C2H4_prf_filepath = joinpath(@__DIR__, "results/econ/prf/C2H4_prf.csv")
    HCOO_prf_filepath = joinpath(@__DIR__, "results/econ/prf/HCOO_prf.csv")
    CO_prf_filepath = joinpath(@__DIR__, "results/econ/prf/CO_prf.csv")

    C2H4_water_rsrc_filepath = joinpath(@__DIR__, "results/resource/water/C2H4_water_rsrc.csv")
    HCOO_water_rsrc_filepath = joinpath(@__DIR__, "results/resource/water/HCOO_water_rsrc.csv")
    CO_water_rsrc_filepath = joinpath(@__DIR__, "results/resource/water/CO_water_rsrc.csv")

    C2H4_power_filepath = joinpath(@__DIR__, "results/resource/power/C2H4_power.csv")
    HCOO_power_filepath = joinpath(@__DIR__, "results/resource/power/HCOO_power.csv")
    CO_power_filepath = joinpath(@__DIR__, "results/resource/power/CO_power.csv")

    C2H4_electrolyzer_filepath = joinpath(@__DIR__, "results/econ/electrolyzer/C2H4_electrolyzer.csv")
    HCOO_electrolyzer_filepath = joinpath(@__DIR__, "results/econ/electrolyzer/HCOO_electrolyzer.csv")
    CO_electrolyzer_filepath = joinpath(@__DIR__, "results/econ/electrolyzer/CO_electrolyzer.csv")

    C2H4_separations_filepath = joinpath(@__DIR__, "results/econ/separations/C2H4_separations.csv")
    HCOO_separations_filepath = joinpath(@__DIR__, "results/econ/separations/HCOO_separations.csv")
    CO_separations_filepath = joinpath(@__DIR__, "results/econ/separations/CO_separations.csv")

    C2H4_battery_filepath = joinpath(@__DIR__, "results/econ/battery/C2H4_battery.csv")
    HCOO_battery_filepath = joinpath(@__DIR__, "results/econ/battery/HCOO_battery.csv")
    CO_battery_filepath = joinpath(@__DIR__, "results/econ/battery/CO_battery.csv")

    C2H4_maintenance_filepath = joinpath(@__DIR__, "results/econ/maintenance/C2H4_maintenance.csv")
    HCOO_maintenance_filepath = joinpath(@__DIR__, "results/econ/maintenance/HCOO_maintenance.csv")
    CO_maintenance_filepath = joinpath(@__DIR__, "results/econ/maintenance/CO_maintenance.csv")

    C2H4_water_filepath = joinpath(@__DIR__, "results/econ/water/C2H4_water.csv")
    HCOO_water_filepath = joinpath(@__DIR__, "results/econ/water/HCOO_water.csv")
    CO_water_filepath = joinpath(@__DIR__, "results/econ/water/CO_water.csv")

    # Initialize flow output array
    flow_results = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    # Populate each output array with JuMP varaible value found in optimizer
    for i = 1:num_source_pts
        for j = 1:num_sink_pts
            for k in 1:num_species
                flow_results[i, j, k] = JuMP.value(flow[i, j, k])
            end
        end
    end

    # Initialize NPV output arrays for each molecule
    results_NPV = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    # Initialize electricity cost output arrays for each molecule
    results_elec = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    # Initialize electricity cost output arrays for each molecule
    results_trans = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    # Initialize electricity cost output arrays for each molecule
    results_ccu = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    # Initialize electricity cost output arrays for each molecule
    results_tea = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    # Initialize electricity cost output arrays for each molecule
    results_prf = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    #Initialize water consumption output arrays for each molecule
    results_water = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    #Initialize power consumption output arrays for each molecule
    results_power = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    #Initialize electrolyzer cost output arrays for each molecule
    results_electrolyzer = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    #Initialize separations cost output arrays for each molecule
    results_separations = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    #Initialize battery cost output arrays for each molecule
    results_battery = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    #Initialize electrolyzer cost output arrays for each molecule
    results_maintenance = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    #Initialize water cost output arrays for each molecule
    results_water_rsrc = zeros(Float64, num_source_pts, num_sink_pts, num_species)

    # Populate each output array with JuMP varaible value found in optimizer
    for i = 1:num_source_pts
        for j = 1:num_sink_pts
            for k = 1:num_species
                results_NPV[i, j, k] = sum(
                (
                - (D[i, j, k] * transport_A + transport_B) * flow_results[i, j, k] # Transportation cost
                -  electrolyzer_TEA(
                    flow_results[i, j, k],
                    product_phase[k],
                    product_mw[k],
                    product_density[k],
                    electrons_product[k],
                    electrons_co2[k],
                    system_potential[k],
                    selectivity[k],
                    hours_year,
                    current_density[k],
                    single_pass_conversion[k],
                    electrolyzer_cost,
                    electricity_price[i],
                    water_price,
                    water_rate[k],
                    crf)  # TEA system cost
                - flow_results[i, j, k] * carbons[k] * C_mw / product_mw[k] * CO2_cap_cost[i] # Carbon capture cost
                - flow_results[i, j, k] / hours_year / 3600 * 1e3 / product_mw[k] * electrons_product[k] * F * system_potential[k] / selectivity[k] / 1e3 * hours_year * electricity_price[i] # Electrochemical electricity price
                + flow_results[i, j, k] * value[k] # Product value
                )
                )

                if flow_results[i, j, k] != 0
                    results_electrolyzer[i, j, k],
                    results_separations[i, j, k],
                    results_battery[i, j, k],
                    results_maintenance[i, j, k],
                    results_water[i, j, k] =
                    electrolyzer_TEA_array(
                        flow_results[i, j, k],
                        product_phase[k],
                        product_mw[k],
                        product_density[k],
                        electrons_product[k],
                        electrons_co2[k],
                        system_potential[k],
                        selectivity[k],
                        hours_year,
                        current_density[k],
                        single_pass_conversion[k],
                        electrolyzer_cost,
                        electricity_price[i],
                        water_price,
                        water_rate[k],
                        crf)  # TEA system cost
                    results_trans[i, j, k] = (D[i, j, k] * transport_A + transport_B) * flow_results[i, j, k] # Transportation cost
                    results_tea[i, j, k] =  electrolyzer_TEA(
                        flow_results[i, j, k],
                        product_phase[k],
                        product_mw[k],
                        product_density[k],
                        electrons_product[k],
                        electrons_co2[k],
                        system_potential[k],
                        selectivity[k],
                        hours_year,
                        current_density[k],
                        single_pass_conversion[k],
                        electrolyzer_cost,
                        electricity_price[i],
                        water_price,
                        water_rate[k],
                        crf)  # TEA system cost
                    results_ccu[i, j, k] =  flow_results[i, j, k] * carbons[k] * C_mw / product_mw[k] * CO2_cap_cost[i] # Carbon capture cost
                    results_elec[i, j, k] =  flow_results[i, j, k] / hours_year / 3600 * 1e3 / product_mw[k] * electrons_product[k] * F * system_potential[k] / selectivity[k] / 1e3 * hours_year * electricity_price[i] # Electrochemical electricity price
                    results_prf[i, j, k] = flow_results[i, j, k] * value[k] # Product value
                    results_water_rsrc[i, j, k] = flow_results[i, j, k] / hours_year / 3600 * 1e3 / product_mw[k] * electrons_product[k] * F / selectivity[k] / H2O_electrons / F * H2O_mw / 1e3 * hours_year * 3600
                    results_power[i, j, k] = flow_results[i, j, k] / hours_year / 3600 * 1e3 / product_mw[k] * electrons_product[k] * F * system_potential[k] / selectivity[k] / 1e3 #kW
                end
            end
        end
    end

    # Write sink/source data to CSV files for data visualization
    # Useful for runs where points are excluded from optimization routine
    source_index_path = joinpath(@__DIR__, "results/index/source_index.csv")
    sink_index_path = joinpath(@__DIR__, "results/index/sink_index.csv")

    CSV.write(source_index_path, DataFrame(source_pts))
    CSV.write(sink_index_path, DataFrame(sink_pts))

    # Write output flow files to CSV
    CSV.write(C2H4_result_filepath, DataFrame(flow_results[:, :, 1]))
    CSV.write(HCOO_result_filepath, DataFrame(flow_results[:, :, 2]))
    CSV.write(CO_result_filepath, DataFrame(flow_results[:, :, 3]))

    # write output NPV dataframes to CSV
    CSV.write(C2H4_NPV_filepath, DataFrame(results_NPV[:, :, 1]))
    CSV.write(HCOO_NPV_filepath, DataFrame(results_NPV[:, :, 2]))
    CSV.write(CO_NPV_filepath, DataFrame(results_NPV[:, :, 3]))

    # write output electricity cost dataframes to CSV
    CSV.write(C2H4_trans_filepath, DataFrame(results_trans[:, :, 1]))
    CSV.write(HCOO_trans_filepath, DataFrame(results_trans[:, :, 2]))
    CSV.write(CO_trans_filepath, DataFrame(results_trans[:, :, 3]))

    # write output electricity cost dataframes to CSV
    CSV.write(C2H4_tea_filepath, DataFrame(results_tea[:, :, 1]))
    CSV.write(HCOO_tea_filepath, DataFrame(results_tea[:, :, 2]))
    CSV.write(CO_tea_filepath, DataFrame(results_tea[:, :, 3]))

    # write output electricity cost dataframes to CSV
    CSV.write(C2H4_ccu_filepath, DataFrame(results_ccu[:, :, 1]))
    CSV.write(HCOO_ccu_filepath, DataFrame(results_ccu[:, :, 2]))
    CSV.write(CO_ccu_filepath, DataFrame(results_ccu[:, :, 3]))

    # write output electricity cost dataframes to CSV
    CSV.write(C2H4_elec_filepath, DataFrame(results_elec[:, :, 1]))
    CSV.write(HCOO_elec_filepath, DataFrame(results_elec[:, :, 2]))
    CSV.write(CO_elec_filepath, DataFrame(results_elec[:, :, 3]))

    # write output electricity cost dataframes to CSV
    CSV.write(C2H4_prf_filepath, DataFrame(results_prf[:, :, 1]))
    CSV.write(HCOO_prf_filepath, DataFrame(results_prf[:, :, 2]))
    CSV.write(CO_prf_filepath, DataFrame(results_prf[:, :, 3]))

    # write output water consumption dataframes to CSV
    CSV.write(C2H4_water_rsrc_filepath, DataFrame(results_water_rsrc[:, :, 1]))
    CSV.write(HCOO_water_rsrc_filepath, DataFrame(results_water_rsrc[:, :, 2]))
    CSV.write(CO_water_rsrc_filepath, DataFrame(results_water_rsrc[:, :, 3]))

    # write output power consumption dataframes to CSV
    CSV.write(C2H4_power_filepath, DataFrame(results_power[:, :, 1]))
    CSV.write(HCOO_power_filepath, DataFrame(results_power[:, :, 2]))
    CSV.write(CO_power_filepath, DataFrame(results_power[:, :, 3]))

    # write output electrolyzer cost dataframes to CSV
    CSV.write(C2H4_electrolyzer_filepath, DataFrame(results_electrolyzer[:, :, 1]))
    CSV.write(HCOO_electrolyzer_filepath, DataFrame(results_electrolyzer[:, :, 2]))
    CSV.write(CO_electrolyzer_filepath, DataFrame(results_electrolyzer[:, :, 3]))

    # write output separations cost dataframes to CSV
    CSV.write(C2H4_separations_filepath, DataFrame(results_separations[:, :, 1]))
    CSV.write(HCOO_separations_filepath, DataFrame(results_separations[:, :, 2]))
    CSV.write(CO_separations_filepath, DataFrame(results_separations[:, :, 3]))

    # write output battery cost dataframes to CSV
    CSV.write(C2H4_battery_filepath, DataFrame(results_battery[:, :, 1]))
    CSV.write(HCOO_battery_filepath, DataFrame(results_battery[:, :, 2]))
    CSV.write(CO_battery_filepath, DataFrame(results_battery[:, :, 3]))

    # write output maintenance cost dataframes to CSV
    CSV.write(C2H4_maintenance_filepath, DataFrame(results_maintenance[:, :, 1]))
    CSV.write(HCOO_maintenance_filepath, DataFrame(results_maintenance[:, :, 2]))
    CSV.write(CO_maintenance_filepath, DataFrame(results_maintenance[:, :, 3]))

    # write output water cost dataframes to CSV
    CSV.write(C2H4_water_filepath, DataFrame(results_water[:, :, 1]))
    CSV.write(HCOO_water_filepath, DataFrame(results_water[:, :, 2]))
    CSV.write(CO_water_filepath, DataFrame(results_water[:, :, 3]))

    println("Optimization Routine Completed!")

end

@time co2_optimizer()


text:
function double_iex_multinomial_edit_accept_reject!(
    S_curr::InteractionSequence{Int},
    S_prop::InteractionSequence{Int},
    posterior::SimPosterior,
    γ_curr::Float64,
    mcmc::SimIexInsertDeleteProportional,
    P::CumCondProbMatrix,
    aux_data::InteractionSequenceSample{Int},
    suff_stat_curr::Float64,
    aux_init_at_prev::Bool
    )

    N = length(S_curr)
    n_tot = sum(length, S_curr)
    m_tot = n_tot
    dist = posterior.dist
    V = posterior.V
    K_inner, K_outer = (posterior.K_inner, posterior.K_outer)
    K_in_lb, K_in_ub = (K_inner.l, K_inner.u)
    data = posterior.data
    γ_prior = posterior.S_prior.γ
    mode_prior = posterior.S_prior.mode

    aux_mcmc = mcmc.aux_mcmc

    δ = rand(1:mcmc.ν_ed)  # Number of edits to enact
    rem_edits = δ # Remaining edits to allocate
    j = 0 # Keeps track how many interaction have been edited
    Z = 1
    log_ratio = 0.0

    for (i,(I_curr,I_prop)) in enumerate(zip(S_curr, S_prop))

        # If at end we just assign all remaining edits to final interaction
        n = length(I_curr)
        if i == N
            δ_tmp = rem_edits
        # Otherwise we sample the number of edits via rescaled Binomial
        else
            log_p_tmp = log(n) - log(n_tot)
            p = exp(log_p_tmp - log(Z))
            δ_tmp = rand(Binomial(rem_edits, p)) # Number of edits to perform on ith interaction
            Z -= exp(log_p_tmp)
        end

        # println("   Index $i getting $δ_tmp edits...")
        # If we sampled zero edits we skip to next iteration
        if δ_tmp == 0
            continue
        else
            j += 1 # Increment j
            # Make edits ....
            @inbounds n = length(I_curr)
            d = rand(0:min(n-K_in_lb, δ_tmp))
            m = n + δ_tmp - 2*d
            m_tot += δ_tmp - 2*d

            # Catch invalid proposals
            if (m > K_in_ub)
                # Here we just reject the proposal
                for (I_prop,I_curr) in zip(S_prop,S_curr)
                    @inbounds copy!(I_prop, I_curr)
                end
                return 0, suff_stat_curr
            end

            ind_del = view(mcmc.ind_del, 1:d)
            ind_add = view(mcmc.ind_add, 1:(δ_tmp-d))
            vals_del = view(mcmc.vals, 1:d)

            # Sample indexing info and new entries (all in-place)
            StatsBase.seqsample_a!(1:n, ind_del)
            StatsBase.seqsample_a!(1:m, ind_add)

            # *** HERE IS DIFFERENT FROM MODEL SAMPLER ***
            # The delete_insert_informed() function does the sampling + editing
            @inbounds log_ratio += delete_insert_informed!(
                I_prop,
                ind_del, ind_add, vals_del,
                P)

            @inbounds mcmc.ind_update[j] = i # Store which interaction was updated

            # Add to log_ratio
            # log_prod_term += log(b - a + 1) - log(ub(m, δ_tmp) - lb(m, δ_tmp, model) +1)
            log_ratio += log(min(n-K_in_lb, δ_tmp)+1) - log(min(m-K_in_lb, δ_tmp)+1) + log(m) - log(n)

        end

        # Update rem_edits
        rem_edits -= δ_tmp

        # If no more left terminate
        if rem_edits == 0
            break
        end

    end

    log_ratio += N*(log(n_tot) - log(m_tot))

    aux_model = SIM(
        S_prop, γ_curr,
        dist,
        V,
        K_inner,
        K_outer
    )

    if aux_init_at_prev
        tmp = deepcopy(aux_data[1])
        draw_sample!(aux_data, aux_mcmc, aux_model, init=tmp)
    else
        draw_sample!(aux_data, aux_mcmc, aux_model)
    end

    aux_log_lik_ratio = -γ_curr * (
        mapreduce(x -> dist(x, S_curr), + , aux_data)
        - mapreduce(x -> dist(x, S_prop), +, aux_data)
    )

    suff_stat_prop = mapreduce(x -> dist(x, S_prop), + , data)
    log_lik_ratio = -γ_curr * (
        suff_stat_prop - suff_stat_curr
    )

    log_prior_ratio = -γ_prior * (
        dist(S_prop, mode_prior) - dist(S_curr, mode_prior)
    )

    log_multinom_term = log_multinomial_ratio(S_curr, S_prop)

    # Log acceptance probability
    log_α = log_lik_ratio + log_prior_ratio + aux_log_lik_ratio + log_ratio + log_multinom_term
    # println("Edit alloc proposing...")
    # println("$(S_prop)")
    # Accept-reject step. Use info in mcmc.ind_update to know which interaction are to be copied over
    if log(rand()) < log_α
        for i in view(mcmc.ind_update, 1:j)
            @inbounds copy!(S_curr[i], S_prop[i])
        end
        return 1, suff_stat_prop
    else
        for i in view(mcmc.ind_update, 1:j)
            @inbounds copy!(S_prop[i], S_curr[i])
        end
        return 0, suff_stat_curr
    end
end


text:
using LinearAlgebra, ForwardDiff, Optim
using DataFrames, Plots, Distributions, Statistics

ProjDir = @__DIR__
cd(ProjDir)

include("./dynamic_pricing/price_demand_models.jl")
include("./demand_curve/select_optimal_demand_curve.jl")
include("./intervals/intervals.jl")
include("./utils/logx.jl")
include("./actual_demand/actual_demand_curve.jl")

# Generate demand hypothesis

a_range = range(-1.2, -0.6, length=4)
b_range = range(60, 70, length=4)
h_vec = generate_demand_hypothesis(a_range, b_range)
println("Price demand hypotheses (see also fig_01.png):")
display(h_vec)
println()

plot(xlim=(0, 60), ylim=(0, 80),
    xlabel="Price", ylabel="Demand",
    title="Price Demand hypotheses with optimal price points."
)
x = range(0, 60, length=101)
for i in 1:size(h_vec, 1)
	plot!(x, h_vec[i, :d].(x), color=:red, leg=false)
end
scatter!(h_vec[:, :p_opt], h_vec[:, :d_opt], m=(:cross))
savefig(joinpath(ProjDir, "fig_01.png"))

# Emperical demand curve

println()
curve = emperical_demand_curve(30, 60, 20)
println("Emperical price demand curve:")
display(curve)

# Intervals for price updates

T = 24 * 1                       # time step is one hour, flash offering for 1 day
m = 4                            # not more than 4 price updates

t_mask = intervals(m, T, 2)      # intervals for price lvels
p = h_vec[1, :p_opt]             # initial price

println("\nSelected intervals:")
display(t_mask')
println("\nInitial price = $p\n")

# Do a first iteration in selecting demand hypothesis

(price, row, hist_d) = select_optimal_demand_curve(t_mask, p)
println()
display((price, row))
println()
display((hist_d))

plot(xlim=(0, 60), ylim=(0, 80), xlabel="Price", ylabel="Demand")
plot!(curve[:, :prices], curve[:, :demands],
    xlabel="price", ylabel="demand", line=(:dot, :gray))
x = range(0, 60, length=101)
plot!(x, h_vec[row, :d].(x), color=:red, leg=false)
savefig(joinpath(ProjDir, "fig_02.png"))

#=
colors = ['#8EA604', '#F5BB00', '#EC9F05', '#D76A03']

def visualize_snapshot(t):
    fig.clear()

    plt.subplot(2, 1, 1)
    plt.xlabel('Price')
    plt.ylabel('Demand')
    plt.title('Price-demand curve fitting')
    plt.plot(curve[0], curve[1], 'k:')
    for i in range(0, t):
        plt.plot(prices, list(map(np.array(history)[i, 2], prices)), 'k-')
    plt.plot(np.array(history)[0:t,0], np.array(history)[0:t, 1], 'rx')

    plt.subplot(2, 1, 2)
    plt.xlabel('Time')
    plt.ylabel('Demand/price')
    plt.title('Realized demand and price')
    plt.plot(range(0, t+1), np.array(history)[0:t+1,0], 'r-')
    bars = plt.bar(range(0, T-1), np.pad(np.array(history)[0:t+1, 1], (0, T-2-t), 'constant'), 0.35)
    for i in range(0, t+1):
        bars[i].set_color(colors[t_mask[i]])

fig = plt.figure(figsize=(10, 8))
plt.subplots_adjust(hspace=0.5)
visualize_snapshot(T-2)               # visualize the state in the end of the simulation
plt.show()
=#


text:
### Definition of Binary trees
# This was mostly copied over from the AbstractTrees example, credits are due
#
# Copied from AbstractTrees.jl Jan. 2021

## most of this is taken from the AbstractTrees example
mutable struct BinaryNode{T}
  data::Union{T, Nothing}
  left::Union{BinaryNode{T}, Nothing}
  right::Union{BinaryNode{T}, Nothing}

  # Root constructor
  BinaryNode{T}() where T = new{T}(nothing, nothing, nothing)
  BinaryNode{T}(data) where T = new{T}(data, nothing, nothing)
  BinaryNode{T}(data, left::BinaryNode{T}, right::BinaryNode{T}) where T = new{T}(data, left, right)
  BinaryNode{T}(left::BinaryNode{T}, right::BinaryNode{T}) where T = new{T}(nothing, left, right)
end
BinaryNode(data) = BinaryNode{typeof(data)}(data)
BinaryNode(left::BinaryNode{T}, right::BinaryNode{T}) where T = BinaryNode{T}(left, right)
BinaryNode(data, left::BinaryNode, right::BinaryNode) = BinaryNode{typeof(data)}(data, left, right)

# some utility functions
isleaf(node::BinaryNode) = isnothing(node.left) && isnothing(node.right)
isbranch(node::BinaryNode) = !isnothing(node.left) && !isnothing(node.right)

# counts the maximum depth of the tree
depth(node::BinaryNode) = _depth(node, 1)
function _depth(node::BinaryNode, level)
  if isleaf(node)
    return level
  elseif !isnothing(node.left)
    if !isnothing(node.right)
      return max(_depth(node.left, level+1), _depth(node.right, level+1))
    else
      return _depth(node.left, level+1)
    end
  else
    return _depth(node.right, level+1)
  end
end

function nleaves(node::BinaryNode)
  if isleaf(node)
    return 1
  elseif !isnothing(node.left)
    if !isnothing(node.right)
      return nleaves(node.left) + nleaves(node.right)
    else
      return nleaves(node.left)
    end
  else
    return nleaves(node.right)
  end
end

function AbstractTrees.children(node::BinaryNode)
  if !isnothing(node.left)
    if !isnothing(node.right)
      return (node.left, node.right)
    end
    return (node.left,)
  end
  !isnothing(node.right) && return (node.right,)
  return ()
end

# routine to check equivalence of tree structures
function compatible(node1::BinaryNode, node2::BinaryNode)
  if isnothing(node1.left) ⊻ isnothing(node2.left) return false end
  if isnothing(node1.right) ⊻ isnothing(node2.right) return false end
  if !isnothing(node1.left) && !compatible(node1.left, node2.left) return false end
  if !isnothing(node1.right) && !compatible(node1.right, node2.right) return false end
  return true
end

==(node1::BinaryNode, node2::BinaryNode) = (node1.data == node2.data) && (node1.left == node2.left) && (node1.right == node2.right)
!=(node1::BinaryNode, node2::BinaryNode) = (node1.data != node2.data) || (node1.left != node2.left) || (node1.right != node2.right)

# remove leaves on the bottom level
function prune_leaves!(node::BinaryNode)
  if isleaf(node.left) && isleaf(node.right)
    node.left = nothing
    node.right = nothing
  else
    if !isleaf(node.left)
      node.left = prune_leaves!(node.left)
    end
    if !isleaf(node.right)
      node.right = prune_leaves!(node.right)
    end
  end
  return node
end

## Optional enhancements
# These next two definitions allow inference of the item type in iteration.
Base.eltype(::Type{BinaryNode{T}}) where T = T
Base.eltype(::Type{<:TreeIterator{BinaryNode{T}}}) where T = BinaryNode{T}
Base.IteratorEltype(::Type{<:TreeIterator{BinaryNode{T}}}) where T = Base.HasEltype()

# routines for displaying trees
AbstractTrees.printnode(io::IO, node::BinaryNode) = print(io, node.data)
Base.show(io::IO, node::BinaryNode) = print(io, "BinaryNode{$(eltype(node))}")

# Implement iteration over the immediate children of a node
function Base.iterate(node::BinaryNode)
  isdefined(node, :left) && return (node.left, false)
  isdefined(node, :right) && return (node.right, true)
  return nothing
end
function Base.iterate(node::BinaryNode, state::Bool)
  state && return nothing
  isdefined(node, :right) && return (node.right, true)
  return nothing
end
Base.IteratorSize(::Type{BinaryNode{T}}) where T = Base.SizeUnknown()
#Base.eltype(::Type{BinaryNode{T}}) where T = BinaryNode{T}

# ## Things we need to define to leverage the native iterator over children
# ## for the purposes of AbstractTrees.
# # Set the traits of this kind of tree
# Base.eltype(::Type{<:TreeIterator{BinaryNode{T}}}) where T = BinaryNode{T}
# Base.IteratorEltype(::Type{<:TreeIterator{BinaryNode{T}}}) where T = Base.HasEltype()
# AbstractTrees.parentlinks(::Type{BinaryNode{T}}) where T = AbstractTrees.StoredParents()
# AbstractTrees.siblinglinks(::Type{BinaryNode{T}}) where T = AbstractTrees.StoredSiblings()
# # Use the native iteration for the children
# AbstractTrees.children(node::BinaryNode) = node

# Base.parent(root::BinaryNode, node::BinaryNode) = isdefined(node, :parent) ? node.parent : nothing

# function AbstractTrees.nextsibling(tree::BinaryNode, child::BinaryNode)
#   isdefined(child, :parent) || return nothing
#   p = child.parent
#   if isdefined(p, :right)
#       child === p.right && return nothing
#       return p.right
#   end
#   return nothing
# end


text:
using LightGraphs
using Nauty
using Test

@testset "Nauty.jl" begin
    # Petersen graph
    g = DenseNautyGraph(10)
    for i=1:5
        add_edge!(g,i,1+(i%5))
        add_edge!(g,i,i+5)
        add_edge!(g,i+5,6+((i+2)%5))
    end
    @test nv(g) == 10
    @test ne(g) == 15
    dict = nauty(g,automgroup=true)
    @test dict[:grpsize] == 120
end


text:
using PlotlyJS

function logistic(r,x)
    r*x*(1.0-x)
end

function logisticmap(;x0=0.21,r=4,n=10)
    results = [(x0,logistic(r,x0))];
    for i in 1:(n-1)
      xold = results[end][2]
      xnew = logistic(r,xold)
      append!(results,[(xold,xold)])
      append!(results,[(xold,xnew)])
  end
  results
end

simres = logisticmap(x0=0.21,r=4,n=100)

xvals = [x for (x, y) in simres]
yvals = [y for (x, y) in simres]

trace = scatter(;x=xvals, y=yvals, mode="lines")
plot(trace)

#this part is for generating many different images for an animation

for x0 in 0:0.01:1.0
    simres = logisticmap(x0=x0,r=4,n=100)
    xvals = [x for (x, y) in simres]
    yvals = [y for (x, y) in simres]
    trace = scatter(;x=xvals, y=yvals, mode="lines")
    plot(trace)
end

#xvals = []
#yvals = []

#for (x, y) in simres
    #append!(xvals, x)
    #append!(yvals, y)
#end


text:
# A test case for wide lines and mitering at joints

@block Contributors ["GLMakie backend tests"] begin
    @cell "Miter Joints for line rendering" [lines] begin
        scene = Scene()

        r = 4
        sep = 4*r
        scatter!(scene, (sep+2*r)*[-1,-1,1,1], (sep+2*r)*[-1,1,-1,1])

        for i=-1:1
            for j=-1:1
                angle = pi/2 + pi/4*i
                x = r*[-cos(angle/2),0,-cos(angle/2)]
                y = r*[-sin(angle/2),0,sin(angle/2)]

                linewidth = 40 * 2.0^j
                lines!(scene, x .+ sep*i, y .+ sep*j, color=RGBAf0(0,0,0,0.5), linewidth=linewidth)
                lines!(scene, x .+ sep*i, y .+ sep*j, color=:red)
            end
        end
        scene
    end
end


text:
# Test routine to compare performance for a number of SOCP Lasso problems
workspace()
include("../../../src/Solver.jl")
include("../solverComparison/Compare.jl")

using OSQP, OSSDP, Base.Test, Compare, JuMP, Mosek

rng = Random.MersenneTwister(12345)

nn = 25
timestamp = Dates.format(now(), "yyddmm_HH-MM")
dataFile = "./SC_" * timestamp * "LassoSOCP.jld"
problemType = "SOCP-Lasso"

sr1 = SolverResult(nn, problemType,"COSMO",timestamp,0,false)
sr2 = SolverResult(nn, problemType,"COSMO",timestamp,0,true)
resData = [sr1;sr2]
ws = 0
for iii =1:1:nn
  # generate problem data
  n = 8
  m = 50*n
  F = rand(m,n)

  vtrue = sprand(n,1,0.1 )
  noise = 0.1*rand(m,1)
  b = F*vtrue + noise
  μMax = norm(F'*b,Inf)
  μ = 0.1*μMax

  pDims = [m;n;m*n]

  # define lasso problem as SOCP

  # create augmented matrices
  Aa = sparse([1 zeros(1,2*n+1) 1 zeros(1,m);
        -1 zeros(1,2*n) 1 zeros(1,m+1);
        zeros(m,1) -2*F zeros(m,n+2) Matrix(1.0I,m,m);
        zeros(n,1) eye(n) -eye(n) zeros(n,m+2);
        zeros(n,1) -eye(n) -eye(n) zeros(n,m+2);
       zeros(1,2*n+1) -1 zeros(1,m+1);
       zeros(1,2*n+2) -1 zeros(1,m);
       zeros(m,2n+3) -eye(m)])

  ba = [1;1;-2*b;zeros(2*n+m+2)]
  q = [1;zeros(n);μ*ones(n,1);zeros(m+2,1)]
  P = spzeros(length(q),length(q))

  # define cone membership
  K = Cone(2+m,2*n,[m+2],[])

setOFF = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,check_termination=1,scaling = 0,eps_abs = 1e-3,eps_rel=1e-3,adaptive_rho=true)
setON = OSSDPSettings(rho=0.1,sigma=1e-6,alpha=1.6,max_iter=1500,verbose=false,check_termination=1,scaling = 10,eps_abs = 1e-3,eps_rel=1e-3,adaptive_rho=true)

  # Solve with OSSDP
  resOSSDP_unscaled,nothing = OSSDP.solve(P,q,Aa,ba,K,setOFF);
  print("\n.")
  resOSSDP_scaled,nothing = OSSDP.solve(P,q,Aa,ba,K,setON);
  print(".")

  model = Model(solver=MosekSolver())
  @variable(model, x[1:n])
  @variable(model, t[1:n])
  @variable(model, y)
  @variable(model, w[1:m+1])
  @objective(model, Min, y + μ*ones(n)'*t)
  @constraint(model, -t .<= x)
  @constraint(model, x .<= t)
  @constraint(model, w[1] == 1-y)
  @constraint(model, w[2:m+1] .== 2*(F*x-b))
  @constraint(model, norm(w) <= 1+y)
  status = JuMP.solve(model)
  println("Objective value: ", getobjectivevalue(model))
  println("x = ", getvalue(x))

  resArray = [resOSSDP_unscaled;resOSSDP_scaled]
  updateResults!(dataFile,resData,resArray,pDims,"SOCP-$(iii)",0.,true)
  printStatus(iii,nn,"SOCP-$(iii)",resData)
end


text:
using Test

@info "CANalyze.Decode tests..."
@testset "bit" begin
    import CANalyze.Utils
    import CANalyze.Frames
    import CANalyze.Signals
    import CANalyze.Messages
    import CANalyze.Decode

    @testset "bit_1" begin
        for start=0:63
            m = Utils.mask(UInt64, 1, start)
            signal = Signals.Bit(start=start)
            frame = Frames.CANFrame(0x1FF, Utils.to_bytes(m))
            @test Decode.decode(signal, frame)
        end
    end

    @testset "bit_2" begin
        for start=1:63
            m = Utils.mask(UInt64, 1, start)
            signal = Signals.Bit(start=start-1)
            frame = Frames.CANFrame(0x1FF, Utils.to_bytes(m))
            @test !Decode.decode(signal, frame)
        end
    end

    @testset "bit_3" begin
        signal = Signals.Bit(start=8)
        frame = Frames.CANFrame(0x1FF, 1)
        @test_throws DomainError Decode.decode(signal, frame)
    end

    @testset "bit_4" begin
        signal = Signals.Bit(start=8)
        frame = Frames.CANFrame(0x1FF, 1)
        @test Decode.decode(signal, frame, nothing) == nothing
    end
end

@testset "unsigned" begin
    import CANalyze.Utils
    import CANalyze.Frames
    import CANalyze.Signals
    import CANalyze.Messages
    import CANalyze.Decode

    @testset "unsigned_1" begin
        for start=0:63
            for len=1:(64-start)
                m = Utils.mask(UInt64, len, start)
                signal = Signals.Unsigned{Float64}(start=start, length=len, factor=2.0,
                                                   offset=1337,
                                                   byte_order=:little_endian)
                frame = Frames.CANFrame(0x1FF, Utils.to_bytes(m))
                decode = Decode.decode(signal, frame)
                value = Utils.mask(UInt64, len) * Signals.factor(signal) + Signals.offset(signal)
                @test decode == value
            end
        end
    end

    @testset "unsigned_2" begin
        signal = Signals.Unsigned{Float64}(start=7, length=8, factor=2.0, offset=1337,
                                           byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, [i for i=1:8])
        decode = Decode.decode(signal, frame)
        value = 1 * Signals.factor(signal) + Signals.offset(signal)
        @test decode == value
    end

    @testset "unsigned_3" begin
        signal = Signals.Unsigned{Float64}(start=7, length=16, factor=2.0, offset=1337,
                                           byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0xAB, 0xCD)
        decode = Decode.decode(signal, frame)
        value = 0xABCD * Signals.factor(signal) + Signals.offset(signal)
        @test decode == value
    end

    @testset "unsigned_4" begin
        signal = Signals.Unsigned{Float64}(start=7, length=24, factor=2.0, offset=1337,
                                           byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0xAB, 0xCD, 0xEF)
        decode = Decode.decode(signal, frame)
        value = 0xABCDEF * Signals.factor(signal) + Signals.offset(signal)
        @test decode == value
    end

    @testset "unsigned_5" begin
        signal = Signals.Unsigned{Float64}(start=8, length=1, factor=2.0, offset=1337,
                                           byte_order=:little_endian)
        frame = Frames.CANFrame(0x1FF, 0x01)
        @test_throws DomainError Decode.decode(signal, frame)
    end

    @testset "unsigned_6" begin
        signal = Signals.Unsigned{Float64}(start=6, length=8, factor=2.0, offset=1337,
                                           byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0x01)
        @test_throws DomainError Decode.decode(signal, frame)
    end
end

@testset "signed" begin
    import CANalyze.Utils
    import CANalyze.Frames
    import CANalyze.Signals
    import CANalyze.Messages
    import CANalyze.Decode
    using Random

    @testset "signed_1" begin
        for start=0:62
            for len=1:(64-start)
                m = Utils.mask(UInt64, len-1, start)
                signal = Signals.Signed{Float64}(start=start, length=len, factor=2.0,
                                                 offset=1337,
                                                 byte_order=:little_endian)
                frame = Frames.CANFrame(0x1FF, Utils.to_bytes(m))
                decode = Decode.decode(signal, frame)
                value = Utils.mask(UInt64, len-1) * Signals.factor(signal) + Signals.offset(signal)
                @test decode == value
            end
        end
    end

    @testset "signed_2" begin
        for start=0:63
            for len=1:(64-start)
                m = Utils.mask(UInt64, len, start)
                signal = Signals.Signed{Float64}(start=start, length=len, factor=2.0,
                                                 offset=1337,
                                                 byte_order=:little_endian)
                frame = Frames.CANFrame(0x1FF, Utils.to_bytes(m))
                decode = Decode.decode(signal, frame)
                value = Utils.mask(Int64, len) + ~Utils.mask(Int64, len)
                value = value * Signals.factor(signal) + Signals.offset(signal)
                @test decode == value
            end
        end
    end

    @testset "signed_3" begin
        for len=1:64
            for choice=1:64-len
                m = Utils.bit_mask(Int64, len-1, rand(0:(len-1), choice)...)
                signal = Signals.Signed{Float64}(start=0, length=len, factor=2.0,
                                                 offset=1337,
                                                 byte_order=:little_endian)
                frame = Frames.CANFrame(0x1FF, Utils.to_bytes(m))
                decode = Decode.decode(signal, frame)
                value = m + ~Utils.mask(Int64, len)
                value = value * Signals.factor(signal) + Signals.offset(signal)
                @test decode == value
            end
        end
    end

    @testset "signed_4" begin
        signal = Signals.Signed{Float64}(start=7, length=8, factor=set=1337,
                                         byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0xFE)
        decode = Decode.decode(signal, frame)
        value = -2 * Signals.factor(signal) + Signals.offset(signal)
        @test decode == value
    end

    @testset "signed_5" begin
        signal = Signals.Signed{Float64}(start=8, length=1, factor=2.0, offset=1337,
                                         byte_order=:little_endian)
        frame = Frames.CANFrame(0x1FF, 0x01)
        @test_throws DomainError Decode.decode(signal, frame)
    end

    @testset "signed_6" begin
        signal = Signals.Signed{Float64}(start=6, length=8, factor=2.0, offset=1337,
                                         byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0x01)
        @test_throws DomainError Decode.decode(signal, frame)
    end
end

@testset "float_signal" begin
    import CANalyze.Utils
    import CANalyze.Frames
    import CANalyze.Signals
    import CANalyze.Messages
    import CANalyze.Decode
    using Random

    @testset "float_signal_1" begin
        for T in [Float16, Float32, Float64]
            data = [i for i=0:(sizeof(T)-1)]
            signal = Signals.FloatSignal{T}(start=0, factor=2.0, offset=1337,
                                            byte_order=:little_endian)
            frame = Frames.CANFrame(0x1FF, data)
            decode = Decode.decode(signal, frame)
            value = reinterpret(T, data)[1] * Signals.factor(signal) + Signals.offset(signal)
            @test decode == value
        end
    end

    @testset "float_signal_2" begin
        for T in [Float16, Float32, Float64]
            data = [i for i=0:(sizeof(T)-1)]
            signal = Signals.FloatSignal{T}(start=7, factor=2.0, offset=1337,
                                            byte_order=:big_endian)
            frame = Frames.CANFrame(0x1FF, data)
            decode = Decode.decode(signal, frame)
            value = reinterpret(T, reverse(data))[1] * Signals.factor(signal)
            value += Signals.offset(signal)
            @test decode == value
        end
    end

    @testset "float_signal_3" begin
        signal = Signals.Float16Signal(start=1, byte_order=:little_endian)
        frame = Frames.CANFrame(0x1FF, 0x01, 0x02)
        @test_throws DomainError Decode.decode(signal, frame)
    end

    @testset "float_signal_4" begin
        signal = Signals.Float16Signal(start=6, byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0x01, 0x02)
        @test_throws DomainError Decode.decode(signal, frame)
    end

    @testset "float_signal_5" begin
        signal = Signals.Float32Signal(start=1, byte_order=:little_endian)
        frame = Frames.CANFrame(0x1FF, 0x01, 0x02, 0x03, 0x04)
        @test_throws DomainError Decode.decode(signal, frame)
    end

    @testset "float_signal_6" begin
        signal = Signals.Float32Signal(start=6, byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0x01, 0x02, 0x03, 0x04)
        @test_throws DomainError Decode.decode(signal, frame)
    end

    @testset "float_signal_7" begin
        signal = Signals.Float64Signal(start=1, byte_order=:little_endian)
        frame = Frames.CANFrame(0x1FF, 0:7)
        @test_throws DomainError Decode.decode(signal, frame)
    end

    @testset "float_signal_8" begin
        signal = Signals.Float64Signal(start=6, byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0:7)
        @test_throws DomainError Decode.decode(signal, frame)
    end
end

@testset "raw" begin
    import CANalyze.Utils
    import CANalyze.Frames
    import CANalyze.Signals
    import CANalyze.Messages
    import CANalyze.Decode

    @testset "raw_1" begin
        for start=0:63
            for len=1:(64-start)
                m = Utils.mask(UInt64, len, start)
                signal = Signals.Raw(start=start, length=len, byte_order=:little_endian)
                frame = Frames.CANFrame(0x1FF, Utils.to_bytes(m))
                decode = Decode.decode(signal, frame)
                value = Utils.mask(UInt64, len)
                @test decode == value
            end
        end
    end

    @testset "raw_2" begin
        signal = Signals.Raw(start=7, length=8, byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, [i for i=1:8])
        decode = Decode.decode(signal, frame)
        @test decode == 1
    end

    @testset "raw_3" begin
        signal = Signals.Raw(start=7, length=16, byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0xAB, 0xCD)
        decode = Decode.decode(signal, frame)
        @test decode == 0xABCD
    end

    @testset "raw_4" begin
        signal = Signals.Raw(start=7, length=24, byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0xAB, 0xCD, 0xEF)
        decode = Decode.decode(signal, frame)
        @test decode == 0xABCDEF
    end

    @testset "raw_5" begin
        signal = Signals.Raw(start=7, length=64, byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, [i for i=1:8])
        decode = Decode.decode(signal, frame)
        @test decode == 0x01_02_03_04_05_06_07_08
    end

    @testset "raw_6" begin
        signal = Signals.Raw(start=3, length=8, byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0x21, 0xAB)
        decode = Decode.decode(signal, frame)
        @test decode == 0x1A
    end

    @testset "raw_7" begin
        signal = Signals.Raw(start=3, length=16, byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0x21, 0xAB, 0xCD)
        decode = Decode.decode(signal, frame)
        @test decode == 0x1ABC
    end

    @testset "raw_8" begin
        signal = Signals.Raw(start=8, length=1, byte_order=:little_endian)
        frame = Frames.CANFrame(0x1FF, 0x01)
        @test_throws DomainError Decode.decode(signal, frame)
    end

    @testset "raw_9" begin
        signal = Signals.Raw(start=6, length=8, byte_order=:big_endian)
        frame = Frames.CANFrame(0x1FF, 0x01)
        @test_throws DomainError Decode.decode(signal, frame)
    end
end

@testset "named_signal" begin
    import CANalyze.Utils
    import CANalyze.Frames
    import CANalyze.Signals
    import CANalyze.Messages
    import CANalyze.Decode

    @testset "named_signal_1" begin
        signal = Signals.Signed{Float64}(start=0,
                                         length=16,
                                         factor=2.0,
                                         offset=1337,
                                         byte_order=:little_endian)
        named_signal = Signals.NamedSignal("SIG", nothing, nothing, signal)
        frame = Frames.CANFrame(0x1FF, 0x01, 0x02)
        @test Decode.decode(signal, frame) == Decode.decode(named_signal, frame)
    end

    @testset "named_signal_2" begin
        signal = Signals.Signed{Float64}(start=1,
                                         length=16,
                                         factor=2.0,
                                         offset=1337,
                                         byte_order=:little_endian)
        named_signal = Signals.NamedSignal("SIG", nothing, nothing, signal)
        frame = Frames.CANFrame(0x1FF, 0x01, 0x02)
        @test_throws DomainError Decode.decode(signal, frame)
        @test Decode.decode(named_signal, frame) == nothing
    end

    @testset "named_signal_3" begin
        signal = Signals.Signed{Float64}(start=1,
                                         length=16,
                                         factor=2.0,
                                         offset=1337,
                                         byte_order=:little_endian)
        named_signal = Signals.NamedSignal("SIG", nothing, 42.0, signal)
        frame = Frames.CANFrame(0x1FF, 0x01, 0x02)
        @test_throws DomainError Decode.decode(signal, frame)
        @test Decode.decode(named_signal, frame) == 42
    end
end

@testset "message" begin
    @testset "message_1" begin
        sig1 = Signals.Signed{Float64}(start=0, length=8, byte_order=:little_endian)
        sig2 = Signals.Signed{Float64}(start=8, length=8, byte_order=:little_endian)
        sig3 = Signals.Signed{Float64}(start=16, length=8, byte_order=:little_endian)

        named_signal_1 = Signals.NamedSignal("A", nothing, nothing, sig1)
        named_signal_2 = Signals.NamedSignal("B", nothing, nothing, sig2)
        named_signal_3 = Signals.NamedSignal("C", nothing, nothing, sig3)
        signals = [named_signal_1, named_signal_2, named_signal_3]

        frame = Frames.CANFrame(0x1FF, 0x01, 0x02, 0x03)
        m = Messages.Message(0x1FF, 8, "M", named_signal_1, named_signal_2, named_signal_3)

        value = Decode.decode(m, frame)
        for signal in signals
            @test value[Signals.name(signal)] == Decode.decode(signal, frame)
        end
    end
end


text:
### utils.jl
#
# Copyright (C) 2016, 2017 Mosè Giordano.
#
# Maintainer: Mosè Giordano <mose AT gnu DOT org>
# Keywords: periodogram, lomb scargle
#
# This file is a part of LombScargle.jl.
#
# License is BSD 3-clause "New" or "Revised".
#
### Commentary
#
# This file contains some utilities for the LombScargle.jl package.
#
### Code:

export power, freq, freqpower, findmaxpower, findmaxfreq, period, periodpower,
    findmaxperiod, prob, probinv, fap, fapinv

"""
    power(p::Periodogram)

Return the power vector of Lomb–Scargle periodogram `p`.
"""
power(p::Periodogram) = p.power

"""
    freq(p::Periodogram)

Return the frequency vector of Lomb–Scargle periodogram `p`.
"""
freq(p::Periodogram) = p.freq

"""
    freqpower(p::Periodogram)

Return the 2-tuple `(freq(p), power(p))`, where `freq(p)` and `power(p)` are the
frequency vector and the power vector of Lomb–Scargle periodogram `p`
respectively.
"""
freqpower(p::Periodogram) = (freq(p), power(p))

"""
    findmaxpower(p::Periodogram)

Return the highest power of the periodogram `p`.
"""
findmaxpower(p::Periodogram) = maximum(power(p))

_findmaxfreq(freq::AbstractVector{<:Real}, power::AbstractVector{<:Real}, threshold::Real) =
    freq[findall(x -> x >= threshold, power)]

"""
    findmaxfreq(p::Periodogram, [interval::AbstractVector{Real}], threshold::Real=findmaxpower(p))

Return the array of frequencies with the highest power in the periodogram `p`.
If a scalar real argument `threshold` is provided, return the frequencies with
power larger than or equal to `threshold`.  If you want to limit the search to a
narrower frequency range, pass as second argument a vector with the extrema of
the interval.
"""
findmaxfreq(p::Periodogram, threshold::Real=findmaxpower(p)) =
    _findmaxfreq(freqpower(p)..., threshold)

function findmaxfreq(p::Periodogram, interval::AbstractVector{<:Real}, threshold::Real=NaN)
    f = freq(p)
    lo, hi = extrema(interval)
    indices = findall(x -> lo <= x <= hi, f)
    pow = power(p)[indices]
    if isnan(threshold)
        threshold = maximum(pow)
    end
    return _findmaxfreq(f[indices], pow, threshold)
end

"""
    period(p::Periodogram)

Return the period vector of Lomb–Scargle periodogram `p`.  It is equal to `1 ./ freq(p)`.
"""
period(p::Periodogram) = 1 ./ freq(p)

"""
    periodpower(p::Periodogram)

Return the 2-tuple `(period(p), power(p))`, where `period(p)` and `power(p)` are
the period vector and the power vector of Lomb–Scargle periodogram `p`
respectively.
"""
periodpower(p::Periodogram) = (period(p), power(p))

"""
    findmaxperiod(p::Periodogram, [interval::AbstractVector{Real}], threshold::Real=findmaxpower(p))

Return the array of periods with the highest power in the periodogram `p`.  If a
scalar real argument `threshold` is provided, return the period with power
larger than or equal to `threshold`.  If you want to limit the search to a
narrower period range, pass as second argument a vector with the extrema of the
interval.
"""
findmaxperiod(p::Periodogram, threshold::Real=findmaxpower(p)) =
    1 ./ findmaxfreq(p, threshold)
findmaxperiod(p::Periodogram, interval::AbstractVector{<:Real},
              threshold::Real=NaN) =
                  1 ./ findmaxfreq(p, 1 ./ interval, threshold)

"""
    autofrequency(times::AbstractVector{Real};
                  samples_per_peak::Integer=5,
                  nyquist_factor::Integer=5,
                  minimum_frequency::Real=NaN,
                  maximum_frequency::Real=NaN)

Determine a suitable frequency grid for the given vector of `times`.

Optional keyword arguments are:

* `samples_per_peak`: the approximate number of desired samples across the
  typical peak
* `nyquist_factor`: the multiple of the average Nyquist frequency used to choose
  the maximum frequency if `maximum_frequency` is not provided
* `minimum_frequency`: if specified, then use this minimum frequency rather than
  one chosen based on the size of the baseline
* `maximum_frequency`: if specified, then use this maximum frequency rather than
  one chosen based on the average Nyquist frequency

This is based on prescription given at
<https://jakevdp.github.io/blog/2015/06/13/lomb-scargle-in-python/> and uses the
same keywords names adopted in Astropy.
"""
function autofrequency(times::AbstractVector{<:Real};
                       samples_per_peak::Integer=5,
                       nyquist_factor::Integer=5,
                       minimum_frequency::Real=NaN,
                       maximum_frequency::Real=NaN)
    T = maximum(times) - minimum(times)
    δf = inv(samples_per_peak * T)
    f_min = isfinite(minimum_frequency) ? minimum_frequency : (δf / 2)
    if isfinite(maximum_frequency)
        return f_min:δf:maximum_frequency
    else
        return f_min:δf:(nyquist_factor * length(times) / 2T)
    end
end

"""
    prob(P::Periodogram, pow::Real)

Return the probability that the periodogram power can exceed the value `pow`.

Its inverse is the `probinv` function.
"""
function prob(P::Periodogram, pow::Real)
    N = length(P.times)
    normalization = P.norm
    if normalization == :standard
        return (1 - pow)^((N - 3) / 2)
    elseif normalization == :Scargle
        return exp(-pow)
    elseif normalization == :HorneBaliunas
        return (1 - 2*pow/(N - 1))^((N - 3) / 2)
    elseif normalization == :Cumming
        return (1 + 2*pow/(N - 3))^((3 - N) / 2)
    else
        error("normalization \"", string(normalization), "\" not supported")
    end
end

"""
    probinv(P::Periodogram, prob::Real)

Return the power value of the periodogram power whose probability is `prob`.

This is the inverse of `prob` function.
"""
function probinv(P::Periodogram, prob::Real)
    N = length(P.times)
    normalization = P.norm
    if normalization == :standard
        return 1 - prob^(2/(N - 3))
    elseif normalization == :Scargle
        return -log(prob)
    elseif normalization == :HorneBaliunas
        return (N - 1) * (1 - prob ^ (2 / (N - 3))) / 2
    elseif normalization == :Cumming
        return (N - 3) * (prob ^ (2 / (3 - N)) - 1) / 2
    else
        error("normalization \"", string(normalization), "\" not supported")
    end
end

"""
    LombScargle.M(P::Periodogram)

Estimates the number of independent frequencies in the periodogram `P`.
"""
function M(P::Periodogram)
    tmin, tmax = extrema(P.times)
    fmin, fmax = extrema(P.freq)
    return (tmax - tmin)*(fmax - fmin)
end

"""
    fap(P::Periodogram, pow::Real)

Return the false-alarm probability for periodogram `P` and power value `pow`.

Its inverse is the [`fapinv`](@ref) function.
"""
fap(P::Periodogram, pow::Real) = 1 - (1 - prob(P, pow))^M(P)

"""
    fapinv(P::Periodogram, prob::Real)

Return the power value of the periodogram whose false-alarm probability is
`prob`.

This is the inverse of [`fap`](@ref) function.
"""
fapinv(P::Periodogram, prob::Real) = probinv(P, 1 - (1 - prob)^(inv(M(P))))

"""
    LombScargle.model(times::AbstractVector{Real},
                      signal::AbstractVector{R2},
                      [errors::AbstractVector{R3},]
                      frequency::Real,
                      [times_fit::AbstractVector{R4}];
                      center_data::Bool=true,
                      fit_mean::Bool=true)

Return the best fitting Lomb–Scargle model for the given signal at the given
frequency.

Mandatory arguments are:

* `times`: the observation times
* `signal`: the signal, sampled at `times` (must have the same length as
  `times`)
* `frequency`: the frequency at which to calculate the model

Optional arguments are:

* `errors`: the vector of uncertainties of the signal.  If provided, it must
  have the same length as `signal` and `times`, and be the third argument.  Like
  for [`lombscargle`](@ref), if the signal has uncertainties, the `signal`
  vector can also be a vector of `Measurement` objects, and this argument should
  be omitted

* `times_fit`: the vector of times at which the model will be calculated.  It
  defaults to `times`.  If provided, it must come after `frequency`

Optional keyword arguments `center_data` and `fit_mean` have the same meaning as
in [`lombscargle`](@ref):

- `fit_mean`: whether to fit for the mean. If this is `false`, like in the
  original Lomb--Scargle periodogram, ``\\mathbf{A}`` does not have the third
  column of ones, ``c_f`` is set to ``0`` and the unknown vector to be determined
  becomes ``x = [a_f, b_f]^\\text{T}``
- `center_data`: whether the data should be pre-centered before solving the
  linear system. This is particularly important if `fit_mean=false`
"""
function model(t::AbstractVector{<:Real}, s::AbstractVector{T},
               errors::AbstractVector{<:Real}, f::Real,
               t_fit::AbstractVector{<:Real}=t;
               center_data::Bool=true, fit_mean::Bool=true) where {T<:Real}
    @assert length(t) == length(s) == length(errors)
    if center_data
        # Compute weights vector
        w = 1 ./ (errors .^ 2)
        m = (w⋅s)/sum(w)
        y = (s .- m) ./ errors
    else
        # We don't want to center the data: the mean is 0 and the signal is left
        # unchanged
        m = zero(T)
        y = s./errors
    end
    # The Lomb–Scargle periodogram looks for the best fitting sinusoidal
    # function
    #     a·cos(ωt) + b·sin(ωt) + c
    # In order to find the coefficients a, b, c for the given frequency we can
    # solve the linear system A·x = y, where A is the matrix with rows:
    # [cos(ωt) sin(ωt) 1]; x is the column vector [a, b, c], and y is the
    # column vector of the signal
    ω = 2 * f
    if fit_mean
        a, b, c = [cospi.(ω .* t) sinpi.(ω .* t)  fill(1, size(t))] ./ errors \ y
        return a .* cospi.(ω .* t_fit) .+ b .* sinpi.(ω .* t_fit) .+ (c + m)
    else
        # If fit_mean==false, the model to be fitted is a·cos(ωt) + b·sin(ωt)
        a, b = [cospi.(ω .* t) sinpi.(ω .* t)] ./ errors \ y
        return a .* cospi.(ω .* t_fit) .+ b .* sinpi.(ω .* t_fit) .+ m
    end
end

# No uncertainties: errors=ones(s)
model(t::AbstractVector{<:Real}, s::AbstractVector{<:Real},
      f::Real, t_fit::AbstractVector{<:Real}=t; kwargs...) =
          model(t, s, fill(1, size(s)), f, t_fit; kwargs...)

# Uncertainties provided via Measurement type
model(t::AbstractVector{<:Real}, s::AbstractVector{<:Measurement},
      f::Real, t_fit::AbstractVector{<:Real}=t; kwargs...) =
          model(t, Measurements.value.(s), Measurements.uncertainty.(s),
                f, t_fit; kwargs...)


text:
mutable struct WalkingPulseEffect{T<:AbstractFloat} <: ConfigurableEffect
    config::EffectConfig
    pulse_length::Int
    pulse::Array{UInt8}
    pulse_shape::Array{T}
    spec::Array{T}
    maxspec::T
    assgn::Array{Int}
    starts::StepRange
    walking_offset::Float64
    function WalkingPulseEffect{T}(leddata, config, ana) where T<:AbstractFloat
        pulse_length = round(Int, 15*config.scaling)
        pulse = zeros(UInt8, pulse_length)
        pulse_shape = zeros(pulse_length)
        spec, maxspec = binFFT(ana, 25)
        walking_offset = 0.0
        numPulses = floor(Int, size(leddata.channels[1],1)/30)
        starts = (1:floor(Int, size(leddata.channels[1],1)/numPulses ):size(leddata.channels[1],1))[1:end-1] + round(Int, walking_offset*size(leddata.channels[1],1)/100)
        assgn = (starts[1]:(starts[1]+size(pulse,1)-1)).%size(leddata.channels[1],1) + 1
        return new(config, pulse_length, pulse, pulse_shape, spec, maxspec, assgn, starts, walking_offset)
    end
end

function update_vars!(leddata::LEDArray, ana::AudioAnalysis, effect::WalkingPulseEffect)
    effect.spec, effect.maxspec = binFFT(ana, 25)
    effect.spec./=fft_scale[]
    effect.spec[isnan.(effect.spec)] = 0.0
    effect.spec .= abs.(effect.spec)
    effect.pulse_length = round(Int, 15*effect.config.scaling)
end

function update!(leddata::LEDArray, ana::AudioAnalysis, effect::WalkingPulseEffect)

    update_vars!(leddata, ana, effect)
    default_color = hsl_to_rgb(effect.config.secondary_color)
    pulse_coef = min(effect.spec[1]+effect.spec[2]+mean(effect.spec[3:5])/3.5, 1)
    pulse_coef = pulse_coef>0.05 ? pulse_coef : 0.0
    effect.pulse_shape = [i%effect.pulse_length/4==0 ? abs(rand()*0.5*pulse_coef) : abs(sin(pi/4*i)*pulse_coef) for i in linspace(0,4,effect.pulse_length)]
    effect.pulse_shape[isnan.(effect.pulse_shape)] = 0.0
    effect.pulse = zeros(UInt8, effect.pulse_length, 3)
    effect.pulse_shape .= (effect.pulse_shape).^0.25

    effect.pulse = vcat([hsl_to_rgb(effect.config.primary_color.h/360, 0.95*effect.pulse_shape[i]+0.05, 0.45*effect.pulse_shape[i]+0.05)' for i in 1:size(effect.pulse,1)]...)
    starts_dict = Dict{Int, StepRange}()
    assgn_dict = Dict{Tuple{Int,Int}, Array{Int,1}}()

    for channel in leddata.channels
        for j in 1:size(channel,1)
            for k in 1:size(channel,2)
                channel[j,k] = default_color[k]
            end
        end
        numPulses = floor(Int, size(channel,1)/30)

        if !(size(channel,1) in keys(starts_dict) && starts_dict[size(channel,1)][1] == 1 + round(Int, effect.walking_offset*size(channel,1)/100))
            effect.starts = (1:floor(Int, size(channel,1)/numPulses ):size(channel,1))[1:end-1] + round(Int, effect.walking_offset*size(channel,1)/100)
        else
            effect.starts = starts_dict[size(channel,1)]
        end
        for k in eachindex(effect.starts)
            if !((size(channel,1), effect.starts[k]) in keys(assgn_dict))
                effect.assgn = (effect.starts[k]:(effect.starts[k]+size(effect.pulse,1)-1)).%size(channel,1) + 1
                assgn_dict[(size(channel,1), effect.starts[k])] = effect.assgn
            else
                effect.assgn = assgn_dict[(size(channel,1), effect.starts[k])]
            end
            for i in eachindex(effect.assgn)
                for j in 1:3
                    channel[effect.assgn[i],j] = effect.pulse[i,j]
                end
            end
        end
        update!(channel)
    end

    if effect.walking_offset != 10000 && effect.walking_offset % 100 == 0
        effect.walking_offset = 10000
    end
    effect.walking_offset += 0.5*effect.config.speed
end