text:
using DashBootstrapComponents, DashHtmlComponents

modal = html_div([
    dbc_button("Open modal", id = "open-dismiss"),
    dbc_modal(
        [
            dbc_modalheader(dbc_modaltitle("Dismissing"), close_button = false),
            dbc_modalbody(
                "This modal has no close button and can't be dismissed by " *
                "pressing ESC. Try clicking on the backdrop or the below " *
                "close button.",
            ),
            dbc_modalfooter(dbc_button("Close", id = "close-dismiss")),
        ],
        id = "modal-dismiss",
        keyboard = false,
        backdrop = "static",
    ),
],)

callback!(
    app,
    Output("modal-dismiss", "is_open"),
    Input("open-dismiss", "n_clicks"),
    Input("close-dismiss", "n_clicks"),
    State("modal-dismiss", "is_open"),
) do n_open, n_close, is_open
    return n_open > 0 || n_close > 0 ? is_open == 0 : is_open
end;

max_stars_repo_name: glsdown/dash-bootstrap-components

max_issues_repo_name: glsdown/dash-bootstrap-components

max_forks_repo_name: glsdown/dash-bootstrap-components


text:
module EquationsOfState

include("prelude.jl")
include("Collections.jl")
include("NonlinearFitting.jl")
include("FiniteStrains.jl")
include("LinearFitting.jl")
include("Find.jl")

end # module

ext: jl


text:
module IDError
struct UnIdentifierError <: Exception
    msg::AbstractString
end
end

using .IDError: UnIdentifierError

function Base.showerror(io::IO, err::UnIdentifierError)
    print(io, err.msg)
end


text:
function clean(phone_number)
    # Main machinery is in the regex expression.
    NXX = "[2-9][\\d]{2}"
    country = "(?:\\+?1)?"
    area = "\\(?\\s*($NXX)\\s*\\)?"
    exchange = "($NXX)"
    subscriber = "([\\d]{4})"
    fillers = "\\s*(?:\\.|-)?\\s*"
    r = Regex("^\\s*$country$fillers$area$fillers$exchange$fillers$subscriber\\s*\$")

    result = match(r, phone_number)

    if result === nothing
        throw(ArgumentError("\"$phone_number\" is not a valid NANP phone number."))
    else
        return Base.string(result.captures...)
    end
end


text:
## general tape transformations

"""
Trace function execution, apply transformation on the resulting tape and
recompile function back. Returns function as is if transformation didn't result
in any changes (as determined by the second return value from `tform`). This way
`transform()` can be safely applied to any function and changes only the ones
for which `tform` is applicable.

Params:

 * tform::Function - transfomrmation function of type (Tape) -> (Tape, Bool)
   where second return value shows whether
"""
function transform(tform::Function, f, args...)
    _, tape = trace(f, args...)
    new_tape, changed = tform(tape)
    return changed ? compile(new_tape; bind=false) : f
end


text:
export get_positions,
    get_velocities,
    get_mass_densities,
    get_kernel_widths,
    get_specific_energies,
    get_pressures,
    step!

abstract type Stepper{N} end

"""
    get_positions(stepper::Stepper{N}) -> Vector{SVectorF{N}}

Return the positions of the particles.
"""
function get_positions end

"""
    get_velocities(stepper::Stepper{N}) -> Vector{SVectorF{N}}

Return the velocities of the particles.
"""
function get_velocities end

"""
    get_mass_densities(stepper::Stepper) -> Vector{Float}

Return the mass densities of the particles.
"""
function get_mass_densities end

"""
    get_kernel_widths(stepper::Stepper) -> Vector{Float}

Return the widths of the particle smoothing kernels.
"""
function get_kernel_widths end

"""
    get_specific_energies(stepper::Stepper) -> Vector{Float}

Return the specific energies of the particles.
"""
function get_specific_energies end

"""
    get_pressures(stepper::Stepper) -> Vector{Float}

Return the pressures of the particles.
"""
function get_pressures end

"""
    step!(stepper::Stepper, Δt::Number)

Advance the state of the fluid with the given time step.
"""
function step! end

include("stepping/euler.jl")

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
module AppDemo

Base.Experimental.@optlevel 1

using ReduxImGui
using ReduxImGui.Redux
using ReduxImGui.CImGui
using ReduxImGui.Buttons: AbstractButtonState, AbstractButtonAction
using ReduxImGui.RadioButtonGroups: AbstractRadioButtonGroupAction
using ReduxImGui.Checkboxes: AbstractCheckboxAction
using ReduxImGui.Menus: AbstractMenu, AbstractMenuAction
using ReduxImGui.MainMenuBars: AbstractMainMenuBar, AbstractMainMenuBarAction

# additional widgets
include("../widgets/Repeaters.jl")
using .Repeaters

include("../utils.jl")

# UIs
include("ui/buttons.jl")
include("ui/checkboxes.jl")
include("ui/combos.jl")
include("ui/repeaters.jl")
include("ui/inputs.jl")
include("ui/drags.jl")
include("ui/sliders.jl")
include("ui/colors.jl")
include("ui/listboxes.jl")

# actions
abstract type AbstractAppDemoAction <: AbstractSyncAction end

# state
struct State <: AbstractImmutableState
    buttons::Dict{String,Buttons.AbstractButtonState}
    checkboxes::Dict{String,Checkboxes.State}
    radio_button_group::RadioButtonGroups.State
    colorful_buttons::Vector{ColorButtons.State}
    repeaters::Dict{String,Repeaters.State}
    combos::Dict{String,Combos.State}
    input_texts::Dict{String,InputTexts.State}
    input_text_with_hints::Dict{String,InputTextWithHints.State}
    input_ints::Dict{String,InputInts.State}
    input_floats::Dict{String,InputFloats.State}
    input_doubles::Dict{String,InputDoubles.State}
    drag_ints::Dict{String,DragInts.State}
    drag_floats::Dict{String,DragFloats.State}
    slider_ints::Dict{String,SliderInts.State}
    slider_floats::Dict{String,SliderFloats.State}
    slider_angles::Dict{String,SliderAngles.State}
    slider_strings::Dict{String,SliderStrings.State}
    color_edits::Dict{String,ColorEdits.State}
    listboxes::Dict{String,ListBoxes.State}
    main_menubar::MainMenuBar
    extra::Dict{String,Bool}
end

# initial states
include("init/buttons.jl")
include("init/checkboxes.jl")
include("init/combos.jl")
include("init/repeaters.jl")
include("init/inputs.jl")
include("init/drags.jl")
include("init/sliders.jl")
include("init/colors.jl")
include("init/listboxes.jl")
include("init/menus.jl")

const IMGUI_DEMO_STATE = AppDemo.State(
    IMGUI_DEMO_BUTTON_STATES,
    IMGUI_DEMO_CHECKBOX_STATES,
    IMGUI_DEMO_RADIO_BUTTON_GROUP_STATES,
    IMGUI_DEMO_COLOR_BUTTON_STATES,
    IMGUI_DEMO_REPEATER_STATES,
    IMGUI_DEMO_COMBO_STATES,
    IMGUI_DEMO_INPUT_TEXT_STATES,
    IMGUI_DEMO_INPUT_TEXT_WITH_HINT_STATES,
    IMGUI_DEMO_INPUT_INT_STATES,
    IMGUI_DEMO_INPUT_FLOAT_STATES,
    IMGUI_DEMO_INPUT_DOUBLE_STATES,
    IMGUI_DEMO_DRAG_INT_STATES,
    IMGUI_DEMO_DRAG_FLOAT_STATES,
    IMGUI_DEMO_SLIDER_INT_STATES,
    IMGUI_DEMO_SLIDER_FLOAT_STATES,
    IMGUI_DEMO_SLIDER_ANGLE_STATES,
    IMGUI_DEMO_SLIDER_STRING_STATES,
    IMGUI_DEMO_COLOR_EDIT_STATES,
    IMGUI_DEMO_LISTBOX_STATES,
    IMGUI_DEMO_MAIN_MENUBAR_STATES,
    Dict(
        "show_app_main_manu_bar"=>true,
    )
)

# reducers
reducer(state::AbstractState, action::AbstractAction) = state
reducer(state::Vector{<:AbstractState}, action::AbstractAction) = state
reducer(state::Dict{String,<:AbstractState}, action::AbstractAction) = state

function reducer(state::Dict{String,Buttons.AbstractButtonState}, action::AbstractButtonAction)
    next_state = Buttons.reducer(state, action)
    next_state = OnOffButtons.reducer(next_state, action)
    next_state = ColorButtons.reducer(next_state, action)
    return next_state
end

reducer(state::Dict{String,Checkboxes.State}, action::AbstractCheckboxAction) =
    Checkboxes.reducer(state, action)

reducer(state::RadioButtonGroups.State, action::AbstractRadioButtonGroupAction) =
    RadioButtonGroups.reducer(state, action)

reducer(state::Vector{ColorButtons.State}, action::AbstractButtonAction) =
    ColorButtons.reducer(state, action)

reducer(state::Dict{String,Menu}, action::AbstractMenuAction) =
    Menus.reducer(state, action)

reducer(state::MainMenuBar, action::AbstractMainMenuBarAction) =
    MainMenuBars.reducer(state, action)

reducer(state::State, action::AbstractAction) =
    State(reducer(state.buttons, action),
          reducer(state.checkboxes, action),
          reducer(state.radio_button_group, action),
          reducer(state.colorful_buttons, action),
          Repeaters.repeater(state.repeaters, action),
          Combos.combo(state.combos, action),
          InputTexts.input_text(state.input_texts, action),
          InputTextWithHints.input_text_with_hint(state.input_text_with_hints, action),
          InputInts.input_int(state.input_ints, action),
          InputFloats.input_float(state.input_floats, action),
          InputDoubles.input_double(state.input_doubles, action),
          DragInts.drag_int(state.drag_ints, action),
          DragFloats.drag_float(state.drag_floats, action),
          SliderInts.slider_int(state.slider_ints, action),
          SliderFloats.slider_float(state.slider_floats, action),
          SliderAngles.slider_angle(state.slider_angles, action),
          SliderStrings.slider_string(state.slider_strings, action),
          ColorEdits.color_edit(state.color_edits, action),
          ListBoxes.listbox(state.listboxes, action),
          reducer(state.main_menubar, action),
          state.extra,
    )

# helper
"""
    ImGui_Demo(store::AbstractStore, get_state=Redux.get_state)
"""
function ImGui_Demo(store::AbstractStore, get_state=Redux.get_state)
    CImGui.ShowDemoWindow(Ref(true))
    CImGui.SetNextWindowPos((650, 20), CImGui.ImGuiCond_FirstUseEver)
    CImGui.SetNextWindowSize((550, 680), CImGui.ImGuiCond_FirstUseEver)
    CImGui.Begin("Demo", Ref(true), CImGui.ImGuiWindowFlags_NoSavedSettings)

    get_state(store).extra["show_app_main_manu_bar"] &&
        get_state(store).main_menubar(store, s->get_state(s).main_menubar)

    ReduxImGui.TreeNode("Basic") do
        basic_button(store, get_state)
        basic_checkbox(store, get_state)
        radio_button_group(store, get_state)
        color_buttons(store, get_state)

        # use AlignTextToFramePadding() to align text baseline to the baseline of framed elements (otherwise a Text+SameLine+Button sequence will have the text a little too high by default)
        CImGui.AlignTextToFramePadding()
        CImGui.Text("Hold to repeat:")
        CImGui.SameLine()
        arrow_buttons_with_repeater(store, get_state)

        CImGui.Text("Hover over me")
        CImGui.IsItemHovered() && CImGui.SetTooltip("I am a tooltip")

        CImGui.SameLine()
        CImGui.Text("- or me")
        CImGui.IsItemHovered() && ReduxImGui.Tooltip() do
            CImGui.Text("I am a fancy tooltip")
            CImGui.PlotLines("Curve", Cfloat[0.6, 0.1, 1.0, 0.5, 0.92, 0.1, 0.2], 7)
        end

        CImGui.Separator()

        CImGui.LabelText("label", "Value")
        naive_combo(store, get_state)
        naive_inputs(store, get_state)
        naive_drags(store, get_state)
        naive_sliders(store, get_state)
        naive_color_edits(store, get_state)
        naive_listboxes(store, get_state)

    end

    ReduxImGui.TreeNode("Trees") do
        ReduxImGui.TreeNode("Basic trees") do
            for i = 0:4
                if CImGui.TreeNode(Ptr{Cvoid}(i), "Child $i")
                    CImGui.Text("blah blah")
                    CImGui.SameLine()
                    CImGui.SmallButton("button") && @info "Trigger `Basic trees`'s button | find me here: $(@__FILE__) at line $(@__LINE__)"
                    CImGui.TreePop()
                end
            end
        end
    end

    ReduxImGui.TreeNode("Advanced, with Selectable nodes") do
        ShowHelpMarker("This is a more standard looking tree with selectable nodes.\nClick to select, CTRL+Click to toggle, click on arrows or double-click to open.")

    end
    CImGui.End()
end

end # module


text:
using Bisquit
using Genie, Genie.Router, Genie.Renderer.Json, Genie.Requests
using HTTP
using UUIDs

import JSON

# initialize an unique id for the new node instance
node_identifier = uuid4( ) |> string

# initialize a genesis block at startup
blockchain = create_genesis( )

# setup route adress on localhost
route("/mine") do
  # define some utiled variables
  last_block = blockchain.chain[ end ]
  last_proof = last_block.proof
  proof = proof_of_work( last_proof )

  # after finding the proof we'll get a reward
  # "0" to identify that the current node has mined a new block
  new_transaction( blockchain, "0", node_identifier, 1 )

  previous_hash = Bisquit.hash( last_block )
  block = Bisquit.new_block

  response = Dict(
    "message"=> "Mined new Block",
    "index"=> block.index,
    "transactions"=> block.transactions,
    "proof"=> block.proof,
    "previous_hash"=> block.previous_hash,
  )
  json( response )
end

# setup route adress on localhost
# view the current chain transaction
route("/chain") do
  response = Dict(
    "chain" => blockchain.chain,
    "length" => length(blockchain.chain)
  )

  # format the output
  response |> JSON.json
end

route("/transactions/new", method = POST) do
  # starting the json request
  data = jsonpayload( )

  # show the current sha256 key for following objects
  print( keys( data ) )

  # add array to the requested keys
  requested_keys = [ "sender", "recipient", "amount" ]

  # check if the current transaction contains the required keys. If not do: give a shoutout error message
  if !( length( keys( data ) ) == length( requested_keys ) && all( k -> in( k, requested_keys ), keys( data ) ) )
    # format the error message in json format
    "Missing valuess" |> json
  end

  # initialize a new transaction
  index = new_transaction( blockchain, data[ "sender" ], data[ "recipipient" ], data[ "amount" ])

  # shout out the success message and format in json
  response = Dict("message" => "Transaction is beeing added to Block $index")
  json( response )
end

# setup route on localhost
# view current chain and send transactions
route("/send") do
  response = HTTP.request("POST", "http://localhost:8000/transactions/new", [("Content-Type", "application/json")], """{"sender":"hello", "recipient":"thanks", "amount":1}""")

  # format body output
  response.body |> String |> JSON
end

# run the server
Genie.startup(async = false)

max_stars_repo_name: timo-cmd/Bisquit


text:
function paralax(gun::Gun,target::AbstractTarget)
    h = gun.hg*(target.position-gun.zerotage)/(target.position*gun.zerotage)
    v = gun.vg*(target.position-gun.zerotage)/(target.position*gun.zerotage)
    return h,v
end


text:
using WebIO
using JSExpr

function hue_app()
    scope = Scope()
    import!(scope, ["//cdnjs.cloudflare.com/ajax/libs/p5.js/0.5.7/p5.js"])

    # Note: we explicitly qualify @js, @var, and @new with the JSExpr module
    # because Blink.jl exports macros with the same name.
    sketch = JSExpr.@js function (p5)
        JSExpr.@var s = function(p)
            JSExpr.@var barWidth = 20
            JSExpr.@var lastBar = -1

            p.setup = function ()
               p.createCanvas(720, 400)
               p.colorMode(p.HSB, p.height, p.height, p.height)
               p.noStroke()
               p.background(0)
            end

            p.draw = function ()
                JSExpr.@var whichBar = p.mouseX / barWidth
                JSExpr.@var barX = whichBar * barWidth;
                if whichBar != lastBar
                   p.fill(p.mouseY, p.height, p.height)
                   p.rect(barX, 0, barWidth, p.height)
                   lastBar = whichBar
               end
            end
        end
        this.dom.querySelector("#p5container").innerText = "";
        JSExpr.@new p5(s, "p5container");
    end
    onimport(scope, sketch)

    scope.dom = dom"div#p5container"("Loading p5...")

    scope
end

# Display in whatever frontend is avalaible
function main()
    if @isdefined(IJulia) || @isdefined(Juno)
        return hue_app()
    elseif @isdefined(Blink)
        win = Window()
        body!(win, hue_app())
        win
    elseif @isdefined(Mux)
        @sync webio_serve(page("/", req -> hue_app()), rand(8000:9000))
    else
        error("do one of using Mux, using Blink before running the
               example, or run it from within IJulia or Juno")
    end
end

main()


text:
@testset "1029.two-city-scheduling.jl" begin
    @test two_city_sched_cost([[10,20],[30,200],[400,50],[30,20]]) == 110
    @test two_city_sched_cost([[259,770],[448,54],[926,667],[184,139],[840,118],[577,469]]) == 1859
    @test two_city_sched_cost([[515,563],[451,713],[537,709],[343,819],[855,779],[457,60],[650,359],[631,42]]) == 3086
end


text:
@testset "Well-known text" begin
    @testset "Writing WKT" begin
	    @test wkt(SVector(1.1, 2.2)) == "POINT (1.1 2.2)"
	    @test wkt(SVector(1.1, 2.2, 3.3)) == "POINT Z (1.1 2.2 3.3)"

	    @test wkt(Line(SVector(1.1, 1.2), SVector(2.1, 2.2))) == "LINESTRING (1.1 1.2, 2.1 2.2)"
	    @test wkt(Line(SVector(1.1, 1.2, 1.3), SVector(2.1, 2.2, 2.3))) == "LINESTRING Z (1.1 1.2 1.3, 2.1 2.2 2.3)"

        @test wkt(LineString(SVector{2,Float64}[])) == "LINESTRING EMPTY"
   	    @test wkt(LineString([SVector(1.1, 1.2), SVector(2.1, 2.2)])) == "LINESTRING (1.1 1.2, 2.1 2.2)"
   	    @test wkt(LineString(SVector{3,Float64}[])) == "LINESTRING Z EMPTY"
	    @test wkt(LineString([SVector(1.1, 1.2, 1.3), SVector(2.1, 2.2, 2.3)])) == "LINESTRING Z (1.1 1.2 1.3, 2.1 2.2 2.3)"

	    @test wkt(Polygon(SVector{2,Float64}[])) == "POLYGON EMPTY"
	    @test wkt(Polygon([SVector(1.0, 1.0), SVector(1.0, 2.0), SVector(2.0, 2.0), SVector(1.0, 1.0)])) ==
	        "POLYGON ((1.0 1.0, 1.0 2.0, 2.0 2.0, 1.0 1.0))"
	    @test wkt(Polygon(SVector{3,Float64}[])) == "POLYGON Z EMPTY"
	    @test wkt(Polygon([SVector(1.0, 1.0, 0.0), SVector(1.0, 2.0, 0.0), SVector(2.0, 2.0, 0.0), SVector(1.0, 1.0, 0.0)])) ==
	        "POLYGON Z ((1.0 1.0 0.0, 1.0 2.0 0.0, 2.0 2.0 0.0, 1.0 1.0 0.0))"

        @test wkt(SVector{2,Float64}[]) == "MULTIPOINT EMPTY"
        @test wkt([SVector(1.1, 2.2)]) == "MULTIPOINT ((1.1 2.2))"
        @test wkt([SVector(1.1, 2.2), SVector(1.1, 2.2)]) == "MULTIPOINT ((1.1 2.2), (1.1 2.2))"
	    @test wkt(SVector{3,Float64}[]) == "MULTIPOINT Z EMPTY"
	    @test wkt([SVector(1.1, 2.2, 3.3)]) == "MULTIPOINT Z ((1.1 2.2 3.3))"
	    @test wkt([SVector(1.1, 2.2, 3.3), SVector(1.1, 2.2, 3.3)]) == "MULTIPOINT Z ((1.1 2.2 3.3), (1.1 2.2 3.3))"

	    @test wkt([Line(SVector(1.1, 1.2), SVector(2.1, 2.2))]) == "MULTILINESTRING ((1.1 1.2, 2.1 2.2))"
	    @test wkt([Line(SVector(1.1, 1.2), SVector(2.1, 2.2)), Line(SVector(1.1, 1.2), SVector(2.1, 2.2))]) ==
	        "MULTILINESTRING ((1.1 1.2, 2.1 2.2), (1.1 1.2, 2.1 2.2))"
	    @test wkt([Line(SVector(1.1, 1.2, 1.3), SVector(2.1, 2.2, 2.3))]) == "MULTILINESTRING Z ((1.1 1.2 1.3, 2.1 2.2 2.3))"
	    @test wkt([Line(SVector(1.1, 1.2, 1.3), SVector(2.1, 2.2, 2.3)), Line(SVector(1.1, 1.2, 1.3), SVector(2.1, 2.2, 2.3))]) ==
	        "MULTILINESTRING Z ((1.1 1.2 1.3, 2.1 2.2 2.3), (1.1 1.2 1.3, 2.1 2.2 2.3))"

   	    @test wkt(LineString{2,Float64,Vector{SVector{2,Float64}}}[]) == "MULTILINESTRING EMPTY"
	    @test wkt([LineString([SVector(1.1, 1.2), SVector(2.1, 2.2)])]) == "MULTILINESTRING ((1.1 1.2, 2.1 2.2))"
   	    @test wkt([LineString([SVector(1.1, 1.2), SVector(2.1, 2.2)]), LineString([SVector(1.1, 1.2), SVector(2.1, 2.2)])]) ==
   	    	"MULTILINESTRING ((1.1 1.2, 2.1 2.2), (1.1 1.2, 2.1 2.2))"
	    @test wkt(LineString{3,Float64,Vector{SVector{3,Float64}}}[]) == "MULTILINESTRING Z EMPTY"
	    @test wkt([LineString([SVector(1.1, 1.2, 1.3), SVector(2.1, 2.2, 2.3)])]) == "MULTILINESTRING Z ((1.1 1.2 1.3, 2.1 2.2 2.3))"
	    @test wkt([LineString([SVector(1.1, 1.2, 1.3), SVector(2.1, 2.2, 2.3)]), LineString([SVector(1.1, 1.2, 1.3), SVector(2.1, 2.2, 2.3)])]) ==
	        "MULTILINESTRING Z ((1.1 1.2 1.3, 2.1 2.2 2.3), (1.1 1.2 1.3, 2.1 2.2 2.3))"

        @test wkt(Polygon{2,Float64,LineString{2,Float64,Vector{SVector{2,Float64}}},Vector{LineString{2,Float64,Vector{SVector{2,Float64}}}}}[]) ==
            "MULTIPOLYGON EMPTY"
	    @test wkt([Polygon([SVector(1.0, 1.0), SVector(1.0, 2.0), SVector(2.0, 2.0), SVector(1.0, 1.0)])]) ==
	        "MULTIPOLYGON (((1.0 1.0, 1.0 2.0, 2.0 2.0, 1.0 1.0)))"
	    @test wkt([Polygon([SVector(1.0, 1.0), SVector(1.0, 2.0), SVector(2.0, 2.0), SVector(1.0, 1.0)]), Polygon([SVector(1.0, 1.0), SVector(1.0, 2.0), SVector(2.0, 2.0), SVector(1.0, 1.0)])]) ==
	        "MULTIPOLYGON (((1.0 1.0, 1.0 2.0, 2.0 2.0, 1.0 1.0)), ((1.0 1.0, 1.0 2.0, 2.0 2.0, 1.0 1.0)))"
	    @test wkt(Polygon{3,Float64,LineString{3,Float64,Vector{SVector{3,Float64}}},Vector{LineString{3,Float64,Vector{SVector{3,Float64}}}}}[]) ==
            "MULTIPOLYGON Z EMPTY"
	    @test wkt([Polygon([SVector(1.0, 1.0, 0.0), SVector(1.0, 2.0, 0.0), SVector(2.0, 2.0, 0.0), SVector(1.0, 1.0, 0.0)])]) ==
	        "MULTIPOLYGON Z (((1.0 1.0 0.0, 1.0 2.0 0.0, 2.0 2.0 0.0, 1.0 1.0 0.0)))"
	    @test wkt([Polygon([SVector(1.0, 1.0, 0.0), SVector(1.0, 2.0, 0.0), SVector(2.0, 2.0, 0.0), SVector(1.0, 1.0, 0.0)]), Polygon([SVector(1.0, 1.0, 0.0), SVector(1.0, 2.0, 0.0), SVector(2.0, 2.0, 0.0), SVector(1.0, 1.0, 0.0)])]) ==
	        "MULTIPOLYGON Z (((1.0 1.0 0.0, 1.0 2.0 0.0, 2.0 2.0 0.0, 1.0 1.0 0.0)), ((1.0 1.0 0.0, 1.0 2.0 0.0, 2.0 2.0 0.0, 1.0 1.0 0.0)))"

        @test wkt(Any[]) == "GEOMETRYCOLLECTION EMPTY"
	    @test wkt(Any[SVector(1.0,2.0)]) == "GEOMETRYCOLLECTION (POINT (1.0 2.0))"
	    @test wkt(Any[SVector(1.0,2.0), Line(SVector(1.1, 1.2), SVector(2.1, 2.2))]) == "GEOMETRYCOLLECTION (POINT (1.0 2.0), LINESTRING (1.1 1.2, 2.1 2.2))"
    end

    @testset "Reading WKT" begin
        @test_broken read_wkt(IOBuffer("POINT EMPTY")) == SVector(NaN, NaN) # Not sure what this should be?
	    @test read_wkt("POINT (1.1 2.2)") == SVector(1.1, 2.2)
	    @test_broken read_wkt(IOBuffer("POINT Z EMPTY")) == SVector(NaN, NaN, NaN) # Not sure what this should be?
	    @test read_wkt("POINT Z (1.1 2.2 3.3)") == SVector(1.1, 2.2, 3.3)

	    @test read_wkt("LINESTRING EMPTY") == LineString(SVector{2,Float64}[])
	    @test read_wkt("LINESTRING (1.1 1.2, 2.1 2.2)") == LineString([SVector(1.1, 1.2), SVector(2.1, 2.2)])
	    @test read_wkt("LINESTRING Z EMPTY") == LineString(SVector{3,Float64}[])
	    @test read_wkt("LINESTRING Z (1.1 1.2 1.3, 2.1 2.2 2.3)") == LineString([SVector(1.1, 1.2, 1.3), SVector(2.1, 2.2, 2.3)])

	    @test read_wkt("POLYGON EMPTY") == Polygon(SVector{2,Float64}[])
	    @test read_wkt("POLYGON ((1.0 1.0, 1.0 2.0, 2.0 2.0, 1.0 1.0))") ==
	        Polygon([SVector(1.0, 1.0), SVector(1.0, 2.0), SVector(2.0, 2.0), SVector(1.0, 1.0)])
	    @test read_wkt("POLYGON Z EMPTY") == Polygon(SVector{3,Float64}[])
	    @test read_wkt("POLYGON Z ((1.0 1.0 0.0, 1.0 2.0 0.0, 2.0 2.0 0.0, 1.0 1.0 0.0))") ==
	        Polygon([SVector(1.0, 1.0, 0.0), SVector(1.0, 2.0, 0.0), SVector(2.0, 2.0, 0.0), SVector(1.0, 1.0, 0.0)])

        @test read_wkt(IOBuffer("MULTIPOINT EMPTY"))::Vector{SVector{2,Float64}} == []
	    @test read_wkt(IOBuffer("MULTIPOINT ((1.1 2.2))")) == [SVector(1.1, 2.2)]
        @test read_wkt(IOBuffer("MULTIPOINT ((1.1 2.2), (1.1 2.2))")) == [SVector(1.1, 2.2), SVector(1.1, 2.2)]
        @test read_wkt(IOBuffer("MULTIPOINT Z EMPTY"))::Vector{SVector{3,Float64}} == []
	    @test read_wkt(IOBuffer("MULTIPOINT Z ((1.1 2.2 3.3))")) == [SVector(1.1, 2.2, 3.3)]
	    @test read_wkt(IOBuffer("MULTIPOINT Z ((1.1 2.2 3.3), (1.1 2.2 3.3))")) == [SVector(1.1, 2.2, 3.3), SVector(1.1, 2.2, 3.3)]

	    @test read_wkt("MULTILINESTRING EMPTY")::Vector{<:LineString{2}} == []
	    @test read_wkt("MULTILINESTRING ((1.1 1.2, 2.1 2.2))") == [LineString([SVector(1.1, 1.2), SVector(2.1, 2.2)])]
   	    @test read_wkt("MULTILINESTRING ((1.1 1.2, 2.1 2.2), (1.1 1.2, 2.1 2.2))") ==
   	        [LineString([SVector(1.1, 1.2), SVector(2.1, 2.2)]), LineString([SVector(1.1, 1.2), SVector(2.1, 2.2)])]
   	    @test read_wkt("MULTILINESTRING Z EMPTY")::Vector{<:LineString{3}} == []
	    @test read_wkt("MULTILINESTRING Z ((1.1 1.2 1.3, 2.1 2.2 2.3))") == [LineString([SVector(1.1, 1.2, 1.3), SVector(2.1, 2.2, 2.3)])]
	    @test read_wkt("MULTILINESTRING Z ((1.1 1.2 1.3, 2.1 2.2 2.3), (1.1 1.2 1.3, 2.1 2.2 2.3))") ==
	        [LineString([SVector(1.1, 1.2, 1.3), SVector(2.1, 2.2, 2.3)]), LineString([SVector(1.1, 1.2, 1.3), SVector(2.1, 2.2, 2.3)])]

	    @test read_wkt("MULTIPOLYGON EMPTY")::Vector{<:Polygon{2}} == []
        @test read_wkt("MULTIPOLYGON (((1.0 1.0, 1.0 2.0, 2.0 2.0, 1.0 1.0)))") ==
	        [Polygon([SVector(1.0, 1.0), SVector(1.0, 2.0), SVector(2.0, 2.0), SVector(1.0, 1.0)])]
	    @test read_wkt("MULTIPOLYGON (((1.0 1.0, 1.0 2.0, 2.0 2.0, 1.0 1.0)), ((1.0 1.0, 1.0 2.0, 2.0 2.0, 1.0 1.0)))") ==
	        [Polygon([SVector(1.0, 1.0), SVector(1.0, 2.0), SVector(2.0, 2.0), SVector(1.0, 1.0)]), Polygon([SVector(1.0, 1.0), SVector(1.0, 2.0), SVector(2.0, 2.0), SVector(1.0, 1.0)])]
	    @test read_wkt("MULTIPOLYGON Z EMPTY")::Vector{<:Polygon{3}} == []
	    @test read_wkt("MULTIPOLYGON Z (((1.0 1.0 0.0, 1.0 2.0 0.0, 2.0 2.0 0.0, 1.0 1.0 0.0)))") ==
	        [Polygon([SVector(1.0, 1.0, 0.0), SVector(1.0, 2.0, 0.0), SVector(2.0, 2.0, 0.0), SVector(1.0, 1.0, 0.0)])]
	    @test read_wkt("MULTIPOLYGON Z (((1.0 1.0 0.0, 1.0 2.0 0.0, 2.0 2.0 0.0, 1.0 1.0 0.0)), ((1.0 1.0 0.0, 1.0 2.0 0.0, 2.0 2.0 0.0, 1.0 1.0 0.0)))") ==
	        [Polygon([SVector(1.0, 1.0, 0.0), SVector(1.0, 2.0, 0.0), SVector(2.0, 2.0, 0.0), SVector(1.0, 1.0, 0.0)]), Polygon([SVector(1.0, 1.0, 0.0), SVector(1.0, 2.0, 0.0), SVector(2.0, 2.0, 0.0), SVector(1.0, 1.0, 0.0)])]

        @test read_wkt("GEOMETRYCOLLECTION EMPTY")::Vector{Any} == []
	    @test read_wkt("GEOMETRYCOLLECTION (POINT (1.0 2.0))") == Any[SVector(1.0,2.0)]
	    @test read_wkt("GEOMETRYCOLLECTION (POINT (1.0 2.0), LINESTRING (1.1 1.2, 2.1 2.2))") ==
	        Any[SVector(1.0,2.0), LineString([SVector(1.1, 1.2), SVector(2.1, 2.2)])]

	    @test_throws RoamesGeometry.WKTParsingError read_wkt(IOBuffer("POINT (1.1 2.2 3.3)"))
	    @test_throws RoamesGeometry.WKTParsingError read_wkt(IOBuffer("POINT Z (1.1 2.2)"))
	    @test_throws RoamesGeometry.WKTParsingError read_wkt(IOBuffer("POINT M (1.1 2.2 0.0)"))
	    @test_throws RoamesGeometry.WKTParsingError read_wkt(IOBuffer("POINT ZM (1.1 2.2 3.3 0.0)"))
	    @test_throws RoamesGeometry.WKTParsingError read_wkt(IOBuffer("PONT (1.1 2.2)"))
	    @test_throws RoamesGeometry.WKTParsingError read_wkt(IOBuffer("POINT 1.1 2.2)"))
	    @test_throws RoamesGeometry.WKTParsingError read_wkt(IOBuffer("POINT (1.1 2.2"))
	    @test_throws RoamesGeometry.WKTParsingError read_wkt(IOBuffer("POINT (1.1 2.2))"))
	    @test_throws RoamesGeometry.WKTParsingError read_wkt(IOBuffer("POINT (1.1 2.2)a"))
    end

    @testset "eWKT" begin
        @test wkt(SVector(1.1, 2.2); srid="32756") == "SRID=32756; POINT (1.1 2.2)"
        @test wkt(SVector(1.1, 2.2); srid="32756,ausgeoid09") == "SRID=32756,ausgeoid09; POINT (1.1 2.2)"

        srid_ref = Ref("")
        @test read_wkt("SRID=32756; POINT (1.1 2.2)", srid=srid_ref) == SVector(1.1, 2.2)
        @test srid_ref[] == "32756"

        srid_ref2 = Ref("")
        @test read_wkt("SRID=32756,ausgeoid09; POINT (1.1 2.2)",srid=srid_ref2) == SVector(1.1, 2.2)
        @test srid_ref2[] == "32756,ausgeoid09"
    end
end


text:
module FEAutodiffTests

using LinearAlgebra
using Gridap.FESpaces
using Gridap.Arrays
using Gridap.Fields
using Gridap.Geometry
using Gridap.TensorValues
using Gridap.CellData

domain = (0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)

V = FESpace(model=model,valuetype=Float64,reffe=:Lagrangian,order=2,conformity=:H1)
U = TrialFESpace(V)

dv = get_cell_basis(V)
du = get_cell_basis(U)
uh = FEFunction(U,rand(num_free_dofs(U)))

trian = Triangulation(model)
quad = CellQuadrature(trian,2)

function user_uh_to_cell_energy(uh)
  cell_e = integrate(0.5*∇(uh)⋅∇(uh),trian,quad)
end

function user_uh_to_cell_residual(uh)
  cell_r = integrate(∇(uh)⋅∇(dv),trian,quad)
end

function user_uh_to_cell_jacobian(uh)
  cell_j = integrate(∇(du)⋅∇(dv),trian,quad)
end

cell_r = user_uh_to_cell_residual(uh)
cell_j = user_uh_to_cell_jacobian(uh)
cell_h = cell_j

cell_r_auto = autodiff_cell_residual_from_energy(user_uh_to_cell_energy,uh)
cell_h_auto = autodiff_cell_jacobian_from_energy(user_uh_to_cell_energy,uh)
cell_j_auto = autodiff_cell_jacobian_from_residual(user_uh_to_cell_residual,uh)

test_array(cell_r_auto,cell_r)
test_array(cell_j_auto,cell_j)
test_array(cell_h_auto,cell_h)

trian_Γ = BoundaryTriangulation(model)
quad_Γ = CellQuadrature(trian_Γ,2)
cell_ids = get_cell_id(trian_Γ)

function user_uh_to_cell_energy_Γ(uh)
  uh_Γ = restrict(uh,trian_Γ)
  cell_e = integrate(0.5*∇(uh_Γ)⋅∇(uh_Γ),trian_Γ,quad_Γ)
end

function user_uh_to_cell_residual_Γ(uh)
  uh_Γ = restrict(uh,trian_Γ)
  dv_Γ = restrict(dv,trian_Γ)
  cell_r = integrate(∇(uh_Γ)⋅∇(dv_Γ),trian_Γ,quad_Γ)
end

function user_uh_to_cell_jacobian_Γ(uh)
  uh_Γ = restrict(uh,trian_Γ)
  dv_Γ = restrict(dv,trian_Γ)
  du_Γ = restrict(du,trian_Γ)
  cell_j = integrate(∇(du_Γ)⋅∇(dv_Γ),trian_Γ,quad_Γ)
end

cell_e_Γ = user_uh_to_cell_energy_Γ(uh)
cell_r_Γ = user_uh_to_cell_residual_Γ(uh)
cell_j_Γ = user_uh_to_cell_jacobian_Γ(uh)
cell_h_Γ = cell_j_Γ

cell_r_Γ_auto = autodiff_cell_residual_from_energy(user_uh_to_cell_energy_Γ,uh,cell_ids)
cell_h_Γ_auto = autodiff_cell_jacobian_from_energy(user_uh_to_cell_energy_Γ,uh,cell_ids)
cell_j_Γ_auto = autodiff_cell_jacobian_from_residual(user_uh_to_cell_residual_Γ,uh,cell_ids)

test_array(cell_r_Γ_auto,cell_r_Γ)
test_array(cell_j_Γ_auto,cell_j_Γ)
test_array(cell_h_Γ_auto,cell_h_Γ)

const p = 3
@law j(∇u) = norm(∇u)^(p-2) * ∇u
@law dj(∇du,∇u) = (p-2)*norm(∇u)^(p-4)*inner(∇u,∇du)*∇u + norm(∇u)^(p-2)*∇du
f(x) = 0

res(u,v) = ∇(v)⋅j(∇(u)) - v*f
jac(u,du,v) = ∇(v)⋅dj(∇(du),∇(u))

function user_uh_to_cell_residual_2(uh)
  cell_r = integrate(res(uh,dv),trian,quad)
end

function user_uh_to_cell_jacobian_2(uh)
  cell_j = integrate(jac(uh,du,dv),trian,quad)
end

cell_j = user_uh_to_cell_jacobian_2(uh)

cell_j_auto = autodiff_cell_jacobian_from_residual(user_uh_to_cell_residual_2,uh)

test_array(cell_j_auto,cell_j)

end # module


text:
using PkgBenchmark, CSV, WeakRefStrings
if !is_windows()
    using DecFP
end

prep(io::IO, ::Type{Int64}) = write(io, "10")
prep(io::IO, ::Type{Float64}) = write(io, "10.0")
prep(io::IO, ::Type{WeakRefString{UInt8}}) = write(io, "hey there sailor")
prep(io::IO, ::Type{String}) = write(io, "hey there sailor")
prep(io::IO, ::Type{Date}) = write(io, "2016-09-28")
prep(io::IO, ::Type{DateTime}) = write(io, "2016-09-28T03:21:00")
if !is_windows()
    prep(io::IO, ::Type{Dec64}) = write(io, "10.0")
end

function prep{T}(::Type{IOBuffer}, ::Type{T})
    io = IOBuffer()
    prep(io, T)
    seekstart(io)
    return io, ()->return
end
function prep{T}(::Type{IOStream}, ::Type{T})
    t = tempname()
    io = open(t, "w")
    prep(io, T)
    close(io)
    io = open(t, "r")
    return io, ()->rm(t)
end

TYPES = !is_windows() ? (Int, Float64, WeakRefString{UInt8}, String, Date, DateTime, Dec64) : (Int, Float64, WeakRefString{UInt8}, String, Date, DateTime)

@benchgroup "CSV" begin
    @benchgroup "CSV.parsefield" begin
        opts = CSV.Options()
        row = col = 1
        state = Ref{CSV.ParsingState}(CSV.None)
        for I in (IOBuffer, IOStream)
            for T in TYPES
                io, f = prep(I, T)
                @bench "$I - $T" CSV.parsefield($io, $T, opts, row, col, state)
                f()
            end
        end
    end
    FILE = joinpath(dirname(@__FILE__), "randoms_small.csv")

    @benchgroup "CSV.read" begin
        @bench "CSV.read" CSV.read(FILE)
    end

    @benchgroup "CSV.write" begin
        df = CSV.read(FILE)
        t = tempname()
        @bench "CSV.write" CSV.write(t, df)
    end

end

struct Select{T, NT, names}
    x::T
end

function Select(x, names...)
    NT = Tables.schema(x)
    return Select{typeof(x), NT, names}(x)
end

Base.length(s::Select) = length(s.x)
Base.eltype(s::Select{T, NT, names}) where {T, NT, names} = select(NT, names)
function select(::Type{NamedTuple{names, types}}, nms) where {names, types}
    typs = []
    for nm in nms
        for (i, n) in enumerate(names)
            nm === n && push!(typs, types.parameters[i])
        end
    end
    return NamedTuple{nms, Tuple{typs...}}
end

struct SelectRow{T}
    x::T
end

@inline function Base.iterate(s::Select, st=())
    state = iterate(s.x, st...)
    state === nothing && return nothing
    row, st = state
    return SelectRow(row), (st,)
end

@inline Base.getproperty(row::SelectRow, name::Symbol) = getproperty(getfield(row, 1), name)

function profile(N)
    dir = "/Users/jacobquinn/.julia/dev/CSV/test/test_files"
    file = joinpath(dir, "pandas_zeros.csv")
    f = CSV.File(file)
    Profile.clear()
    for i = 1:N
        seek(f.io, f.rowpositions[1])
        @profile f |> columntable
    end
    return
end

# generate single column files w/ 1M rows for each type
using WeakRefStrings

val = "hey"
for i in (1001, 100.1, WeakRefString{UInt8}(pointer(val), 3, 0), Date(2008, 1, 3), DateTime(2008, 3, 4))
    open("/Users/jacobquinn/Downloads/randoms_$(typeof(i)).csv", "w") do f
        for j = 1:1_000_000
            write(f, string(i))
            write(f, "\n")
        end
    end
end

using CSV, TextParse
for T in (Int, Float64, WeakRefStrings.WeakRefString{UInt8}, Date, DateTime)
    println("comparing for T = $T...")
    # T == WeakRefStrings.WeakRefString{UInt8} && continue
    @time CSV.read("/Users/jacobquinn/Downloads/randoms_$(T).csv");
    # @time TextParse.csvread("/Users/jacobquinn/Downloads/randoms_$T.csv");
end

for T in (Int, Float64, WeakRefStrings.WeakRefString{UInt8}, Date, DateTime)
    println("comparing for T = $T...")
    # T == WeakRefStrings.WeakRefString{UInt8} && continue
    # @time CSV.read("/Users/jacobquinn/Downloads/randoms_$(T).csv");
    @time TextParse.csvread("/Users/jacobquinn/Downloads/randoms_$T.csv");
end

@time CSV.read("/Users/jacobquinn/Downloads/yellow_tripdata_2015-01.csv");
@time TextParse.csvread("/Users/jacobquinn/Downloads/yellow_tripdata_2015-01.csv");

for T in ('Int64', 'Float64', 'WeakRefString{UInt8}', 'Date', 'DateTime'):
    start = time.time()
    delim = ','
    table = pandas.read_csv("/Users/jacobquinn/Downloads/randoms_" + T + ".csv", delimiter=delim)
    end = time.time()
    print(end - start)

start = time.time()
delim = ','
table = pandas.read_csv("/Users/jacobquinn/Downloads/yellow_tripdata_2015-01.csv", delimiter=delim)
end = time.time()
print(end - start)
@time df = CSV.read("/Users/jacobquinn/Downloads/file.txt"; delim=' ');
@time TextParse.csvread("/Users/jacobquinn/Downloads/randoms_$(T).csv")
# julia> for T in (Int, Float64, WeakRefStrings.WeakRefString{UInt8}, Date, DateTime)
#            println("comparing for T = $T...")
#            @time CSV.read("/Users/jacobquinn/Downloads/randoms_$(T).csv");
#            @time TextParse.csvread("/Users/jacobquinn/Downloads/randoms_$(T).csv");
#        end
# comparing for T = Int64...
# pre-allocating DataFrame w/ rows = 999999
#   0.043684 seconds (1.00 M allocations: 22.929 MiB, 31.61% gc time)
#   0.045556 seconds (460 allocations: 15.575 MiB, 3.20% gc time)
# comparing for T = Float64...
# pre-allocating DataFrame w/ rows = 999999
#   0.080026 seconds (1.00 M allocations: 22.974 MiB, 23.80% gc time)
#   0.082530 seconds (457 allocations: 16.528 MiB)
# comparing for T = WeakRefString{UInt8}...
# pre-allocating DataFrame w/ rows = 999999
#   0.058446 seconds (1.89 k allocations: 22.986 MiB, 8.53% gc time)
#   0.069034 seconds (595 allocations: 5.188 MiB)
# comparing for T = Date...
# pre-allocating DataFrame w/ rows = 999999
#   0.125229 seconds (2.00 M allocations: 53.504 MiB, 20.94% gc time)
#   0.120472 seconds (1.00 M allocations: 51.846 MiB, 6.73% gc time)
# comparing for T = DateTime...
# pre-allocating DataFrame w/ rows = 999999
#   0.175855 seconds (2.00 M allocations: 53.504 MiB, 23.30% gc time)
#   0.187619 seconds (1.00 M allocations: 60.516 MiB, 4.40% gc time)

T = Int64
@time source = CSV.Source("/Users/jacobquinn/Downloads/randoms_$(T).csv";)
@time source = CSV.Source("/Users/jacobquinn/Downloads/randoms_small.csv"; allowmissing=:auto)
@time source = CSV.Source("/Users/jacobquinn/Downloads/randoms_small.csv"; allowmissing=:none)
# source.schema = DataStreams.Data.Schema(DataStreams.Data.header(source.schema), (Int, String, String, Float64, Float64, Date, DateTime), 9)
# @time df = CSV.read(source, NamedTuple);
sink = Si = NamedTuple
transforms = Dict{Int,Function}(1=>x->x-1)
append = false
args = kwargs = ()
source_schema = DataStreams.Data.schema(source)
sink_schema, transforms2 = DataStreams.Data.transform(source_schema, transforms, true);
sinkstreamtype = DataStreams.Data.Field
sink = Si(sink_schema, sinkstreamtype, append, args...; kwargs...);
columns = []
filter = x->true
@code_warntype DataStreams.Data.stream!(source, sinkstreamtype, sink, source_schema, sink_schema, transforms2, filter, columns)
@time DataStreams.Data.stream!(source, sinkstreamtype, sink, source_schema, sink_schema, transforms2, filter, columns)

function testt(t)
    a = getfield(t, 1)
    b = getfield(t, 2)
    c = getfield(t, 3)
    d = getfield(t, 4)
    e = getfield(t, 5)
    f = getfield(t, 6)
    g = getfield(t, 7)
    return (a, b, c, d, e, f, g)
end
@code_warntype testt((i1=(?Int)[], i2=(?String)[], i3=(?String)[], i4=(?Float64)[], i5=(?Float64)[], i6=(?Date)[], i7=(?DateTime)[]))

@code_llvm DataStreams.Data.stream!(source, sinkstreamtype, sink, source_schema, sink_schema, transforms2, filter, columns)
@time DataStreams.Data.stream!(source, sinkstreamtype, sink, source_schema, sink_schema, transforms2, filter, columns)

@code_warntype @time CSV.parsefield(IOBuffer(), ?Int, CSV.Options(), 0, 0, CSV.STATE)

t = Vector{Int}(1000000)

# having CSV.parsefield(io, T) where T !>: Missing decreases allocations by 1.00M
# inlining CSV.parsefield also dropped allocations
# making CSV.Options not have a type parameter also sped things up
#

using BenchmarkTools

g(x) = x < 5 ? x : -1
A = [i for i = 1:10]
function get_then_set(A)
    @simd for i = 1:10
        @inbounds A[i] = g(i)
    end
    return A
end
@code_warntype g(1)
@code_warntype get_then_set(A)
@benchmark get_then_set(A) # 20ns

@inline g3(x) = g2(x)
@inline function g2(x)
    if x < 20
        return x * 20
    end

    if x < 15
        return nothing
    end

    if x < 12
        return 2x
    end

    if x * 20 / 4 % 2 == 0
        return 1
    end

    if x < 0
        return nothing
    end
    return nothing
end

A = Union{Int, Void}[i for i = 1:10]
@inline function get_then_set2(A)
    @simd for i = 1:10
        # Base.arrayset(A, g2(i), i)
        val = g3(i)
        if val isa Void
            @inbounds A[i] = val#::Union{Int, Void}
        else
            @inbounds A[i] = val#::Union{Int, Void}
        end
    end
    return A
end
function run_lots(N)
    A = Union{Int, Void}[i for i = 1:10]
    for i = 1:N
        get_then_set2(A)
    end
    return
end

@code_warntype g2(1)
@code_warntype get_then_set2(A)
@code_llvm get_then_set2(A)
@benchmark get_then_set2(A) # 155ns

g4(x::Int) = 1
g4(x::Void) = 0

A = [i for i = 1:10]
function get_sum(A)
    s = 0
    for a in A
        s += g4(a)
    end
    return s
end
@code_warntype get_sum(A)
@code_llvm get_sum(A)
@benchmark get_sum(A) # 24ns

A = Union{Int, Void}[i for i = 1:10]
A[[3, 5, 7]] = nothing
function get_sum2(A)
    s = 0
    for a in A
        s += g4(a)
    end
    return s
end
@code_warntype get_sum2(A)
@code_llvm get_sum(A)
@benchmark get_sum2(A) # 100ns

function getstatic{T}(t::T)
    return t[1]
end

function Base.getproperty(csvrow::Row{F}, ::Type{T}, col::Int, name::Symbol) where {T, F <: File{NT, transpose}} where {NT, transpose}
    col === 0 && return missing
    f = getfield(csvrow, 1)
    row = getfield(csvrow, 2)
    if transpose
        @inbounds Parsers.fastseek!(f.io, f.positions[col])
    else
        lastparsed = f.lastparsedcol[]
        # print("$name: lastparsed=$lastparsed, col=$col, ")
        if col === lastparsed + 1
            # print("reading the next sequential column; lastnewline=$(newline(f.lastparsedcode[])), ")
            if newline(f.lastparsedcode[])
                # println("result=missing")
                f.lastparsedcol[] = col
                return missing
            end
        elseif col > lastparsed + 1
            # print("lastnewline=$(newline(f.lastparsedcode[])), ")
            if newline(f.lastparsedcode[])
                f.lastparsedcol[] = col
                return missing
            end
            sk = skipcells(f, col - (lastparsed + 1))
            # print("skipping cells forward, skipped $(col - (lastparsed + 1)) cells=$sk, ")
            # skipping cells
            if !sk
                # println("result=missing")
                f.lastparsedcol[] = col
                return missing
            end
        else
            @inbounds Parsers.fastseek!(f.io, f.positions[row])
            sk = skipcells(f, col - 1)
            # print("reverse skipping cells, skipped $(col - 1) cells=$sk, ")
            # randomly seeking within row
            if !sk
                # println("result=missing")
                f.lastparsedcol[] = col
                return missing
            end
        end
    end
    # @show position(f.io)
    r = parsefield(f, parsingtype(T), row, col, f.strict)
    # println("result=$r")
    if transpose
        @inbounds f.positions[col] = position(f.io)
    else
        f.lastparsedcol[] = col
    end
    return r
end

# code for generating various-sized csv files
# used to determine the row vs. column-access threshold
function gencsv(rows, cols)
    df = DataFrame([round.(rand(rows), digits=4) for _ ∈ 1:cols], Symbol.(["col$i" for i ∈ 1:cols]))
    CSV.write("random_$(rows)_$(cols).csv", df)
end

function go(compile=true)
    # for cols in (10, 25, 50, 75, 100, 250)
    #     for rows in (10, 50, 100, 1000, 5000, 10000, 50000, 100000)
    #         println("generating $rows by $cols csv file...")
    #         @time gencsv(rows, cols)
    #     end
    # end
    rt = NamedTuple{(:compiled, :rows, :cols, :time), Tuple{Int, Int, Float64}}[]
    for cols in (10, 25, 50, 75, 100, 250)
        for rows in (10, 50, 100, 1000, 5000, 10000, 50000, 100000)
            println("reading $rows by $cols csv file...")
            e = @elapsed begin
                df = CSV.File("random_$(rows)_$(cols).csv") |> DataFrame;
            end
            push!(rt, (compiled=compile, rows=rows, cols=cols, time=e))
        end
    end
    CSV.write("results_$compile.csv", rt)
end


text:
using ForwardDiff

# Define Test Problems
# ====================

tests = []

# [1] some random 2D problem
f1(x) = [x[1] + x[2] + x[1]^2, x[2] + x[1]*x[2]]
df1(x) = ForwardDiff.jacobian(f1, x)
init1() = [1.0,1.0]
randinit1() = rand(2)
push!(tests, (f1, df1, init1, randinit1, "toy2d"))

# [2] Chandrasekhar H-equation; cf http://www4.ncsu.edu/~ctk/newton/Chapter3/heq.m
let c = 0.9, n = 100
   gr = ((1:n) .- 0.5) / n
   A_heq = gr * ones(n)'
   A_heq = (0.5 * c / n) * A_heq ./ (A_heq+A_heq')
   f2(x) = x - ones(n) ./ (ones(n) - (A_heq * x))
   df2(x) = ForwardDiff.jacobian(f2, x)
   init2() = ones(n)
   randinit2() = ones(n) + 0.1 * (rand(n) .- 0.5)
   push!(tests, (f2, df2, init2, randinit2, "H-equation"))
end

# [3] MCP example by Josephy - stolen from NLsolve.jl
f3(x) = [3*x[1]^2+2*x[1]*x[2]+2*x[2]^2+x[3]+3*x[4]-6,
         2*x[1]^2+x[1]+x[2]^2+3*x[3]+2*x[4]-2,
         3*x[1]^2+x[1]*x[2]+2*x[2]^2+2*x[3]+3*x[4]-1,
         x[1]^2+3*x[2]^2+2*x[3]+3*x[4]-3 ]
df3(x) = ForwardDiff.jacobian(f3, x)
init3() = zeros(4)
randinit3() =  [1.22474487, 0., 0., 0.5 ] + (rand(4) .-0.5)
push!(tests, (f3, df3, init3, randinit3, "Josephy"))


text:
"""
Helper function for ReweighingWrapper and ReweighingSamplingWrapper.
grps is an array of values of protected attribute. y is an array of ground truth values.
Array of (Frequency) weights are returned from this function.
"""
function _calculateWeights(grps, y)
    grps_classes = levels(grps)
    n_classes = length(grps_classes)
    n = length(y)

    fav = y.==levels(y)[2] # Boolean
    unfav = y.==levels(y)[1]
    n₊ = sum(fav) # Total Number of favourable outcomes

    weights = zeros(Float64, n)
    for i in 1:n_classes
        class = grps.==grps_classes[i]
        weights[class .& fav] .= n₊*sum(class)/(n*sum(class .& fav))
        weights[class .& unfav] .= (n-n₊)*sum(class)/(n*sum(class .& unfav))
    end
    return weights
end

"""
    ReweighingWrapper

ReweighingWrapper is a preprocessing algorithm wrapper in which Weights for each group-label combination is calculated.
These calculated weights are then passed to the classifier model which further uses it to make training fair.
"""
struct ReweighingWrapper{M<:MLJBase.Model} <: DeterministicComposite
    grp::Symbol
    classifier::M
end

"""
    ReweighingWrapper(classifier=nothing, grp=:class)

Instantiates a ReweighingWrapper which wrapper the `classifier` with the Reweighing fairness algorithm.
The sensitive attribute can be specified by the parameter `grp`.
If `classifier` doesn't support weights while training, an error is thrown.
"""
function ReweighingWrapper(; classifier::MLJBase.Model=nothing, grp::Symbol=:class)
    model = ReweighingWrapper(grp, classifier)
    message = MLJBase.clean!(model)
    isempty(message) || @warn message
    return model
end

function MLJBase.clean!(model::ReweighingWrapper)
    warning = ""
    model.classifier!=nothing || (warning *= "No classifier specified in model\n")
    target_scitype(model) <: AbstractVector{<:Finite} || (warning *= "Only Binary Classifiers are supported\n")
    supports_weights(model.classifier) || (warning *= "Classifier provided does not support weights\n")
    return warning
end

function MLJBase.fit(model::ReweighingWrapper,
    verbosity::Int, X, y)

    grps = X[:, model.grp]

    weights = _calculateWeights(grps, y)

    Xs = source(X)
    ys = source(y)
    Ws = source(weights)

    classifier = model.classifier
    mach1 = machine(classifier, Xs, ys, Ws)
    ŷ = MMI.predict(mach1, Xs)

    mach2 = machine(Deterministic(), Xs, ys; predict=ŷ)
    return!(mach2, model, verbosity)
end

MMI.input_scitype(::Type{<:ReweighingWrapper{M}}) where M = input_scitype(M)
MMI.target_scitype(::Type{<:ReweighingWrapper{M}}) where M = AbstractVector{<:Finite{2}}

"""
    ReweighingSamplingWrapper

ReweighingSamplingWrapper is a preprocessing algorithm wrapper in which Weights for each group-label combination is calculated.
Using the calculated weights, rows are sampled uniformly. The weight is used to sample uniformly.
The number of datapoints used to train after sampling from the reweighed dataset can be controlled by `factor`.
"""
struct ReweighingSamplingWrapper{M<:MLJBase.Model} <: DeterministicComposite
    grp::Symbol
    classifier::M
    factor::Float64
    rng::Union{Int,AbstractRNG}
end

"""
    ReweighingSamplingWrapper(classifier=nothing, grp=:class, factor=1, rng=Random.GLOBAL_RNG)

Instantiates a ReweighingSamplingWrapper which wrapper the classifier with the Reweighing fairness algorithm together with sampling.
The sensitive attribute can be specified by the parameter `grp`.
`factor`*number_of_samples_in_original_data datapoints are sampled using calculated weights and then used to train after sampling from the reweighed dataset.
A negative or no value value for `factor` parameter instructs the algorithm to use the same number of datapoints as in original sample.
"""
function ReweighingSamplingWrapper(; classifier::MLJBase.Model=nothing, grp::Symbol=:class, factor::Float64=1.0, rng=nothing)
    if rng isa Integer
        rng = MersenneTwister(rng)
    end
    if rng == nothing
        rng = Random.GLOBAL_RNG
    end
    model = ReweighingSamplingWrapper(grp, classifier, factor, rng)
    message = MLJBase.clean!(model)
    isempty(message) || @warn message
    return model
end

function MLJBase.clean!(model::ReweighingSamplingWrapper)
    warning = ""
    model.classifier!=nothing || (warning *= "No classifier specified in model\n")
    target_scitype(model) <: AbstractVector{<:Finite} || (warning *= "Only Binary Classifiers are supported\n")
    return warning
end

function MLJBase.fit(model::ReweighingSamplingWrapper,
    verbosity::Int, X, y)
    grps = X[:, model.grp]
    n = model.factor<0 ? length(y) : Int(round(model.factor*length(y)))
    weights = _calculateWeights(grps, y)
    indices = sample(model.rng, 1:length(y), FrequencyWeights(weights), n; replace=true, ordered=false)
    Xnew = X[indices, :]
    ynew = y[indices]
    Xs = source(Xnew)
    ys = source(ynew)
    mach1 = machine(model.classifier, Xs, ys)
    ŷ = MMI.predict(mach1, Xs)

    mach2 = machine(Deterministic(), Xs, ys; predict=ŷ)
    return!(mach2, model, verbosity)
end

MMI.input_scitype(::Type{<:ReweighingSamplingWrapper{M}}) where M = input_scitype(M)
MMI.target_scitype(::Type{<:ReweighingSamplingWrapper{M}}) where M = AbstractVector{<:Finite{2}}


text:
using Test

@testset "Column" begin

    include(string(@__DIR__) * "/Column/" * "ColumnTest1.jl")
    include(string(@__DIR__) * "/Column/" * "ColumnTest2.jl")

end


text:
######################
### Mutable states ###
######################

mutable struct DAState{T<:AbstractFloat}
    m     :: Int
    ϵ     :: T
    μ     :: T
    x_bar :: T
    H_bar :: T
end

function DAState(ϵ::AbstractFloat)
    μ = computeμ(ϵ)
    return DAState(0, ϵ, μ, 0.0, 0.0)
end

function computeμ(ϵ::AbstractFloat)
    return log(10 * ϵ) # see NUTS paper sec 3.2.1
end

function reset!(dastate::DAState{T}) where {T<:AbstractFloat}
    dastate.μ = computeμ(dastate.ϵ)
    dastate.m = 0
    dastate.x_bar = zero(T)
    dastate.H_bar = zero(T)
end

mutable struct MSSState{T<:AbstractFloat}
    ϵ :: T
end

################
### Adaptors ###
################

abstract type StepSizeAdaptor <: AbstractAdaptor end

# finalize!(adaptor::T) where {T<:StepSizeAdaptor} = nothing

getϵ(ss::StepSizeAdaptor) = ss.state.ϵ

struct FixedStepSize{T<:AbstractFloat} <: StepSizeAdaptor
    ϵ :: T
end
Base.show(io::IO, a::FixedStepSize) = print(io, "FixedStepSize($(a.ϵ))")

getϵ(fss::FixedStepSize) = fss.ϵ

"""
An implementation of the Nesterov dual averaging algorithm to tune step size.

References

Hoffman, M. D., & Gelman, A. (2014). The No-U-Turn Sampler: adaptively setting path lengths in Hamiltonian Monte Carlo. Journal of Machine Learning Research, 15(1), 1593-1623.
Nesterov, Y. (2009). Primal-dual subgradient methods for convex problems. Mathematical programming, 120(1), 221-259.
"""
struct NesterovDualAveraging{T<:AbstractFloat} <: StepSizeAdaptor
  γ     :: T
  t_0   :: T
  κ     :: T
  δ     :: T
  state :: DAState{T}
end
Base.show(io::IO, a::NesterovDualAveraging) = print(io, "NesterovDualAveraging(γ=$(a.γ), t_0=$(a.t_0), κ=$(a.κ), δ=$(a.δ), state.ϵ=$(getϵ(a)))")

NesterovDualAveraging(γ::T, t_0::T, κ::T, δ::T, ϵ::T) where {T<:AbstractFloat} = NesterovDualAveraging(γ, t_0, κ, δ, DAState(ϵ))
NesterovDualAveraging(δ::T, ϵ::T) where {T<:AbstractFloat} = NesterovDualAveraging(0.05, 10.0, 0.75, δ, ϵ)

struct ManualSSAdaptor{T<:AbstractFloat} <:StepSizeAdaptor
    state :: MSSState{T}
end
Base.show(io::IO, a::ManualSSAdaptor) = print(io, "ManualSSAdaptor()")

ManualSSAdaptor(initϵ::T) where {T<:AbstractFloat} = ManualSSAdaptor{T}(MSSState(initϵ))

# Ref: https://github.com/stan-dev/stan/blob/develop/src/stan/mcmc/stepsize_adaptation.hpp
# Note: This function is not merged with `adapt!` to empahsize the fact that
#       step size adaptation is not dependent on `θ`.
function adapt_stepsize!(da::NesterovDualAveraging{T}, α::T) where {T <: AbstractFloat}
    DEBUG && @debug "Adapting step size..." α

    # Clip average MH acceptance probability
    α = α > 1 ? one(T) : α

    @unpack state, γ, t_0, κ, δ = da
    @unpack μ, m, x_bar, H_bar = state

    m = m + 1

    η_H = one(T) / (m + t_0)
    H_bar = (one(T) - η_H) * H_bar + η_H * (δ - α)

    x = μ - H_bar * sqrt(m) / γ     # x ≡ logϵ
    η_x = m^(-κ)
    x_bar = (one(T) - η_x) * x_bar + η_x * x

    ϵ = exp(x)
    DEBUG && @debug "Adapting step size..." "new ϵ = $ϵ" "old ϵ = $(da.state.ϵ)"

    # TODO: we might want to remove this when all other numerical issues are correctly handelled
    if isnan(ϵ) || isinf(ϵ)
        @warn "Incorrect ϵ = $ϵ; ϵ_previous = $(da.state.ϵ) is used instead."
        @unpack m, ϵ, x_bar, H_bar = state

    end

    @pack! state = m, ϵ, x_bar, H_bar
end

adapt!(da::NesterovDualAveraging, θ::AbstractVector{<:AbstractFloat}, α::AbstractFloat) = adapt_stepsize!(da, α)

reset!(da::NesterovDualAveraging) = reset!(da.state)

function finalize!(da::NesterovDualAveraging)
    da.state.ϵ = exp(da.state.x_bar)
end


text:
module locate

end

max_stars_repo_name: tmaone/xsh

max_issues_repo_name: turing-complete/xsh

max_forks_repo_name: turing-complete/xsh


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля/
#РЕЗУЛЬТАТ: Робот - в исходном положении, и клетки поля промакированы так:
#нижний ряд - полностью, следующий - весь, за исключением одной последней клетки на Востоке и т.д.
function mark_diag(r)
    x,y = num_sudost!(r)
    mark_diag!(r)
    return_sudost!(r,x,y)
end
function num_sudost!(r)
    x=0
    y=0
    while isborder(r,Sud)==0
        move!(r,Sud)
        x+=1
    end
    while isborder(r,Ost)==0
        move!(r,Ost)
        y+=1
    end
    return x,y
end
function mark_diag!(r)
    while isborder(r,Nord)==0
        putmarker!(r)
        move!(r,Nord)
        move!(r,West)
    end
    while ismarker(r)==0
        for side in [West,Ost]
            while isborder(r,side)==0 && ismarker(r)==0
                putmarker!(r)
                move!(r,side)
            end
            putmarker!(r)
            if isborder(r,Sud)==0
                move!(r,Sud)
            end
        end
    end
    while isborder(r,Ost)==0
        move!(r,Ost)
    end
end
function return_sudost!(r,x,y)
    for i in 1:x
        move!(r,Nord)
    end
    for i in 1:y
        move!(r,West)
    end
end
mark_diag(r)


text:
module Green
include("common.jl")
using LinearAlgebra, Combinatorics
export Model, GreenN
export density, Heisenberg, thermalavg
export Gn, Gnc
# include("hilbert.jl")
# using .Hilbert

"""
    Model(H::Operator, c⁺::Vector{Operator}, c⁻::Vector{Operator})

Construct a model.  All input operators are in the Fock space. But all members of Model struct is in the eigenspace.
Operator should be a matrix type.

# Arguments:
- β: inverse temperature
- H: Hamiltonian in the Fock space
- c⁺, c⁻: creation/anniliation operators for each orbital, all in the Fock space
"""
struct Model{N,No}
    isfermi::Bool
    β::Float
    dim::Int # dim of Hilbert space
    Norbital::Int # number of orbitals
    E::SVector{N,Float} # eigen energy
    Z::Float # partition sum
    H::Operator # diagnoalized Hamitlonian
    c⁺::SVector{No,Operator} # creation operator in the eigenspace
    c⁻::SVector{No,Operator} # anniliation operator in the eigenspace
    n::SVector{No,Operator}  # density operator in the eigenspace

    function Model(β, H, _c⁺, isfermi=true)
        dim = size(H)[1]
        @assert size(H) == size(_c⁺[1])
        @assert size(H) == (dim, dim)
        Norbital = length(_c⁺)

        E, U = eigen(Float64.(Matrix(H)))  # U'*H*U will diagonalize the Hamiltonian
        Z = sum(exp.(-β * E))
        E = sort(E)

        Hdiag = zeros(Float, (dim, dim))
        Hdiag[diagind(Hdiag)] = E
        c⁺ = [U' * o * U for o in _c⁺]
        c⁻ = [U' * o' * U for o in _c⁺]
        n = [c⁺[i] * c⁻[i] for i in 1:Norbital]

        return new{dim,Norbital}(isfermi, β, dim, Norbital, E, Z, Hdiag, c⁺, c⁻, n)
    end
end

function thermalavg(o::Operator, E, β, Z)
    if !(size(o) == (length(E), length(E)))
        throw(AssertionError("Dimension of Operator[$(o.m), $(o.n)] doesn't match with $(length(E))"))
    end
    avg = 0.0
    for di in 1:length(E)
        avg += o[di, di] * exp(-β * E[di])
    end
    return avg / Z
end
"""
    Heisenberg(o::Operator, E, τ)

   Transform operator o into Heisenberg picture
```math
    \\exp(Hτ)\\cdot o \\cdot \\exp(-Hτ)
```
"""
function Heisenberg(o::Operator, E, τ)
    if !(size(o) == (length(E), length(E)))
        throw(AssertionError("Dimension of Operator[$(o.m), $(o.n)] doesn't match with $(length(E))"))
    end
    if abs(τ) < 1.0e-10
        return o  # if τ≈0, no evoluation is needed
    end
    oh = copy(o)
    for ri in 1:length(E)
        oh[ri, :] *= exp(E[ri] * τ)
    end
    for ci in 1:length(E)
        oh[:, ci] *= exp(-E[ci] * τ)
    end
    return oh
end

function parity(p)
    """
    calculate the parity of a given permutation of the array [1, 2, 3, ...]
    """
    n = length(p)
    not_seen = Set{Int}(1:n)
    seen = Set{Int}()
    cycles = Array{Int,1}[]
    while !isempty(not_seen)
        cycle = Int[]
        x = pop!(not_seen)
        while !in(x, seen)
            push!(cycle, x)
            push!(seen, x)
            x = p[x]
            pop!(not_seen, x, 0)
        end
        push!(cycles, cycle)
    end
    cycle_lengths = map(length, cycles)
    even_cycles = filter(i -> i % 2 == 0, cycle_lengths)
    length(even_cycles) % 2 == 0 ? 1 : -1
end

struct Partition
    """
    2-partition of 2N-legs of Green's functions
    """
    l::Vector{Vector{Int}}
    r::Vector{Vector{Int}}
    sign::Vector{Int}
    function Partition(order)
        N = 2^order - 2
        function addincoming(l, r)
            # map 1, 2, 3, ... to order+1, order+2, ...
            l = l .+ order
            r = r .+ order
            nl, nr = length(l), length(r)
            inl = [i for i in 1:nl]
            inr = [i for i in nl+1:order]
            append!(inl, l)
            append!(inr, r)
            return inl, inr
        end

        left = Vector{Vector{Int}}(undef, N)
        right = similar(left)
        sign = Vector{Int}(undef, N)
        for (si, s) in enumerate(partitions(1:order, 2))
            # generate 2-partition of the outgoing legs
            # all incoming legs are kept the same
            s1, s2 = addincoming(s[1], s[2])
            left[2si-1], right[2si-1] = s1, s2
            sign[2si-1] = parity(vcat(s1, s2))

            q1, q2 = addincoming(s[2], s[1])
            left[2si], right[2si] = q1, q2
            sign[2si] = parity(vcat(q1, q2))
            # println("($s1, $s2) -> $p1, ($q1, $q2) -> $p2")
        end
        return new(left, right, sign)
    end
end

"""
    GreenN(m::Model, τ, orbital, isfermi=true)

Construct struct to store the variables to evaulate N-body Green's functions.
The leg index is assumed to be [1, 2, 3, ...,2N], where the incoming legs are 1:N, and the outgoing legs are N+1:2N
The full Green's function is defined as,
```math
G_{2N}=\\left<\\mathcal{T} c_{2N}c_{2N-1}...c_{N+1}c^+_{N}...c^+_{2}c^+_{1}\\right>
```

e.g.,

    1->------->-3
        | G4 |
    2->------->-4

All other Green's function are derived from the above full Green's function

#Arguments
- m: Model struct
- τ: array stores the imaginary-time of legs
- orbital: array stores the orbital/spin of legs
- N: particle number (Note that the leg numebr is 2N)
"""
struct GreenN
    N::Int # n-body or 2n-point
    τ::Vector{Float} # 1:2n, array of imaginary-time
    orbital::Vector{Int} # 1:2n, array of orbitals
    hop::Vector{Operator}
    partition::Vector{Partition}
    function GreenN(m::Model, τ, orbital, N=Int(length(τ) / 2))
        # @assert length(τ) == length(orbital)
        # N = Int(length(τ) / 2)

        hop = Vector{Operator}(undef, 2N)
        for idx in 1:N
            hop[idx] = Heisenberg(m.c⁺[orbital[idx]], m.E, τ[idx])
        end
        for idx in N+1:2N
            hop[idx] = Heisenberg(m.c⁻[orbital[idx]], m.E, τ[idx])
        end

        partition = [Partition(o) for o in 2:N]
        return new(N, τ, orbital, hop, partition)
    end
end

function density(m::Model, orbital)
    n = 0.0
    for (ei, e) in enumerate(m.E)
        n += exp(-m.β * e) * m.n[orbital][ei, ei]
    end
    return n / m.Z
end

"""
    @inline function Gn(m::Model, g::GreenN)

    Evaluate the full N-body Green's function.

```math
G_{2N}=\\left<\\mathcal{T} c_{2N}c_{2N-1}...c_{N+1}c^+_{N}...c^+_{2}c^+_{1}\\right>
```

e.g.,

    1->------->-3
        | G4 |
    2->------->-4

"""
@inline function Gn(m::Model, g::GreenN)
    return Gn(m, g, [i for i in 1:g.N*2])
end

function Gn(m::Model, g::GreenN, idx, level=1)
    # printstyled("...."^level, "start $level Gn-$idx\n", color=:green)
    if g.N == 1 # fast treatment of one-body Green's function
        i, o = idx[1], idx[2]
        # println(g.hop[o])
        # println(g.hop[i])
        if g.τ[i] < g.τ[o]
            return thermalavg(g.hop[o] * g.hop[i], m.E, m.β, m.Z)
        else
            G = thermalavg(g.hop[i] * g.hop[o], m.E, m.β, m.Z)
            m.isfermi && (G *= -1.0)
            # printstyled("...."^level, "get $level Gn-$idx -> $(idx[perm])*$(parity(perm))   $G\n", color=:green)
            return G
        end
    else
        perm = sortperm(g.τ[idx])
        # if length(idx) == 4
        #     println(g.τ[idx], "  and  ", g.orbital[idx], " --> ", perm)
        # end
        # M=c(τ_{2N})...c(2)c(1)
        M = prod(g.hop[reverse(idx[perm])])  # idx[perm] gives the idx array with the correct time ordering
        # note that the storage order of hop operators are the opposite of the phyical order of the operators
        G = thermalavg(M, m.E, m.β, m.Z) # G=tr(exp(-βH) M)/Z
        m.isfermi && (G *= parity(perm))
        # printstyled("...."^level, "get $level Gn-$idx -> $(idx[perm])*$(parity(perm))   $G\n", color=:green)
        return G
    end
end

"""
    function Gnc(m::Model, g::GreenN)

    Evaluate the connected N-body Green's function.

e.g., G4c(43;21) =
    1->------->-4       1->------->-4     1->---   --->-4
        | G4 |      -                  +         X
    2->------->-3       2->------->-3     2->---   --->-3
"""
@inline function Gnc(m::Model, g::GreenN)
    return Gnc(m, g, [i for i in 1:g.N*2])
end

function Gnc(m::Model, g::GreenN, idx, level=1)
    # printstyled("...."^level, "start $level Gc-$idx\n", color=:red)
    if length(idx) == 2 # case of G2
        return Gn(m, g, idx, level)
    end

    G = Gn(m, g, idx, level)

    order = Int(length(idx) / 2)
    p = g.partition[order-1]
    for (li, l) in enumerate(p.l)
        # printstyled("...."^level, "begin $(level): ", idx, " => ", idx[l], ", ", idx[p.r[li]], ": parity: ", 1, " --> ", p.sign[li], "\n", color=:yellow)
        disconnected = Gnc(m, g, idx[l], level + 1) * Gn(m, g, idx[p.r[li]], level + 1)
        m.isfermi && (disconnected *= p.sign[li])
        G -= disconnected
        # printstyled("...."^level, "end $(level)\n", color=:yellow)
    end
    # printstyled("...."^level, "end $level Gc-$idx $G\n", color=:red)
    return G
end

function G2(m::Model, g::GreenN)
    if g.N != 1
        throw(AssertionError("G2 function is for n=2"))
    end
    τi, τo = g.τ[1], g.τ[2]
    β, E = m.β, m.E
    c⁺, c⁻ = m.c⁺[g.orbital[1]], m.c⁻[g.orbital[2]] # in is creation operator, out is anniliation operator

    G = 0.0
    if (τi < τo)
        for j in 1:m.dim
            for k in 1:m.dim
                G += exp(-(β - τo + τi) * E[j] - (τo - τi) * E[k]) * c⁻[j, k] * c⁺[k, j]
            end
        end
    else # equal-time means 0^-
        for j in 1:m.dim
            for k in 1:m.dim
                G += -exp(-(β - τi + τo) * E[j] - (τi - τo) * E[k]) * c⁺[j, k] * c⁻[k, j]
            end
        end
    end
    return G / m.Z
end

end


text:
for i in 1:2:9
     println(i)
end #> 1 3 5 7 9

search("Julia","uli") #> 2:4
typeof(1:1000) #> UnitRange{Int64}

a = split("A,B,C,D",",")
typeof(a) #> Array{SubString{ASCIIString},1}
show(a) #> SubString{ASCIIString}["A","B","C","D"]

arr = [100, 25, 37]
arra = Any[100, 25, "ABC"] #> element type Any
arr[1] #> 100
arr[end] #> 37
# arr[6] #> BoundsError
println(eltype(arr)) #> Int64
println(length(arr)) #> 3
println(ndims(arr)) #> 1
println(size(arr,1)) #> 3
println(size(arr)) #> (3,)

arr2 = Array(Int64,5)
show(arr2) #> [0,0,0,0,0]
sizehint(arr2, 10^5)

arr3 = Float64[] #> 0-element Array{Float64,1}
push!(arr3, 1.0) #> 1-element Array{Float64,1}
arr4 = [1:7] #> 7-element Array{Int64,1}: [1,2,3,4,5,6,7]
join(arr4, ", ") #> "1, 2, 3, 4, 5, 6, 7"
arr4[1:3] # => [1, 2, 3]
arr4[2:end] # => [2, 3, 4, 5, 6, 7]  , b[2:] is deprecated

b = [1:7]
c = [100,200,300]
append!(b,c) # Now b is [1, 2, 3, 4, 5, 6, 7, 100, 200, 300]
pop!(b) #> 300, b is now [1, 2, 3, 4, 5, 6, 7, 100, 200]
push!(b, 42) # b is now [1, 2, 3, 4, 5, 6, 7, 100, 200, 42]
shift!(b) #> 1, b is now [2, 3, 4, 5, 6, 7, 100, 200, 42]
unshift!(b, 42) # b is now [42, 2, 3, 4, 5, 6, 7, 100, 200, 42]
splice!(b,8) #> 100, b is now [42, 2, 3, 4, 5, 6, 7, 200, 42]
in(42, b) #> true
in(43, b) #> false
println()
println("sorting:")
sort(b) #> [2,3,4,5,6,7,42,42,200], b is not changed
println(b) #> [42,2,3,4,5,6,7,200,42]
sort!(b) #> b is now changed to [2,3,4,5,6,7,42,42,200]
println(b) #> [2,3,4,5,6,7,42,42,200]
println()
arr = [1, 5, 3]
# looping
for e in arr
    print("$e ")
end # prints 1 5 3
arr + 2 #> [3, 7, 5]
arr * 2 #> [2, 10, 6]

a1 = [1, 2, 3]
a2 = [4, 5, 6]
a1 .* a2 #> [4, 10, 18]
dot(a1, a2) #> 32

arr = [1,2,3,4,5]
arr[2:4] = [8,9,10]
println(arr) #> 1 8 9 10 5

zeros(5)
ones(4)
ones(3, 2)
eqa = linspace(0, 10, 5)
println()
show(eqa) #> [0.0,2.5,5.0,7.5,10.0]
fill!(arr, 42) #> [42,42,42]
println()
println(cell(4)) #> {#undef,#undef,#undef,#undef}

v1 = rand(Int32, 5)
println(v1) #> Int32[1735276173,972339632,1303377282,1493859467,-788555652]
int64(v1) #> converts elements to 64 bit integers

a = [1,2,4,6]
a1 = a
show(a1)
a[4] = 0
show(a) #> [1,2,4,0]
show(a1) #> [1,2,4,0]
b = copy(a)
b = deepcopy(a)

a = [1,2,3]
function change_array(arr)
	arr[2] = 25
end
change_array(a)
println(a) #>[ 1, 25, 3]

# the splat operator:
arr = ['a', 'b', 'c']
show(join(arr)) #> "abc"
show(utf32(arr)) #> "abc"
show(string(arr)) #> "['a', 'b', 'c']"
show(string(arr...)) #> "abc"

max_stars_repo_path: Module 1/Chapter02/arrays.jl

max_issues_repo_path: Module 1/Chapter02/arrays.jl

max_forks_repo_path: Module 1/Chapter02/arrays.jl


text:
using FractionalDiffEq

K = 1
fdorder = 1.9
dx = pi/20
dt = 0.01
n = 2
xStart = 0
xEnd = pi
x = collect(0:dx:xEnd)
t = collect(0:dt:n)
u0t = 0
uendt = 0
u0(x) = sin(x)

using Plots
plotlyjs()
U=solve(fdorder, dx, dt, xStart, xEnd, n, K, u0t, uendt, u0, FiniteDiffEx())
plt=plot(x, t, U, st=:surface)


text:
mutable struct BernoulliDist
    _mu::Float64
    _beta::BetaDist
    # if use BetaDist, _bayes == true
    _bayes::Bool
    function BernoulliDist(mu::AbstractArray{Float64,1})
        _mu::Float64 = sum(mu) / size(mu)[1]
        new(_mu, BetaDist(0, 0), false)
    end
    function BernoulliDist(beta::BetaDist)
        new(0, beta, true)
    end
end

function pdf(dist::BernoulliDist, x::Int64)
    return dist._mu^x * (1.0 - dist._mu)^(1 - x)
end

function fitting(dist::BernoulliDist, X::AbstractArray{Float64,1})
    if dist._bayes == false
        return fitting_ml(dist, X)
    else
        return fitting_map(dist, X)
    end
end

function fitting_ml(dist::BernoulliDist, X::AbstractArray{Float64,1})
    # X is the AbstractArray of 0 or 1
    # X .== 0 is the AbstractArray indicating if X[i] == 0
    n_zeros = sum(X .== 0.0)
    # X .== 1 is the AbstractArray indicating if X[i] == 1
    n_ones = sum(X .== 1.0)
    dist._mu = n_ones * 1.0 / (n_ones + n_zeros)
end

function fitting_map(dist::BernoulliDist, X::AbstractArray{Float64,1})
    # X is the AbstractArray of 0 or 1
    # X .== 0 is the AbstractArray indicating if X[i] == 0
    n_zeros = sum(X .== 0.0)
    # X .== 1 is the AbstractArray indicating if X[i] == 1
    n_ones = sum(X .== 1.0)
    dist._beta._n_ones += n_ones
    dist._beta._n_zeros += n_zeros
end

function draw(bern::BernoulliDist, n_trials::Int64)
    trials = rand(Uniform(0, 1), n_trials)

    mu = 0
    if (bern._bayes == true)
        # sample from beta
        mu = bern._beta._n_ones / (bern._beta._n_ones + bern._beta._n_zeros)
    else
        # sample ml
        mu = bern._mu
    end

    return sum(mu .> trials)
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
@testset "396.rotate-function.jl" begin
    @test max_rotate_function([4, 3, 2, 6]) == 26
end


text:
### A Pluto.jl notebook ###
# v0.12.10

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 15a4ba3e-f0d1-11ea-2ef1-5ff1dee8795f
using Pkg

# ╔═╡ 21e744b8-f0d1-11ea-2e09-7ffbcdf43c37
begin
	Pkg.activate(mktempdir())

	Pkg.add("Gadfly")
	Pkg.add("Compose")
	Pkg.add("Statistics")
	Pkg.add("Hyperscript")
	Pkg.add("Colors")
	Pkg.add("Images")
	Pkg.add("ImageMagick")
	Pkg.add("ImageFiltering")

	using Gadfly
	using Images
	using Compose
	using Hyperscript
	using Colors
	using Statistics
	using PlutoUI
	using ImageMagick
	using ImageFiltering
end

# ╔═╡ 10f850fc-f0d1-11ea-2a58-2326a9ea1e2a
set_default_plot_size(12cm, 12cm)

# ╔═╡ 7b4d5270-f0d3-11ea-0b48-79005f20602c
function convolve(M, kernel)
    height, width = size(kernel)

    half_height = height ÷ 2
    half_width = width ÷ 2

    new_image = similar(M)

    # (i, j) loop over the original image
	m, n = size(M)
    @inbounds for i in 1:m
        for j in 1:n
            # (k, l) loop over the neighbouring pixels
			accumulator = 0 * M[1, 1]
			for k in -half_height:-half_height + height - 1
				for l in -half_width:-half_width + width - 1
					Mi = i - k
					Mj = j - l
					# First index into M
					if Mi < 1
						Mi = 1
					elseif Mi > m
						Mi = m
					end
					# Second index into M
					if Mj < 1
						Mj = 1
					elseif Mj > n
						Mj = n
					end

					accumulator += kernel[k, l] * M[Mi, Mj]
				end
			end
			new_image[i, j] = accumulator
        end
    end

    return new_image
end

# ╔═╡ 6fd3b7a4-f0d3-11ea-1f26-fb9740cd16e0
function disc(n, r1=0.8, r2=0.8)
	white = RGB{Float64}(1,1,1)
	blue = RGB{Float64}(colorant"#4EC0E3")
	convolve(
		[(i-n/2)^2 + (j-n/2)^2 <= (n/2-5)^2 ? white : blue for i=1:n, j=1:n],
		Kernel.gaussian((1,1))
	)
end

# ╔═╡ fe3559e0-f13b-11ea-06c8-a314e44c20d6
brightness(c) = 0.3 * c.r + 0.59 * c.g + 0.11 * c.b

# ╔═╡ 0ccf76e4-f0d9-11ea-07c9-0159e3d4d733
@bind img_select Radio(["disc", "mario"], default="disc")

# ╔═╡ 236dab08-f13d-11ea-1922-a3b82cfc7f51
begin
	url = "http://files.softicons.com/download/game-icons/super-mario-icons-by-sandro-pereira/png/32/Retro%20Mario.png"
	img = Dict(
		"disc" => disc(25),
		"mario" => load(download(url))
	)[img_select]
end

# ╔═╡ 03434682-f13b-11ea-2b6e-11ad781e9a51
md"""Show $G_x$ $(@bind Gx CheckBox())

     Show $G_y$ $(@bind Gy CheckBox())"""

# ╔═╡ ca13597a-f168-11ea-1a2c-ff7b98b7b2c7
function partial_derivatives(img)
	Sy,Sx = Kernel.sobel()
	∇x, ∇y = zeros(size(img)), zeros(size(img))

	if Gx
		∇x = convolve(brightness.(img), Sx)
	end
	if Gy
		∇y = convolve(brightness.(img), Sy)
	end
	return ∇x, ∇y
end

# ╔═╡ b369584c-f183-11ea-260a-35dc797e63ad

# ╔═╡ b2cbe058-f183-11ea-39dc-23d4a5b92796

# ╔═╡ 9d9cccb2-f118-11ea-1638-c76682e636b2
function arrowhead(θ)
	eq_triangle = [(0, 1/sqrt(3)),
		           (-1/3, -2/(2 * sqrt(3))),
		           (1/3, -2/(2 * sqrt(3)))]

	compose(context(units=UnitBox(-1,-1,2,2), rotation=Rotation(θ, 0, 0)),
				polygon(eq_triangle))
end

# ╔═╡ b7ea8a28-f0d7-11ea-3e98-7b19a1f58304
function quiver(points, vecs)
	xmin = minimum(first.(points))
	ymin = minimum(last.(points))
	xmax = maximum(first.(points))
	ymax = maximum(last.(points))
	hs = map(x->hypot(x...), vecs)
	hs = hs / maximum(hs)

	vector(p, v, h) = all(iszero, v) ? context() :
		(context(),
		    (context((p.+v.*6 .- .2)..., .4,.4),
				arrowhead(atan(v[2], v[1]) - pi/2)),
		stroke(RGBA(90/255,39/255,41/255,h)),
		fill(RGBA(90/255,39/255,41/255,h)),
		line([p, p.+v.*8]))

	compose(context(units=UnitBox(xmin,ymin,xmax,ymax)),
         vector.(points, vecs, hs)...)
end

# ╔═╡ c821b906-f0d8-11ea-2df0-8f2d06964aa2
function sobel_quiver(img, ∇x, ∇y)
    quiver([(j-1,i-1) for i=1:size(img,1), j=1:size(img,2)],
           [(∇x[i,j], ∇y[i,j]) for i=1:size(img,1), j=1:size(img,2)])
end

# ╔═╡ 6da3fdfe-f0dd-11ea-2407-7b85217b35cc
# render an Image using squares in Compose
function compimg(img)
	xmax, ymax = size(img)
	xmin, ymin = 0, 0
	arr = [(j-1, i-1) for i=1:ymax, j=1:xmax]

	compose(context(units=UnitBox(xmin, ymin, xmax, ymax)),
		fill(vec(img)),
		rectangle(
			first.(arr),
			last.(arr),
			fill(1.0, length(arr)),
			fill(1.0, length(arr))))
end

# ╔═╡ f22aa34e-f0df-11ea-3053-3dcdc070ec2f
let
	∇x, ∇y = partial_derivatives(img)

	compose(context(),
		sobel_quiver(img, ∇x, ∇y),
		compimg(img))
end

# ╔═╡ 885ec336-f146-11ea-00c4-c1d1ab4c0001
	function show_colored_array(array)
		pos_color = RGB(0.36, 0.82, 0.8)
		neg_color = RGB(0.99, 0.18, 0.13)
		to_rgb(x) = max(x, 0) * pos_color + max(-x, 0) * neg_color
		to_rgb.(array) / maximum(abs.(array))
	end

# ╔═╡ 9232dcc8-f188-11ea-08fe-b787ea93c598
begin
	Sy, Sx = Kernel.sobel()
	show_colored_array(Sx)
	Sx
end

# ╔═╡ 7864bd00-f146-11ea-0020-7fccb3913d8b
let
	∇x, ∇y = partial_derivatives(img)

	to_show = (x -> RGB(0, 0, 0)).(zeros(size(img)))
	if Gx && Gy
		edged = sqrt.(∇x.^2 + ∇y.^2)
		to_show = Gray.(edged) / maximum(edged)
	elseif Gx
		to_show = show_colored_array(∇x)
	elseif Gy
		to_show = show_colored_array(∇y)
	end
	compose(
		context(),
		compimg(to_show)
	)
end

# ╔═╡ Cell order:
# ╠═15a4ba3e-f0d1-11ea-2ef1-5ff1dee8795f
# ╟─21e744b8-f0d1-11ea-2e09-7ffbcdf43c37
# ╠═10f850fc-f0d1-11ea-2a58-2326a9ea1e2a
# ╟─7b4d5270-f0d3-11ea-0b48-79005f20602c
# ╠═6fd3b7a4-f0d3-11ea-1f26-fb9740cd16e0
# ╟─fe3559e0-f13b-11ea-06c8-a314e44c20d6
# ╟─b7ea8a28-f0d7-11ea-3e98-7b19a1f58304
# ╟─0ccf76e4-f0d9-11ea-07c9-0159e3d4d733
# ╟─236dab08-f13d-11ea-1922-a3b82cfc7f51
# ╟─03434682-f13b-11ea-2b6e-11ad781e9a51
# ╟─ca13597a-f168-11ea-1a2c-ff7b98b7b2c7
# ╟─f22aa34e-f0df-11ea-3053-3dcdc070ec2f
# ╟─9232dcc8-f188-11ea-08fe-b787ea93c598
# ╠═7864bd00-f146-11ea-0020-7fccb3913d8b
# ╠═b369584c-f183-11ea-260a-35dc797e63ad
# ╠═b2cbe058-f183-11ea-39dc-23d4a5b92796
# ╟─9d9cccb2-f118-11ea-1638-c76682e636b2
# ╟─c821b906-f0d8-11ea-2df0-8f2d06964aa2
# ╟─6da3fdfe-f0dd-11ea-2407-7b85217b35cc
# ╠═885ec336-f146-11ea-00c4-c1d1ab4c0001

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
using RunTests
using Base.Test

@testmodule LinearRegressionTest begin

  function test_linearRegression()
  #  @test false
  end

end

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -8.232182974914569 <= q <= 6.070970235703473)

                     add_NL_constraint(m, :(gelu(-0.17705543456567296 + 0.5634748900178774*gelu(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])) + 0.6106029608527352*gelu(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])) + 0.5384388086098779*gelu(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))) + gelu(0.08053767290922664 + 0.911477397877237*gelu(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])) + 0.7134526118957916*gelu(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])) + -0.34097520379571256*gelu(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(x[3]))) + gelu(0.5992621446321471 + 0.15304126986204158*gelu(-0.9628319588335761 + 0.5826659544327617*$(x[1]) + -0.6099922924008303*$(x[2]) + 0.9101904633653426*$(x[3])) + -0.37134199071135*gelu(0.06629760484384972 + -0.7340864295918403*$(x[1]) + 0.2309908240081593*$(x[2]) + 0.7867070707759196*$(x[3])) + 0.8747239793986465*gelu(-0.07370832879656497 + -0.32073040453400337*$(x[1]) + -0.07649911788086428*$(x[2]) + 0.03499506906851524*$(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


text:
# ============================================================================
# Composition operation for one step of an approximation model with one step
# backward of the same model
# ============================================================================

"""
    StepIntersect{DM<:AbstractApproximationModel} <: AbstractApproximationModel

Approximation model that composes (intersecting) one step Forward of a given
model with one step backward of the same model.

### Fields

- `model`  -- approximation model
- `setops` -- set operations method

### Notes

Let ``x' = Ax`` with ``x(0) ∈ X₀``. This methods consists of:

- Compute the discretized system with step-size ``δ`` obtaining ``Ω0`` and the given
  approxatmion model `method`.

- Compute the (lazy) linear map ``ΦX₀``. This set constains of the (exact) reachable
  states at the time point ``δ``.

- Apply the approximation model `method` with initial condition ``ΦX₀`` one step
  backward in time, with the state transition matrix ``-A``, call this reach-set
  ``Ω₀₊``

- Intersect ``Ω₀₋`` and ``Ω₀₊`` and return such set. The intersection is done either
  lazily or concretely depending on the specified `setops` field.
"""
struct StepIntersect{DM<:AbstractApproximationModel, SO} <: AbstractApproximationModel
    model::DM
    setops::SO
end

hasbackend(alg::StepIntersect) = !isnothing(alg.model)

# convenience constructor using symbols
function StepIntersect(alg=Forward(); setops=alg.setops)
    return StepIntersect(alg, _alias(setops))
end

function Base.show(io::IO, alg::StepIntersect)
    print(io, "`StepIntersect` approximation model with: \n")
    print(io, "    - model: $(alg.model) \n")
    print(io, "    - set operations method: $(alg.setops)\n")
end

Base.show(io::IO, m::MIME"text/plain", alg::StepIntersect) = print(io, alg)

# ------------------------------------------------------------
# Homogeneous case
# ------------------------------------------------------------

# x' = Ax, x(0) ∈ X₀, x ∈ X
function discretize(ivp::IVP{<:CLCS, <:LazySet}, δ, alg::StepIntersect)
    A = state_matrix(ivp)
    X0 = initial_state(ivp)
    X = stateset(ivp)

    ivpd = discretize(ivp, δ, alg.model)
    Φ = state_matrix(ivpd)
    Ω0₊ = initial_state(ivpd)

    Aneg = -A
    ΦX0 = Φ * X0
    ivpneg = IVP(CLCS(Aneg, X), ΦX0)
    ivpnegd = discretize(ivpneg, δ, alg.model)
    Ω0₋ = initial_state(ivpnegd)

    Ω0 = _apply_setops(Ω0₊ ∩ Ω0₋, alg.setops)
    return IVP(CLDS(Φ, X), Ω0)
end

# ------------------------------------------------------------
# Inhomogeneous case
# ------------------------------------------------------------


text:
############################################################
# Utilities ################################################
############################################################

struct joUtilsException <: Exception
    msg :: String
end

############################################################
## jo_eye/speye/full #######################################
export jo_eye, jo_speye, jo_full
"""
return identity array

    jo_eye(m::Integer,n::Integer=m)
    jo_eye(DT::DataType,m::Integer,n::Integer=m)

"""
@inline jo_eye(m::Integer,n::Integer=m) = Matrix{joFloat}(LinearAlgebra.I,m,n)
@inline jo_eye(DT::DataType,m::Integer,n::Integer=m) = Matrix{DT}(LinearAlgebra.I,m,n)
"""
return sparse identity array

    jo_speye(m::Integer,n::Integer=m)
    jo_speye(DT::DataType,m::Integer,n::Integer=m)

"""
@inline jo_speye(m::Integer,n::Integer=m) = sparse(one(joFloat)*LinearAlgebra.I,m,n)
@inline jo_speye(DT::DataType,m::Integer,n::Integer=m) = sparse(one(DT)*LinearAlgebra.I,m,n)
"""
return full array

    jo_full(A::AbstractArray)

"""
@inline jo_full(A::AbstractArray) = Array(A)

############################################################
## jo_method_error #########################################
function jo_method_error(O::joAbstractOperator,s::String)
    op_type=typeof(O)
    op_name=O.name
    @info "JOLI method error for operator:" op_name op_type
    throw(joUtilsException(s))
end
function jo_method_error(L,R,s::String)
    left_type=typeof(L)
    left_name=typeof(L)<:joAbstractOperator ? L.name : "non-JOLI"
    right_type=typeof(R)
    right_name=typeof(R)<:joAbstractOperator ? R.name : "non-JOLI"
    @info "JOLI method error for combination:" left_name left_type right_name right_type
    throw(joUtilsException(s))
end
function jo_method_error(L,M,R,s::String)
    left_type=typeof(L)
    left_name=typeof(L)<:joAbstractOperator ? L.name : "non-JOLI"
    mid_type=typeof(M)
    mid_name=typeof(M)<:joAbstractOperator ? M.name : "non-JOLI"
    right_type=typeof(R)
    right_name=typeof(R)<:joAbstractOperator ? R.name : "non-JOLI"
    @info "JOLI method error for combination:" left_name left_type mid_name mid_type right_name right_type
    throw(joUtilsException(s))
end

############################################################
## type tree ###############################################
export jo_type_tree
"""
Show type tree of JOLI operators, or any type with editional Type argument.

    jo_type_tree()
    jo_type_tree(Number)
"""
function jo_type_tree(top::Type=joAbstractOperator;bl::String="* ",in::String="  ",super::Bool=true)
    super ? println(bl,top," <: ",supertype(top)) : println(bl,top)
    ts = subtypes(top)
    if length(ts) > 0
        for t in ts
            T=Base.unwrap_unionall(t)
            jo_type_tree(T;bl=in*bl,in=in,super=false)
        end
    end
end

############################################################
## default types ###########################################
export joInt, joFloat, joComplex
global joInt=Int64
global joFloat=Float64
global joComplex=ComplexF64
export jo_joInt_set, jo_joFloat_set, jo_joComplex_set, jo_jo32bit_set, jo_jo64bit_set, jo_joTypes_get
"""
set default integer type joInt

    function jo_joInt_set(DT::DataType=joInt)

"""
function jo_joInt_set(DT::DataType=joInt)
    global joInt=DT
    return joInt
end
"""
set default float type joFloat

    function jo_joFloat_set(DT::DataType=joFloat)

"""
function jo_joFloat_set(DT::DataType=joFloat)
    global joFloat=DT
    return joFloat
end
"""
set default complex type joComplex

    function jo_joComplex_set(DT::DataType=joComplex)

"""
function jo_joComplex_set(DT::DataType=joComplex)
    global joComplex=DT
    return joComplex
end
"""
set default type joInt, joFloat, joComplex to 32 bit

    function jo_jo32bit_set()

"""
function jo_jo32bit_set()
    global joInt=Int32
    global joFloat=Float32
    global joComplex=ComplexF32
    global joTol=sqrt(eps(Float32))
    return joInt, joFloat, joComplex
end
"""
set default type joInt, joFloat, joComplex to 64 bit

    function jo_jo64bit_set()

"""
function jo_jo64bit_set()
    global joInt=Int64
    global joFloat=Float64
    global joComplex=ComplexF64
    global joTol=sqrt(eps(Float64))
    return joInt, joFloat, joComplex
end

"""
get default types joInt, joFloat, joComplex

    function jo_joTypes_get()

"""
jo_joTypes_get() = (joInt, joFloat, joComplex)

############################################################
## macros ##################################################
export @joNF
"""
Nullable{Function} macro for null function

    @joNF
"""
macro joNF()
    return :(Nullable{Function}())
end
"""
Nullable{Function} macro for given function

    @joNF ... | @joNF(...)
"""
macro joNF(fun::Expr)
    return :(Nullable{Function}($fun))
end

############################################################
## dafault iterative solver for square operator ############
global jo_iterative_solver4square = (A,v)->gmres(A,v)
export jo_iterative_solver4square_set
"""
Set default iterative solver for \\(jo,vec) and square jo

    jo_iterative_solver4square_set(f::Function)

Where f must take two arguments (jo,vec) and return vec.

# Example (using IterativeSolvers)
- jo_iterative_solver4square_set((A,v)->gmres(A,v))

"""
function jo_iterative_solver4square_set(f::Function)
    global jo_iterative_solver4square
    jo_iterative_solver4square = (A,b)->f(A,b)
end

############################################################
## dafault iterative solver for tall operator ##############
global jo_iterative_solver4tall = @joNF
export jo_iterative_solver4tall_set
"""
Set default iterative solver for \\(jo,vec) and tall jo

    jo_iterative_solver4tall_set(f::Function)

Where f must take two arguments (jo,vec) and return vec.

# Example
- jo_iterative_solver4tall_set((A,v)->tall_solve(A,v))

"""
function jo_iterative_solver4tall_set(f::Function)
    global jo_iterative_solver4tall
    jo_iterative_solver4tall = (A,b)->f(A,b)
end

############################################################
## dafault iterative solver for wide operator ##############
global jo_iterative_solver4wide = @joNF
export jo_iterative_solver4wide_set
"""
Set default iterative solver for \\(jo,vec) and wide jo

    jo_iterative_solver4wide_set(f::Function)

Where f must take two arguments (jo,vec) and return vec.

# Example
- jo_iterative_solver4wide_set((A,v)->wide_solve(A,v))

"""
function jo_iterative_solver4wide_set(f::Function)
    global jo_iterative_solver4wide
    jo_iterative_solver4wide = (A,b)->f(A,b)
end

############################################################
## smart precision type ####################################
export jo_precision_type
"""
Type of the real number or element type of complex number.

# Example
- jo_precision_type(1.)
- jo_precision_type(1+im*3.)
"""
jo_precision_type(x::Tx) where {ITx<:Number, Tx<:Union{Complex{ITx}, ITx}} = ITx

############################################################
## complex precision type ##################################
export jo_complex_eltype
"""
Type of element of complex scalar

    jo_complex_eltype(a::Complex)

# Example
- jo_complex_eltype(1.+im*1.)
- jo_complex_eltype(zero(ComplexF64))

"""
jo_complex_eltype(a::Complex{T}) where {T} = T
"""
Type of element of complex data type

    jo_complex_eltype(DT::DataType)

# Example
- jo_complex_eltype(ComplexF32)

"""
function jo_complex_eltype(DT::DataType)
    DT<:Complex || throw(joUtilsException("jo_complex_eltype: Input type must be Complex"))
    a=zero(DT)
    return jo_complex_eltype(a)
end

############################################################
## type checks #############################################
export jo_check_type_match, jo_type_mismatch_error_set
global jo_type_mismatch_warn=false
global jo_type_mismatch_error=true
"""
Toggle between warning and error for type mismatch

    jo_type_mismatch_error_set(flag::Bool)

# Examples
- jo_type_mismatch_error_set(false) turns on warnings instead of errors
- jo_type_mismatch_error_set(true) reverts to errors

"""
function jo_type_mismatch_error_set(flag::Bool)
    global jo_type_mismatch_warn
    global jo_type_mismatch_error
    if flag
        jo_type_mismatch_warn=false
        jo_type_mismatch_error=true
    else
        jo_type_mismatch_warn=true
        jo_type_mismatch_error=false
    end
    return
end
function jo_type_mismatch_warn_set(flag::Bool)
    global jo_type_mismatch_warn
    @warn "Very, very bad idea! You are a sneaky fellow."
    @warn "Function jo_type_mismatch_warn_set is deprecated and will be removed in the future."
    jo_type_mismatch_warn=flag
    return
end
"""
Check type match

    jo_check_type_match(DT1::DataType,DT2::DataType,where::String)

The bahaviour of the function while types do not match depends on
values of jo_type_mismatch_warn and jo_type_mismatch_error flags.
Use jo_type_mismatch_error_set to toggle those flags from warning
mode to error mode.

# EXAMPLE
- jo_check_type_match(Float32,Float64,"my session")

"""
function jo_check_type_match(DT1::DataType,DT2::DataType,where::String)
    if !(DT1==DT2)
        jo_type_mismatch_error && throw(joUtilsException("type mismatch: expected $DT1 instead of $DT2 in $where"))
        jo_type_mismatch_warn && @warn "type mismatch: $DT1 vs. $DT2 in $where"
    end
    return
end

############################################################
## type conversion utlis ###################################
export jo_convert, jo_convert_warn_set
global jo_convert_warn=true
"""
Set warning mode for jo_convert

    jo_convert_warn_set(flag::Bool)

# Example
- jo_convert_warn_set(false) turns of the warnings

"""
function jo_convert_warn_set(flag::Bool)
    global jo_convert_warn
    jo_convert_warn=flag
    return
end
"""
Convert vector to new type

    jo_convert(DT::DataType,v::AbstractArray,warning::Bool=true)

# Limitations
- converting integer array to shorter representation will throw an error
- converting float/complex array to integer will throw an error
- converting from complex to float drops immaginary part and issues warning;
  use jo_convert_warn_set(false) to turn off the warning

# Example
- jo_convert(ComplexF32,rand(3))

"""
function jo_convert(DT::DataType,vin::AbstractArray{VT},warning::Bool=true) where {VT<:Integer}
    DT==VT && return vin
    #println("jo_convert{VT<:Integer}")
    if DT<:Integer
        if typemax(DT)>typemax(VT)
            vout=convert(AbstractArray{DT},vin)
        else
            throw(joUtilsException("jo_convert: Refused conversion from $VT to $DT."))
        end
    else
        vout=convert(AbstractArray{DT},vin)
    end
    return vout
end
function jo_convert(DT::DataType,vin::AbstractArray{VT},warning::Bool=true) where {VT<:AbstractFloat}
    DT==VT && return vin
    #println("jo_convert{VT<:AbstractFloat}")
    if !(DT<:Integer)
        vout=convert(AbstractArray{DT},vin)
    else
        throw(joUtilsException("jo_convert: Refused conversion from $VT to $DT."))
    end
    return vout
end
function jo_convert(DT::DataType,vin::AbstractArray{VT},warning::Bool=true) where {VT<:Complex}
    DT==VT && return vin
    #println("jo_convert{VT<:Complex}")
    if DT<:Complex
        vout=convert(AbstractArray{DT},vin)
    elseif DT<:AbstractFloat
        (warning && jo_convert_warn) && @warn "jo_convert: Inexact conversion from $VT to $DT. Dropping imaginary part."
        vout=convert(AbstractArray{DT},real(vin))
    else
        throw(joUtilsException("jo_convert: Refused conversion from $VT to $DT."))
    end
    return vout
end
"""
Convert number to new type

    jo_convert(DT::DataType,n::Number,warning::Bool=true)

# Limitations
- converting integer number to shorter representation will throw an error
- converting float/complex number to integer will throw an error
- converting from complex to float drops immaginary part and issues warning;
  use jo_convert_warn_set(false) to turn off the warning

# Example
- jo_convert(ComplexF32,rand())

"""
function jo_convert(DT::DataType,nin::NT,warning::Bool=true) where {NT<:Integer}
    DT==NT && return nin
    #println("jo_convert{NT<:Integer}")
    if DT<:Integer
        if typemax(DT)>typemax(NT)
            nout=convert(DT,nin)
        else
            throw(joUtilsException("jo_convert: Refused conversion from $NT to $DT."))
        end
    else
        nout=convert(DT,nin)
    end
    return nout
end
function jo_convert(DT::DataType,nin::NT,warning::Bool=true) where {NT<:AbstractFloat}
    DT==NT && return nin
    #println("jo_convert{NT<:AbstractFloat}")
    if !(DT<:Integer)
        nout=convert(DT,nin)
    else
        throw(joUtilsException("jo_convert: Refused conversion from $NT to $DT."))
    end
    return nout
end
function jo_convert(DT::DataType,nin::NT,warning::Bool=true) where {NT<:Complex}
    DT==NT && return nin
    #println("jo_convert{NT<:Complex}")
    if DT<:Complex
        nout=convert(DT,nin)
    elseif DT<:AbstractFloat
        (warning && jo_convert_warn) && @warn "jo_convert: Inexact conversion from $NT to $DT. Dropping imaginary part."
        nout=convert(DT,real(nin))
    else
        throw(joUtilsException("jo_convert: Refused conversion from $NT to $DT."))
    end
    return nout
end


text:
const path = joinpath(Dolo.pkg_path, "examples", "models", "dynare", "EA_QUEST3.mod")
const sm = Dolo.load_modfile(path)
const ff = Dolang.FunctionFactory(sm, :dynare)

# functions we will benchmark
regex_create_symbolic_model(path) = Dolo.load_modfile(path)
parse_stringify(sm) = Dolang.stringify(sm.equations[:dynare])
build_function_factory(sm) = Dolang.FunctionFactory(sm, :dynare)
build_levels_func_body(ff) = Dolang.func_body(ff, Dolang.Der{0})
compute_jacobian_derivatives(ff) = Dolang._jacobian_expr_mat(ff)
build_jacobian_func_body(ff) = Dolang.func_body(ff, Dolang.Der{1})
compute_hessian_derivatives(ff) = Dolang._hessian_exprs(ff)
build_hessian_func_body(ff) = Dolang.func_body(ff, Dolang.Der{2})
solve_triangular_system(sm) = Dolo.solve_triangular_system(sm.calibration)
dynare_import(path) = Dolo.dynare_import(path)

# build benchmark group
suite["eaquest"] = BenchmarkGroup(["eaquest"])

for (f, arg) in [(regex_create_symbolic_model, path),
                 (parse_stringify, sm),
                 (build_function_factory, sm),
                 (build_levels_func_body, ff),
                 (compute_jacobian_derivatives, ff),
                 (build_jacobian_func_body, ff),
                 (compute_hessian_derivatives, ff),
                 (build_hessian_func_body, ff),
                 (solve_triangular_system, sm),
                 (dynare_import, path),
                 ]
    # remove `DoloBenchmarks` module prefix before each function name
    k = Symbol(split(string(f), ".")[2])
    suite["eaquest"][k]= @benchmarkable $(f)($(arg)) seconds=10
end


text:
import Base: isnan
import DataArrays: DataArray, isna, dropna, NA, NAtype
if VERSION < v"0.6.0"
  import DataArrays: array
else
  import DataArrays: isnan
end
import DataFrames: DataFrame, ncol, convert
export GLRM

import Missings: Missing, missing, ismissing

# TODO: identify categoricals automatically from PooledDataArray columns

function GLRM(df::DataFrame, k::Int;
              losses = Loss[], rx = QuadReg(.01), ry = QuadReg(.01),
              offset = true, scale = false,
              prob_scale = true, NaNs_to_Missing = true)
    if NaNs_to_Missing
        df = copy(df)
        NaNs_to_Missing!(df)
    end
    if losses == Loss[] # if losses not specified, identify ordinal, boolean and real columns
        # change the wal get_reals, etc work.

        #reals, real_losses = get_reals(df)
        #bools, bool_losses = get_bools(df)
        #ordinals, ordinal_losses = get_ordinals(df)
        #easier to use just one function for this usecase.
        reals, real_losses, bools, bool_losses, ordinals, ordinal_losses = get_loss_types(df)
        A = [df[reals] df[bools] df[ordinals]]
        labels = [names(df)[reals]; names(df)[bools]; names(df)[ordinals]]
        losses = [real_losses; bool_losses; ordinal_losses]
    else # otherwise require one loss function per column
        A = df
        ncol(df)==length(losses) ? labels = names(df) : error("please input one loss per column of dataframe")
    end
    # identify which entries in data frame have been observed (ie are not N/A)
    obs = observations(A)
    # initialize X and Y
    X = randn(k,size(A,1))
    Y = randn(k,size(A,2))

    # form model
    rys = Array{Regularizer}(length(losses))
    for i=1:length(losses)
        if isa(losses[i].domain, OrdinalDomain) && embedding_dim(losses[i])>1 #losses[i], MultinomialOrdinalLoss) || isa(losses[i], OrdisticLoss)
            rys[i] = OrdinalReg(copy(ry))
        else
            rys[i] = copy(ry)
        end
    end
    glrm = GLRM(A, losses, rx, rys, k, obs=obs, X=X, Y=Y, offset=offset, scale=scale)

    # scale model so it really computes the MAP estimator of the parameters
    if prob_scale
        prob_scale!(glrm)
    end
    return glrm, labels
end

function get_loss_types(df::DataFrame)
    m,n = size(df)
    reals = fill(false,n)
    bools = fill(false,n)
    ordinals = fill(false,n)

    for j in 1:n
        # assuming there are no columns with *all* values missing. (which would make it a redundant column)
        t = eltype(collect(skipmissing(df[j]))[1])
        if(t == Float64)
            reals[j] = true
        elseif t == Bool
            bools[j] = true
        elseif t == (Int || Int32 || Int64)
            ordinals[j] = true
        end
    end

    n1 = sum(reals)
    real_losses = Array{Loss}(n1)
    for i=1:n1
        real_losses[i] = QuadLoss()
    end

    n2 = sum(bools)
    bool_losses = Array{Loss}(n2)
    for i in 1:n2
        bool_losses[i] = HingeLoss()
    end

    n3 = sum(ordinals)
    ord_idx = (1:size(df,2))[ordinals]
    maxs = zeros(n3,1)
    mins = zeros(n3,1)
    for j in 1:n3
        col = df[ord_idx[j]]
        try
            maxs[j] = maximum(skipmissing(col))
            mins[j] = minimum(skipmissing(col))
        end
    end

    # set losses and regularizers
    ord_losses = Array{Loss}(n3)
    for j=1:n3
        ord_losses[i] = OrdinalHinge(mins[i],maxs[i])
    end
    return reals,real_losses,bools,bool_losses,ordinals,ord_losses
end

# leaving the functions below untouched for now.
function get_reals(df::DataFrame)
    m,n = size(df)
    reals = [typeof(df[i])<:AbstractArray{Float64,1} for i in 1:n]
    n1 = sum(reals)
    losses = Array{Loss}(n1)
    for i=1:n1
        losses[i] = QuadLoss()
    end
    return reals, losses
end

function get_bools(df::DataFrame)
    m,n = size(df)
    bools = [isa(df[i], AbstractArray{Bool,1}) for i in 1:n]
    n1 = sum(bools)
    losses = Array{Loss}(n1)
    for i=1:n1
        losses[i] = HingeLoss()
    end
    return bools, losses
end

function get_ordinals(df::DataFrame)
    m,n = size(df)
    # there must be a better way to check types...
    ordinals = [(isa(df[i], AbstractArray{Int,1}) ||
                 isa(df[i], AbstractArray{Int32,1}) ||
                 isa(df[i], AbstractArray{Int64,1})) for i in 1:n]
    nord = sum(ordinals)
    ord_idx = (1:size(df,2))[ordinals]
    maxs = zeros(nord,1)
    mins = zeros(nord,1)
    for i in 1:nord
        col = df[ord_idx[i]]
        try
            maxs[i] = maximum(dropna(col))
            mins[i] = minimum(dropna(col))
        end
    end

    # set losses and regularizers
    losses = Array{Loss}(nord)
    for i=1:nord
        losses[i] = OrdinalHinge(mins[i],maxs[i])
    end
    return ordinals, losses
end


text:
using Sockets

include("i3ipc.jl")

s = connect(read(`i3 --get-socketpath`, String) |> chomp)

write(s, Message(RUN_COMMAND, "exec notify-send Success"))


text:
# This file was generated, do not modify it. # hide
boston2 = coerce(boston, :Chas => OrderedFactor);