text:
B_data = [ -2.5   0.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   5.26957  -2.45701   -0.812567   0.0        0.0         0.0         0.0         0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
1.0  -3.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   5.26957  -2.45701   -0.812567   0.0        0.0         0.0         0.0         0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.0       0.0        0.0        0.0        0.0         0.0         0.0         0.0        0.0         0.728864   -0.283275   0.0544112    0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0  -2.5  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       2.0        0.0        0.0        0.0         0.0         0.0         0.0        0.0         1.05443    -0.549114  -0.00531726   0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        2.0        0.0        0.0         0.0         0.0         0.0        0.0         0.564472   -0.211292   0.14682      0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        2.0        0.0         0.0         0.0         0.0        0.0         0.0         0.0        0.0          0.728864   -0.283275    0.0544112    0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0  -2.5  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        2.0         0.0         0.0         0.0        0.0         0.0         0.0        0.0          1.05443    -0.549114   -0.00531726   0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         2.0         0.0         0.0        0.0         0.0         0.0        0.0          0.564472   -0.211292    0.14682      0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         2.0         0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          0.728864  -0.283275   0.0544112    0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         2.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          1.05443   -0.549114  -0.00531726   0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         0.0        2.0         0.0         0.0        0.0          0.0         0.0         0.0          0.564472  -0.211292   0.14682      0.0   0.0
0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   2.91546  -1.1331     0.217645   0.0        0.0         0.0         0.0         0.0        0.0         0.5         0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   4.21772  -2.19645   -0.021269   0.0        0.0         0.0         0.0         0.0        0.0         0.0         0.5        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0   2.25789  -0.845166   0.587279   0.0        0.0         0.0         0.0         0.0        0.0         0.0         0.0        0.5          0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        2.91546   -1.1331      0.217645    0.0         0.0        0.0         0.0         0.0        0.0          0.5         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5  -0.0   0.0   0.0   0.0   0.0       0.0        0.0        4.21772   -2.19645    -0.021269    0.0         0.0        0.0         0.0         0.0        0.0          0.0         0.5         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0   0.0       0.0        0.0        2.25789   -0.845166    0.587279    0.0         0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.5          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         2.91546    -1.1331     0.217645    0.0         0.0        0.0          0.0         0.0         0.0          0.5        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5  -0.0   0.0       0.0        0.0        0.0        0.0         0.0         4.21772    -2.19645   -0.021269    0.0         0.0        0.0          0.0         0.0         0.0          0.0        0.5        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -2.5   0.0       0.0        0.0        0.0        0.0         0.0         2.25789    -0.845166   0.587279    0.0         0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.5          0.0   0.0
0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -3.50473   0.0        0.0        3.96463   -1.84856    -0.611345    0.0         0.0        0.0         1.45773    -0.56655    0.108822     0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0      -3.50473   -1.55879    8.07172   -3.76354    -1.24466     0.0         0.0        0.0         2.10886    -1.09823   -0.0106345    0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       1.55879   -3.50473   -0.142449   0.0664188   0.0219656   0.0         0.0        0.0         1.12894    -0.422583   0.29364      0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0       -6.51418    0.0         0.0        11.8939     -5.54569   -1.83404     0.0         0.0        0.0          1.45773    -0.56655     0.108822     0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0       -6.51418    -4.67637    24.2152    -11.2906    -3.73397     0.0         0.0        0.0          2.10886    -1.09823    -0.0106345    0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        4.67637    -6.51418    -0.427348    0.199257   0.0658969   0.0         0.0        0.0          1.12894    -0.422583    0.29364      0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0        -2.56427     0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          2.94447   -1.25976   -0.120432     0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0        -2.56427   -0.584547    0.0         0.0        0.0          0.0         0.0         0.0          5.13576   -2.50956   -0.477381     0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         0.584547  -2.56427     0.0         0.0        0.0          0.0         0.0         0.0          1.07552   -0.397676   0.301877    -0.0  -0.0
0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   8.86241  -3.90594   -0.699373   0.0        0.0         0.0         0.0         0.0        0.0        -5.25709     0.0        0.0          0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  16.3253   -7.84177   -1.88825    0.0        0.0         0.0         0.0         0.0        0.0         0.0        -5.25709   -2.33819      0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
-0.0  -0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0   2.04421  -0.745538   0.620228   0.0        0.0         0.0         0.0         0.0        0.0         0.0         2.33819   -5.25709      0.0         0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        2.91546   -1.1331      0.217645    0.0         0.0        0.0        17.8409     -8.31853   -2.75105     -9.77127     0.0         0.0          0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        4.21772   -2.19645    -0.021269    0.0         0.0        0.0        36.3227    -16.9359    -5.60095      0.0        -9.77127    -7.01456      0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0   0.0       0.0        0.0        2.25789   -0.845166    0.587279    0.0         0.0        0.0        -0.641022    0.298885   0.0988454    0.0         7.01456    -9.77127      0.0        0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         2.91546    -1.1331     0.217645    0.0         0.0        0.0          2.23011    -1.03982    -0.343882    -3.84641    0.0        0.0          0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         4.21772    -2.19645   -0.021269    0.0         0.0        0.0          4.54034    -2.11699    -0.700119     0.0       -3.84641   -0.87682      0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0  -0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   1.0   0.0       0.0        0.0        0.0        0.0         0.0         2.25789    -0.845166   0.587279    0.0         0.0        0.0         -0.0801277   0.0373606   0.0123557    0.0        0.87682   -3.84641      0.0   0.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          1.31739   -0.614252  -0.203142    -2.5   2.0
0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0   0.0       0.0        0.0        0.0        0.0         0.0         0.0         0.0        0.0         0.0         0.0        0.0          0.0         0.0         0.0          2.63479   -1.2285    -0.406283     1.0  -2.0]


text:
using Test
import Pluto: update_run!, WorkspaceManager, ClientSession, ServerSession, Notebook, Cell

@testset "Rich output" begin
    @testset "Interactive inspector" begin
        🍭 = ServerSession()
        🍭.options.evaluation.workspace_use_distributed = false
        fakeclient = ClientSession(:fake, nothing)
        🍭.connected_clients[fakeclient.id] = fakeclient

        notebook = Notebook([
                Cell("[1,1,[1]]"),
                Cell("Dict(:a => [:b, :c])"),
                Cell("[3, Dict()]"),
                Cell("[4,[3, Dict()]]"),
                Cell("[5, missing, 5]"),
                Cell("[]"),
                Cell("(7,7)"),
                Cell("(a=8,b=[8])"),
                Cell("Ref(9)"),
                Cell("struct Ten x end"),
                Cell("Ten(11)"),
                Cell("Vector(undef, 12)"),
            ])
            fakeclient.connected_notebook = notebook

        update_run!(🍭, notebook, notebook.cells)

        @test notebook.cells[1].repr_mime isa MIME"application/vnd.pluto.tree+xml"
        @test notebook.cells[2].repr_mime isa MIME"application/vnd.pluto.tree+xml"
        @test notebook.cells[3].repr_mime isa MIME"application/vnd.pluto.tree+xml"
        @test notebook.cells[4].repr_mime isa MIME"application/vnd.pluto.tree+xml"
        @test notebook.cells[5].repr_mime isa MIME"application/vnd.pluto.tree+xml"
        @test notebook.cells[6].repr_mime isa MIME"application/vnd.pluto.tree+xml"
        @test notebook.cells[7].repr_mime isa MIME"application/vnd.pluto.tree+xml"
        @test notebook.cells[8].repr_mime isa MIME"application/vnd.pluto.tree+xml"
        @test notebook.cells[9].repr_mime isa MIME"application/vnd.pluto.tree+xml"
        @test notebook.cells[11].repr_mime isa MIME"application/vnd.pluto.tree+xml"
        @test notebook.cells[12].repr_mime isa MIME"application/vnd.pluto.tree+xml"

        WorkspaceManager.unmake_workspace((🍭, notebook))

    end
end

ext: jl


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function add_chp_fuel_burn_constraints(m, p; _n="")
    # Fuel burn slope and intercept
    fuel_burn_full_load = 1.0 / p.s.chp.elec_effic_full_load  # [kWt/kWe]
    fuel_burn_half_load = 0.5 / p.s.chp.elec_effic_half_load  # [kWt/kWe]
    fuel_burn_slope = (fuel_burn_full_load - fuel_burn_half_load) / (1.0 - 0.5)  # [kWt/kWe]
    fuel_burn_intercept = fuel_burn_full_load - fuel_burn_slope * 1.0  # [kWt/kWe_rated]

    # Fuel cost
    fuel_cost_per_kwh = p.s.chp.fuel_cost_per_mmbtu / MMBTU_TO_KWH
    fuel_cost_series = per_hour_value_to_time_series(fuel_cost_per_kwh, p.s.settings.time_steps_per_hour,
                                                     "CHP.fuel_cost_per_mmbtu")
    m[:TotalCHPFuelCosts] = @expression(m, p.pwf_fuel["CHP"] *
        sum(m[:dvFuelUsage]["CHP", ts] * fuel_cost_series[ts] for ts in p.time_steps)
    )
    # Conditionally add dvFuelBurnYIntercept if coefficient p.FuelBurnYIntRate is greater than ~zero
    if fuel_burn_intercept > 1.0E-7
        dv = "dvFuelBurnYIntercept"*_n
        m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv, lower_bound=0)

        #Constraint (1c1): Total Fuel burn for CHP **with** y-intercept fuel burn and supplementary firing
        @constraint(m, CHPFuelBurnCon[t in p.techs.chp, ts in p.time_steps],
            m[Symbol("dvFuelUsage"*_n)][t,ts]  == p.hours_per_timestep * (
                m[Symbol("dvFuelBurnYIntercept"*_n)][t,ts] +
                p.production_factor[t,ts] * fuel_burn_slope * m[Symbol("dvRatedProduction"*_n)][t,ts] +
                m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] / p.s.chp.supplementary_firing_efficiency
            )
        )
        #Constraint (1d): Y-intercept fuel burn for CHP
        @constraint(m, CHPFuelBurnYIntCon[t in p.techs.chp, ts in p.time_steps],
                    fuel_burn_intercept * m[Symbol("dvSize"*_n)][t] - p.s.chp.max_kw *
                    (1-m[Symbol("binCHPIsOnInTS"*_n)][t,ts])  <= m[Symbol("dvFuelBurnYIntercept"*_n)][t,ts]
                    )
    else
        #Constraint (1c2): Total Fuel burn for CHP **without** y-intercept fuel burn
        @constraint(m, CHPFuelBurnConLinear[t in p.techs.chp, ts in p.time_steps],
            m[Symbol("dvFuelUsage"*_n)][t,ts]  == p.hours_per_timestep * (
                p.production_factor[t,ts] * fuel_burn_slope * m[Symbol("dvRatedProduction"*_n)][t,ts] +
                m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] / p.s.chp.supplementary_firing_efficiency
            )
        )
    end
end

function add_chp_thermal_production_constraints(m, p; _n="")
    # Thermal production slope and intercept
    thermal_prod_full_load = 1.0 / p.s.chp.elec_effic_full_load * p.s.chp.thermal_effic_full_load  # [kWt/kWe]
    thermal_prod_half_load = 0.5 / p.s.chp.elec_effic_half_load * p.s.chp.thermal_effic_half_load   # [kWt/kWe]
    thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5)  # [kWt/kWe]
    thermal_prod_intercept = thermal_prod_full_load - thermal_prod_slope * 1.0  # [kWt/kWe_rated

    dv = "dvProductionToWaste"*_n
    m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv, lower_bound=0)

    # Conditionally add dvThermalProductionYIntercept if coefficient p.s.chpThermalProdIntercept is greater than ~zero
    if thermal_prod_intercept > 1.0E-7
        dv = "dvThermalProductionYIntercept"*_n
        m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv, lower_bound=0)

        #Constraint (2a-1): Upper Bounds on Thermal Production Y-Intercept
        @constraint(m, CHPYInt2a1Con[t in p.techs.chp, ts in p.time_steps],
            m[Symbol("dvThermalProductionYIntercept"*_n)][t,ts] <= thermal_prod_intercept * m[Symbol("dvSize"*_n)][t]
        )
        # Constraint (2a-2): Upper Bounds on Thermal Production Y-Intercept
        @constraint(m, CHPYInt2a2Con[t in p.techs.chp, ts in p.time_steps],
            m[Symbol("dvThermalProductionYIntercept"*_n)][t,ts] <= thermal_prod_intercept * p.s.chp.max_kw
            * m[Symbol("binCHPIsOnInTS"*_n)][t,ts]
        )
        #Constraint (2b): Lower Bounds on Thermal Production Y-Intercept
        @constraint(m, CHPYInt2bCon[t in p.techs.chp, ts in p.time_steps],
            m[Symbol("dvThermalProductionYIntercept"*_n)][t,ts] >= thermal_prod_intercept * m[Symbol("dvSize"*_n)][t]
            - thermal_prod_intercept * p.s.chp.max_kw * (1 - m[Symbol("binCHPIsOnInTS"*_n)][t,ts])
        )
        # Constraint (2c): Thermal Production of CHP
        # Note: p.HotWaterAmbientFactor[t,ts] * p.HotWaterThermalFactor[t,ts] removed from this but present in math
        @constraint(m, CHPThermalProductionCon[t in p.techs.chp, ts in p.time_steps],
            m[Symbol("dvThermalProduction"*_n)][t,ts] ==
            thermal_prod_slope * p.production_factor[t,ts] * m[Symbol("dvRatedProduction"*_n)][t,ts]
            + m[Symbol("dvThermalProductionYIntercept"*_n)][t,ts] +
            m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts]
        )
    else
        @constraint(m, CHPThermalProductionConLinear[t in p.techs.chp, ts in p.time_steps],
            m[Symbol("dvThermalProduction"*_n)][t,ts] ==
            thermal_prod_slope * p.production_factor[t,ts] * m[Symbol("dvRatedProduction"*_n)][t,ts] +
            m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts]
        )
    end

end

"""
    add_chp_supplementary_firing_constraints(m, p; _n="")

Used by add_chp_constraints to add supplementary firing constraints if
    p.s.chp.supplementary_firing_max_steam_ratio > 1.0 to add CHP supplementary firing operating constraints.
    Else, the supplementary firing dispatch and size decision variables are set to zero.
"""
function add_chp_supplementary_firing_constraints(m, p; _n="")
    thermal_prod_full_load = 1.0 / p.s.chp.elec_effic_full_load * p.s.chp.thermal_effic_full_load  # [kWt/kWe]
    thermal_prod_half_load = 0.5 / p.s.chp.elec_effic_half_load * p.s.chp.thermal_effic_half_load   # [kWt/kWe]
    thermal_prod_slope = (thermal_prod_full_load - thermal_prod_half_load) / (1.0 - 0.5)  # [kWt/kWe]

    # Constrain upper limit of dvSupplementaryThermalProduction, using auxiliary variable for (size * useSupplementaryFiring)
    @constraint(m, CHPSupplementaryFireCon[t in p.techs.chp, ts in p.time_steps],
                m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <=
                (p.s.chp.supplementary_firing_max_steam_ratio - 1.0) * p.production_factor[t,ts] * (thermal_prod_slope * m[Symbol("dvSupplementaryFiringSize"*_n)][t] + m[Symbol("dvThermalProductionYIntercept"*_n)][t,ts])
                )
    # Constrain lower limit of 0 if CHP tech is off
    @constraint(m, NoCHPSupplementaryFireOffCon[t in p.techs.chp, ts in p.time_steps],
                !m[Symbol("binCHPIsOnInTS"*_n)][t,ts] => {m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts] <= 0.0}
                )
end

function add_binCHPIsOnInTS_constraints(m, p; _n="")
    # Note, min_turn_down_pct for CHP is only enforced in p.time_steps_with_grid
    @constraint(m, [t in p.techs.chp, ts in p.time_steps_with_grid],
        m[Symbol("dvRatedProduction"*_n)][t, ts] <= p.s.chp.max_kw * m[Symbol("binCHPIsOnInTS"*_n)][t, ts]
    )
    @constraint(m, [t in p.techs.chp, ts in p.time_steps_with_grid],
        p.s.chp.min_turn_down_pct * m[Symbol("dvSize"*_n)][t] - m[Symbol("dvRatedProduction"*_n)][t, ts] <=
        p.s.chp.max_kw * (1 - m[Symbol("binCHPIsOnInTS"*_n)][t, ts])
    )
end

function add_chp_rated_prod_constraint(m, p; _n="")
    @constraint(m, [t in p.techs.chp, ts in p.time_steps],
        m[Symbol("dvSize"*_n)][t] >= m[Symbol("dvRatedProduction"*_n)][t, ts]
    )
end

"""
    add_chp_hourly_om_charges(m, p; _n="")

- add decision variable "dvOMByHourBySizeCHP"*_n for the hourly CHP operations and maintenance costs
- add the cost to TotalPerUnitHourOMCosts
"""
function add_chp_hourly_om_charges(m, p; _n="")
    dv = "dvOMByHourBySizeCHP"*_n
    m[Symbol(dv)] = @variable(m, [p.techs.chp, p.time_steps], base_name=dv, lower_bound=0)

    #Constraint CHP-hourly-om-a: om per hour, per time step >= per_unit_size_cost * size for when on, >= zero when off
	@constraint(m, CHPHourlyOMBySizeA[t in p.techs.chp, ts in p.time_steps],
        p.s.chp.om_cost_per_hr_per_kw_rated * m[Symbol("dvSize"*_n)][t] -
        p.max_size[t] * p.s.chp.om_cost_per_hr_per_kw_rated * (1-m[Symbol("binCHPIsOnInTS"*_n)][t,ts])
            <= m[Symbol("dvOMByHourBySizeCHP"*_n)][t, ts]
    )
	#Constraint CHP-hourly-om-b: om per hour, per time step <= per_unit_size_cost * size for each hour
	@constraint(m, CHPHourlyOMBySizeB[t in p.techs.chp, ts in p.time_steps],
        p.s.chp.om_cost_per_hr_per_kw_rated * m[Symbol("dvSize"*_n)][t]
            >= m[Symbol("dvOMByHourBySizeCHP"*_n)][t, ts]
    )
	#Constraint CHP-hourly-om-c: om per hour, per time step <= zero when off, <= per_unit_size_cost*max_size
	@constraint(m, CHPHourlyOMBySizeC[t in p.techs.chp, ts in p.time_steps],
        p.max_size[t] * p.s.chp.om_cost_per_hr_per_kw_rated * m[Symbol("binCHPIsOnInTS"*_n)][t,ts]
            >= m[Symbol("dvOMByHourBySizeCHP"*_n)][t, ts]
    )

    m[:TotalHourlyCHPOMCosts] = @expression(m, p.third_party_factor * p.pwf_om *
    sum(m[Symbol(dv)][t, ts] * p.hours_per_timestep for t in p.techs.chp, ts in p.time_steps))
    nothing
end

"""
    add_chp_constraints(m, p; _n="")

Used in src/reopt.jl to add_chp_constraints if !isempty(p.techs.chp) to add CHP operating constraints and
cost expressions.
"""
function add_chp_constraints(m, p; _n="")
    # TODO if chp.min_turn_down_pct is 0.0, and there is no fuel burn or thermal y-intercept, we don't need the binary below
    @warn """Adding binary variable to model CHP.
                Some solvers are very slow with integer variables"""
    @variables m begin
        binCHPIsOnInTS[p.techs.chp, p.time_steps], Bin  # 1 If technology t is operating in time step; 0 otherwise
    end

    m[:TotalHourlyCHPOMCosts] = 0
    m[:TotalCHPFuelCosts] = 0
    m[:TotalCHPPerUnitProdOMCosts] = @expression(m, p.third_party_factor * p.pwf_om *
        sum(p.s.chp.om_cost_per_kwh * p.hours_per_timestep *
        m[:dvRatedProduction][t, ts] for t in p.techs.chp, ts in p.time_steps)
    )

    if p.s.chp.om_cost_per_hr_per_kw_rated > 1.0E-7
        add_chp_hourly_om_charges(m, p)
    end

    add_chp_fuel_burn_constraints(m, p; _n=_n)
    add_chp_thermal_production_constraints(m, p; _n=_n)
    add_binCHPIsOnInTS_constraints(m, p; _n=_n)
    add_chp_rated_prod_constraint(m, p; _n=_n)

    if p.s.chp.supplementary_firing_max_steam_ratio > 1.0
        add_chp_supplementary_firing_constraints(m,p; _n=_n)
    else
        for t in p.techs.chp
            fix(m[Symbol("dvSupplementaryFiringSize"*_n)][t], 0.0, force=true)
            for ts in p.time_steps
                fix(m[Symbol("dvSupplementaryThermalProduction"*_n)][t,ts], 0.0, force=true)
            end
        end
    end
end


text:
######### MachO.jl - An implementation of the MachO File format ################
#
## General design philosophy
#
# All methods should operate on seekable IO object. In particular this means
# that the basic interface buffers as little as possible (e.g. the MachOHandle
# object does not contain the mach_header).
#
################################################################################

__precompile__()
module MachO

# For printing
import Base: show, print, unsafe_string, showcompact, readuntil

# For endianness-handling
using StructIO
import StructIO: unpack

# This package implements the ObjFileBase interface
import ObjFileBase
import ObjFileBase: sectionsize, sectionoffset, readheader, ObjectHandle, readmeta,
    strtab_lookup, debugsections, endianness, load_strtab, deref, sectionname,
    sectionaddress, handle, symbolnum, printfield, symname, symbolvalue, isundef,
    intptr, isBSS

# For MachO datatypes (e.g. fixed size string)
import Base: ==, *

# Reexports from ObjFileBase
export sectionsize, sectionoffset, readheader, readmeta,
    debugsections

export LoadCmds, symname, segname

# Tree Interface for visualization
using AbstractTrees

import AbstractTrees: children

############################ Data Structures ###################################

#
# Contains all constants defined in the MachO standard with values taken
# from the appropriate headers on MacOS 10.9
#
include("constants.jl")

#
# Represents the actual MachO file
#
immutable MachOHandle{T<:IO} <: ObjectHandle
    # The IO object. This field is speciallized on to avoid dispatch performance
    # hits, especially when operating on an IOBuffer, which is an important
    # usecase for in-memory files
    io::T
    # position(io) of the start of the file in the io stream.
    start::Int
    # Whether or not the data is bswap'ed in memory (i.e. has different
    # endianness)
    bswapped::Bool
    # Whether or not the file is 64bit
    is64::Bool
end
__init__() = push!(ObjFileBase.ObjHandles, MachOHandle)
ObjFileBase.handle(handle::MachOHandle) = handle
Base.eof(handle::MachOHandle) = eof(handle.io)

endianness(oh::MachOHandle) = oh.bswapped ? :SwappedEndian : :NativeEndian
intptr(oh::MachOHandle) = oh.is64 ? UInt64 : UInt32

function show(io::IO,h::MachOHandle)
    print(io,"MachO handle (")
    print(io,h.is64?"64-bit":"32-bit")
    h.bswapped && print(io,",swapped")
    print(io,")")
end

############################ Data Structures ###################################
#
# This section contains data structures as defined by the MachO specification.
# They are used below to acutally read in the data and may be ocassionally
# referenced from interface structs where this is convenient and no other
# interface exists
#
################################################################################

abstract MachOLC

# Dummy lc that we use to return when we don't know what a certain load command is
immutable dummy_lc <: MachOLC
end

@io immutable mach_header
    magic::UInt32
    cputype::UInt32
    cpusubtype::UInt32
    filetype::UInt32
    ncmds::UInt32
    sizeofcmds::UInt32
    flags::UInt32
end

@io immutable mach_header_64
    magic::UInt32
    cputype::UInt32
    cpusubtype::UInt32
    filetype::UInt32
    ncmds::UInt32
    sizeofcmds::UInt32
    flags::UInt32
    reserved::UInt32
end
mach_header_64(magic,cputype,cpusubtype,filetype,ncmds,sizeofcmds,flags) =
    mach_header_64(magic,cputype,cpusubtype,filetype,ncmds,sizeofcmds,flags,0)

@io immutable load_command
    cmd::UInt32
    cmdsize::UInt32
end

# A 16 byte string, represented as a UInt128, but shown as a string
@io immutable small_fixed_string
    string::UInt128
end

@io immutable uuid_command <: MachOLC
    uuid::UInt128
end

immutable thread_command <: MachOLC
    flavor::UInt32
    count::UInt32
    data::Vector{UInt}
end

@io immutable entry_point_command <: MachOLC
    cmdsize::UInt32
    entryoff::UInt64
    stacksize::UInt64
end

@io immutable segment_command <: MachOLC
    segname::small_fixed_string
    vmaddr::UInt32
    vmsize::UInt32
    fileoff::UInt32
    filesize::UInt32
    maxprot::UInt32
    initprot::UInt32
    nsects::UInt32
    flags::UInt32
end

@io immutable segment_command_64 <: MachOLC
    segname::small_fixed_string
    vmaddr::UInt64
    vmsize::UInt64
    fileoff::UInt64
    filesize::UInt64
    maxprot::Cint
    initprot::Cint
    nsects::UInt32
    flags::UInt32
end

@io immutable section <: ObjFileBase.Section{MachOHandle}
    sectname::small_fixed_string
    segname::small_fixed_string
    addr::UInt32
    size::UInt32
    offset::UInt32
    align::UInt32
    reloff::UInt32
    nreloc::UInt32
    flags::UInt32
    reserved1::UInt32
    reserved2::UInt32
end

@io immutable section_64 <: ObjFileBase.Section{MachOHandle}
    sectname::small_fixed_string
    segname::small_fixed_string
    addr::UInt64
    size::UInt64
    offset::UInt32
    align::UInt32
    reloff::UInt32
    nreloc::UInt32
    flags::UInt32
    reserved1::UInt32
    reserved2::UInt32
end
isBSS(sec::Union{section, section_64}) = (sec.flags & SECTION_TYPE) == S_ZEROFILL

@io immutable relocation_info <: ObjFileBase.Relocation{MachOHandle}
    address::Int32
    target::UInt32
end

@io immutable symtab_command <: MachOLC
    symoff::UInt32
    nsyms::UInt32
    stroff::UInt32
    strsize::UInt32
end
symtab_command() = symtab_command(0,0,0,0)

@io immutable dysymtab_command <: MachOLC
    ilocalsym::UInt32
    nlocalsym::UInt32
    iextdefsym::UInt32
    nextdefsym::UInt32
    iundefsym::UInt32
    nundefsym::UInt32
    tocoff::UInt32
    ntoc::UInt32
    modtaboff::UInt32
    nmodtab::UInt32
    extrefsymoff::UInt32
    nextrefsyms::UInt32
    indirectsymoff::UInt32
    nindirectsyms::UInt32
    extreloff::UInt32
    nextrel::UInt32
    locreloff::UInt32
    nlocrel::UInt32
end

@io immutable version_min_macosx_command <: MachOLC
    version::UInt32
    sdk::UInt32
end

immutable dylib_command <: MachOLC
    offset::UInt32
    timestamp::UInt32
    current_version::UInt32
    compatibilty::UInt32

    # Read in automatically, when possible, via offset
    name::AbstractString
end

immutable dylinker_command <: MachOLC
    name::AbstractString
end

@io immutable routines_command_64 <: MachOLC
    init_address::UInt64
    init_module::UInt64
    reserverd::NTuple{6, UInt64}
end

immutable sub_client_command  <: MachOLC
    name::AbstractString
end

# Read in a C string, until we reach the end of the string or max out at max_len
function read_unsafe_string(io, max_len)
    str = UInt8[]
    idx = 0
    c = read(io, UInt8)
    while c != 0x00 && idx < max_len
        push!(str, c)
        c = read(io, UInt8)
        idx += 1
    end
    return String(str)
end

function unpack_lcstr{ioT<:IO}(h::MachOHandle{ioT}, offset, min_offset, max_offset)
    # Perform sanity checking on offset; if it is too small or too large,
    # don't try to read the lc_str, just assign it "<lc_str offset corrupt>"
    if offset >= min_offset && offset < max_offset
        # Seek to the previously extracted offset, minus minlen
        skip(h.io, offset - min_offset)

        # Read in the cstring
        lc_str = read_unsafe_string(h.io, max_offset - offset)
    else
        # If we are outside the bounds, (either the string begins in the middle
        # of the rest of the structure, or it begins outside of this load command)
        # do not attempt to automatically read it
        lc_str = "<lc_str offset corrupt>"
    end
end

function unpack{ioT<:IO}(h::MachOHandle{ioT},::Type{dylib_command},cmdsize::UInt32)
    # Get the offset
    offset = unpack(h, UInt32)

    # Now get timestamp, current_version and compatibilty
    timestamp = unpack(h, UInt32)
    current_version = unpack(h, UInt32)
    compatibilty = unpack(h, UInt32)

    # Grab our name, if we can (e.g. if offset is within bounds)
    name = unpack_lcstr(h, offset, 6*sizeof(UInt32), cmdsize)
    return dylib_command(offset, timestamp, current_version, compatibilty, name)
end

function unpack{ioT<:IO}(h::MachOHandle{ioT},
    T::Union{Type{dylinker_command},Type{sub_client_command}},cmdsize::UInt32)
    offset = unpack(h, UInt32)
    return T(unpack_lcstr(h, offset, 6*sizeof(UInt32), cmdsize))
end

@io immutable dyld_info_command <: MachOLC
    rebase_off::UInt32
    rebase_size::UInt32
    bind_off::UInt32
    bind_size::UInt32
    weak_bind_off::UInt32
    weak_bind_size::UInt32
    lazy_bind_off::UInt32
    lazy_bind_size::UInt32
    export_off::UInt32
    export_size::UInt32
end

@io immutable source_version_command <: MachOLC
    version::UInt64
end

@io immutable linkedit_data_command <: MachOLC
    offset::UInt32
    size::UInt32
end

immutable sub_framework_command <: MachOLC
    offset::UInt32

    # Read in automatically, when possible, via offset
    umbrella::AbstractString
end

immutable rpath_command <: MachOLC
    offset::UInt32

    # Read in automatically, when possible, via offset
    path::AbstractString
end

for T in [sub_framework_command, rpath_command]
    @eval function unpack{ioT<:IO}(h::MachOHandle{ioT},::Type{$T},cmdsize::UInt32)
        # Get the offset
        offset = unpack(h, UInt32)

        # Grab our path if we can (e.g. if offset is within bounds)
        path = unpack_lcstr(h, offset, 3*sizeof(UInt32), cmdsize)
        return $T(offset, path)
    end
end

@io immutable nlist <: ObjFileBase.SymtabEntry{MachOHandle}
    n_strx::UInt32
    n_type::UInt8
    n_sect::UInt8
    n_desc::UInt16
    n_value::UInt32
end

@io immutable nlist_64 <: ObjFileBase.SymtabEntry{MachOHandle}
    n_strx::UInt32
    n_type::UInt8
    n_sect::UInt8
    n_desc::UInt16
    n_value::UInt64
end

@io immutable fat_header
    nfat_arch::UInt32
end

@io immutable fat_arch
    cputype::UInt32
    cpusubtype::UInt32
    offset::UInt32
    size::UInt32
    align::UInt32
end

########################### Printing Data Structures ###########################
#
# This prints the basic Mach-O data structures above. Where there is no good
# reason not to, the output matches that of otool.
#
################################################################################

function unsafe_string(x::small_fixed_string)
    a8 = reinterpret(UInt8,[x.string])
    z = findfirst(a8,0)
    String(a8[1:(z == 0 ? length(a8) : z-1)])
end
show(io::IO,x::small_fixed_string) = show(io,unsafe_string(x))
print(io::IO,x::small_fixed_string) = print(io,unsafe_string(x))

# These can all be made a lot faster if they ever show up in a profile
Base.isempty(x::small_fixed_string) = isempty(unsafe_string(x))
Base.length(x::small_fixed_string) = length(unsafe_string(x))

==(x::small_fixed_string,y::AbstractString) = unsafe_string(x) == y
==(x::AbstractString,y::small_fixed_string) = y==x

*(a::String,b::small_fixed_string) = a*unsafe_string(b)

# TODO: Implement
printflags(io,flags) = nothing

function show(io::IO,h::Union{mach_header,mach_header_64})
    print(io,"Mach header (",isa(h,mach_header_64)?"64-bit":"32-bit",")\n")
    print(io,"")
    print(io,"""
          magic  cputype cpusubtype     filetype ncmds sizeofcmds    flags
    """)
    printfield(io,MAGICS[h.magic],11)
    printfield(io,CPUTYPES[h.cputype],9)
    print(io," "^11)
    # TODO: cpusubtype printing
    #printfield(io,CPUSUBTYPES[h.cpusubtype],8)
    # otool has caps == h.cpusubtype & CPU_SUBTYPE_MASK here, I don't consider
    # it all that useful. If somebody has a strong case why it's useful to have
    # that extra field, prlease submit a PR!
    printfield(io,FILETYPES[h.filetype],13)
    printfield(io,dec(h.ncmds),6)
    printfield(io,dec(h.sizeofcmds),11)
    printflags(io,h.flags)
end

const VM_PROT_READ      =   0x01 # read permission
const VM_PROT_WRITE     =   0x02 # write permission
const VM_PROT_EXECUTE   =   0x04 # execute permission

protstr(prot) = string(
    (prot & VM_PROT_READ)    > 0 ? "r" : "",
    (prot & VM_PROT_WRITE)   > 0 ? "w" : "",
    (prot & VM_PROT_EXECUTE) > 0 ? "x" : "")

function segflags(flags)
    if flags == 0
        return "(none)"
    else
        return "(unimplemented)"
    end
end

#     printentry(io,"offset",dec(l.offset))
#    printentry(io,"align","2^",l.align," (",2^l.align,")")

printentry(io::IO,header,values...) = (printfield(io,header,15);println(io," ",values...))
# Showing Load commands
function show(io::IO,l::Union{segment_command_64,segment_command})
    println(io,"Load Command (",isa(l,segment_command_64)?"SEGMENT_64":
        "SEGMENT","):")

    printentry(io,"name",l.segname)
    printentry(io,"addr",l.vmaddr)
    printentry(io,"size",l.vmsize)
    printentry(io,"fileoff",dec(l.fileoff))
    printentry(io,"filesize",dec(l.filesize))
    printentry(io,"maxprot",protstr(l.maxprot))
    printentry(io,"initprot",protstr(l.initprot))
    printentry(io,"nsects",dec(l.nsects))
    printentry(io,"flags",segflags(l.flags))
end

function sattrs(attributes)
    strings = String[]
    for (k,v) in SECATTRS
        if attributes & k > 0
            push!(strings, v)
        end
    end
    join(strings,",")
end

function show(io::IO, s::Union{section,section_64})
    #function body
    println(io,"  Section:")

    printentry(io,"sectname",s.sectname)
    printentry(io,"segname",s.segname)
    printentry(io,"addr","0x",hex(s.addr,2*sizeof(s.addr)))
    printentry(io,"size","0x",hex(s.size,2*sizeof(s.addr)))
    printentry(io,"offset",dec(s.offset))
    printentry(io,"align","2^",s.align," (",2^s.align,")")
    printentry(io,"reloff",s.reloff)
    printentry(io,"nreloc",s.nreloc)
    printentry(io,"type",SECTYPES[s.flags&SECTION_TYPE])
    printentry(io,"attributes",sattrs(s.flags&SECTION_ATTRIBUTES))
end

function show(io::IO,l::symtab_command)
    println(io,"Load Command (SYMTAB):")
    printentry(io,"symoff",dec(l.symoff))
    printentry(io,"nsyms",dec(l.nsyms))
    printentry(io,"stroff",dec(l.stroff))
    printentry(io,"strsize",dec(l.strsize))
end

function show(io::IO,l::dysymtab_command)
    println(io,"Load Command (DYSYMTAB):")
    printentry(io,"ilocalsym",dec(l.ilocalsym))
    printentry(io,"nlocalsym",dec(l.nlocalsym))
    printentry(io,"iextdefsym",dec(l.iextdefsym))
    printentry(io,"nextdefsym",dec(l.nextdefsym))
    printentry(io,"iundefsym",dec(l.iundefsym))
    printentry(io,"nundefsym",dec(l.nundefsym))
    printentry(io,"tocoff",dec(l.tocoff))
    printentry(io,"ntoc",dec(l.ntoc))
    printentry(io,"modtaboff",dec(l.modtaboff))
    printentry(io,"nmodtab",dec(l.nmodtab))
    printentry(io,"extrefsymoff",dec(l.extrefsymoff))
    printentry(io,"nextrefsyms",dec(l.nextrefsyms))
    printentry(io,"indirectsymoff",dec(l.indirectsymoff))
    printentry(io,"nindirectsyms",dec(l.nindirectsyms))
    printentry(io,"extreloff",dec(l.extreloff))
    printentry(io,"nextrel",dec(l.nextrel))
    printentry(io,"locreloff",dec(l.locreloff))
    printentry(io,"nlocrel",dec(l.nlocrel))
end

function decodeversion(io::IO, version::UInt32)
    print(io,version>>16)
    print(io,'.')
    print(io,(version>>8)&0xFF)
    print(io,'.')
    print(io,version&0xFF)
end

function show(io::IO,l::version_min_macosx_command)
    println(io,"Load Command (VERSION_MIN_MACOSX):")
    printentry(io,"version",sprint(decodeversion,l.version))
    printentry(io,"sdk",sprint(decodeversion,l.sdk))
end

################################ Interface #####################################

import Base: show,
    # IO methods
    read, write, seek, seekstart, position,
    # Iteration
    start, next, done,
    # Indexing
    length, getindex

#
# Note that this function is different from ObjFileBase.readmeta
# Constructs and initializes the MachOHandle object
#
function readmeta(io::IO,::Type{MachOHandle})
    start = position(io)
    magic = read(io,UInt32)
    if magic == MH_MAGIC
        return MachOHandle(io,start,false,false)
    elseif magic == MH_CIGAM
        return MachOHandle(io,start,true,false)
    elseif magic == MH_MAGIC_64
        return MachOHandle(io,start,false,true)
    elseif magic == MH_CIGAM_64
        return MachOHandle(io,start,true,true)
    elseif magic == FAT_CIGAM || magic == FAT_MAGIC
        return FatMachOHandle(io,start)
    else
        throw(ObjFileBase.MagicMismatch("Invalid Magic ($(hex(magic)))!"))
    end
end

function readloadcmd(h::MachOHandle)
    cmd = unpack(h,load_command)
    ccmd = cmd.cmd & ~LC_REQ_DYLD
    if ccmd == LC_UUID
        return (cmd,unpack(h, uuid_command))
    elseif ccmd == LC_SEGMENT
        return (cmd,unpack(h, segment_command))
    elseif ccmd == LC_SEGMENT_64
        return (cmd,unpack(h, segment_command_64))
    elseif ccmd == LC_SYMTAB
        return (cmd,unpack(h, symtab_command))
    elseif ccmd == LC_DYSYMTAB
        return (cmd,unpack(h, dysymtab_command))
    elseif ccmd == LC_VERSION_MIN_MACOSX
        return (cmd,unpack(h, version_min_macosx_command))
    elseif ccmd == LC_ID_DYLIB || ccmd == LC_LOAD_DYLIB ||
            ccmd == LC_REEXPORT_DYLIB || ccmd == LC_LOAD_UPWARD_DYLIB ||
            ccmd == LC_LOAD_WEAK_DYLIB
        return (cmd,unpack(h, dylib_command, cmd.cmdsize))
    elseif ccmd == LC_ID_DYLINKER || ccmd == LC_LOAD_DYLINKER ||
            ccmd == LC_DYLD_ENVIRONMENT
        return (cmd,unpack(h, dylinker_command, cmd.cmdsize))
    elseif ccmd == LC_DYLD_INFO
        return (cmd,unpack(h, dyld_info_command))
    elseif ccmd == LC_SOURCE_VERSION
        return (cmd,unpack(h, source_version_command))
    elseif ccmd == LC_CODE_SIGNATURE || ccmd == LC_SEGMENT_SPLIT_INFO ||
            ccmd == LC_FUNCTION_STARTS || ccmd == LC_DATA_IN_CODE ||
            ccmd == LC_DYLIB_CODE_SIGN_DRS
        return (cmd,unpack(h, linkedit_data_command))
    elseif ccmd == LC_SUB_FRAMEWORK
        return (cmd,unpack(h, sub_framework_command, cmd.cmdsize))
    elseif ccmd == LC_RPATH
        return (cmd,unpack(h, rpath_command, cmd.cmdsize))
    elseif ccmd == LC_UNIXTHREAD
        flavor = read(h, UInt32)
        count = read(h, UInt32)
        data = read(h, UInt, div(cmd.cmdsize - 4sizeof(UInt32),sizeof(UInt)))
        return (cmd,thread_command(flavor, count, data))
    elseif ccmd == LC_MAIN
        return (cmd,unpack(h, entry_point_command))
    elseif ccmd == LC_ROUTINES_64
        return (cmd,unpack(h, routines_command_64))
    elseif ccmd == LC_SUB_CLIENT
        return (cmd,unpack(h, sub_client_command, cmd.cmdsize))
    else
        info("Unimplemented load command $(LCTYPES[ccmd]) (0x$(hex(ccmd)))")
        return (cmd,dummy_lc())
    end
end

## Iteration

import Base: eltype

# Iterate over the load commands of a file
# If the header has
# already been read, it may be passed to the constructor. Optionally
# the start address of the load commands may also be passed in
# passed in. If a header is passed, it is assumed that the start address
# is the current position of the stream. Otherwise, it will be determined when
# reading the header.
immutable LoadCmds
    h::MachOHandle
    start::UInt64
    ncmds::UInt32
    sizeofcmds::UInt32
end

immutable LoadCmd{T<:MachOLC}
    h::MachOHandle
    off::UInt64
    cmd_id::UInt32
    cmd::T
end
deref(lc::LoadCmd) = lc.cmd
handle(lc::LoadCmd) = lc.h

function children(lc::LoadCmd)
    isa(lc.cmd, segment_commands) ? Sections(lc) : ()
end

show{T}(io::IO, x::LoadCmd{T}) = (print(io,"0x",hex(x.off,8),":\n "); show(io,x.cmd); print(io,'\n'))
eltype{T<:MachOLC}(::LoadCmd{T}) = T

function LoadCmds(h::MachOHandle, header = nothing, start = -1)
    if header === nothing
        header = readheader(h)
    end
    if start == -1
        start = position(h)
    end
    LoadCmds(h,start,header.ncmds,header.sizeofcmds)
end
children(h::MachOHandle) = LoadCmds(h)

# A tuple of the position before the current load command,
# the number of the load command and the size of the current load
# command. I.e. the next load command will be found at state[1]+state[3]
start(l::LoadCmds) = (l.start,0,0)
seek(l::LoadCmds,state) = seek(l.h,state[1]+state[3])
function next(l::LoadCmds,state)
    seek(l,state)
    cmdh,cmd = readloadcmd(l.h)
    (LoadCmd(l.h,state[1]+state[3],cmdh.cmd,cmd),(state[1]+state[3],state[2]+1,cmdh.cmdsize))
end
length(l::LoadCmds) = l.ncmds
done(l::LoadCmds,state) = state[2] >= length(l)

# Access to sections

const segment_commands = Union{segment_command_64,segment_command}

immutable Sections <: ObjFileBase.Sections{MachOHandle}
    h::MachOHandle
    command::segment_commands
    start::Int
    function Sections(h::MachOHandle, segment::segment_commands, start=-1)
        if start==-1
            start = position(h)
        end
        new(h,segment,start)
    end
    function Sections(segment::LoadCmd)
        if !(eltype(segment) <: segment_commands)
            error("Load command is not a segment")
        end
        start = segment.off + sizeof(eltype(segment)) + sizeof(load_command)
        new(segment.h,segment.cmd,start)
    end
end
ObjFileBase.handle(sects::Sections) = sects.h
ObjFileBase.Sections(segment::LoadCmd) = Sections(segment)
ObjFileBase.Sections(h::MachOHandle, args...) = Sections(h,args...)
ObjFileBase.mangle_sname(h::MachOHandle, name) = string("__", name)

# Covenience function for the (common) case where the object
# only has one segment
function Sections(h::MachOHandle)
    seglcs = collect(filter(LoadCmds(h)) do lc
        isa(lc.cmd, segment_commands)
    end)
    # @assert length(seglcs) == 1
    Sections(seglcs[])
end

immutable SectionRef{H<:MachOHandle, T<:Union{section,section_64}} <: ObjFileBase.SectionRef{MachOHandle}
    handle::MachOHandle
    header::T
end
@Base.pure ObjFileBase.SectionRef{H<:MachOHandle}(::Type{H}) = SectionRef{H}
deref(x::SectionRef) = x.header
handle(x::SectionRef) = x.handle
sectionname(x::SectionRef) = sectionname(deref(x))
sectionaddress(x::SectionRef) = sectionaddress(deref(x))

length(s::Sections) = s.command.nsects
function getindex(s::Sections,n)
    if n < 1 || n > length(s)
        throw(BoundsError())
    end
    sT = isa(s.command,segment_command_64) ? section_64 : section
    seek(s.h,s.start + (n-1)*sizeof(sT))
    SectionRef{typeof(s.h),sT}(s.h,unpack(s.h, sT))
end

start(s::Sections) = 1
done(s::Sections,n) = n > length(s)
next(s::Sections,n) = (s[n],n+1)

for f in (:readuntil,:write)
    @eval $(f){T<:IO}(io::MachOHandle{T},args...) = $(f)(io.io,args...)
end
read{T<:IO}(io::MachOHandle{T},num::Integer) = read(io.io,num)

seek{T<:IO}(io::MachOHandle{T},pos::Integer) = seek(io.io,io.start+pos)
seekstart(io::MachOHandle) = seek(io.io,io.start)
position{T<:IO}(io::MachOHandle{T}) = position(io.io)-io.start

unpack{T,ioT<:IO}(h::MachOHandle{ioT},::Type{T}) =
    unpack(h.io,T,h.bswapped ? :SwappedEndian : :NativeEndian)

# Access to relocations
immutable Relocations
    sec::SectionRef
end

immutable RelocationRef <: ObjFileBase.RelocationRef{MachOHandle}
    h::MachOHandle
    reloc::relocation_info
end

deref(x::RelocationRef) = x.reloc

entrysize(s::Relocations) = sizeof(relocation_info)
endof(s::Relocations) = s.sec.header.nreloc
length(r::Relocations) = endof(r)
function getindex(s::Relocations,n)
    if n < 1 || n > length(s)
        throw(BoundsError())
    end
    offset = s.sec.header.reloff +
        (n-1)*entrysize(s)
    seek(s.sec.handle,offset)
    RelocationRef(s.sec.handle,unpack(s.sec.handle, relocation_info))
end

start(s::Relocations) = 1
done(s::Relocations,n) = n > length(s)
next(s::Relocations,n) = (x=s[n];(x,n+1))

function unpack{ioT<:IO}(h::MachOHandle{ioT}, ::Type{UInt32})
    return unpack(h.io, UInt32, endianness(h))
end

pack{T,ioT<:IO}(h::MachOHandle{ioT},::Type{T}) =
    pack(h.io,T,h.bswapped ? :SwappedEndian : :NativeEndian)

function readheader(h::MachOHandle)
    seekstart(h)
    unpack(h,h.is64 ? mach_header_64 : mach_header)
end

ObjFileBase.isrelocatable(handle::MachOHandle) =
    readheader(handle).filetype == MachO.MH_OBJECT

sectionsize(sect::Union{section,section_64}) = sect.size
# In non-relocatable files, address and offset are the same and offset is 0
sectionoffset(sect::Union{section,section_64}) = sect.offset == 0 ? sect.addr : sect.offset

### Access to Symbols

immutable Symbols <: ObjFileBase.Symbols{MachOHandle}
    lc::LoadCmd{symtab_command}
end
ObjFileBase.Symbols(lc::LoadCmd) = Symbols(lc)
ObjFileBase.Symbols(h::MachOHandle) = Symbols(h)
symname(syms::Symbols, sym) = symname(syms.lc, sym)
ObjFileBase.StrTab(s::Symbols) = ObjFileBase.StrTab(s.lc)

immutable SymbolRef <: ObjFileBase.SymbolRef{MachOHandle}
    symbols::Symbols
    num::UInt32
    entry::ObjFileBase.SymtabEntry{MachOHandle}
end
deref(x::SymbolRef) = x.entry
symbolnum(x::SymbolRef) = x.num
symbolvalue(x::SymbolRef, args...) = deref(x).n_value
handle(x::SymbolRef) = x.handle

isglobal(x) = (x.n_type & N_EXT) != 0
islocal(x) = !isglobal(x) && ((x.n_type & N_UNDF) == 0)
isweak(x) = (x.n_desc & (N_WEAK_REF | N_WEAK_DEF)) != 0
isdebug(x) = (x.n_type & N_TYPE) == N_STAB
isundef(x::ObjFileBase.SymtabEntry{MachOHandle}) = (x.n_type == N_UNDF) || ((x.n_type & N_EXT) != 0 && (x.n_sect == NO_SECT))

# Symbol printing stuff
function showcompact(io::IO, x::SymbolRef)
    print(io,'[')
    printfield(io,dec(symbolnum(x)),5)
    print(io,"] ")
    showcompact(io, x.entry; syms = x.symbols, sections = Sections(handle(x.symbols.lc)))
end

# Try to follow the same format as llvm-objdump
function showcompact(io::IO,x::ObjFileBase.SymtabEntry{MachOHandle};
        syms = nothing, sections = nothing)
    # Value
    print(io,string("0x",hex(x.n_value,2*sizeof(x.n_value))))
    print(io," ")

    # Symbol flags
    print(io, isglobal(x) ? "g" : islocal(x) ? "l" : "-")
    print(io, isweak(x) ? "w" : "-")
    print(io, "-"^3) # Unsupported
    print(io, isdebug(x) ? "d" : "-")

    # Skip Symbol type (TODO)
    print(io, " ")

    print(io, " ")
    if (x.n_type & N_TYPE) == N_UNDF && isglobal(x)
            printfield(io,
                x.n_value != 0 ? "*COM*" : "*UND*",
                20; align = :left)
    elseif (x.n_type & N_TYPE) == N_ABS
        printfield(io,"*ABS*",20; align = :left)
    elseif sections !== nothing
        printfield(io, sectionname(sections[x.n_sect]), 20; align = :left)
    else
        printfield(io, "Section #$(x.n_sect)", 20; align = :left)
    end
    print(io, " ")

    if syms !== nothing
        print(io,symname(syms, x))
    else
        print(io,"@",x.n_strx)
    end
end

function Symbols(oh::MachOHandle)
    symlcs = [lc for lc in LoadCmds(oh) if isa(lc.cmd, symtab_command)]
    @assert length(symlcs) == 1
    return Symbols(symlcs[1])
end

start(s::Symbols) = 1
done(s::Symbols,n) = n > length(s)
next(s::Symbols,n) = (s[n],n+1)

length(s::Symbols) = s.lc.cmd.nsyms
function getindex(s::Symbols,n)
    if n < 1 || n > length(s)
        throw(BoundsError())
    end
    sT = s.lc.h.is64 ? nlist_64 : nlist
    seek(s.lc.h,s.lc.cmd.symoff + (n-1)*sizeof(sT))
    SymbolRef(s, n-1, unpack(s.lc.h, sT))
end

immutable SectionStrTab <: ObjFileBase.StrTab
    strtab::SectionRef
end
immutable OffsetStrTab <: ObjFileBase.StrTab
    h::MachOHandle
    offset::UInt
    size::UInt
end
ObjFileBase.StrTab(s::SectionRef) = SectionStrTab(s)
ObjFileBase.StrTab(s::LoadCmd{symtab_command}) =
    OffsetStrTab(handle(s),deref(s).stroff,deref(s).strsize)

function strtab_lookup(s::SectionStrTab,index)
    seek(s.strtab.handle,sectionoffset(s.strtab)+index)
    strip(readuntil(s.strtab.handle,'\0'),'\0')
end
function strtab_lookup(s::OffsetStrTab,index)
    @assert index < s.size
    seek(s.h,s.offset+index)
    strip(readuntil(s.h,'\0'),'\0')
end
load_strtab(strtab::SectionRef) = SectionStrTab(strtab)

function strtable_lookup(io::MachOHandle,command::symtab_command,offset)
    seek(io,command.stroff+offset)
    strip(readuntil(io,'\0'),'\0')
end

symname(io::MachOHandle,command::symtab_command,sym) = strtable_lookup(io, command, sym.n_strx)
symname(sym::SymbolRef; strtab = ObjFileBase.StrTab(sym.symbols), kwargs...) = strtab_lookup(strtab, deref(sym).n_strx)
symname(x::LoadCmd{symtab_command}, sym) = symname(x.h, x.cmd, sym)
segname(x::Union{segment_command_64,section_64}) = x.segname
segname(x::LoadCmd{segment_command_64}) = segname(x.cmd)
sectionname(x::section_64) = x.sectname
sectionaddress(x::section_64) = x.addr

### Fat Handle
immutable FatMachOHandle
    io::IO
    start::UInt64
    header::fat_header
    archs::Vector{fat_arch}
end

function FatMachOHandle(io,start)
    header = unpack(io,fat_header,:BigEndian)
    archs = Array(fat_arch,header.nfat_arch)
    for i in 1:header.nfat_arch
        archs[i] = unpack(io,fat_arch,:BigEndian)
    end
    FatMachOHandle(io,start,header,archs)
end

function show(io::IO,h::fat_arch)
    printentry(io,"cputype",CPUTYPES[h.cputype])
    # TODO: subtype printing
    #printentry(io,"cputype",CPUSUBTYPES[h.cputype])
    printentry(io,"offset","0x",hex(h.offset,2sizeof(UInt32)))
    printentry(io,"size","0x",hex(h.size,2sizeof(UInt32)))
    printentry(io,"align",dec(h.align))
end

function show(io::IO,h::FatMachOHandle)
    print(io, "Fat Mach Handle (",length(h.archs), " architectures)")
    for (i,arch) in enumerate(h.archs)
        println(io)
        println(io, "architecture ", i)
        print(io, arch)
    end
end

function getindex(h::FatMachOHandle,i)
    seek(h.io,h.start + h.archs[i].offset)
    readmeta(h.io,MachOHandle)
end

### Compact Unwind Support

include("compact_unwind.jl")

### Relocation Support

include("relocate.jl")

### DWARF support

function debugsections{T<:segment_commands}(seg::LoadCmd{T})
    sects = collect(Sections(seg))
    snames = map(sectionname,sects)
    sections = Dict{String,SectionRef}()
    for i in 1:length(snames)
        # remove leading "__"
        ind = findfirst(ObjFileBase.DEBUG_SECTIONS,unsafe_string(snames[i])[3:end])
        if ind != 0
            sections[ObjFileBase.DEBUG_SECTIONS[ind]] = sects[i]
        end
    end
    sections
end

function debugsections(oh::MachOHandle)
    segs = collect(filter(LoadCmds(oh)) do lc
        isa(lc.cmd, segment_commands) &&
            (segname(lc) == "__DWARF" || isempty(segname(lc)))
    end)
    isempty(segs) && error("No debug sections present")
    @assert length(segs) == 1
    ObjFileBase.DebugSections(oh,debugsections(first(segs)))
end

include("precompile.jl")

end # module


text:
using PiGPIO

red_pin = 18 #change this accordingly with your GPIO pin number
p=Pi()
set_mode(p, red_pin, PiGPIO.OUTPUT)
try
    for i in 1:10
        PiGPIO.write(p, red_pin, PiGPIO.HIGH)
        sleep(0.5)
        PiGPIO.write(p, red_pin, PiGPIO.LOW)
        sleep(0.5)
    end
finally
    println("Cleaning up!")
    set_mode(p, red_pin, PiGPIO.INPUT)
end


text:
# Make sure filtration building is also benchmarked.
function _ripserer(data, complex, complex_kwargs, ripserer_kwargs)
    return ripserer(complex(data; complex_kwargs...); ripserer_kwargs...)
end

function _ripserer(filename::String, complex, complex_kwargs, ripserer_kwargs)
    data = load_data(filename)
    return ripserer(complex(data; complex_kwargs...); ripserer_kwargs...)
end

function ripserer_benchmark(
    b::Benchmark; extras=false, from_file=false, sparse=nothing, kwargs...
)
    if !isnothing(b.threshold)
        complex_kwargs = (;threshold=b.threshold)
    else
        complex_kwargs = (;)
    end
    if extras
        complex_kwargs = (;complex_kwargs..., b.extra_kwargs...)
    end
    if !isnothing(sparse)
        complex_kwargs = (;complex_kwargs..., sparse=sparse)
    end
    if from_file
        input = b.filename
    else
        input = b.data
    end
    kwargs = (;kwargs..., b.kwargs...)
    return @benchmarkable _ripserer(
        $(input), $(b.complex), $(complex_kwargs), $(kwargs)
    ) samples=5 seconds=1000 gcsample=true
end


text:
using Pkg
cd(@__DIR__)
Pkg.activate("../"); Pkg.instantiate()

using Plots; gr()
using LaTeXStrings
using DataFrames
using CSV
using OMJulia
using OrdinaryDiffEq
using LinearAlgebra
using ControlSystems
using DataInterpolations
using BlackBoxOptim
using BenchmarkTools

# Plot path
figpath = "../figs/";
# Linewidths and styles
LW1 = 2.5
LW2 = 1.5
LS1 = :solid
LS2 = :dot
LS3 = :dash
LC1 = :blue
LC2 = :red
LC3 = :green
LC4 = :black
LA1 = 1.0
LA2 = 0.7
LA3 = 0.3
;

# Function for helicopter model
function helicopter_ode(dx,x,Ps,t,P0,u)
# P0 = (rc=rc,h=h,m=m,Jbc=Jbc,Js=Js,g=g,K_th_th=K_th_th,K_th_psi=K_th_psi,K_psi_th=K_psi_th,K_psi_psi=K_psi_psi,D_th=D_th,D_psi=D_psi);
    th = x[1]
    psi = x[2]
    w_th = x[3]
    w_psi = x[4]
    #
    M11 = Ps[4]*P0.Jbc + Ps[3]*P0.m*((Ps[1]*P0.rc)^2+(Ps[2]*P0.h)^2)
    M22 = Ps[4]*P0.Jbc*cos(th)^2 + Ps[3]*P0.m*(Ps[1]*P0.rc*cos(th)-Ps[2]*P0.h*sin(th))^2 + Ps[5]*P0.Js
    #
    dLdth = -w_psi^2*(((Ps[4]*P0.Jbc)+(Ps[3]*P0.m)*((Ps[1]*P0.rc)^2-(Ps[2]*P0.h)^2))*cos(th)*sin(th)-(Ps[1]*P0.rc)*(Ps[2]*P0.h)*(cos(th)^2-sin(th)^2)) - (Ps[3]*P0.m)*(Ps[6]*P0.g)*((Ps[1]*P0.rc)*cos(th) - (Ps[2]*P0.h)*sin(th))
    #
    dp_psi_dth = -2*(((Ps[4]*P0.Jbc)+(Ps[3]*P0.m)*((Ps[1]*P0.rc)^2-(Ps[2]*P0.h)^2))*cos(th)*sin(th) + (Ps[3]*P0.m)*(Ps[1]*P0.rc)*(Ps[2]*P0.h)*(cos(th)^2-sin(th)^2))*w_psi
    #
    F_th = Ps[7]*P0.K_th_th*u.th(t) - Ps[8]*P0.K_th_psi*u.psi(t) - P0.D_th*w_th
    F_psi = Ps[9]*P0.K_psi_th*u.th(t) - Ps[10]*P0.K_psi_psi*u.psi(t) - P0.D_psi*w_psi
    #
    dx[1] = w_th
    dx[2] = w_psi
    dx[3] = (dLdth + F_th)/M11
    dx[4] = (-dp_psi_dth*w_th + F_psi)/M22
end
#
# Converting from radians to degrees in plots
#
rad2deg(t,u) = (t,u*180/pi)
;

# Model parameters
rc = 1.5e-2 # m
h = 0.5e-2 # m
m = 0.5 # kg
Jbc = 15e-3 # kg.m2
Js = 5e-3 # kg.m2
g = 9.81 # m/s2
K_th_th = 55e-3 # Nm/V
K_th_psi = 5e-3 # Nm/V
K_psi_th = 0.15 # Nm/V
K_psi_psi = 0.2 # Nm/V
D_th = 1e-2 # Nm/(rad/s)
D_psi = 8e-2 # Nm/(rad/s)
#
P0 = (rc=rc,h=h,m=m,Jbc=Jbc,Js=Js,g=g,K_th_th=K_th_th,K_th_psi=K_th_psi,K_psi_th=K_psi_th,K_psi_psi=K_psi_psi,D_th=D_th,D_psi=D_psi)
Ps = fill(1.0,length(P0))
;

# Set-up for simulating helicopter
#
th0 = pi/2
psi0 = 0.0
w_th0 = -0.1
w_psi0 = 0.1
#
x0 = [th0,psi0,w_th0,w_psi0]
u_c = (th=t->0.0, psi=t->0.0)
#
helicopter! = (dx,x,Ps,t) -> helicopter_ode(dx,x,Ps,t,P0,u_c)
#
tspan = (0.,20)
prob = ODEProblem(helicopter!,x0,tspan,Ps)
#
sol = solve(prob,AutoTsit5(Rosenbrock23()),reltol=1e-6)
#
plot(sol,vars=(rad2deg,0,1),lw=LW1,lc=LC1,label=L"$\theta$ [deg]")
plot!(sol,vars=(rad2deg,0,2),lw=LW1,lc=LC2,label=L"$\psi$ [deg]")
plot!(xlim = tspan,xlabel=L"time $t$ [s]")
figname="Helicopter_theoretical_simulation.svg"
savefig(figpath*figname)

df=CSV.read("../data/Lab-Helicopter_Experimental-data.csv")
ENV["COLUMNS"]=100 # change the number of IJulia columns from default 80 to 100
df[1:5,:]

u_th = df[:,1]
u_psi = df[:,2]
theta = df[:,3]
psi = df[:,4]
dt = df[1,5]
tm = 0:dt:(length(psi)-1)*dt
;

data = [theta';psi']
Nd = length(tm);

plot(tm,theta*180/pi,lw=LW2,lc=LC1,label=L"\theta",legend=:topleft)
plot!(twinx(),tm,psi*180/pi,lw=LW2,lc=LC2,label=L"\psi")
plot!(title="Experimental helicopter outputs: measured angles", ylabel=L"angles [$^\circ$]", xlabel=L"time $t$ [s]", box=true)
figname="Helicopter_experimental_angles.svg"
savefig(figpath*figname)

plot(tm,u_th,lw=LW1,lc=LC1,label=L"u_\theta",legend=:topleft)
plot!(twinx(),tm,u_psi,lw=LW1,lc=LC2,label=L"u_\psi",legend=:bottomright)
plot!(title="Experimental helicopter inputs: motor voltages", ylabel=L"$u$ [V]", xlabel=L"time $t$ [s]",box=true)
figname="Helicopter_experimental_voltages.svg"
savefig(figpath*figname)

U_th = ConstantInterpolation(u_th,tm)
U_psi = ConstantInterpolation(u_psi,tm)
Theta = ConstantInterpolation(theta,tm)
Psi = ConstantInterpolation(psi,tm);

# Set-up for simulating helicopter
#
th0 = Theta(0.0)
psi0 = Psi(0.0)
w_th0 = 0.0
w_psi0 = 0.0
#
x0 = [th0,psi0,w_th0,w_psi0]
u_t = (th=U_th, psi=U_psi)
#
helicopter! = (dx,x,Ps,t) -> helicopter_ode(dx,x,Ps,t,P0,u_t)
#
tspan = (tm[1],tm[end])
#
prob = ODEProblem(helicopter!,x0,tspan,Ps)
#
sol = solve(prob,AutoTsit5(Rosenbrock23()),reltol=1e-6)
#
plot(sol,vars=(rad2deg,0,1),lw=LW1,lc=LC1,label=L"$\theta$ [deg]")
plot!(sol,vars=(rad2deg,0,2),lw=LW1,lc=LC2,label=L"$\psi$ [deg]")
plot!(xlim = tspan,xlabel=L"time $t$ [s]", ylabel="angles")
plot!(title="Angles: experimental inputs, nominal parameters",box=true)
figname="Helicopter_simulated_angles_experimental_inputs_nominal.svg"
savefig(figpath*figname)

function cost(p,Ps,pidx,x0,xidx,data,Nd,tspan)
    #
    xx0 = copy(x0)
    xx0[xidx] = p[xidx].*x0[xidx]
    Ps[pidx] .= p[length(xidx)+1:end]
    #
    prob = ODEProblem(helicopter!,xx0,tspan,Ps)
    sol = solve(prob, AutoTsit5(Rosenbrock23()), saveat=dt)
    #
    return norm(sol[1:2,:]-data,2)^2/Nd
    #
end
#
;

#
# (:rc, :h, :m, :Jbc, :Js, :g, :K_th_th, :K_th_psi, :K_psi_th, :K_psi_psi, :D_th, :D_psi)
Ps_keys = [:rc, :h, :m, :Jbc, :Js, :K_th_th,:K_th_psi,:K_psi_th,:K_psi_psi,:D_th, :D_psi]
pidx = [findall(Ps_keys[i] .== keys(P0))[1] for i in 1:length(Ps_keys)]
#
xidx = collect(3:4)
#
p = fill(1.0,length(xidx)+length(pidx))
#
p_lo = 0.2*p
p_hi = 2*p
p_lo_hi = collect(zip(p_lo,p_hi))
#
# Initial value of loss function
Ps = fill(1.0,length(P0))
x0 = [Theta(0),Psi(0),1.0,1.0]
#
loss = (p) -> cost(p,Ps,pidx,x0,xidx,data,Nd,tspan)
;

loss(p)

@time res = bboptimize(loss; SearchRange = p_lo_hi, NumDimensions = length(p));

best_fitness(res)

p_best = best_candidate(res)

loss(p_best)

x0_best = copy(x0)
x0_best[xidx] = p_best[xidx].*x0[xidx]
Ps_best=fill(1.0,length(P0))
Ps_best[pidx] .= p_best[length(xidx)+1:end]
#
prob_best = ODEProblem(helicopter!,x0_best,tspan,Ps_best)
#
sol_best = solve(prob_best,AutoTsit5(Rosenbrock23()),reltol=1e-6)
#
plot(sol_best,vars=(rad2deg,0,1),lw=LW1,lc=LC1,label=L"$\theta$ [deg]")
plot!(tm,theta*180/pi,lw=LW1,ls=LS2,lc=LC2,label=L"\theta^\mathrm{d}")
plot!(xlim = tspan,xlabel=L"time $t$ [s]", ylabel=L"$\theta$ [$^\circ$]")
plot!(title="Pitch angle: model (blue, solid) vs. data (red, dotted)",box=true)
figname="Helicopter_pitch-angle_model_fit.svg"
savefig(figpath*figname)

plot(sol_best,vars=(rad2deg,0,2),lw=LW1,lc=LC1,label=L"$\psi$ [deg]")
plot!(tm,psi*180/pi,lw=LW1,ls=LS2,lc=LC2,label=L"\psi^\mathrm{d}")
plot!(xlim = tspan,xlabel=L"time $t$ [s]", ylabel=L"$\psi$ [$^\circ$]")
plot!(title="Yaw angle: model (blue, solid) vs. data (red, dotted)",box=true)
figname="Helicopter_yaw-angle_model_fit.svg"
savefig(figpath*figname)

# Estimated initial states -- states 3 and 4 have been estimated
x0_best[xidx] = p_best[xidx].*x0[xidx]
show(x0_best)

# Estimated model parameters -- all parameters except gravity g have been assumed uncertain
Ps_best=fill(1.0,length(P0))
Ps_best[pidx] .= p_best[length(xidx)+1:end]
collect(zip(keys(P0),Ps_best.*values(P0)))

# Nominal parameters, taken from laboratory task description (and rounded somewhat)
P0

open("../optimization_results/globaloptimizationresults.jl", "w") do f
    write(f,"P0 = $P0 \n")
    write(f,"Ps_best = $Ps_best \n")
    write(f,"x0_best = $x0_best \n")
    write(f,"tspan = $tspan \n")
end


text:
using Sims
using NineML
using Winston

al = NineML.parse_ComponentClass("IaF.xml")
uls = NineML.parse_Components("IaF_UL.xml")

tf = 100.0
dt = 0.025

ul= uls[1]
NineML.eval_Dynamics(al.dynamics, ul.name, ul.parameterValues, ul.initialValues, al.parameters, al.states, al.ports)
model = ul.name

iaf = eval(quote $model () end)
iaf_s = create_sim(iaf)
iaf_ptr = setup_sunsim(iaf_s)

iaf_yout = sunsim(iaf_ptr, tf, int(tf/dt))
f = figure()
pl = plot (iaf_yout.y[:,1], iaf_yout.y[:,2])
display(pl)


text:
using Distributed
import Nanosoldier, GitHub

nodes = Dict(Any => addprocs(1))
@everywhere import Nanosoldier

auth = GitHub.authenticate(ENV["GITHUB_AUTH"])
secret = ENV["GITHUB_SECRET"]

config = Nanosoldier.Config(ENV["USER"], nodes, auth, secret;
                            workdir = joinpath(homedir(), "test_workdir"),
                            trackrepo = "vtjnash/julia",
                            reportrepo = "vtjnash/NanosoldierReports",
                            admin = "vtjnash",
                            testmode = true)

server = Nanosoldier.Server(config)


text:
struct ConditionalModel{A,B,M,Argvals,Obs} <: AbstractModel{A,B,M,Argvals,Obs}
    model :: Model{A,B,M}
    argvals :: Argvals
    obs :: Obs
end

function Base.show(io::IO, cm::ConditionalModel)
    println(io, "ConditionalModel given")
    println(io, "    arguments    ", keys(argvals(cm)))
    println(io, "    observations ", keys(observations(cm)))
    println(io, Model(cm))
end

export argvals
argvals(c::ConditionalModel) = c.argvals
argvals(m::Model) = NamedTuple()

export observations
observations(c::ConditionalModel) = c.obs

export observed
function observed(cm::ConditionalModel{A,B,M,Argvals,Obs}) where {A,B,M,Argvals,Obs}
    keys(schema(Obs))
end

Model(c::ConditionalModel) = c.model

Model(::Type{<:ConditionalModel{A,B,M}}) where {A,B,M} = type2model(Model{A,B,M})

ConditionalModel(m::Model) = ConditionalModel(m,NamedTuple(), NamedTuple())

(m::Model)(nt::NamedTuple) = ConditionalModel(m)(nt)

(cm::ConditionalModel)(nt::NamedTuple) = ConditionalModel(cm.model, merge(cm.argvals, nt), cm.obs)

(m::Model)(;argvals...)= m((;argvals...))

(m::Model)(args...) = m(NamedTuple{Tuple(m.args)}(args))

import Base

Base.:|(m::Model, nt::NamedTuple) = ConditionalModel(m) | nt

Base.:|(cm::ConditionalModel, nt::NamedTuple) = ConditionalModel(cm.model, cm.argvals, merge(cm.obs, nt))


text:
if Compat.Sys.iswindows()
	include("installer.jl")
end
function build_app_bundle(script::String;
				resources = String[],
				libraries = String[],
				builddir = "builddir",
				appname = "nothing",
				create_installer = false)

	# Create build directory
	script = abspath(script)
	base_path = dirname(script)
	builddir = joinpath(base_path, builddir, appname)
	@info "Building at path $builddir"
	mkpath(builddir)

	core_path = joinpath(builddir, "core")
	lib_path = joinpath(builddir, "lib")
	res_path = joinpath(builddir, "res")

	# Create resources and libraries dirctories
	mkpath(core_path)
	mkpath(lib_path)
	mkpath(res_path)

	delim = Compat.Sys.iswindows() ? '\\' : '/'

	@info "Copying resources:"
	for res in resources
		print("Copying $res...")
		Compat.cp(res, joinpath(res_path, split(res, delim)[end]), force = true)
		println("... done.")
	end

	@info "Copying libraries"
	for lib in libraries
		print("Copying $lib...")
		Compat.cp(lib, joinpath(lib_path, split(lib, delim)[end]), force = true)
		println("... done.")
	end

	build_executable(script, builddir = core_path)

    (create_installer && Compat.Sys.islinux()) && throw(error("Cannot create installer on Linux"))

    if Compat.Sys.iswindows()
        create_installer && installer(builddir, name = appname)
    end

end


text:
""" Normalize spectrum based on blaze model (also tryies to include black body function) from FITS file. """
function blaze_normalize_spectrum!(spectrum::ST) where { ST<:AbstractSpectra }
    @assert haskey(spectrum.metadata,:continuum)
    if spectrum.metadata[:normalization] == :raw
        blaze_smoothed = smooth_blaze(blaze)
        spectrum.flux ./= blaze_smoothed
        spectrum.var ./= (blaze_smoothed ) .^2
        spectrum.metadata[:normalization] = :blaze
    elseif spectrum.metadata[:normalization] == :continuum
        spectrum.flux .*= spectrum.metadata[:continuum]
        spectrum.var .*= (spectrum.metadata[:continuum] ) .^2
        spectrum.metadata[:normalization] = :blaze
    elseif spectrum.metadata[:normalization] == :blaze
        # do nothing
    else
        @error "Normalizing from " * string(spectrum.metadata.normalization) * " to blaze isn't implemented yet"
    end
    return spectrum
end

""" Normalize each spectrum based on blaze model (also tryies to include black body function) from FITS files. """
function blaze_normalize_spectra!(spectra::AS) where { ST<:AbstractSpectra, AS<:AbstractArray{ST} }
    for spectrum in spectra
        blaze_normalize_spectrum!(spectrum)
    end
    return spectra
end

""" Normalize spectrum based on continuum model from FITS file. """
function continuum_normalize_spectrum!(spectrum::ST) where { ST<:AbstractSpectra }
    @assert haskey(spectrum.metadata,:continuum)
    if spectrum.metadata[:normalization] == :raw
        blaze_smoothed = smooth_blaze(blaze)
        spectrum.flux ./= spectrum.metadata[:continuum] .* blaze_smoothed
        spectrum.var ./= (spectrum.metadata[:continuum].* blaze_smoothed ) .^2
        spectrum.metadata[:normalization] = :continuum
    elseif spectrum.metadata[:normalization] == :blaze
        spectrum.flux ./= spectrum.metadata[:continuum]
        spectrum.var ./= (spectrum.metadata[:continuum] ) .^2
        spectrum.metadata[:normalization] = :continuum
    elseif spectrum.metadata[:normalization] == :continuum
        # do nothing
    else
        @error "Normalizing from " * string(spectrum.metadata.normalization) * " to continuum isn't implemented yet"
    end
    return spectrum
end

""" Normalize each spectrum based on continuum model from FITS files. """
function continuum_normalize_spectra!(spectra::AS) where { ST<:AbstractSpectra, AS<:AbstractArray{ST} }
    for spectrum in spectra
        continuum_normalize_spectrum!(spectrum)
    end
    return spectra
end

""" `filter_line_list( linelist, inst ; λmin, λmax )`
Return DataFrame based on linelist, but filtered to only include lines between λmin and λmax.
Defaults for λmin and λmax specified are in traits for EXPRES.
"""
function filter_line_list(df::DataFrame, inst::IT ;
                            λmin::Real = default_filter_line_list_λmin, λmax::Real = default_filter_line_list_λmax ) where {
                                   IT<:EXPRES.AnyEXPRES }
   df |> @filter(λmin <= _.lambda <= λmax) |>
    #    @filter( _.lambda < 6000.0 ) |>                       # Avoid tellurics at redder wavelengths
    #    @filter( _.lambda >6157 || _.lambda < 6155  ) |>   # Avoid "line" w/ large variability
    DataFrame
end

""" `find_worst_telluric_in_each_chunk( chunk_list_timeseries, spectra )`
Return an array of the worst (smallest) telluric for each chunk at any time for the provided ChunkListTimeseries.
"""
function find_worst_telluric_in_each_chunk( clt::AbstractChunkListTimeseries, data::AbstractArray{AS,1} )  where { T1<:Real, T2<:Real, T3<:Real, A1<:AbstractArray{T1,2}, A2<:AbstractArray{T2,2}, A3<:AbstractArray{T3,2}, IT<:AnyEXPRES, AS<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT}  }
   num_lines = num_chunks(clt)
   num_obs = length(clt)
   min_telluric_model_one_obs = ones(num_lines, num_obs )
   min_telluric_model_all_obs  = ones(num_lines)
   for ch_idx in 1:num_chunks(clt)
       for t_idx in 1:num_obs
          view_indices = clt.chunk_list[t_idx].data[ch_idx].λ.indices
          cols = view_indices[1]
          order = view_indices[2]
          min_telluric_model_one_obs[ch_idx, t_idx] = minimum(view(data[t_idx].metadata[:tellurics], cols, order))
      end # times
      min_telluric_model_all_obs[ch_idx] = minimum(min_telluric_model_one_obs[ch_idx, :])
  end # lines
  return min_telluric_model_all_obs
end

"""
   `make_clean_line_list_from_tellurics(line_list, expres_data; Δv_to_avoid_tellurics )`
Returns a new line list that excludes lines with telluric contamination.
Inputs:
- `line_list`:  Dataframe containing field lambda
- `expres_data`:  Array of spectra
- `Δv_to_avoid_tellurics`:  in m/s
Outputs:
- `line_list_without_tellurics`:   DataFrame with fields: `lambda`, `weight`, `lambda_lo`, and `lambda_hi`.
Warning: Currently, assumes a tellurics value in metadata for each spectra, such as is provided by EXPRES.
"""
function make_clean_line_list_from_tellurics(line_list::DataFrame, expres_data::DT; Δv_to_avoid_tellurics::Real = default_Δv_to_avoid_tellurics,
            v_center_to_avoid_tellurics::Real = 0.0, telluric_threshold::Real = 1
               ) where { T1<:Real, A1<:AbstractArray{T1}, T2<:Real, A2<:AbstractArray{T2}, T3<:Real, A3<:AbstractArray{T3}, IT<:EXPRES.AnyEXPRES, ST<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT}, DT<:AbstractArray{ST,1} }

   @assert 0.5*RvSpectMLBase.max_bc_earth_rotation <= Δv_to_avoid_tellurics <= 4*RvSpectMLBase.max_bc
   line_list_to_search_for_tellurics = copy(line_list)
   # TODO: Check sign convention for v_center_to_avoid_tellurics
   line_list_to_search_for_tellurics.lambda_lo = line_list_to_search_for_tellurics.lambda./calc_doppler_factor(Δv_to_avoid_tellurics).*calc_doppler_factor(v_center_to_avoid_tellurics)
   line_list_to_search_for_tellurics.lambda_hi = line_list_to_search_for_tellurics.lambda.*calc_doppler_factor(Δv_to_avoid_tellurics).*calc_doppler_factor(v_center_to_avoid_tellurics)
   chunk_list_timeseries = RvSpectMLBase.make_chunk_list_timeseries_around_lines(expres_data,line_list_to_search_for_tellurics)
   line_list_to_search_for_tellurics.min_telluric_model_all_obs = find_worst_telluric_in_each_chunk( chunk_list_timeseries, expres_data)
   line_list_no_tellurics_df = line_list_to_search_for_tellurics |> @filter(_.min_telluric_model_all_obs >= telluric_threshold ) |> DataFrame
end

default_min_Δv_clean = 8000.0

function find_ranges_with_tellurics_in_order(spectrum::ST, order::Integer; telluric_threshold::Real = 1, min_Δv_clean::Real = default_min_Δv_clean, verbose::Bool = false) where { T1<:Real, A1<:AbstractArray{T1}, T2<:Real, A2<:AbstractArray{T2}, T3<:Real, A3<:AbstractArray{T3}, IT<:EXPRES.AnyEXPRES, ST<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT}  }

    @assert haskey(spectrum.metadata,:tellurics)
    telluric_ranges = Vector{Tuple{Int64,Int64}}(undef,0)

    c = RvSpectMLBase.speed_of_light_mps
    all_wavelengths = range(first(spectrum.λ[:,order]),stop=last(spectrum.λ[:,order]),length=size(spectrum.λ,1) )
    tellurics = spectrum.metadata[:tellurics]
    start = findfirst(x->!isnan(x),view(tellurics,:,order) )
    stop = findlast(x->!isnan(x),view(tellurics,:,order) )
    if isnothing(start) || isnothing(stop)
        if verbose   println("# The entire order ", order, " was NaNs!?!")   end
        return Vector{Tuple{eltype(spectrum.λ), eltype(spectrum.λ)} }(undef,0)
    end
    @assert stop >= start +1
    in_telluric = tellurics[start,order] < telluric_threshold
    idx_start_this_telluric = in_telluric ? start : 0
    idx_stop_this_telluric = 0
    for i in start+1:stop
        if in_telluric
            if !(tellurics[i,order] < telluric_threshold)  # exitted telluric
                idx_stop_this_telluric = i-1
                #println("# Order = " , order, "  Adding pixels = ", idx_start_this_telluric, " - ", idx_stop_this_telluric, " λ = ", spectrum.λ[idx_start_this_telluric,order], " - ", spectrum.λ[idx_stop_this_telluric,order] )
                push!(telluric_ranges, (idx_start_this_telluric,idx_stop_this_telluric) )
                #push!(telluric_ranges,(spectrum.λ[idx_start_this_telluric,order], spectrum.λ[idx_stop_this_telluric,order]) )
                in_telluric = false
            end
        else
            if (tellurics[i,order] < telluric_threshold)  # entered telluric
                idx_start_this_telluric = i
                if length(telluric_ranges) >= 1
                    idx_last_start = last(telluric_ranges)[1]
                    idx_last_stop = last(telluric_ranges)[2]
                    if spectrum.λ[idx_start_this_telluric,order] - spectrum.λ[idx_last_stop,order] <= min_Δv_clean/c * spectrum.λ[idx_start_this_telluric,order]
                        idx_start_this_telluric = first(last(telluric_ranges))
                        pop!(telluric_ranges)
                    end
                end
                in_telluric = true
            end
        end
    end
    if in_telluric
        idx_stop_this_telluric = stop
        #println("# Order = " , order, "  Adding pixels = ", idx_start_this_telluric, " - ", idx_stop_this_telluric, " λ = ", spectrum.λ[idx_start_this_telluric,order], " - ", spectrum.λ[idx_stop_this_telluric,order] )
        push!(telluric_ranges,(idx_start_this_telluric,idx_stop_this_telluric) )
    end
    #lambda_range = map(r->(spectrum.λ[first(r),order], spectrum.λ[last(r),order] ), telluric_ranges)
    DataFrame(:lambda_lo=>map(r->spectrum.λ[first(r),order], telluric_ranges), :lambda_hi=>map(r->spectrum.λ[last(r),order], telluric_ranges) )
end

import EchelleInstruments.merge_sorted_wavelength_ranges

function find_ranges_with_tellurics(spectrum::ST; min_order::Integer = 1, max_order::Integer = size(spectrum.λ,2), min_Δv_clean::Real = default_min_Δv_clean, telluric_threshold::Real = 1 ) where {
            T1<:Real, A1<:AbstractArray{T1}, T2<:Real, A2<:AbstractArray{T2}, T3<:Real, A3<:AbstractArray{T3}, IT<:EXPRES.AnyEXPRES, ST<:Spectra2DBasic{T1,T2,T3,A1,A2,A3,IT}  }
    @assert haskey(spectrum.metadata,:tellurics)
    c = RvSpectMLBase.speed_of_light_mps

    #list_of_ranges = Vector{Tuple{eltype(spectrum.λ),eltype(spectrum.λ)}}(undef,0)
    list_of_ranges = DataFrame()
    for ord in min_order:max_order
        ranges_for_order = find_ranges_with_tellurics_in_order(spectrum,ord, min_Δv_clean=min_Δv_clean, telluric_threshold=telluric_threshold)
        if size(ranges_for_order,1) >= 1
            append!(list_of_ranges,ranges_for_order)
        end
    end

    if size(list_of_ranges,1)<=1   return list_of_ranges   end
    sort!(list_of_ranges, :lambda_lo )
    non_overlapping_ranges = merge_sorted_wavelength_ranges(list_of_ranges,min_Δv_clean=min_Δv_clean)
    return non_overlapping_ranges
end

#import EchelleInstruments.choose_obs_idx_for_init_guess

function choose_obs_idx_for_init_guess(df::DataFrame, inst::IT ) where { IT<:EXPRES.AnyEXPRES }
   @assert size(df,1) >= 1
   @assert hasproperty(df,:snr_prelim)
   idx = findmax(df.snr_prelim)[2]
   @assert 1 <= idx <= size(df,1)
   return idx
end

max_stars_repo_name: RvSpectML/EchelleInstruments

max_issues_repo_name: RvSpectML/EchelleInstruments

max_forks_repo_name: RvSpectML/EchelleInstruments


text:
using CirruSepal
using Base.Test

f = open("../examples/demo.cirru", "r")

ast = CirruSepal.compile(readall(f))

dump(ast)

eval(ast)


text:
# Second-order Soliton model from J. Lightwave Tech. Vol. 25 No. 12 Dec. 2007
# Johan Hult
T_FWHM = 0.1ps
T_T0 = T_FWHM/(2*log(1+sqrt(2)))
β₂ = -0.01ps^2/m
γnl = 0.01W/m
P_pulse = (abs(β₂/γnl)/(T_T0))^2
soliton_period = π*T_T0^2/(2*abs(β₂))
soliton_fundamental = ToyModel(;betacoeff=[0.0/m, 0.0ps/m, β₂],
                                nonlinearcoeff=γnl,
                                linearloss=0.0/m,
                                power_in=P_pulse,
                                pulsetime=T_FWHM,
                                length=soliton_period)

P_pulse_secondorder = (abs(β₂/γnl)/(T_T0)*2)^2
soliton_secondorder = ToyModel(;ω0=200.0THz, betacoeff=[0.0/m, 0.0ps/m, β₂],
                                nonlinearcoeff=γnl,
                                linearloss=0.0/m,
                                power_in=P_pulse_secondorder,
                                pulsetime=T_FWHM,
                                length=soliton_period)

sol_fundamental = solve(soliton_fundamental, DynamicNLSE();
                        tpoints=2^12, time_window=25.0ps)
sol_secondorder = solve(soliton_secondorder, DynamicNLSE();
                        tpoints=2^12, time_window=25.0ps)

sq_diff(x, y) = abs(sum(abs2.(x) - abs2.(y)))
@test sq_diff(sol_fundamental(0m), sol_fundamental(soliton_period)) < 15
@test sq_diff(sol_secondorder(soliton_period), sol_secondorder(0m)) < 15
@test sq_diff(sol_fundamental(0m), sol_secondorder(0m)) > 5
@test sq_diff(FT(sol_fundamental, 0m), FT(sol_fundamental, soliton_period)) < 15/ps^2
@test sq_diff(sol_fundamental(soliton_period/2), sol_fundamental(0m)) < 15


text:
struct material{R<:Real}
    e1::R
    e2::R
    g12::R
    nu12::R
    rho::R
    xt::R
    xc::R
    yt::R
    yc::R
    s::R
    t::R
end

getQ(mat::material, theta::Real=0.0) = getQ(mat.e1, mat.e2, mat.g12, mat.nu12,
    theta)

getQ(mat::Array{<:material,1}, lam::Laminate) = getQ.(mat, lam.theta)


text:
module LagrangianParticleTracking

export LagrangianParticles, update_particle_properties!

using Adapt
using KernelAbstractions
using StructArrays

using Oceananigans.Grids
using Oceananigans.Architectures: device
using Oceananigans.Fields: interpolate, datatuple, compute!, location

import Base: size, length, show

abstract type AbstractParticle end

struct Particle{T} <: AbstractParticle
    x :: T
    y :: T
    z :: T
end

struct LagrangianParticles{P, R, T, D, Π}
        properties :: P
       restitution :: R
    tracked_fields :: T
          dynamics :: D
        parameters :: Π
end

@inline no_dynamics(args...) = nothing

"""
    LagrangianParticles(; x, y, z, restitution=1.0, dynamics=no_dynamics, parameters=nothing)

Construct some `LagrangianParticles` that can be passed to a model. The particles will have initial locations
`x`, `y`, and `z`. The coefficient of restitution for particle-wall collisions is specified by `restitution`.

`dynamics` is a function of `(lagrangian_particles, model, Δt)` that is called prior to advecting particles.
`parameters` can be accessed inside the `dynamics` function.
"""
function LagrangianParticles(; x, y, z, restitution=1.0, dynamics=no_dynamics, parameters=nothing)
    size(x) == size(y) == size(z) ||
        throw(ArgumentError("x, y, z must all have the same size!"))

    (ndims(x) == 1 && ndims(y) == 1 && ndims(z) == 1) ||
        throw(ArgumentError("x, y, z must have dimension 1 but ndims=($(ndims(x)), $(ndims(y)), $(ndims(z)))"))

    particles = StructArray{Particle}((x, y, z))

    return LagrangianParticles(particles; restitution, dynamics, parameters)
end

"""
    LagrangianParticles(particles::StructArray; restitution=1.0, tracked_fields::NamedTuple=NamedTuple(), dynamics=no_dynamics)

Construct some `LagrangianParticles` that can be passed to a model. The `particles` should be a `StructArray`
and can contain custom fields. The coefficient of restitution for particle-wall collisions is specified by `restitution`.

A number of `tracked_fields` may be passed in as a `NamedTuple` of fields. Each particle will track the value of each
field. Each tracked field must have a corresponding particle property. So if `T` is a tracked field, then `T` must also
be a custom particle property.

`dynamics` is a function of `(lagrangian_particles, model, Δt)` that is called prior to advecting particles.
`parameters` can be accessed inside the `dynamics` function.
"""
function LagrangianParticles(particles::StructArray; restitution=1.0, tracked_fields::NamedTuple=NamedTuple(),
                             dynamics=no_dynamics, parameters=nothing)

    for (field_name, tracked_field) in pairs(tracked_fields)
        field_name in propertynames(particles) ||
            throw(ArgumentError("$field_name is a tracked field but $(eltype(particles)) has no $field_name field! " *
                                "You might have to define your own particle type."))
    end

    return LagrangianParticles(particles, restitution, tracked_fields, dynamics, parameters)
end

size(lagrangian_particles::LagrangianParticles) = size(lagrangian_particles.properties)
length(lagrangian_particles::LagrangianParticles) = length(lagrangian_particles.properties)

function Base.show(io::IO, lagrangian_particles::LagrangianParticles)
    particles = lagrangian_particles.properties
    properties = propertynames(particles)
    fields = lagrangian_particles.tracked_fields
    print(io, "$(length(particles)) Lagrangian particles with\n",
        "├── $(length(properties)) properties: $properties\n",
        "└── $(length(fields)) tracked fields: $(propertynames(fields))")
end

include("update_particle_properties.jl")

end # module


text:
# TODO more specific
const CashQuantity = Unitful.Quantity
const CashType = Unitful.FreeUnits

####
#### Brokerage
####

abstract type Brokerage end

####
#### Asset
####

abstract type Asset end

fullsymbol(c::Asset) = "$(string(nameof(typeof(c)))):$(_fullsymbol(c))"

struct Stock <: Asset
    symbol::String
    exchange::String
    currency::CashType
end

_fullsymbol(c::Stock) = "$(c.exchange):$(c.symbol)"
currency(c::Stock) = c.currency

####
#### Position
####

struct Position{C<:Asset}
    contract::C
    position::Float64
    market_price::CashQuantity
    market_value::CashQuantity
    average_cost::CashQuantity
end

function Position(contract, position, market_price, market_value, average_cost::Real)
    # TODO is avg_cost in units of contract or account currency?
    Position(contract, position, market_price, market_value, average_cost * currency(contract))
end

####
#### Account
####

abstract type Account end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.7.0-DEV
#     language: julia
#     name: julia-1.7
# ---

# %% [markdown]
# See https://twitter.com/taketo1024/status/1400990265364549634
#
# A straightforward translation in Swift to Julia of
# * https://gist.github.com/taketo1024/e1e1b92b0c5c9c441bd870476356204e

# %%
step = 10

function forAnySmall(p)
    for i in 0:step-1
        ε = (step - i) / step
        if !p(ε)
            return false
        end
    end
    return true
end

function existsSmall(p)
    for i in 0:step-1
        δ = (step - i) / step
        if p(δ)
            return true
        end
    end
    return false
end

function forAnyInRange((from, to), p)
    for _ in 0:step-1
        x = from + rand() * (to - from)
        if !p(x)
            return false
        end
    end
    return true
end

function isContinuous(f, at, a)
    forAnySmall(ε ->
        existsSmall(δ ->
            forAnyInRange((a - δ, a + δ), x ->
                abs(f(x) - f(a)) < ε
            )
        )
    )
end

println(isContinuous(x -> x^2, :at, 0))
# true: x^2 is continuous at x = 0

println(isContinuous(floor, :at, 1))
# false: floor func is not continuous at x = 1

println(isContinuous(floor, :at, 1/2))
# true: floor func is continuous at x = 1/2

# %%
macro for_any_positive(ε, P)
    :(forAnySmall($(esc(ε)) -> $(esc(P))))
end

macro there_exists_some_positive(δ, P)
    :(existsSmall($(esc(δ)) -> $(esc(P))))
end

macro for_any(x_in_ran, P)
    x = x_in_ran.args[2]
    ran = x_in_ran.args[3]
    :(forAnyInRange($(esc(ran)), $(esc(x)) -> $(esc(P))))
end

macro p(x)
    :(println($(esc(x))))
end

# %%
a = 0
@p @for_any_positive ε @there_exists_some_positive δ @for_any x ∈ (a - δ, a + δ) abs(x^2 - a^2) < ε

# %%
a = 1
@p @for_any_positive ε @there_exists_some_positive δ @for_any x ∈ (a - δ, a + δ) abs(floor(x) - floor(a)) < ε

# %%
a = 0.5
@p @for_any_positive ε @there_exists_some_positive δ @for_any x ∈ (a - δ, a + δ) abs(floor(x) - floor(a)) < ε

# %%
a = 0
@p @for_any_positive ε (
    @there_exists_some_positive δ (
        @for_any x ∈ (a - δ, a + δ) abs(x^2 - a^2) < ε))

# %%
a = 1
@p @for_any_positive ε (
    @there_exists_some_positive δ (
        @for_any x ∈ (a - δ, a + δ) abs(floor(x) - floor(a)) < ε))

# %%
a = 0.5
@p @for_any_positive ε (
    @there_exists_some_positive δ (
        @for_any x ∈ (a - δ, a + δ) abs(floor(x) - floor(a)) < ε))

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public


text:
"""
    setcallbacks!(learner, callbacks)

Set `learner`'s callbacks to `callbacks`, removing all current callbacks.
"""
function setcallbacks!(learner, callbacks)
    learner.callbacks = Callbacks(callbacks)
end

"""
    addcallback!(learner, callback)

Adds `callback` to `learner` and updates the dependency graph.
"""
function addcallback!(learner, callback::AbstractCallback)
    learner.callbacks = Callbacks(vcat(learner.callbacks.cbs, callback))
    init!(callback, learner)
    return
end

"""
    getcallback(learner, C)

Find callback of type `C` in `learner`'s callbacks and return it.
If there is none, return `nothing`.
"""
function getcallback(learner, C::Type{<:FluxTraining.Callback})
    cbidx = findfirst(isa.(learner.callbacks.cbs, C))
    return isnothing(cbidx) ? nothing : learner.callbacks.cbs[cbidx]
end

"""
    replacecallback!(learner, callback::C)

Replace existing callback of type `C` on learner with `callback`.
Return the replaced callback.

If `learner` doesn't have a callback of type `C`, add `callback` and
return `nothing`.
"""
function replacecallback!(learner, callback::C) where {C <: FluxTraining.Callback}
    cbidx = findfirst(isa.(learner.callbacks.cbs, C))
    if isnothing(cbidx)
        FluxTraining.addcallback!(learner, callback)
        return nothing
    else
        oldcb = learner.callbacks.cbs[cbidx]
        learner.callbacks.cbs[cbidx] = callback
        FluxTraining.setcallbacks!(learner, learner.callbacks.cbs)
        return oldcb
    end
end

"""
    removecallback!(learner, C)

Remove the first callback of type `C` from `learner` and return it.
If there is none, return `nothing`.
"""
function removecallback!(learner, C::Type{<:FluxTraining.Callback})
    cbidx = findfirst(isa.(learner.callbacks.cbs, C))
    if isnothing(cbidx)
        return nothing
    end
    cb = popat!(learner.callbacks.cbs, cbidx)
    learner.callbacks = Callbacks(learner.callbacks.cbs)
    return cb
end


text:
abstract AbstractShrinker
immutable NoShrinker <: AbstractShrinker end

"""
Type for encoding a method of generating examples.

All strategies should implement `shrink`, `arbitrary` and `returntype`
methods
"""
abstract AbstractStrategy

"""
Takes an example (`x`) and a strategy and returns an iterable of subexamples
"""
shrink(x, strategy::AbstractStrategy) = error("shrink not implemented for $(typeof(strategy))")
"""
Generates a random example using `startegy`
"""
arbitrary(strategy::AbstractStrategy) = error("arbitrary not implemented for $(typeof(strategy))")
"""
Return type of the `arbitrary` function
"""
returntype(strategy::Type{AbstractStrategy}) = error("returntype not implemented for type $(typeof(strategy))")

immutable CounterexampleError{T} <: Exception
    n_tests::Int
    n_shrinks::Int
    counterexample::T
end

function Base.showerror(io::IO, ex::CounterexampleError)
    print(io, """
    Found counterexample after $(ex.n_tests) tests, $(ex.n_shrinks) shrinks:

    counterexample: $(ex.counterexample)

    """)
end

"""
Search method to find a counterexample in the event of propositions
failure.

Finds smaller counterexamples by recursively shrinking bad counterexamples
until no smaller counterexample can be found.
"""
function find_minimal_counterexample(prop, args::Tuple, strategy::AbstractStrategy)
    cand_feas = false
    args_old = deepcopy(args)
    n_shrinks = 0
    while true
        shrink_strategies = shrink(args_old, strategy)
        for args_new in shrink_strategies
            if !prop(args_new...)
                cand_feas = true
                args_old = args_new
                break
            end
        end

        if !cand_feas
            return n_shrinks, args_old
        end
        cand_feas = false
        n_shrinks += 1
    end
end

"""
Generate `n` test samples using `strategy` on a proposition (`prop`).

If the proposition holds for all test samples then the test passes
"""
function forAll(prop, strategies::Tuple;
                n::Int=100)
    strategy = tuples(strategies...)
    for i=1:n
        args = arbitrary(strategy)
        if !prop(args...)
            n_shrinks, example = find_minimal_counterexample(prop, args, strategy)
            throw(CounterexampleError(i, n_shrinks, example))
        end
    end
    true
end

max_stars_repo_name: cpritcha/PropertyCheck

max_issues_repo_name: cpritcha/PropertyCheck

max_forks_repo_name: cpritcha/PropertyCheck


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -0.27492067550218746 <= q <= 2.084201684402027)

                     add_NL_constraint(m, :(1/(1 + exp(-(-0.0416569430659357 + 0.13675627334278007*1/(1 + exp(-(0.1752604319892206 + -0.8557289632032381*1/(1 + exp(-(-0.028885055816546323 + 0.15091441453709065*1/(1 + exp(-(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])))) + -0.8760277243880852*1/(1 + exp(-(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))))) + 0.21149394336428795*1/(1 + exp(-(-0.29077418744555983 + -0.11119024401630329*1/(1 + exp(-(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])))) + 0.4817613487326513*1/(1 + exp(-(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3])))))))))) + 0.703085236566185*1/(1 + exp(-(-0.4914467550189556 + 0.5078124686195609*1/(1 + exp(-(-0.028885055816546323 + 0.15091441453709065*1/(1 + exp(-(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])))) + -0.8760277243880852*1/(1 + exp(-(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))))) + 0.7523099505151869*1/(1 + exp(-(-0.29077418744555983 + -0.11119024401630329*1/(1 + exp(-(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])))) + 0.4817613487326513*1/(1 + exp(-(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))))))))))) + 1/(1 + exp(-(0.6991595157542867 + 0.5132010992296534*1/(1 + exp(-(0.1752604319892206 + -0.8557289632032381*1/(1 + exp(-(-0.028885055816546323 + 0.15091441453709065*1/(1 + exp(-(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])))) + -0.8760277243880852*1/(1 + exp(-(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))))) + 0.21149394336428795*1/(1 + exp(-(-0.29077418744555983 + -0.11119024401630329*1/(1 + exp(-(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])))) + 0.4817613487326513*1/(1 + exp(-(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3])))))))))) + 0.6682334830395109*1/(1 + exp(-(-0.4914467550189556 + 0.5078124686195609*1/(1 + exp(-(-0.028885055816546323 + 0.15091441453709065*1/(1 + exp(-(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])))) + -0.8760277243880852*1/(1 + exp(-(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))))) + 0.7523099505151869*1/(1 + exp(-(-0.29077418744555983 + -0.11119024401630329*1/(1 + exp(-(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3])))) + 0.4817613487326513*1/(1 + exp(-(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions