text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* File containg user models of columns
#*
#*
#* The default nomenclature is:
#*		Type_Column_reboilertype_condensertyper
#*
#* where:
#*	Type = refluxed or reboiled or section
#*	Column = Stripping, Absorption, Rectifier, Distillation
#*	Reboiler type (if exists) = kettle or thermosyphon
#*	Condenser type (if exists) = with subccoling or without subcooling
#*
#*-----------------------------------------------------------------------
#* Author: Based on Models written by Paula B. Staudt
#* $Id$
#*---------------------------------------------------------------------
# The complete documentation for these models needs to be updated !!!
#----------------------------------------------------------------------
#* Model of a  basic column section with:
#*	- NumberOfTrays=number of trays.
#*
#*---------------------------------------------------------------------
type User_Section_Column
	User_Section_Column()=begin
		new(
			User_Section_ColumnBasic(),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour Outlet in the section",
				:PosX=>1,
				:PosY=>0.35,
				:Protected=>true
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid Outlet in the section",
				:PosX=>1,
				:PosY=>0.65,
				:Protected=>true
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Liquid Inlet in the section",
				:PosX=>0.80,
				:PosY=>0
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour Outlet in the section",
				:PosX=>0.30,
				:PosY=>0
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Vapour Inlet in the section",
				:PosX=>0.30,
				:PosY=>1
			)),
			liquid_stream(Dict{Symbol,Any}(
				:Brief=>"Liquid Outlet in the section",
				:PosX=>0.80,
				:PosY=>1
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Liquid connection at the middle trays",
				:PosX=>0.75,
				:PosY=>1,
				:Hidden=>true
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Vapour connection at the middle trays",
				:PosX=>0.55,
				:PosY=>0,
				:Hidden=>true
			)),
			[
				:(LiquidConnector.F= LiquidInlet.F),
				:(LiquidConnector.T = LiquidInlet.T),
				:(LiquidConnector.P = LiquidInlet.P),
				:(LiquidConnector.z = LiquidInlet.z),
				:(LiquidConnector.v = LiquidInlet.v),
				:(LiquidConnector.h = LiquidInlet.h),
				:(VapourConnector.F= VapourInlet.F),
				:(VapourConnector.T = VapourInlet.T),
				:(VapourConnector.P = VapourInlet.P),
				:(VapourConnector.z = VapourInlet.z),
				:(VapourConnector.v = VapourInlet.v),
				:(VapourConnector.h = VapourInlet.h),
				:(LiquidOutlet.F= trays(_base_1.NumberOfTrays).OutletL.F),
				:(LiquidOutlet.T = trays(_base_1.NumberOfTrays).OutletL.T),
				:(LiquidOutlet.P = trays(_base_1.NumberOfTrays).OutletL.P),
				:(LiquidOutlet.z = trays(_base_1.NumberOfTrays).OutletL.z),
				:(VapourOutlet.F= trays(1).OutletV.F),
				:(VapourOutlet.T = trays(1).OutletV.T),
				:(VapourOutlet.P = trays(1).OutletV.P),
				:(VapourOutlet.z = trays(1).OutletV.z),
				:(VapourDrawOff.F*_base_1.VapSideTrayIndex= trays.VapourSideStream.F),
				:(VapourDrawOff.T = trays(_base_1.VapourSideStreamLocation).VapourSideStream.T),
				:(VapourDrawOff.P = trays(_base_1.VapourSideStreamLocation).VapourSideStream.P),
				:(VapourDrawOff.z = trays(_base_1.VapourSideStreamLocation).VapourSideStream.z),
				:(LiquidDrawOff.F*_base_1.LiqSideTrayIndex= trays.LiquidSideStream.F),
				:(LiquidDrawOff.T = trays(_base_1.LiquidSideStreamLocation).LiquidSideStream.T),
				:(LiquidDrawOff.P = trays(_base_1.LiquidSideStreamLocation).LiquidSideStream.P),
				:(LiquidDrawOff.z = trays(_base_1.LiquidSideStreamLocation).LiquidSideStream.z),
				:(VapourDrawOffFlow = VapourDrawOff.F),
				:(LiquidDrawOffFlow = LiquidDrawOff.F),
			],
			[
				"","","","","","","","","","","","","","","","","","","","","","","","","","","","","","",
			],
			[:VapourDrawOff,:LiquidDrawOff,:LiquidInlet,:VapourOutlet,:VapourInlet,:LiquidOutlet,:LiquidConnector,:VapourConnector,]
		)
	end
	_base_1::User_Section_ColumnBasic
	VapourDrawOff::vapour_stream
	LiquidDrawOff::liquid_stream
	LiquidInlet::stream
	VapourOutlet::vapour_stream
	VapourInlet::stream
	LiquidOutlet::liquid_stream
	LiquidConnector::stream
	VapourConnector::stream
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export User_Section_Column
function setEquationFlow(in::User_Section_Column)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
	addEquation(15)
	addEquation(16)
	addEquation(17)
	addEquation(18)
	addEquation(19)
	addEquation(20)
	addEquation(21)
	addEquation(22)
	addEquation(23)
	addEquation(24)
	addEquation(25)
	addEquation(26)
	addEquation(27)
	addEquation(28)
	addEquation(29)
	addEquation(30)
end
function atributes(in::User_Section_Column,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/SectionColumn"
	fields[:Brief]="Model of a column section."
	fields[:Info]="== Model of a column section containing ==
* NumberOfTrays trays.

== Specify ==
* the feed stream of each tray (Inlet);
* the Murphree eficiency for each tray Emv;
* the InletL stream of the top tray;
* the InletV stream of the bottom tray.

== Initial Conditions ==
* the trays temperature (OutletL.T);
* the trays liquid level (Level) OR the trays liquid flow (OutletL.F);
* (NoComps - 1) OutletL (OR OutletV) compositions for each tray.
"
	drive!(fields,_)
	return fields
end
User_Section_Column(_::Dict{Symbol,Any})=begin
	newModel=User_Section_Column()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
export DependencyTree
export set_arc, del_arc, has_head, arc_present

# Root node will have id == 0, others nodes id depend from their position in sentence
mutable struct DependencyTree
  root::TreeNode
  nodes::Vector{TreeNode}
  length::Integer

  DependencyTree(root::TreeNode, nodes::Vector{TreeNode}, length::Integer) = new(root, nodes, length)
  function DependencyTree(sentence::Sentence)
    nodes = map(collect(1:sentence.length)) do index
      TreeNode(index, sentence.tokens[index].name)
    end

    root = TreeNode(0)

    new(root, nodes, length(nodes))
  end
end

function set_arc(tree::DependencyTree, word_id::Integer, head_id::Integer, label::String)
  if word_id < 1 || word_id > length(tree.nodes)
    return
  end

  if head_id < 0 || head_id > length(tree.nodes)
    return
  end

  dependent_node = tree.nodes[word_id]
  dependent_node.head_id = head_id
  dependent_node.label = label
end

function del_arc(tree::DependencyTree, word_id::Integer)
  if word_id < 1 || word_id > length(tree.nodes)
    return
  end

  dependent_node = tree.nodes[word_id]
  dependent_node.head_id = EMPTY_NODE
  dependent_node.label = EMPTY_LABEL
end

function has_head(tree::DependencyTree, word_id::Integer)
  if word_id == 0
    return false
  end

  node = tree.nodes[word_id]

  node.head_id != EMPTY_NODE
end

function arc_present(tree::DependencyTree, head_id::Integer, word_id::Integer)
  if 1 <= word_id <= tree.length
    return tree.nodes[word_id].head_id == head_id
  end

  false
end

function get_left_child(tree::DependencyTree, word_id::Integer; child_number = 1)
   (word_id < 0 || word_id > tree.length) && return NONEXIST_TOKEN

   num = 0

   for child_id = 1:(word_id - 1)
    tree.nodes[child_id].head_id != word_id && continue

    num += 1
    num == child_number && return child_id
   end

   NONEXIST_TOKEN
end

function get_right_child(tree::DependencyTree, word_id::Integer; child_number = 1)
  (word_id < 0 || word_id > tree.length) && return NONEXIST_TOKEN

  num = 0

  for child_id = (word_id + 1):tree.length
   tree.nodes[child_id].head_id != word_id && continue

   num += 1
   num == child_number && return child_id
  end

  NONEXIST_TOKEN
end

function convert_to_string(tree::DependencyTree)
  lines = map(tree.nodes) do node
    head_word = if node.head_id == 0
      tree.root.token
    elseif node.head_id > 0
      tree.nodes[node.head_id].token
    else
      NO_HEAD
    end

    "$(node.label)($(head_word)-$(node.head_id), $(node.token)-$(node.word_id))"
  end

  join(lines, "\n")
end

ext: jl

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
function dijkstra(g::AbstractGraph, start::Vertex)
    # clear vertex marks
    g.mark = g.mark * 0

    # init distance vector
    distance = fill(Inf, countvertices(g))
    distance[start] = 0

    for i = 1:countvertices(g)
        v = minimumvertex!(g, distance)
        (distance[v] == Inf) && throw(ErrorException("Unreachable vertices."))
        setmark!(g, v, 1)
        let w = first(g, v)
            while w <= countvertices(g)
                newdistance = distance[v] + getweight(g, v, w)
                if distance[w] > newdistance
                    distance[w] = newdistance
                end
                w = next(g, v, w)
            end
        end
    end

    # clear vertex marks
    g.mark = g.mark * 0
    return distance
end


text:
struct Homography{T₁ <: AbstractMatrix} <: ProjectiveEntity
    𝐇::AbstractMatrix
end

function to_matrix(entity::Homography)
    entity.𝐇
end


text:
module Repo # Octo

using ..DBMS
using ..Backends
using ..AdapterBase
using ..Queryable: Structured, SubQuery, FromItem
using ..Octo: Raw, SQLKeyword
using ..Schema # Schema.validates
using ..Pretty: show

struct NeedsConnectError <: Exception
    msg
end

"""
    Connection
"""
struct Connection
    multiple::Bool
    dbname::String
    loader::Module
    adapter::Module
    conn
end

const ExecuteResult = Union{Nothing, NamedTuple}

@enum RepoLogLevel::Int32 begin
    LogLevelDebugSQL = -1
    LogLevelInfo = 0
end

const current = Dict{Symbol, Union{Nothing, Connection, RepoLogLevel}}(
    :connection => nothing,
    :log_level => LogLevelInfo,
)
current_connection() = current[:connection]

const color_params = (; color=:yellow,)
const color_multiple_database = (; color=:blue, bold=true)

function set_log_level(level::RepoLogLevel)
    current[:log_level] = level
end

"""
    Repo.debug_sql(debug::Bool = true)
"""
function debug_sql(debug::Bool = true)
    set_log_level(debug ? LogLevelDebugSQL : LogLevelInfo)
end

function print_debug_sql_params(io, params)
    print(io, repeat(' ', 4))
    for (idx, x) in enumerate(params)
        printstyled(io, repr(x); color_params...)
        length(params) != idx && printstyled(io, ", ")
    end
end

function print_debug_sql(db::Connection, stmt::Structured, params = nothing)
    if current[:log_level] <= LogLevelDebugSQL
        buf = IOBuffer()
        io = IOContext(buf, :color=>true)
        if db.multiple
            printstyled(io, string('(', db.dbname, ')'); color_multiple_database...)
            print(io, ' ')
        end
        show(io, MIME"text/plain"(), stmt)
        !(params isa Nothing) && print_debug_sql_params(io, params)
        println(io)
        @info String(take!(buf))
    end
end

# Repo.connect
"""
    Repo.connect(; adapter::Module, database::Union{Nothing,Type{D} where {D <: DBMS.AbstractDatabase}}=nothing, multiple::Bool=false, kwargs...)::Connection
"""
function connect(; adapter::Module, database::Union{Nothing,Type{D} where {D <: DBMS.AbstractDatabase}}=nothing, multiple::Bool=false, kwargs...)::Connection
    loader = Backends.backend(adapter)
    conn = Base.invokelatest(loader.db_connect; kwargs...)
    dbname = Base.invokelatest(loader.db_dbname, (; kwargs...))
    connection = Connection(multiple, dbname, loader, adapter, conn)
    if !multiple
        current[:connection] = connection
    end
    connection
end

# Repo.disconnect
"""
    Repo.disconnect(; db::Union{Nothing, Connection}=nothing)
"""
function disconnect(; db::Union{Nothing, Connection}=nothing)
    if db === nothing
        connection = current_connection()
        disconnected = connection.loader.db_disconnect(connection.conn)
        current[:connection] = nothing
        disconnected
    else
        db.loader.db_disconnect(db.conn)
    end
end

function get_primary_key(db, M)::Union{Nothing,Symbol,Vector{Symbol}}
    Tname = Base.typename(M)
    tbl = Schema.tables[Tname]
    Base.get(tbl, :primary_key, nothing)
end

function _field_for_primary_key(db, M) # throw Schema.PrimaryKeyError
    primary_key = get_primary_key(db, M)
    if primary_key === nothing
        throw(Schema.PrimaryKeyError(""))
    else
        a = db.adapter
        table = a.from(M)
        a.Field(table, primary_key)
    end
end

function _field_for_primary_key(db, M, nt::NamedTuple) # throw Schema.PrimaryKeyError
    primary_key = get_primary_key(db, M)
    if primary_key === nothing
        throw(Schema.PrimaryKeyError(""))
    else
        key = _field_for_primary_key(db, M)
        pk = getfield(nt, primary_key)
        (key, pk)
    end
end

# Repo.query

"""
    Repo.query(stmt::Structured; db::Connection=current_connection())
"""
function query(stmt::Structured; db::Connection=current_connection())
    a = db.adapter
    sql = a.to_sql(stmt)
    print_debug_sql(db, stmt)
    db.loader.query(db.conn, sql)
end

"""
    Repo.query(M::Type; db::Connection=current_connection())
"""
function query(M::Type; db::Connection=current_connection())
    a = db.adapter
    table = a.from(M)
    query([a.SELECT * a.FROM table]; db=db)
end

"""
    Repo.query(from::FromItem; db::Connection=current_connection())
"""
function query(from::FromItem; db::Connection=current_connection())
    query(from.__octo_model; db=db)
end

"""
    Repo.query(subquery::SubQuery; db::Connection=current_connection())
"""
function query(subquery::SubQuery; db::Connection=current_connection())
    query(subquery.__octo_query; db=db)
end

"""
    Repo.query(rawquery::Octo.Raw; db::Connection=current_connection())
"""
function query(rawquery::Raw; db::Connection=current_connection())
    query([rawquery]; db=db)
end

"""
    Repo.query(str::AbstractString; db::Connection=current_connection())
"""
function query(str::AbstractString; db::Connection=current_connection())
    query(Raw(str); db=db)
end

### Repo.query - vals::Vector
"""
    Repo.query(stmt::Structured, vals::Vector; db::Connection=current_connection())
"""
function query(stmt::Structured, vals::Vector; db::Connection=current_connection()) # throw Backends.UnsupportedError
    a = db.adapter
    prepared = a.to_sql(stmt)
    print_debug_sql(db, stmt, vals)
    db.loader.query(db.conn, prepared, vals)
end

### Repo.query - pk
"""
    Repo.query(M::Type, pk::Union{Int, String}; db::Connection=current_connection())
"""
function query(M::Type, pk::Union{Int, String}; db::Connection=current_connection()) # throw Schema.PrimaryKeyError
    a = db.adapter
    table = a.from(M)
    key = _field_for_primary_key(db, M)
    query([a.SELECT * a.FROM table a.WHERE key == pk]; db=db)
end

"""
    Repo.query(from::FromItem, pk::Union{Int, String}; db::Connection=current_connection())
"""
function query(from::FromItem, pk::Union{Int, String}; db::Connection=current_connection()) # throw Schema.PrimaryKeyError
    query(from.__octo_model, pk; db=db)
end

### Repo.query - pk_range
"""
    Repo.query(M::Type, pk_range::UnitRange{Int64}; db::Connection=current_connection())
"""
function query(M::Type, pk_range::UnitRange{Int64}; db::Connection=current_connection()) # throw Schema.PrimaryKeyError
    a = db.adapter
    table = a.from(M)
    key = _field_for_primary_key(db, M)
    query([a.SELECT * a.FROM table a.WHERE key a.BETWEEN pk_range.start a.AND pk_range.stop]; db=db)
end

"""
    Repo.query(from::FromItem, pk_range::UnitRange{Int64}; db::Connection=current_connection())
"""
function query(from::FromItem, pk_range::UnitRange{Int64}; db::Connection=current_connection()) # throw Schema.PrimaryKeyError
    query(from.__octo_model, pk_range; db=db)
end

### Repo.query - nt::NamedTuple
"""
    Repo.query(stmt::Structured, nt::NamedTuple; db::Connection=current_connection())
"""
function query(stmt::Structured, nt::NamedTuple; db::Connection=current_connection())
    a = db.adapter
    v = vec(stmt)
    for (idx, kv) in enumerate(pairs(nt))
        key = a.Field(nothing, kv.first)
        push!(v, key == kv.second)
        idx != length(nt) && push!(v, a.AND)
    end
    query(v; db=db)
end

"""
    Repo.query(M::Type, nt::NamedTuple; db::Connection=current_connection())
"""
function query(M::Type, nt::NamedTuple; db::Connection=current_connection())
    a = db.adapter
    table = a.from(M)
    v = [a.SELECT, *, a.FROM, table, a.WHERE]
    for (idx, kv) in enumerate(pairs(nt))
        key = Base.getproperty(table, kv.first)
        push!(v, key == kv.second)
        idx != length(nt) && push!(v, a.AND)
    end
    query(v; db=db)
end

"""
    Repo.query(from::FromItem, nt::NamedTuple; db::Connection=current_connection())
"""
function query(from::FromItem, nt::NamedTuple; db::Connection=current_connection())
    query(from.__octo_model, nt; db=db)
end

"""
    Repo.query(subquery::SubQuery, nt::NamedTuple; db::Connection=current_connection())
"""
function query(subquery::SubQuery, nt::NamedTuple; db::Connection=current_connection())
    query(subquery.__octo_query, nt; db=db)
end

"""
    Repo.query(rawquery::Octo.Raw, nt::NamedTuple; db::Connection=current_connection())
"""
function query(rawquery::Raw, nt::NamedTuple; db::Connection=current_connection())
    query([rawquery], nt; db=db)
end

# Repo.get
"""
    Repo.get(M::Type, pk::Union{Int, String}; db::Connection=current_connection())
"""
function get(M::Type, pk::Union{Int, String}; db::Connection=current_connection()) # throw Schema.PrimaryKeyError
    query(M, pk; db=db)
end

"""
    Repo.get(M::Type, pk_range::UnitRange{Int64}; db::Connection=current_connection())
"""
function get(M::Type, pk_range::UnitRange{Int64}; db::Connection=current_connection()) # throw Schema.PrimaryKeyError
    query(M, pk_range; db=db)
end

"""
    Repo.get(M::Type, nt::NamedTuple; db::Connection=current_connection())
"""
function get(M::Type, nt::NamedTuple; db::Connection=current_connection())
    query(M, nt; db=db)
end

# Repo.execute
"""
    Repo.execute(stmt::Structured; db::Connection=current_connection())
"""
function execute(stmt::Structured; db::Connection=current_connection())
    a = db.adapter
    sql = a.to_sql(stmt)
    print_debug_sql(db, stmt)
    db.loader.execute(db.conn, sql)
end

"""
    Repo.execute(stmt::Structured, vals::Vector; db::Connection=current_connection())
"""
function execute(stmt::Structured, vals::Vector; db::Connection=current_connection())
    a = db.adapter
    prepared = a.to_sql(stmt)
    print_debug_sql(db, stmt, vals)
    db.loader.execute(db.conn, prepared, vals)
end

"""
    Repo.execute(stmt::Structured, nts::Vector{<:NamedTuple}; db::Connection=current_connection())
"""
function execute(stmt::Structured, nts::Vector{<:NamedTuple}; db::Connection=current_connection())
    a = db.adapter
    prepared = a.to_sql(stmt)
    print_debug_sql(db, stmt, nts)
    db.loader.execute(db.conn, prepared, nts)
end

execute(raw::Raw; db::Connection=current_connection())                            = execute([raw]; db=db)
execute(raw::Raw, nt::NamedTuple; db::Connection=current_connection())            = execute([raw], [nt]; db=db)
execute(raw::Raw, nts::Vector{<:NamedTuple}; db::Connection=current_connection()) = execute([raw], nts; db=db)
execute(raw::Raw, vals::Vector; db::Connection=current_connection())              = execute([raw], vals; db=db)

execute(str::AbstractString; db::Connection=current_connection())                 = execute(Raw(str); db=db)

# Repo.insert!

function do_insert(block, a, returning::Union{Nothing,Symbol,Vector}, db::Connection)
    extra = []
    if a.DatabaseID === DBMS.PostgreSQL && returning !== nothing
        extra = vcat(a.RETURNING, returning isa Symbol ? returning : tuple(Symbol.(returning)...))
    end
    result = block(extra)
    if a.DatabaseID === DBMS.PostgreSQL
        result
    else
        nt = execute_result(a.INSERT; db=db)
        if result === nothing
            nt
        else
            merge(nt, result)
        end
    end
end

"""
    Repo.insert!(M::Type, nts::Vector{<:NamedTuple}; db::Connection=current_connection(), returning::Union{Nothing,Symbol,Vector}=get_primary_key(db, M))
"""
function insert!(M::Type, nts::Vector{<:NamedTuple}; db::Connection=current_connection(), returning::Union{Nothing,Symbol,Vector}=get_primary_key(db, M))
    if !isempty(nts)
        Schema.validates.(M, nts) # throw InvalidChangesetError
        a = db.adapter
        table = a.from(M)
        nt = first(nts)
        fieldnames = a.Enclosed(collect(keys(nt)))
        values = a.placeholders(length(nt))
        do_insert(a, returning, db) do extra
            execute([a.INSERT a.INTO table fieldnames a.VALUES values extra...], nts; db=db)
        end
    else
        nothing
    end
end

"""
    Repo.insert!(M::Type, vals::Vector{<:Tuple}; db::Connection=current_connection(), returning::Union{Nothing,Symbol,Vector}=get_primary_key(db, M))
"""
function insert!(M::Type, vals::Vector{<:Tuple}; db::Connection=current_connection(), returning::Union{Nothing,Symbol,Vector}=get_primary_key(db, M))
    if !isempty(vals)
        a = db.adapter
        table = a.from(M)
        do_insert(a, returning, db) do extra
            execute([a.INSERT a.INTO table a.VALUES a.VectorOfTuples(vals) extra...]; db=db)
        end
    else
        nothing
    end
end

"""
    Repo.insert!(M, nt::NamedTuple; kwargs...)
"""
function insert!(M, nt::NamedTuple; kwargs...)
    insert!(M, [nt]; kwargs...)
end

# Repo.execute_result
function execute_result(command::SQLKeyword; db::Connection=current_connection())::NamedTuple
    db.loader.execute_result(db.conn, command)
end

# Repo.update!
"""
    Repo.update!(M::Type, nt::NamedTuple; db::Connection=current_connection())
"""
function update!(M::Type, nt::NamedTuple; db::Connection=current_connection()) # throw Schema.PrimaryKeyError
    Schema.validates(M, nt) # throw InvalidChangesetError
    a = db.adapter
    (key, pk) = _field_for_primary_key(db, M, nt)
    table = a.from(M)
    rest = filter(kv -> kv.first != key.name, pairs(nt))
    v = Any[a.UPDATE, table, a.SET]

    # tup[1] : idx
    # tup[2] : kv
    push!(v, tuple(map(tup -> a.Field(table, tup[2].first) == a.placeholder(tup[1]), enumerate(rest))...))

    push!(v, a.WHERE)
    push!(v, key == pk)
    execute(v, collect(values(rest)); db=db)
end

function vecjoin(elements::Array{E, N}, delim::D)::Vector{Union{E, D}} where  {E, N, D}
    first = true
    result = Vector{Union{E, D}}()
    for el in elements
        if first
            first = false
        else
            push!(result, delim)
        end
        push!(result, el)
    end
    result
end

# Repo.delete!

"""
    Repo.delete!(M::Type, nt::NamedTuple; db::Connection=current_connection())
"""
function delete!(M::Type, nt::NamedTuple; db::Connection=current_connection())
    a = db.adapter
    table = a.from(M)
    execute(hcat([a.DELETE a.FROM table a.WHERE], vecjoin([a.Field(table, k) == v for (k, v) in pairs(nt)], a.AND)...); db=db)
end

"""
    Repo.delete!(M::Type, pk_range::UnitRange{Int64}; db::Connection=current_connection())
"""
function delete!(M::Type, pk_range::UnitRange{Int64}; db::Connection=current_connection()) # throw Schema.PrimaryKeyError
    a = db.adapter
    table = a.from(M)
    key = _field_for_primary_key(db, M)
    execute([a.DELETE a.FROM table a.WHERE key a.BETWEEN pk_range.start a.AND pk_range.stop]; db=db)
end

function sql_startswith_insert_update_delete(sql::String)::Bool
    lowercase(first(split(lstrip(sql), ' '))) in ("insert", "update", "delete")
end

end # module Octo.Repo


text:
using Qsosed
using Test, QuadGK

@testset "Disk radiation" begin
    @testset "Blackbody radiation" begin
        @testset "Test spectral radiance" begin
            bb = BlackBody(1e5)
            @test spectral_radiance_frequency(bb, 100) ≈ 3.07235837e-28 rtol = 1e-3
            @test spectral_radiance_frequency(bb, 1e15) ≈ 0.02393854 rtol = 1e-4
            @test spectral_radiance_frequency(bb, 1e18) ≈ 0.0 rtol = 1e-6 atol = 1e-15
            bb = BlackBody(5e3)
            @test spectral_radiance_frequency(bb, 100) ≈ 1.53617919e-29 rtol = 1e-4
            @test spectral_radiance_frequency(bb, 1e15) ≈ 1.00024067e-06
            @test spectral_radiance_frequency(bb, 1e18) ≈ 0
        end
        @testset "Test integrating energy band" begin
            bb = BlackBody(1e4)
            @test spectral_band_radiance_frequency(bb, 100e12, 1000e12) ≈ 1.30413e11 rtol =
                1e-2
            @test spectral_band_radiance(bb, 0.00041357, 0.0041357) ≈ 1.30413e11 rtol = 1e-2
        end
        @testset "Integrating spectral radiance recovers SB law" begin
            bb = BlackBody(1e4)
            @test spectral_band_radiance(bb, 1e-5, 1e2) ≈ SIGMA_SB * bb.T^4 / π rtol = 1e-6
        end
    end
    @testset "Disk radiation" begin
        @testset "Test relativistic AD correction" begin
            bh = BlackHole(1e8, 0.5, 0.6)
            @test disk_nt_rel_factors(bh, 1) ≈ 0.0
            @test disk_nt_rel_factors(bh, 10) ≈ 0.2615157665830733
            @test disk_nt_rel_factors(bh, 50) ≈ 0.6086578218842061
            @test disk_nt_rel_factors(bh, 200) ≈ 0.7849369908504347
            bh = BlackHole(1e9, 0.5, 0.99)
            @test disk_nt_rel_factors(bh, 10) ≈ 0.49104425779640487
            @test disk_nt_rel_factors(bh, 50) ≈ 0.7072886586539658
            @test disk_nt_rel_factors(bh, 200) ≈ 0.8339124564014411
        end
        @testset "Test disc radiated power" begin
            bh = BlackHole(1e8 * M_SUN, 0.5, 0)
            @test disk_flux(bh, 10) / 1.17 ≈ 6839115834162428.0 rtol = 1e-6
            @test disk_flux(bh, 100) / 1.17 ≈ 3.924336168219069e13 rtol = 1e-6
            bh = BlackHole(1e9 * M_SUN, 0.5, 0.9)
            @test disk_flux(bh, 10) / 1.17 ≈ 889138044477300.2 rtol = 1e-6
            @test disk_flux(bh, 100) / 1.17 ≈ 1662600341587.0435 rtol = 1e-6
        end

        @testset "Integrate spectral radial radiance" begin
            bh = BlackHole(1e8 * M_SUN, 0.5, 0)
            energy_range = 10 .^ range(-3, -1, length = 500)
            radius = 50
            flux = compute_disk_flux_at_radius(bh, radius)
            temp = compute_disk_temperature(bh, radius)
            exp = SIGMA_SB * temp^4
            @test flux ≈ exp rtol = 1e-3
            integ_flux, _ = quadgk(
                e -> compute_disk_spectral_flux_at_radius(bh, radius, e),
                1e-5,
                1e2,
                atol = 0,
                rtol = 1e-4,
            )
            @test integ_flux ≈ exp rtol = 1e-2

            @testset "Recover luminosity" begin
                function kernel(bh, radius)
                    integ_flux, _ = quadgk(
                        e -> compute_disk_spectral_flux_at_radius(bh, radius, e),
                        1e-5,
                        1e2,
                        atol = 0,
                        rtol = 1e-4,
                    )
                    return integ_flux
                end
                lumin_integ, _ = quadgk(
                    r -> 2 * 2 * π * r * kernel(bh, r),
                    6.0,
                    1500.0,
                    atol = 0,
                    rtol = 1e-4,
                )
                lumin_integ *= bh.Rg^2
                lbol = compute_bolometric_luminosity(bh)
                @test lumin_integ ≈ lbol rtol = 1e-1
            end

            @testset "Disk sed" begin
                disk_sed = compute_disk_spectral_luminosity(bh, energy_range)
                integ_lumin = sum(disk_sed[2:end] .* diff(energy_range))
                @test integ_lumin ≈ compute_bolometric_luminosity(bh) rtol = 1e-1
            end

        end
    end
end

bh = BlackHole(1e8 * M_SUN, 0.5, 0.0)
@test gravity_radius(bh) ≈ 1580 rtol = 0.01


text:
###################################################################
#Description
###################################################################

###################################################################
##Packages
###################################################################
using LinearAlgebra, SpecialFunctions, NLsolve

###################################################################
##Parameters
###################################################################
σ = 3 #Elasticity of substitution, Source : Footnote 45, EK2002
θ = 8.28 #Technology dispersion parameter(Comparative Advantage), Source : Table VI, EK2002
β = 0.21 #Labor share in costs, Source : Table VIII, EK2002

###################################################################
#Generate data
#Wages, Trade costs, Technology shift parameter
###################################################################
N=43 #Number of countries

#Generate wages costs
ϵ_temp = randn(N)
w = 100*fill(1.0,N,1) + 100*ϵ_temp.^2

#Generate trade costs
ϵ_temp = randn(N,N)
D = fill(1.0,N,N) + ϵ_temp.^2
D=D-Diagonal(D)+I # Make diagonal entries equal 1 (No internal resistance)
#D = fill(1.0,N,N)

#Generate productivity shift parameter
ϵ_temp = randn(N)
T = 10*fill(1.0,N,1) + ϵ_temp.^2

###################################################################
#Model
###################################################################

############### Step 1 : Fixed point for prices ###################
function P_fixedpoint(P_telde,D,T,w,θ,β,σ,N)
    γ = (gamma((θ+1-σ)/θ))^(1/(1-σ)) #Gamma function
    Part1 = γ^(-θ)*T.*(w.^(-θ*β)) #create vector with product of gamma, wages^(θβ) and shift parameter
    Part2 = repeat(transpose(Part1),N) #create matrix which contain copies of above vector in each row
    G = (D.^(-θ)).*Part2 # point wise multiplication of distance^(-θ) and above matrix
    P_tilde = G*P_telde.^(1-β)
    return P_tilde
end
f(P_tilde) =  P_fixedpoint(P_tilde,D,T,w,θ,β,σ,N) #redefine function with single argument

init=ones(N,1) #initial value
sol = fixedpoint(f, init)
P = sol.zero.^(-1/θ)

############### Step 2 : Trade Shares ###################

γ = (gamma((θ+1-σ)/θ))^(1/(1-σ)) #Gamma function
Part1 = γ^(-θ)*T.*(w.^(-θ*β)).*(P.^(-θ*(1-β))) #create vector with product of gamma, wages^(θβ), P^(-θ(1-β)) and shift parameter
Part2 = repeat(transpose(Part1),N) #create matrix which contain copies of above vector in each row
Part3=transpose(repeat(transpose(P.^(-θ)),N)) #create matrix which contain copies of price vector in each row
Π = (D.^(-θ)).*Part2./Part3 # point wise multiplication of distance^(-θ) and above matrixs

#Correlation check
corr = cor(T,diag(Π))
println(corr)

############### Step 3 : Counterfactual ###################
# RTA between two countries : reduce distance between single pair and look at change in trade share for all countries

############### Step 4 : Graphics ###################
#Trade flows in circle


text:
"""
    NormalCanon(η, λ)

Canonical parametrisation of the Normal distribution with canonical parameters `η` and `λ`.

The two *canonical parameters* of a normal distribution ``\\mathcal{N}(\\mu, \\sigma^2)`` with mean ``\\mu`` and
standard deviation ``\\sigma`` are ``\\eta = \\sigma^{-2} \\mu`` and ``\\lambda = \\sigma^{-2}``.
"""
struct NormalCanon{T<:Real} <: ContinuousUnivariateDistribution
    η::T       # σ^(-2) * μ
    λ::T       # σ^(-2)
    μ::T       # μ

    function NormalCanon{T}(η, λ; check_args::Bool=true) where T
        check_args && @check_args(NormalCanon, λ > zero(λ))
        new{T}(η, λ, η / λ)
    end
end

NormalCanon(η::T, λ::T; check_args::Bool=true) where {T<:Real} = NormalCanon{typeof(η/λ)}(η, λ; check_args=check_args)
NormalCanon(η::Real, λ::Real; check_args::Bool=true) = NormalCanon(promote(η, λ)...; check_args=check_args)
NormalCanon(η::Integer, λ::Integer; check_args::Bool=true) = NormalCanon(float(η), float(λ); check_args=check_args)
NormalCanon() = NormalCanon{Float64}(0.0, 1.0; check_args=false)

@distr_support NormalCanon -Inf Inf

#### Type Conversions
convert(::Type{NormalCanon{T}}, η::S, λ::S) where {T <: Real, S <: Real} = NormalCanon(T(η), T(λ))
convert(::Type{NormalCanon{T}}, d::NormalCanon{S}) where {T <: Real, S <: Real} = NormalCanon(T(d.η), T(d.λ); check_args=false)

## conversion between Normal and NormalCanon

convert(::Type{Normal}, d::NormalCanon) = Normal(d.μ, 1 / sqrt(d.λ))
convert(::Type{NormalCanon}, d::Normal) = (λ = 1 / d.σ^2; NormalCanon(λ * d.μ, λ))
meanform(d::NormalCanon) = convert(Normal, d)
canonform(d::Normal) = convert(NormalCanon, d)

#### Parameters

params(d::NormalCanon) = (d.η, d.λ)
@inline partype(d::NormalCanon{T}) where {T<:Real} = T

#### Statistics

mean(d::NormalCanon) = d.μ
median(d::NormalCanon) = mean(d)
mode(d::NormalCanon) = mean(d)

skewness(d::NormalCanon{T}) where {T<:Real} = zero(T)
kurtosis(d::NormalCanon{T}) where {T<:Real} = zero(T)

var(d::NormalCanon) = 1 / d.λ
std(d::NormalCanon) = sqrt(var(d))

entropy(d::NormalCanon) = (-log(d.λ) + log2π + 1) / 2

location(d::NormalCanon) = mean(d)
scale(d::NormalCanon) = std(d)

#### Evaluation

pdf(d::NormalCanon, x::Real) = (sqrt(d.λ) / sqrt2π) * exp(-d.λ * abs2(x - d.μ)/2)
logpdf(d::NormalCanon, x::Real) = (log(d.λ) - log2π - d.λ * abs2(x - d.μ))/2

zval(d::NormalCanon, x::Real) = (x - d.μ) * sqrt(d.λ)
xval(d::NormalCanon, z::Real) = d.μ + z / sqrt(d.λ)

cdf(d::NormalCanon, x::Real) = normcdf(zval(d,x))
ccdf(d::NormalCanon, x::Real) = normccdf(zval(d,x))
logcdf(d::NormalCanon, x::Real) = normlogcdf(zval(d,x))
logccdf(d::NormalCanon, x::Real) = normlogccdf(zval(d,x))

quantile(d::NormalCanon, p::Real) = xval(d, norminvcdf(p))
cquantile(d::NormalCanon, p::Real) = xval(d, norminvccdf(p))
invlogcdf(d::NormalCanon, lp::Real) = xval(d, norminvlogcdf(lp))
invlogccdf(d::NormalCanon, lp::Real) = xval(d, norminvlogccdf(lp))

#### Sampling

rand(rng::AbstractRNG, cf::NormalCanon) = cf.μ + randn(rng) / sqrt(cf.λ)

#### Affine transformations

function Base.:+(d::NormalCanon, c::Real)
    η, λ = params(d)
    return NormalCanon(η + c * λ, λ)
end
Base.:*(c::Real, d::NormalCanon) = NormalCanon(d.η / c, d.λ / c^2)


text:
function LToggle(parent::Scene; bbox = nothing, kwargs...)

    default_attrs = default_attributes(LToggle, parent).attributes
    theme_attrs = subtheme(parent, :LToggle)
    attrs = merge!(merge!(Attributes(kwargs), theme_attrs), default_attrs)

    @extract attrs (halign, valign, cornersegments, framecolor_inactive,
        framecolor_active, buttoncolor, active, toggleduration, rimfraction)

    decorations = Dict{Symbol, Any}()

    layoutobservables = LayoutObservables{LToggle}(attrs.width, attrs.height, attrs.tellwidth, attrs.tellheight,
        halign, valign, attrs.alignmode; suggestedbbox = bbox)

    markersize = lift(layoutobservables.computedbbox) do bbox
        min(width(bbox), height(bbox))
    end

    button_endpoint_inactive = lift(markersize) do ms
        bbox = layoutobservables.computedbbox[]
        Point2f0(left(bbox) + ms / 2, bottom(bbox) + ms / 2)
    end

    button_endpoint_active = lift(markersize) do ms
        bbox = layoutobservables.computedbbox[]
        Point2f0(right(bbox) - ms / 2, bottom(bbox) + ms / 2)
    end

    buttonvertices = lift(markersize, cornersegments) do ms, cs
        roundedrectvertices(layoutobservables.computedbbox[], ms * 0.499, cs)
    end

    # trigger bbox
    layoutobservables.suggestedbbox[] = layoutobservables.suggestedbbox[]

    framecolor = Node{Any}(active[] ? framecolor_active[] : framecolor_inactive[])
    frame = poly!(parent, buttonvertices, color = framecolor, raw = true)[end]
    decorations[:frame] = frame

    animating = Node(false)
    buttonpos = Node(active[] ? [button_endpoint_active[]] : [button_endpoint_inactive[]])

    # make the button stay in the correct place (and start there)
    on(layoutobservables.computedbbox) do bbox
        if !animating[]
            buttonpos[] = active[] ? [button_endpoint_active[]] : [button_endpoint_inactive[]]
        end
    end

    buttonsize = lift(markersize, rimfraction) do ms, rf
        ms * (1 - rf)
    end

    button = scatter!(parent, buttonpos, markersize = buttonsize, color = buttoncolor, raw = true)[end]
    decorations[:button] = button

    buttonstate = addmousestate!(parent, button, frame)

    onmouseleftclick(buttonstate) do state
        if animating[]
            return
        end
        animating[] = true

        tstart = time()

        anim_posfrac = Animations.Animation(
            [0, toggleduration[]],
            active[] ? [1.0, 0.0] : [0.0, 1.0],
            Animations.sineio())
        coloranim = Animations.Animation(
            [0, toggleduration[]],
            active[] ? [framecolor_active[], framecolor_inactive[]] : [framecolor_inactive[], framecolor_active[]],
            Animations.sineio())

        active[] = !active[]
        @async while true
            t = time() - tstart
            # request endpoint values in every frame if the layout changes during
            # the animation
            buttonpos[] = [Animations.linear_interpolate(anim_posfrac(t),
                button_endpoint_inactive[], button_endpoint_active[])]
            framecolor[] = coloranim(t)
            if t >= toggleduration[]
                animating[] = false
                break
            end
            sleep(1/FPS[])
        end
    end

    LToggle(parent, layoutobservables, attrs, decorations)
end


text:
import Base: -, +, *, /, \, ∘, zero, conj, real, imag, abs, abs2
import LinearAlgebra: transpose!, transpose

export tensorproduct!

# Identity operator

struct Identity end

(*)(::Identity,s::GridData) = s

### On scalar grid data ####

function (-)(p_in::GridData)
  Base.broadcast(-,p_in)
end

for op in (:+, :-, :*)
    @eval function $op(p1::T,p2::T) where {T <: GridData}
       Base.broadcast($op,p1,p2)
    end
    @eval function $op(p1::Number,p2::GridData)
       Base.broadcast($op,p1,p2)
    end
    @eval function $op(p1::GridData,p2::Number)
       Base.broadcast($op,p1,p2)
    end
end

for op in (:/,)
    @eval function $op(p1::GridData,p2::Number)
       Base.broadcast($op,p1,p2)
    end
end

"""
    product!(out::GridData,p::GridData,q::GridData)

Compute the Hadamard (i.e. element by element) product of grid data
`p` and `q` (of the same type) and return the result in `out`.

# Example

```jldoctest
julia> q = Edges(Dual,(8,6));

julia> out = p = deepcopy(q);

julia> q.u[3,2] = 0.3;

julia> p.u[3,2] = 0.2;

julia> product!(out,p,q)
Edges{Dual,8,6,Float64} data
u (in grid orientation)
6×7 Array{Float64,2}:
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.06  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
v (in grid orientation)
5×8 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
@inline product!(out::T, p::T, q::T) where {T <: GridData} = (out .= p.*q)

for f in (:Nodes, :XEdges, :YEdges, :Edges)
  @eval @inline product!(out::$f{C,NX,NY},p::$f{C,NX,NY},q::$f{C,NX,NY}) where {C<:CellType,NX,NY} = (out .= p.*q)
end

@inline product!(out::EdgeGradient{C,D,NX,NY},p::EdgeGradient{C,D,NX,NY},q::EdgeGradient{C,D,NX,NY}) where {C<:CellType,D<:CellType,NX,NY} = (out .= p.*q)

"""
    product(p::Edges/Nodes,q::Edges/Nodes) --> Edges/Nodes

Compute the Hadamard product of edge or nodal (primal or dual) data `p` and `q` and return
the result. This operation can also be carried out with the `∘` operator:

# Example

```jldoctest
julia> q = Edges(Dual,(8,6));

julia> p = deepcopy(q);

julia> q.u[3,2] = 0.3;

julia> p.u[3,2] = 0.2;

julia> p∘q
Edges{Dual,8,6,Float64} data
u (in grid orientation)
6×7 Array{Float64,2}:
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
 0.0  0.0  0.06  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0
v (in grid orientation)
5×8 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
@inline product(p::T, q::T) where {T <: GridData} = product!(T(), p, q)

@inline function (∘)(p::GridData{T,NX,NY}, q::GridData) where {T,NX,NY}
  product!(typeof(p)(), p, q)
end

"""
    tensorproduct!(q::EdgeGradient,u::Edges,v::Edges,ut::EdgeGradient,vt::EdgeGradient)

In-place tensor product of `u` and `v`, with the result returned in `q`. The
`ut` and `vt` are supplied as temporary storage.
"""
function tensorproduct!(q::EdgeGradient{C,D,NX,NY},u::Edges{C,NX,NY},v::Edges{C,NX,NY},
                        ut::EdgeGradient{C,D,NX,NY},vt::EdgeGradient{C,D,NX,NY}) where {C<:CellType,D<:CellType,NX,NY}
    grid_interpolate!(ut,u)
    grid_interpolate!(vt,v)
    q.dudx .= ut.dudx ∘ vt.dudx
    q.dudy .= ut.dvdx ∘ vt.dudy
    q.dvdx .= ut.dudy ∘ vt.dvdx
    q.dvdy .= ut.dvdy ∘ vt.dvdy
    return q
end

tensorproduct!(q::EdgeGradient{C,D,NX,NY},u::Edges{C,NX,NY},v::Edges{C,NX,NY}) where {C<:CellType,D<:CellType,NX,NY} =
        tensorproduct!(q,u,v,zero(q),zero(q))

(*)(u::Edges{C,NX,NY},v::Edges{C,NX,NY}) where {C<:CellType,NX,NY} = tensorproduct!(EdgeGradient(C,u),u,v)

"""
    transpose!(pt::EdgeGradient,p::EdgeGradient)

In-place element-by-element transpose of `EdgeGradient` `p`.
"""
function transpose!(qt::EdgeGradient{C,D,NX,NY},q::EdgeGradient{C,D,NX,NY}) where {C<:CellType,D<:CellType,NX,NY}
    qt.dudx .= q.dudx
    qt.dudy .= q.dvdx
    qt.dvdx .= q.dudy
    qt.dvdy .= q.dvdy
    return qt
end

"""
    transpose(p::EdgeGradient) -> EdgeGradient

Element-by-element transpose of `EdgeGradient` `p`
"""
transpose(q::EdgeGradient{C,D,NX,NY}) where {C<:CellType,D<:CellType,NX,NY} =
    transpose!(zero(q),q)

"""
    zero(p::GridData)

Return data of the same type as `p`, filled with zeros.
"""
zero(::Type{T}) where {T <: GridData} = T()
zero(::T) where {T <: GridData} = T()

#### ON COMPLEX GRID DATA

for f in (:conj,)
    @eval function $f(A::GridData{NX,NY,T}) where {NX,NY,T <: ComplexF64}
        Acopy = similar(A)
        Acopy.data .= broadcast($f,A.data)
        return Acopy
    end
end

for f in (:real, :imag, :abs, :abs2)
  @eval function $f(A::GridData{NX,NY,T}) where {NX,NY,T <: ComplexF64}
      Acopy = similar(A,element_type=Float64)
      Acopy.data .= broadcast($f,A.data)
      return Acopy
  end
end

## BROADCASTING

# This enables fused in-place broadcasting of all grid data
Base.BroadcastStyle(::Type{<:GridData}) = Broadcast.ArrayStyle{GridData}()

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GridData}},::Type{T}) where {T}
    similar(unpack(bc),element_type=T)
end

function Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GridData}})
    similar(unpack(bc))
end

function Base.copyto!(dest::T,bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{GridData}}) where {T <: GridData}
    copyto!(dest.data,unpack_data(bc, nothing))
    dest
end

max_stars_repo_name: JuliaIBPM/CartesianGrids

max_issues_repo_name: JuliaIBPM/CartesianGrids

max_forks_repo_name: JuliaIBPM/CartesianGrids


max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
function setup_hdf5dataset_test()
    datasets = [
        ("d1", rand(2, 10)),
        ("g1/d1", rand(10)),
        ("g1/d2", string.('a':'j')),
        ("g2/g1/d1", "test")
    ]

    fid = h5open("test.h5", "w")
    for (path, data) in datasets
        fid[path] = data
    end
    close(fid)

    return first.(datasets), last.(datasets)
end
cleanup_hdf5dataset_test() = rm("test.h5")

@testset "HDF5Dataset" begin
    paths, datas = setup_hdf5dataset_test()
    @testset "Single path" begin
        dataset = HDF5Dataset("test.h5", "d1")
        for i in 1:10
            @test getobs(dataset, i) == getobs(datas[1], i)
        end
        @test numobs(dataset) == 10
        close(dataset)
    end
    @testset "Multiple paths" begin
        dataset = HDF5Dataset("test.h5", paths)
        for i in 1:10
            data = Tuple(map(x -> (x isa String) ? x : getobs(x, i), datas))
            @test @inferred(Tuple, getobs(dataset, i)) == data
        end
        @test numobs(dataset) == 10
        close(dataset)
    end
    cleanup_hdf5dataset_test()
end


text:
function forward_step(body::Body, dt)
	acc = accel(body)
	body.pos += body.vel * dt
	body.vel += acc * dt
end

function leapfrog_step(body::Body, dt)
	body.vel += accel(body) * 0.5 * dt
	body.pos += body.vel * dt
	body.vel += accel(body) * 0.5 * dt
end

function rk2_step(body::Body, dt)
	old_pos = body.pos
	half_vel = body.vel + accel(body)*0.5*dt
	body.pos += body.vel * 0.5 * dt
	body.vel += accel(body) * dt
	body.pos = old_pos + half_vel*dt
end

function rk4_step(body::Body, dt)
	old_pos = body.pos
	a0 = accel(body)
	body.pos = old_pos + body.vel*0.5*dt + a0*0.125*dt*dt
	a1 = accel(body)
	body.pos = old_pos + body.vel*dt + a1*0.5*dt*dt
	a2 = accel(body)
	body.pos = old_pos + body.vel*dt + (a0 + a1*2)*(1/6.0)*dt*dt
	body.vel += (a0 + a1*4 + a2) * (1/6.0) * dt
end

function yo6_step(body::Body, dt)
	d = [0.784513610477560e0, 0.235573213359357e0,
				-1.17767998417887e0, 1.31518632068391e0 ]

	for i in 1:3; leapfrog_step(body, dt*d[i]); end;
	leapfrog_step(body, dt*d[4])
	for i in 0:2; leapfrog_step(body, dt*d[3-i]); end;
end

function yo8_step(body::Body, dt)
	d = [0.104242620869991e1, 0.182020630970714e1, 0.157739928123617e0,
				0.244002732616735e1, -0.716989419708120e-2, -0.244699182370524e1,
				-0.161582374150097e1, -0.17808286265894516e1]

	for i in 1:7; leapfrog_step(body, dt*d[i]); end;
	leapfrog_step(body, dt*d[8])
	for i in 0:6; leapfrog_step(body, dt*d[7-i]); end;
end

######## Multistep ########

function ms2_step(body::Body, dt)
	if body.step == 0
		set_prev_accel!(body, 1, accel(body))
		rk2_step(body, dt)
	else
		old_acc = accel(body)
		jdt = old_acc - prev_accel(body, 1)
		body.pos += body.vel*dt + 0.5*old_acc*dt*dt
		body.vel += old_acc*dt + 0.5*jdt*dt
		set_prev_accel!(body, 1, old_acc)
	end
end

function ms4_step(body::Body, dt)
	if body.step <= 2
		body.prev_accel[:,3-body.step] = accel(body)
		rk4_step(body, dt)
	else
		ap0 = accel(body)
		jdt = ap0*(11.0/6.0) - 3*body.prev_accel[:,1] + 1.5*body.prev_accel[:,2] - body.prev_accel[:,3]/3
		sdt2 = 2*ap0 - 5*body.prev_accel[:,1] + 4*body.prev_accel[:,2] - body.prev_accel[:,3]
		cdt3 = ap0 - 3*body.prev_accel[:,1] + 3*body.prev_accel[:,2] - body.prev_accel[:,3]
		body.pos += body.vel*dt + (ap0/2.0 + jdt/6.0 + sdt2/24.0)*dt*dt
		body.vel += ap0*dt + (jdt/2.0 + sdt2/6.0 + cdt3/24.0)*dt
		# Update prev acc
		body.prev_accel[:,3] = body.prev_accel[:,2]
		body.prev_accel[:,2] = body.prev_accel[:,1]
		body.prev_accel[:,1] = ap0
	end
end

function ms4pc_step(body::Body, dt)
	if body.step == 0
		set_prev_accel!(body, 4, accel(body))
		rk4_step(body, dt)
	elseif body.step == 1
		set_prev_accel!(body, 3	, accel(body))
		rk4_step(body, dt)
	elseif body.step == 2
		set_prev_accel!(body, 2, accel(body))
		rk4_step(body, dt)
		set_prev_accel!(body, 1, accel(body))
	else
		old_pos = body.pos
		old_vel = body.vel
		# ap0 = accel(body)
		jdt = prev_accel(body, 1)*(11.0/6.0) - 3*prev_accel(body, 2) + 1.5*prev_accel(body, 3) - prev_accel(body, 4)/3
		sdt2 = 2*prev_accel(body, 1) - 5*prev_accel(body, 2) + 4*prev_accel(body, 3) - prev_accel(body, 4)
		cdt3 = prev_accel(body, 1) - 3*prev_accel(body, 2) + 3*prev_accel(body, 3) - prev_accel(body, 4)
		body.pos += body.vel*dt + (prev_accel(body, 1)/2.0 + jdt/6.0 + sdt2/24.0)*dt*dt

		set_prev_accel!(body, 4, prev_accel(body, 3))
		set_prev_accel!(body, 3, prev_accel(body, 2))
		set_prev_accel!(body, 2, prev_accel(body, 1))
		set_prev_accel!(body, 1, accel(body))

		jdt = prev_accel(body, 1)*(11.0/6.0) - 3*prev_accel(body, 2) + 1.5*prev_accel(body, 3) - prev_accel(body, 4)/3
		sdt2 = 2*prev_accel(body, 1) - 5*prev_accel(body, 2) + 4*prev_accel(body, 3) - prev_accel(body, 4)
		cdt3 = prev_accel(body, 1) - 3*prev_accel(body, 2) + 3*prev_accel(body, 3) - prev_accel(body, 4)

		body.vel = old_vel + prev_accel(body, 1)*dt + (-jdt/2.0 + sdt2/6.0 - cdt3/24.0)*dt
		body.pos = old_pos + body.vel*dt + (-prev_accel(body, 1)/2.0 + jdt/6.0 - sdt2/24.0)*dt*dt
	end
end

function hermite_step(body::Body, dt)
	old_pos = body.pos
	old_vel = body.vel
	old_acc = accel(body)
	old_jerk = jerk(body)

	body.pos += body.vel*dt + old_acc*(dt*dt/2.0) + old_jerk*(dt*dt*dt/6.0)
	body.vel += old_acc*dt + old_jerk*(dt*dt/2.0)
	body.vel = old_vel + (old_acc + accel(body))*(dt/2.0) + (old_jerk - jerk(body))*(dt*dt/12.0)
	body.pos = old_pos + (old_vel + body.vel)*(dt/2.0) + (old_acc - accel(body))*(dt*dt/12.0)
end


text:
# represents a rectangle
struct ℛ{T}
  left::T
  bottom::T
  width::T
  height::T
end

left(r::ℛ) = r.left
right(r::ℛ{<:Int}) = r.left + r.width - 1
right(r::ℛ{<:AbstractFloat}) = r.left + r.width
width(r::ℛ) = r.width
bottom(r::ℛ) = r.bottom
top(r::ℛ{<:Int}) = r.bottom + r.height - 1
top(r::ℛ{<:AbstractFloat}) = r.bottom + r.height
height(r::ℛ) = r.height
aspectratio(r::ℛ) = r.width / r.height

left(r::Node) = left(r[])
right(r::Node) = right(r[])
width(r::Node) = width(r[])
bottom(r::Node) = bottom(r[])
top(r::Node) = top(r[])
height(r::Node) = height(r[])
aspectratio(r::Node) = aspectratio(r[])

function ℛ(r::ℛ; left=missing, bottom=missing, width=missing, height=missing)
  left === missing && (left = r.left)
  bottom === missing && (bottom = r.bottom)
  width === missing && (width = r.width)
  height === missing && (height = r.height)
  ℛ(left, bottom, width, height)
end

ℛ(r::Node; kwargs...) = ℛ(r[]; kwargs...)

abstract type Canvas end

Base.@kwdef struct DataSource
  generate!::Function
  xyrect::ℛ{Float64}       # default extents in data coordinates
  data = missing
  rubberband::Bool = true
  clim::Tuple = (0.0, 1.0) # default color limits
  xmin::Float64 = -Inf64
  xmax::Float64 = Inf64
  ymin::Float64 = -Inf64
  ymax::Float64 = Inf64
  cmin::Float64 = -Inf64
  cmax::Float64 = Inf64
  xpan::Bool = true
  ypan::Bool = true
  cpan::Bool = true
  xzoom::Bool = true
  yzoom::Bool = true
  czoom::Bool = true
  xylock::Bool = false
  minwidth::Float64 = 0
  maxwidth::Float64 = Inf64
  minheight::Float64 = 0
  maxheight::Float64 = Inf64
end

Base.@kwdef struct Viz
  scene::Scene
  ijrect::Node{ℛ{Int}}       # size of window in pixels
  selrect::Node{ℛ{Int}}      # last selected rectangle
  children::Vector{Canvas}
end

Base.@kwdef struct HeatmapCanvas <: Canvas
  parent::Viz
  ijhome::Node{ℛ{Int}}       # home location for canvas in pixels
  ijrect::Node{ℛ{Int}}       # current location for canvas in pixels
  xyrect::Node{ℛ{Float64}}   # canvas extents in data coordinates
  clim::Node{Tuple}
  buf::Node{Matrix{Float32}}
  datasrc::DataSource
  dirty::Node{Bool}
  task::Ref{Task}
end

Base.@kwdef struct ScatterCanvas <: Canvas
  parent::Viz
  ijhome::Node{ℛ{Int}}       # home location for canvas in pixels
  ijrect::Node{ℛ{Int}}       # current location for canvas in pixels
  xyrect::Node{ℛ{Float64}}   # canvas extents in data coordinates
  buf::Node{Vector{Point2f0}}
  datasrc::DataSource
  dirty::Node{Bool}
  task::Ref{Task}
end

Base.@kwdef struct LineCanvas <: Canvas
  parent::Viz
  ijhome::Node{ℛ{Int}}       # home location for canvas in pixels
  ijrect::Node{ℛ{Int}}       # current location for canvas in pixels
  xyrect::Node{ℛ{Float64}}   # canvas extents in data coordinates
  buf::Node{Vector{Point2f0}}
  datasrc::DataSource
  dirty::Node{Bool}
  task::Ref{Task}
end

Base.show(io::IO, x::Viz) = println(io, typeof(x))
Base.show(io::IO, x::Canvas) = println(io, typeof(x))


text:
using Reproduce
using JLD2
# using Config

function arg_parse_settings(as::ArgParseSettings = ArgParseSettings())

    @add_arg_table as begin
        "--opt1"
        arg_type=Int64
        "--opt2"
        arg_type=Int64
        "--steps"
        arg_type=Int64
    end
    return as
end

function main_experiment(args::Vector{String}, extra_arg)
    arg_settings = arg_parse_settings()
    parsed = parse_args(args, arg_settings)
    main_experiment(parsed, extra_arg)
end

function main_experiment(parsed::Dict, extra_arg = nothing)

    create_info!(parsed, parsed["save_dir"])

    j = 0
    if parsed["opt1"] == 2
        throw("Oh No!!!")
    end

    @save joinpath(parsed["_SAVE"], "data.jld2") parsed

    return j
end


text:
export Directory, getDirectoryEntries, gidsAllUniquelyOwned

# methods and docs based straight off Epetra_Directory to match Comm

"""
A base type as an interface to allow BlockMap objects to reference non-local
elements.

All subtypes must have the following methods, with DirectoryImpl standing in for
the subtype:

getDirectoryEntries(directory::DirectoryImpl, map::BlockMap, globalEntries::AbstractArray{GID},
        high_rank_sharing_procs::Bool)::Tuple{AbstractArray{PID}, AbstractArray{LID}}
        where GID <: Integer where PID <: Integer where LID <:Integer
    - Returns processor and local id infor for non-local map entries.  Returns a tuple
        containing
            1 - an Array of processors owning the global ID's in question
            2 - an Array of local IDs of the global on the owning processor

gidsAllUniquelyOwned(directory::DirectoryImpl)
    - Returns true if all GIDs appear on just one processor
"""
abstract type Directory{GID <: Integer, PID <:Integer, LID <: Integer}
end

#### required methods documentation stubs ####

"""
    getDirectoryEntries(directory, map::BlockMap{GID, PID, LID}, globalEntries::AbstractArray{GID}, high_rank_sharing_procs::Bool)::Tuple{AbstractArray{PID}, AbstractArray{LID}}

Returns processor and local id information for non-local map entries.  Returns a tuple containing
1. an Array of processors owning the global ID's in question
2. an Array of local IDs of the global on the owning processor
"""
function getDirectoryEntries end

"""
    gidsAllUniquelyOwned(directory)

Returns true if all GIDs appear on just one processor
"""
function gidsAllUniquelyOwned end


text:
function load_input(path)
  map(line -> parse(Int32, line), readlines(path))
end

function solve_1(data)
  values = BitSet()
  for val in data
    complement = 2020 - val
    if complement in values
      return complement * val
    end

    push!(values, val)
  end
end

function solve_2(data)
  sort!(data, alg=InsertionSort)
  len = length(data)
  for val in data
    l = 1
    h = len
    while h > l
      sum = val + data[l] + data[h]
      if sum == 2020
        return val * data[l] * data[h]
      elseif sum > 2020
        h -= 1
      else
        l += 1
      end
    end
  end
end

function solve(data)
  println(@time solve_1(data))
  println(@time solve_2(data))
end

data = load_input("input.txt")
solve(data)

max_stars_repo_name: wbthomason/advent

max_issues_repo_name: wbthomason/advent

max_forks_repo_name: wbthomason/advent


text:
function prepcouplingalgorithm(N_spins::Int)
    # for Hamiltonian.
    Id = getsingleId()
    Ix = getsingleIx()
    Iy_no_im = getsingleIynoim()
    Iz = getsingleIz()

    # for testing or coherence.
    Ip = getsingleI⁺()
    Im = getsingleI⁻()
    Iy = getsingleIy()

    Ims = collect( Im for j = 1:N_spins )
    Im_full = Kronecker.kroneckersum(Ims...)

    Izs = collect( Iz for j = 1:N_spins )
    Iz_full = Kronecker.kroneckersum(Izs...)

    Ixs = collect( Ix for j = 1:N_spins )
    Ix_full = Kronecker.kroneckersum(Ixs...)

    Iys_no_im = collect( Iy_no_im for j = 1:N_spins )
    Iys_no_im_full = Kronecker.kroneckersum(Iys_no_im...)

    return Id, Ix, Iy_no_im, Iz, Ip, Im, Iy,
        Im_full, Iz_full, Ix_full, Iys_no_im_full
end

function prepcouplingalgorithm(N_protons_group::Vector{Int})
    N_groups = length(N_protons_group)

    out_array = Vector{Any}(undef, N_groups)
    for i = 1:N_groups
        out_array[i] = prepcouplingalgorithm(N_protons_group[i])
    end

    out = tuple(out_array...)
    return out
end

function setupcompoundspectrum!( css_sys::Vector{Vector{T}},
    J_vals_sys,
    J_inds_sys,
    intermediates_sys,
    ppm2hzfunc,
    cs_singlets,
    N_spins_singlet::Vector{Int},
    fs::T,
    SW::T;
    tol_coherence = 1e-2) where T <: Real

    @assert length(J_vals_sys) == length(J_inds_sys) == length(css_sys)
    #
    N_sys = length(css_sys)
    N_singlets = length(cs_singlets)

    N_groups = N_sys + N_singlets

    # non-singlet objects.
    states_sys = Vector{Vector{Int}}(undef, N_sys)
    coherence_state_pairs_sys = Vector{Vector{Tuple{Int,Int}}}(undef, N_sys)
    eig_vals_sys = Vector{Vector{T}}(undef, N_sys)
    Q_sys = Vector{Vector{Vector{T}}}(undef, N_sys)
    H_sys = Vector{Matrix{T}}(undef, N_sys)
    coherence_mat_sys = Vector{Matrix{T}}(undef, N_sys)
    ms_sys = Vector{Vector{Vector{T}}}(undef, N_sys)
    M_sys = Vector{Vector{T}}(undef, N_sys)

    # output buffers.
    αs = Vector{Vector{T}}(undef, N_groups)
    Ωs = Vector{Vector{T}}(undef, N_groups)

    # Spin systems.
    for i = 1:N_sys

        αs[i], Ωs[i], coherence_mat_sys[i], eig_vals_sys[i], Q_sys[i], coherence_state_pairs_sys[i],
        H_sys[i], H1, H2, M_sys[i] = evalSCalgorithm!(css_sys[i],
        J_vals_sys[i], J_inds_sys[i],
        intermediates_sys[i],
        ppm2hzfunc;
        tol_coherence = tol_coherence)

        # normalize intensity according to number of spins.
        N_spins = length(css_sys[i])
        normalizetoNspins!(αs[i], N_spins)

        Id = getsingleId()
        ms_sys[i] = computequantumnumbers(Q_sys[i], Id)

        tmp1 = collect( coherence_state_pairs_sys[i][j][1] for j = 1:length(coherence_state_pairs_sys[i]))
        tmp2 = collect( coherence_state_pairs_sys[i][j][2] for j = 1:length(coherence_state_pairs_sys[i]))
        states_sys[i] = unique([tmp1; tmp2])
    end

    # singlets. evalRayleighproxy!() updates singlets for ΩS.
    for i = 1:N_singlets
        αs[i+N_sys] = [ N_spins_singlet[i] ]
        Ωs[i+N_sys] = [ ppm2hzfunc.(cs_singlets[i]) .* (2*π) ]
    end

    return αs, Ωs, coherence_mat_sys, eig_vals_sys, Q_sys,
    coherence_state_pairs_sys, H_sys, states_sys, ms_sys, M_sys
end

function setupapproximateΩmolecule!( css_sys::Vector{Vector{T}},
                                    p_cs_sys::Vector{Vector{T}},
                                    J_vals_sys,
                                    J_IDs_sys,
                                    intermediates_sys,
                                    cs_LUT,
                                    ppm2hzfunc,
                                    cs_singlets_compact,
                                    N_spins_singlet::Vector{Int},
                                    N_spins_sys::Vector{Int},
                                    fs::T,
                                    SW::T;
                                    tol_coherence = 1e-2,
                                    prune_a_tol_factor = 0.05,
                                    fuse_radius = ppm2hzfunc(0.001) - ppm2hzfunc(0.0),
                                    Minkowski_parameter = 3.5,
                                    update_α_flag = false) where T <: Real

    #
    N_sys = length(p_cs_sys)
    N_singlets = length(cs_singlets_compact)

    N_groups = N_sys + N_singlets

    # objects from full solve.
    updateΩFSfuncs = Vector{Function}(undef, N_sys)

    # common objects.
    statess = Vector{Vector{Int}}(undef,N_sys)
    state_labelss = Vector{Vector{Tuple{Int,Int}}}(undef,N_sys)
    eig_vals = Vector{Vector{T}}(undef, N_groups)

    # output buffers.
    αs = Vector{Vector{T}}(undef, N_groups)
    Ωs = Vector{Vector{T}}(undef, N_groups)

    # Spin systems.
    for i = 1:N_sys

        # sanity-check.
        @assert N_spins_sys[i] == length(css_sys[i])

        a0, F0, statess[i], state_labelss[i], H, H1, H2, eig_vals0,
            gs = setupanchorspectrum(p_cs_sys[i], css_sys[i],
                                J_vals_sys[i], J_IDs_sys[i],
                                intermediates_sys[i],
                                cs_LUT[i], ppm2hzfunc;
                                tol_coherence = tol_coherence,
                                prune_a_tol_factor = prune_a_tol_factor,
                                fuse_radius = fuse_radius,
                                Minkowski_parameter = Minkowski_parameter)
        # set up buffers.
        eig_vals[i] = copy(eig_vals0)
        Ωs[i] = copy(F0)
        αs[i] = copy(a0)

        # fullsolve.
        updateΩFSfuncs[i] = setupapproximateΩfullsolve!(eig_vals[i],
            αs[i], Ωs[i], statess[i], state_labelss[i],
            H, gs, fs, SW; update_α_flag = update_α_flag)

    end

    # singlets. evalRayleighproxy!() updates singlets for ΩS.
    for i = 1:N_singlets
        αs[i+N_sys] = [ N_spins_singlet[i] ]
        Ωs[i+N_sys] = [ ppm2hzfunc.(cs_singlets_compact[i]) .* (2*π) ]
    end

    return αs, eig_vals, updateΩFSfuncs, Ωs,
            statess, state_labelss
end

# A version of setupapproximateΩrayleigh() that:
#   - J remains constant, taken from pivot J-coupling values.
#   - fixed coherence using eigenvectors from pivot p_cs.
#   - use full eigensolve to get Ω.
#   - option to use updated eigenvectors to update α.
function setupanchorspectrum( p_cs::Vector{T},
                            css::Vector{T},
                            J_vals::Vector{T},
                            J_IDs::Vector{Tuple{Int,Int}},
                            intermediates,
                            cs_LUT::Vector{Vector{Int}},
                            ppm2hzfunc;
                            tol_coherence = 1e-2,
                            prune_a_tol_factor = 0.05,
                            fuse_radius = ppm2hzfunc(0.001) - ppm2hzfunc(0.0),
                            Minkowski_parameter = 3.5) where T <: Real

    #
    N_spins = length(css)

    F0, a0, state_labels, H, H1, H2, eig_vals, gs,
       _ = getreducedcomponentsforspinsystem(p_cs, css,
                     J_vals, J_IDs, intermediates,
                     cs_LUT, ppm2hzfunc;
                     tol_coherence = tol_coherence,
                     prune_a_tol_factor = prune_a_tol_factor,
                     fuse_radius = fuse_radius,
                     Minkowski_parameter = Minkowski_parameter)

    # prepare intermediate objects and buffers.
    tmp1 = collect( state_labels[i][1] for i = 1:length(state_labels))
    tmp2 = collect( state_labels[i][2] for i = 1:length(state_labels))
    states::Vector{Int} = unique([tmp1; tmp2])

    return a0, F0, states, state_labels, H, H1, H2, eig_vals, gs
end

function getreducedcomponentsforspinsystem(p_cs::Vector{T},
            css::Vector{T},
            J_vals::Vector{T},
            J_IDs::Vector{Tuple{Int,Int}},
            intermediates,
            cs_LUT::Vector{Vector{Int}},
            ppm2hzfunc;
            tol_coherence = 1e-2,
            prune_a_tol_factor = 0.05,
            fuse_radius = ppm2hzfunc(0.001) - ppm2hzfunc(0.0),
            Minkowski_parameter = 3.5) where T <: Real

    #
    a0, F0, coherence_mat, eig_vals, Q, state_labels0,
    H, H1, H2, M_array = evalSCalgorithm!(p_cs, css,
    J_vals, J_IDs, intermediates,
    cs_LUT, ppm2hzfunc;
    tol_coherence = tol_coherence)

    # default: no pruning nor merging.
    F3 = F0
    a2 = a0
    state_labels = state_labels0

    if isfinite(prune_a_tol_factor)

        ### prune.
        prune_zero_tol = maximum(a0) * prune_a_tol_factor
        keep_flags = abs.(a0) .> prune_zero_tol
        a_pruned = a0[keep_flags]
        F_pruned = F0[keep_flags]
        state_labels_pruned = state_labels0[keep_flags]

        # # debug.
        # println("fuse_radius = ", fuse_radius)

        # # drop.
        # a2, F2 = dropcomponents(F, F, a; radius = fuse_radius*2*π,
        #                            Minkowski_parameter = Minkowski_parameter)

        ### merge.
        F2, a2 = approxΩsαs(F_pruned, F_pruned, a_pruned; radius = fuse_radius*2*π,
                                   Minkowski_parameter = Minkowski_parameter)

        # state labels.
        state_labels = findapproxstates(F_pruned, state_labels_pruned, F2)

        F3 = getΩfroms(eig_vals, state_labels)
    end

    # normalize intensity according to number of spins.
    N_spins = length(css)
    normalizetoNspins!(a2, N_spins)

    #TODO old possibly irrelevant note: normalize by # of hydrogens css_sys. hunt down the code for cs_singlet (probably at runtime), and do the same.

    return F3, a2, state_labels, H, H1, H2, eig_vals, Q, state_labels0, a0, F0, coherence_mat, F2
end

# #### proxy to the strong coupling.
# parsefunc converts low-dim p to dim(cs).
# - It is molecule and spin group specific.
# uses p_cs to modify cs.
# TODO: add ! to function name since this mutates cs.
function evalSCalgorithm!(  cs::Vector{T},
                            J_vals::Vector{T},
                            J_inds::Vector{Tuple{Int,Int}},
                            intermediates,
                            ppm2hzfunc;
                            tol_coherence::T = 1e-3) where T <: Real

    # set up.
    #pcstocs!(cs, p_cs, cs_LUT)

    Id, Ix, Iy_no_im, Iz, Ip, Im, Iy, Im_full, Iz_full, Ix_full,
    Iys_no_im_full = intermediates # see prepcouplingalgorithm(Int) for details.

    # strong coupling algorithm.
    ω0 = ppm2hzfunc.(cs) .* 2*π
    #H = getgenericHamiltonian(Id, Ix, Iy_no_im, Iz, ω0, J_vals, J_inds)
    H, H1, H2 = getgenericHamiltoniantest(Id, Ix, Iy_no_im, Iz, ω0, J_vals, J_inds)

    a, F, p, s, Q, coherence_labels, M_array = getaΩ(Iz_full, H,
        Iys_no_im_full, Ix_full; tol = tol_coherence)
    #println("hi")

    return a, F, p, s, Q, coherence_labels, H, H1, H2, M_array
end

# for now, assign p is Δcs, and hard code the conversion of
#   Δcs to ΔΩ0 as unit conversion for fullsolve.
# Here, updateΩfuncs takes Δcs instead of cs.
## TODO later: generic version where p reads into p_cs, then converted to css_sys.
#   the units of p_cs needs to be defined before calling this func.
function evalfullsolveproxy!(css_sys_buffer::Vector{Vector{T}},
                            cs_singlets_compact0::Vector{T},
                            #N_singlets::Int,
                            Ωs::Vector{Vector{T}},
                            p::Vector{T},
                            st_ind::Int,
                            cs_LUT::Vector{Vector{Vector{Int}}},
                            updateΩfuncs::Vector{Function},
                            fs::T,
                            SW::T,
                            ppm2hzfunc) where T <: Real
    ### checks.
    N_sys = length(css_sys_buffer)
    @assert length(updateΩfuncs) == N_sys == length(cs_LUT)

    N_singlets = length(cs_singlets_compact0)

    N_groups = N_sys + N_singlets
    @assert length(Ωs) == N_groups

    ### set up.
    fin_ind = st_ind -1

    ### multi-spin systems.
    for i = 1:N_sys

        ## parse to css_sys.
        N_cs_vars = length(cs_LUT[i])

        st_ind = fin_ind +1
        fin_ind = st_ind + N_cs_vars -1
        p_cs = p[st_ind:fin_ind]
        pcstocs!(css_sys_buffer[i], p_cs, cs_LUT[i])

        ## convert units.
        convertΔcstoΔω0!(css_sys_buffer[i], fs, SW)

        #updaterfuncs(ω0)
        updateΩfuncs[i](css_sys_buffer[i])
    end

    ### singlets.
    for i = 1:N_singlets

        ## parse.
        fin_ind += 1

        # ## convert units and update.
        # Ωs[i+N_sys][1] = p[fin_ind]*2*π*fs/SW

        # p is Δcs, and new cs = old cs + Δcs.
        tmp = cs_singlets_compact0[i] + p[fin_ind]

        ## update.
        Ωs[i+N_sys][1] = ppm2hzfunc(tmp) * 2*π
    end

    return fin_ind
end

function convertΔcstoΔω0!(x::Vector{T}, fs::T, SW::T) where T
    for i = 1:length(x)
        x[i] = x[i]*2*π*fs/SW
    end
end

function convertΔcstoΔω0(x::T, fs::T, SW::T)::T where T
    return x*2*π*fs/SW
end

### from eigen vectors to Zeeman states. See the lower case m_j^{(r)}
#   in sec. 18.2 (pag 468), Spin Dynamics.
function computequantumnumbers(basis::Vector{Vector{T}}, Id) where T

    N_spins = round(Int, log(2,length(basis)))
    I_jz_set = collect( getmultiIz(Id, j, N_spins) for j = 1:N_spins )

    m_basis = Vector{Vector{T}}(undef, length(basis))
    for r = 1:length(basis)
        m_basis[r] = collect( getzangularmomentum(basis[r], I_jz_set[j]) for j = 1:N_spins )
    end

    return m_basis
end


text:
"""
  Package: FishEBM
  File: FishEBM.jl
  Justin Angevaare, Devin Rose
  Module definition for an enviro-based model
  Created: May 2015
"""

module FishEBM
  # Version number
  VersionNumber = "0.2.3"

  # Packages used by this package
  using DataArrays, DataFrames, Distributions, Gadfly, ProgressMeter

  export
    # agents.jl functions
    AgentDB,
    findCurrentStage,
    getAge,
    getAgeSpecificPop,
    getCohortNumber,
    getStagePopulation,
    injectAgents!,
    removeEmptyClass!,

    # environment.jl functions
    hashEnvironment!,
    initEnvironment,
    isEmpty,
    pad_environment!,

    # FileIO.jl functions
    ageData,
    aliveData,
    createDir,
    dateDir,
    getDirChar,
    getYearlyData,
    getYearlyPop,
    harvestData,
    killedData,
    resultsDir,
    runDir,
    setProjPath,
    simDir,
    simReadme,
    spawnData,

    # mortality.jl
    harvest!,
    kill!,
    killAgeSpecific!,

    # move.jl
    move!,

    # simulate.jl functions
    simulate,

    # spawn.jl
    spawn!,

    # summary.jl
    simSummary,

    # types.jl
    AdultAssumptions,
    AgentAssumptions,
    ClassPopulation,
    EnviroAgent,
    EnvironmentAssumptions,

    # utilities.jl functions
    initPopulationDensity,
    updatePopulationDensity!

    # include types.jl in the module first, types defined in this function are used in other files
    include("types.jl")

    # include any source files with a function used in the module
    include("agents.jl")
    include("environment.jl")
    include("FileIO.jl")
    include("mortality.jl")
    include("move.jl")
    include("simulate.jl")
    include("spawn.jl")
    include("summary.jl")
    include("utilities.jl")
end


text:
using ArgParse

using DifferentialEquations
using DiffEqCallbacks
using Random
using SparseArrays
using StatsBase
using Distributions

# Piecewise linear functions

# Piecewise linear interpolation at multiple positions (specialized for Float64)
function interp1d(x::Array{Float64,1}, x0::Array{Float64,1}, y0::Array{Float64,1})
    # x - position to evaluate function
    # x0 - positions of discontinuities in function
    # y0 - values of function on each region: y = y0[i] for x0[i] < x < x0[i+1]

    y = zeros(Float64, size(x))
    jj = 1
    for ii in 1:length(y)
        if x0[1]<=x[ii]<x0[end]
            jj = findnext(val->val>x[ii],x0,jj)-1
            y[ii] = y0[jj] + (y0[jj+1]-y0[jj])/(x0[jj+1]-x0[jj])*(x[ii]-x0[jj])
        end
    end
    return y
end

# Piecewise linear interpolation at a single position
function interp1d(x::Number, x0::Array, y0::Array)
    # x - position to evaluate function
    # x0 - positions of discontinuities in function slope
    # y0 - values of function at each discontinuity. Function takes the value 0 outside the range [x0[1], x0[length(x0)]]

    jj = 1
    y = 0.0
    if x0[1]<=x<x0[end]
        jj = findnext(val->val>x,x0,jj)-1
        y = y0[jj] + (y0[jj+1]-y0[jj])/(x0[jj+1]-x0[jj])*(x-x0[jj])
    end
    return y
end

# Piecewise constant interpolation
function piecewise(x::Number, x0::Array, y0::Array)
    # x - position to evaluate function
    # x0 - positions of discontinuities in function
    # y0 - values of function on each region: y = y0[i] for x0[i] < x < x0[i+1]

    jj = 1
    y = 0.0
    if x0[1]<=x<x0[end]
        jj = findnext(val->val>x,x0,jj)-1
        y = y0[jj]
    end
    return y
end

# Inverse of a piecewise constant density function
function piecewise_density_inv(x0::Array{T,1}, y0::Array{T,1}) where {T <: Number}
    csum = zeros(T, size(x0))
    for i in 2:length(x0)
        csum[i] = (x0[i] - x0[i-1])*y0[i-1] + csum[i-1]
    end
    csum = csum / csum[end]
    csum, x0
end

# Inverse of a piecewise linear density function
function linear_density_inverse(x0::Array{T,1}, y0::Array{T,1}) where {T <: Number}
    csum = zeros(T, size(x0))
    for	i in 2:length(x0)
    	csum[i] = (x0[i] - x0[i-1])*(y0[i-1] + y0[i])/2	+ csum[i-1]
    end
    csum_t = csum[end]
    function inverse(x::T)
        z = x*csum_t
        jj = findnext(val->val>z, csum, 1)-1
        z = z - csum[jj]
        y = y0[jj]
        dy = (y0[jj+1] - y0[jj])/(x0[jj+1]-x0[jj])
        if abs(dy)>1e-12
            return x0[jj]+((sqrt(y*y+2*z*dy)-y)/dy)
        elseif abs(y)>1e-12
            return x0[jj]+z/y
        else
            return x0[jj]
        end
    end
    inverse
end

# L^2 / sum-squared vector norm
function norm(x)
    sqrt(sum(abs2,x))
end

# RI HEI10 escape rate function (\beta(C)*C in equation (2) )
# Modified escape rate to explore importance of the precise functional form of this term

function betaC(C, K, n_c, n_alpha)
    C ./ (1.0 .+ exp( (C ./ K) .* n_alpha) .* n_c)
end

# RI HEI10 escape rate function \beta(C)*C (evaluated in place)
function betaC!(dC, C, K, n_c, n_alpha)
#    @. dC = C .* exp( -((C ./ K) .^ n_alpha ) * n_c)
     @. dC = C ./ (1.0 .+ exp( (C ./ K) .* n_alpha) .* n_c)

end

# Derivative of (\beta(C)*C) wrt C
function d_betaC(C::Float64, K::Float64, n_c::Float64, n_alpha::Float64)::Float64
#     exp( -((C ./ K).^ n_alpha) * n_c) *(1 .- n_c .* n_alpha .* (C .^(n_alpha-1))/ (K .^ n_alpha))
      1.0 ./ (1.0 .+ exp( (C ./ K) .* n_alpha ) .* n_c) -
      C .* n_c .* (n_alpha ./ K) .* exp((C ./ K) .* n_alpha) ./ (1.0 .+ exp( (C ./ K) .* n_alpha ) .* n_c).^2

end

# Derivative of (\beta(C)*C) wrt C (evaluated in place)
function d_betaC!(dbC, C, K, n_c, n_alpha)
    @. dbC =  1.0 ./ (1.0 .+ exp( (C ./ K) .* n_alpha ) .* n_c) -
     C .* n_c .* (n_alpha ./ K) .* exp((C ./ K) .* n_alpha) ./ (1.0 .+ exp( (C ./ K) .* n_alpha ) .* n_c).^2
end

# Function to calculate the RHS of the ODE system
# Result returned in-place as first argument
function basic_ode!(dyy::Array{Float64,1}, yy::Array{Float64,1}, p, t::Float64)
    m,           # Number of SC spatial compartments for discretization
    N,           # Number of RIs
    K,           # Hill threshold
    nodes,       # Indices of SC spatial compartments containing RIs
    alpha,       # Rate of HEI10 production on SC - here set to be zero
    beta,        # Rate of HEI10 escape from SC - here set to be zero
    D,           # RI Diffusion coefficient on SC (D_c)
    dL,          # Length of each spatial compartment
    betaC!,      # beta(C)*C function (passed to ODE for efficiency)
    d_betaC,     # derivative of beta(C)*C function wrt C
    n_c,         # Hill coefficient (\gamma)
    n_alpha,
    a_nodes,     # RI HEI10 absorption rate from SC (\alpha)
    b_nodes,     # RI HEI10 escape rate back onto SC (\beta)
    b2_nodes,    # RI HEI10 escape rate into nucleoplasm  (not used)
    b3_nodes,    # RI HEI10 rate  (not used)
    bC,          # Placeholder to re-use vector
    A = p

    # Extract variables from ODE State vector
    @inbounds y = @view yy[1:m] # HEI10 concentration on the SC (c_i, 1<=i<=m)
    @inbounds dy = @view dyy[1:m]
    @inbounds C = @view yy[m+1:m+N] # HEI10 amounts at the RH (C_j 1<=j<=N)
    @inbounds dC = @view dyy[m+1:m+N]

    h = D/dL/dL

    # Diffusion (production and decay - not used) of HEI10 on the SC.
    # Standard second-order finite difference approximation of 1D laplacian
    @inbounds dyy[1] = alpha - beta * y[1] + h * (y[2] - y[1])
    @inbounds @simd for i = 2:(m-1)
        dyy[i] = alpha - beta * y[i] +
                h * (y[i+1] - 2 * y[i] + y[i-1])
    end

    @inbounds dyy[m] = alpha - beta * y[m] + h * (y[m-1] - y[m])

    # Calculate escape of HEI10 back onto SC.
    betaC!(bC, max.(C,0.0), K, n_c, n_alpha)

    # Calculate net fluxes from SC to RIs, and update ODE RHS for appropriate SC compartment
    @inbounds for j = 1:N
        A[j] = a_nodes * y[nodes[j]] - b_nodes * bC[j]
        dyy[nodes[j]] = dyy[nodes[j]] - A[j] / dL
    end

    # Rate of change of HEI10 amounts at the RIs
    @. dC = A .- b2_nodes .* bC .- b3_nodes .* C
end

# Function to calculate the Jacobian of the ODE system RHS - helps implicit integration method be more efficient
# Result returned in-place as first argument
function jac!(J, yy, p, t)
    m, # For parameter interpretation please see basic ODE
    N,
    K,
    nodes,
    alpha,
    beta,
    D,
    dL,
    betaC!,
    d_betaC,
    n_c,
    n_alpha,
    a_nodes,
    b_nodes,
    b2_nodes,
    b3_nodes,
    bC,
    A = p
    r = D / dL / dL
    C = @view yy[m+1:m+N]
    #@show typeof(J)
    J[1, 1] = -beta - r
    J[1, 2] = r
    @inbounds for j = 2:m-1
        J[j, j-1] = r
        J[j, j] = -beta - 2 * r
        J[j, j+1] = r
    end
    @inbounds J[m, m-1] = r
    @inbounds J[m, m] = -beta - r
    dbC::Float64 = 0.0
    idx::Int64 = 0
    @inbounds for j = 1:N
        dbC = d_betaC(max(C[j],0.0), K, n_c, n_alpha)
        J[m+j, m+j] = -(b_nodes + b2_nodes) * dbC - b3_nodes
        J[nodes[j], m+j] = b_nodes * dbC / dL
        J[m+j, nodes[j]] = a_nodes
        J[nodes[j], nodes[j]] -= a_nodes / dL
    end
end

function simulate(idx, args)

    # Initialize random number generator with seed for this simulation
    mt = MersenneTwister(idx)

    # Generate SC length, sampled from Normal distribution clipped to +/- 3 S.D.
    L::Float64 = Float64(args["L"])+Float64(args["Ls"])*clamp(randn(mt), -3, 3)
    # Calculate appropriate number of RIs for this SC length
    N::Int64 = Int64(round(L*args["density"]))
    # Number of intervals for spatial discretization of SC
    m::Int64 = Int64(args["m"])

    # Generate random RI positions
    x0 = sort(rand(mt, N))
    x::Array{Float64,1} = L * x0

    # Length of each interval in  SC discretization
    dL::Float64 = L / m
    # Calculate indices of SC intervals containing each RI
    nodes = map(a -> ceil(Int, a), x / dL)

    # Midpoints of each SC interval
    xm = (0.5:1.0:m) * dL

    n_c::Float64 = Float64(args["n_c"]) # Hill coefficient (gamma)
    n_alpha::Float64 = Float64(args["n_alpha"]) # Hill coefficient (gamma)
    K::Float64 = Float64(args["K"]) # Hill threshold (K_C)

    alpha = 0.0 # HEI10 production rate on SC (Not used here)
    beta = 0.0 # HEI10 decay rate on SC (Not used here)
    a_nodes::Float64 = Float64(args["a_nodes"]) # RI HEI10 absorption rate from SC (\alpha)
    b_nodes::Float64 = Float64(args["b_nodes"]) # RI HEI10 escape rate back onto SC (\beta)
    b2_nodes = 0.0 # RI escape rate (loss to nuceloplasm) - not used here
    b3_nodes = 0.0 # RI constant escape rate (loss to nuceloplasm) - not used here
    D::Float64 = Float64(args["D"]) # RI Diffusion coefficient on SC (D_c)

    bC = zeros(Float64, N) # Buffer to store values of \beta(C)*C in ODE calculation avoiding allocating memory
    A = zeros(Float64, N) # Buffer to store fluxes of HEI10 from RI to SC in ODE calculation

    # Initial uniform HEI10 concentration (c) on SC
    u0::Array{Float64,1} = fill(Float64(args["u0"]), m)

    # Initialize truncated normal distribution for RI HEI10 loading
    # Mean C0_noise (C_0)
    # Standard deviation C0_noise (\sigma)
    # Truncated at +/- 3 s.d.
    d = Truncated(Normal(Float64(args["C0"]), Float64(args["C0_noise"])), max(0, Float64(args["C0"])-3*Float64(args["C0_noise"])), Float64(args["C0"])+3*Float64(args["C0_noise"]))

    # Generate initial RI HEI10 amounts
    C0::Array{Float64,1} = rand(mt, d, N)
    C0[C0.<0] .= 0.0 # Double-check that initial RI HEI10 concentrations are non-negative

    C0_t = 1.0
    C0_t2::Float64 = Float64(args["t_C0_ratio2"]) # ($f_e$)
    t_L::Float64 = Float64(args["t_L"])

    # Additional HEI10 loading for RI near ends of SCs (telomeres)
    if !args["t_exp"]
        # Piecewise linear function ($f$)
        for i in 1:N
    	    p = x[i]/L
            if p < t_L
               r = p/t_L
            C0[i] *= C0_t*(r) + C0_t2*(1-r)
            elseif p> 1-t_L
               r = (1-p)/t_L
               C0[i] *= C0_t*(r) + C0_t2*(1-r)
       	     end
	 end
    else
        # Exponentially decaying RI end weights
	for i in 1:N
    	    p = x[i]/L
	    C0[i] *= (C0_t2-C0_t)*(exp(-2*p/t_L) + exp(-2*(1-p)/t_L))+ C0_t
	 end
    end

    # Construct initial ODE state vector
    y0::Array{Float64,1} = vcat(u0, C0)

    # ODE system parameters
    p = (
        m,
        N,
        K,
        nodes,
        alpha,
        beta,
        D,
        dL,
        betaC!,
        d_betaC,
        n_c,
        n_alpha,
        a_nodes,
        b_nodes,
        b2_nodes,
        b3_nodes,
        bC,
        A,
    )

    # Initialize ODE problem with jacobian
    j = spzeros(m + N, m + N)
    g2 = ODEFunction(basic_ode!; jac=jac!, jac_prototype = j)
    prob_g = ODEProblem(g2, y0, (0.0, Float64(args["n_ts"]) * Float64(args["dt"])), p)
    cb = PositiveDomain(g2) # Additional constraint to ensure non-negativity of concentration values during integration

    # Solve ODE system using implicit Rodas5 solver, with tight relative error tolerance
    sol = solve(prob_g,Rodas5(autodiff=:false),cb=cb, retol=1e-12, save_everystep=false)
    # Return final state
    (L, x, sol.u[end])
end

function simulate_all(args, io)

    # Parallel threaded loop to run simulations for n consecutive seeds, beginning at "start"
    m::Int64 = args["m"]
    solution_data = Vector{Vector{}}()
    for i in 1:Threads.nthreads()
        push!(solution_data,String[])
    end
    Threads.@threads for i in args["start"]:args["start"]+(args["n"]-1)

        L, x, u = simulate(i, args)
        push!(solution_data[Threads.threadid()], string(L) *'\n'*join(map(string, x), ',')*'\n'*join(map(string, u[m+1:end]), ',')*'\n')
    end

    # Combine data from all threads and write to output file

    data = vcat(solution_data...)
    for i in 1:length(data)
        println(io, data[i])
    end
end

function parse_commandline()

# Parse command line options and parameters
# Note that the driver file (run_simulations.py) overrides most of these parameter values

    s = ArgParseSettings()

    @add_arg_table! s begin
        "--filename_base"
            default="test"
            help="base filename"
        "--start"
            arg_type=Int
            default=1
            help="first seed to use for the simulation"
        "--n"
            arg_type=Int
            default=100 # (No unit)
            help="number of seeds" # Number of different SCs to simulate, using a different RNG seed for each one
        "--n_ts"
            arg_type=Int
            default=100 # (no unit)
            help="number of timesteps" # Number of integration periods (again, timestep for numerical integration is shorter)
	                               # In actual fact, for this simulation there is a single integration period of length
				       # n_ts * dt, as only final state is of interest.
        "--L"
	    arg_type=Float64
	    default=30.0 #(um)
            help="length of chromosome" # Mean length of this SC
       "--Ls"
            arg_type=Float64
	    default=0.0 # (um)
            help="std dev in length of chromosome"  # Variation in length of this SC
        "--dt"
     	    arg_type=Float64
	    default=360.0 # (s)
            help="time interval" # Length of each integration period (timestep is shorter, controlled automatically by error estimate)
        "--m"
	    arg_type=Int
	    default=2000 # (no unit)
            help="number of elements" # (M) Number of intervals for spatial discretization of SC
        "--K"
     	    arg_type=Float64
            default=1.0 # (a.u.)
            help="threshold for RI unbinding" # (K_C) Hill function thresold for escape of HEI10 from RI
        "--n_c"
    	    arg_type=Float64
    	    default=1.25 # (no units)
            help="Hill coefficient for RI unbinding" # (\gamma) Hill coefficient for escape of HEI10 from RI
        "--n_alpha"
            arg_type=Float64
            default=1.0 # (no units)
            help="Hill coefficient for RI unbinding" # (\gamma) Hill coefficient for escape of HEI10 from RI
        "--a_nodes"
            arg_type=Float64
            default=2.1 # (um s^{-1})
            help="binding rate at RIs" # (\alpha) RI HEI10 absorption rate from SC
        "--b_nodes"
            arg_type=Float64
            default=0.5 # (s^{-1})
            help="unbinding rate at RIs" # (\beta) RI HEI10 escape rate back onto SC
        "--u0"
            arg_type=Float64
            default=0.16 # (a.u.)
            help="initial HEI10 concentration" # (c_0) HEI10 initial loading on SC
        "--C0"
            arg_type=Float64
            default=0.93 # (a.u.)
            help="initial RI HEI10" # (C_0) HEI10 initial RI loading
        "--C0_noise"
            arg_type=Float64
            default=0.30 # *(a.u.)
            help="initial RI HEI10 noise" # (\sigma) Noise in HEI10 initial RI loading
        "--D"
            arg_type=Float64
            default=1.1 # (um^2 s^{-1})
            help="HEI10 diffusion constant on SC"  # (D_c) RI Diffusion coefficient on SC
        "--t_C0_ratio2"
            arg_type=Float64
            default=1.5 # (no units)
	    help="initial hei10 ratio for RIs in telomeres"  # (f_e) Increased RI loading at SC telomeres
	"--t_exp"
	    help="exponential functions for telomeres"
	    action = :store_true
        "--t_L"
            arg_type=Float64
            default=0.1 # (no units)
            help="telomere rel length" # (x_e) Telomere length as fraction of bivalent
        "--density"
            arg_type=Float64
            default=0.5         # (um^{-1})
            help="RI density"  # (\rho) SC RI density
    end
    return parse_args(s)
end

function main()
    args = parse_commandline()
    println("Parsed args:")
    for (arg,val) in args
        println("  $arg  =>  $val")
    end

    # Generate filename depending on parameter values
    filename=[]
    for (arg,val) in args
    	if length(arg)<8 || arg[1:8]!="filename"
	   push!(filename, "$val")
	end
    end
    # Add filename prefix
    if !isnothing(args["filename_base"])
        filename = args["filename_base"] * join(filename, ",") * ".dat"
    else
	filename = join(filename, ",") * ".dat"
    end

    println("filename", filename)
    open(filename,"w") do io
        # Write parameter values as first line of output file
        println(io, "#Namespace(" * join(map(x-> string(x[1]) *"="*string(x[2]), collect(args)), ",") * ")")

        # Run (group of) simulations
        simulate_all(args, io)
    end
    println("done")
end

#fac = 0.137; plot_obj([1.1, 6.8*fac, 2.2*fac, 1.2*fac, 1.5, 0.1, 2.1, 0.5, 1.25, 1.0]; n=500, m=1000)

main()


text:
# Use the object Tile and other functions for Tiles.
include("Tile.jl")

# Create function to get user input and recursively validate the user input
function get_input(text_to_print, in_same_line, allowed_options, required_params)

    # If we want to print on the same line, use `print()`. Else, `use println()`
    in_same_line ? print(text_to_print) : println(text_to_print)

    # Get input and split it into an array
    txtin = split(lowercase(readline()));

    # If the required amount of parameters is not met or if the command is not in the list of allowed options,
    # ask for input again
    if length(txtin) < required_params || !(txtin[1] in allowed_options)
        println("Please enter a valid command (see README)")
        get_input(text_to_print, in_same_line, allowed_options, required_params)
    else
        # If input is valid, return the input as a string.
        return join(txtin, " ")
    end
end

# Ask if user wants to play again
function play_again()

    ans = get_input("Play again? (y/n): ", true, ["y", "n"], 1)

    # If user wants to play again, run the game again. Else, exit()
    if ans == "y"
        run(`julia $PROGRAM_FILE`)
        exit()
    else
        exit()
    end
end

# Print the current game state
function print_board(board)

    # Gets basic variables required to iterate over array
    rows = length(board)
    cols = length(board[1])
    arr = []

    # If number in column label is more than one digit, just print 1 space instead of 2 so
    # that grid is aligned.
    for num in 1:cols
        n_spaces = num > 9 ? " " : "  "
        push!(arr, num, n_spaces)
    end

    # Join the column labels
    labels_cols = join(arr)

    # Print the column labels
    println("   | "*labels_cols)
    println("___|"*join(["_" for x in 1:length(labels_cols)]))

    # For each row, print the row labels such that the digits align properly by using spaces
    for i in 1:rows
        boxes = join([to_string(board[i][x]) for x in 1:cols], "  ")
        if length("$i") < 2
            println("$i  | "*boxes)
        else
            println("$i | "*boxes)
        end
    end
end

# If a tile is hit and it has 0 mine neighbours, then expose itself, and non-mine neighbours recursively.
function flood_fill(board, cell)
    # Iterate over all the cells to its left-top diagonal neighbour to its right-bottom diagonal neighbour.
    for x_offset in -1:1
        for y_offset in -1:1
            # Get absolute index of the current cell instead of relative index
            i = cell.x + x_offset
            j = cell.y + y_offset

            # If cell exists in the grid (the left-top corner makes the function check outside the array)
            if (i in 1:length(board)) && (j in 1:length(board[1]))
                neighbour = board[i][j]

                # If the tile is not a mine and is not shown, show it.
                # show() also calls flood_fill() if its neighbour has 0 mine neighbours, making this recursive
                if (!(neighbour.is_mine) && !neighbour.shown)
                    show(board, neighbour)
                end
            end
        end
    end
end

# Show a tile, and call flood_fill() if required.
function show(board, tile)
    tile.shown = true

    # If tile has 0 mine-neighbours, call flood_fill() on it so it exposes the required neighbours
    if tile.num_mines == 0
        flood_fill(board, tile)
    end
end

# If game is over, show the leftover mines.
function game_over(board, mine_pos_arr)
    for pos in mine_pos_arr
        (board[pos[1]][pos[2]]).shown = true
    end
end

# Create board at the start of the game
function init_board(r, c, n)

    # Create an array with r rows. We will insert our cells in an array here.
    board = Array{Any}(nothing, r)

    # Initialize an array that holds the positions of the mines.
    mine_pos_arr = []

    # Initialize a flat array with all possible positions
    options = [[(x, y) for x in 1:r, y in 1:c]...]

    # For every mine, pick up a random position in all possible positions, and remove it
    # from the possible positions so the other mines can not be placed in the same place
    # accidentally.
    for mine_pos in 1:n
        x = rand(1:length(options))
        push!(mine_pos_arr, options[x])
        deleteat!(options, x)
    end

    # For every row, create an array of cells, and place the array in its corresponding row.
    for i in 1:r
        temp = []
        for j in 1:c
            tile_mine = (i, j) in mine_pos_arr ? true : false
            push!(temp, Tile(false, tile_mine, i, j, 0, false))
        end
        board[i] = temp
    end

    # Count the number of mine neighbours for each cell.
    # If the cell is a mine, it is set to -1.
    for i in 1:r
        for j in 1:c
            count_neighbours(board, board[i][j])
        end
    end

    # Return game info
    return [board, options, mine_pos_arr]
end

# Print the difficulty selector of the game
function start_game()
    difficulty = """
    Select difficulty:
    - `easy`: 8x8 grid, with 9 mines
    - `medium`: 16x16 grid, with 20 mines (fullscreen recommended)
    - `custom`: You specify the parameters!
    """
    diff_selected = get_input(difficulty, true, ["easy", "medium", "custom"], 1)

    # Based on the selected difficulty, initialize game.
    if diff_selected == "easy"
        board = init_board(8, 8, 9)
        return board
    elseif diff_selected == "medium"
        board = init_board(16, 16, 20)
        return board
    else
        params = split(get_input("Enter the number of rows, columns and mines seperated by spaces (1-99): ", true, string.(collect(1:99)), 3))
        board = init_board(parse(Int, params[1]), parse(Int, params[2]), parse(Int, params[3]))
        return board
    end
end

# Validate a flag or mine command, ancd check if the input is inside bounds.
function validate_command(board, input)
    if (parse(Int, input[2]) > length(board)) || (parse(Int, input[3]) > length(board[1]))
        return false
    else
        return true
    end
end

# Beginning of game, print the welcome message.
println("Welcome to Minesweeper! Please have a look at the README before playing. Enjoy!")

# Get game info by initializing grid, and mines
game = start_game()

# Save info for game-result validation.
grid = game[1]
options = game[2]
mine_pos_arr = game[3]

# Print the starting state of game
print_board(grid)

# Get input at every step using a while true loop, and breaking out of loop on game end.
while true
    # Get command and parameters
    txt = split(get_input("Enter your command: ", true, ["f", "m", "exit"], 3))

    # If the user wishes to mine, and the parameters are valid, mine that cell.
    if txt[1] == "m" && (validate_command(grid, txt))

        # If the mined cell is a mine, game over.
        if (grid[parse(Int, txt[2])][parse(Int, txt[3])]).is_mine == true
            game_over(grid, mine_pos_arr)
            print_board(grid)
            println("You lose. Better luck next time!")
            break
        end

        # Show the cell the user clicked on, and show the updated state
        show(grid, grid[parse(Int, txt[2])][parse(Int, txt[3])])
        print_board(grid)

        # Check win
        win = true

        # If all non-mine locations are shown, then user wins
        for pos in setdiff(options, mine_pos_arr)
            if (grid[pos[1]][pos[2]]).shown == false
                win = false
                break
            end
        end

        # If user wins, print win message.
        if (win)
            println("Congratulations, you win!")
            break
        # If user hasn't won yet, advance to next loop.
        else
            continue
        end

    # If the user wishes to flag and parameters are valid, flag the cell
    elseif txt[1] == "f"  && (validate_command(grid, txt))

        (grid[parse(Int, txt[2])][parse(Int, txt[3])]).is_flagged = true

        # After flagging, print updated game state
        print_board(grid)

    # If command is to exit, exit game.
    elseif txt[1] == "exit"
        exit()
    # If command is invalid (out of bounds), ask user to type command again
    else
        println("Please enter valid cell coordinates.")
    end
end

# On end of game, ank if user wants to play again.
play_again()

max_stars_repo_name: soumitradev/Minesweeper

max_issues_repo_name: soumitradev/Minesweeper

max_forks_repo_name: soumitradev/Minesweeper