text:
# This is a convenience function for libraries like MaBLAS and PaddedMatrices to use explicitly, and for others to force more type inference in precompilation.

const GEMMLOOPSET = loopset(:(
  for n ∈ 1:N, m ∈ 1:M
    Cₘₙ = zero(eltype(C))
    for k ∈ 1:K
      Cₘₙ += A[m, k] * B[k, n]
    end
    C[m, n] += Cₘₙ
  end
));

# function matmul_params(rs::Int, rc::Int, cls::Int)
#     set_hw!(GEMMLOOPSET, rs, rc, cls, Int(cache_size(StaticInt(1))), Int(cache_size(StaticInt(2))), Int(cache_size(StaticInt(3))))
#     order = choose_order(GEMMLOOPSET)
#     order[5], last(order)
# end
function matmul_params(
  rs::Int,
  rc::Int,
  cls::Int;
  M = nothing,
  K = nothing,
  N = nothing,
  W = 0,
)
  set_hw!(
    GEMMLOOPSET,
    rs,
    rc,
    cls,
    Int(cache_size(StaticInt(1))),
    Int(cache_size(StaticInt(2))),
    Int(cache_size(StaticInt(3))),
  )
  if N ≢ nothing
    nloop = GEMMLOOPSET.loops[1]
    GEMMLOOPSET.loops[1] =
      Loop(:n, MaybeKnown(1), MaybeKnown(N), MaybeKnown(1), nloop.rangesym, nloop.lensym)
  end
  if M ≢ nothing
    mloop = GEMMLOOPSET.loops[2]
    GEMMLOOPSET.loops[2] =
      Loop(:m, MaybeKnown(1), MaybeKnown(M), MaybeKnown(1), mloop.rangesym, mloop.lensym)
  end
  if K ≢ nothing
    kloop = GEMMLOOPSET.loops[3]
    GEMMLOOPSET.loops[3] =
      Loop(:k, MaybeKnown(1), MaybeKnown(K), MaybeKnown(1), kloop.rangesym, kloop.lensym)
  end
  GEMMLOOPSET.vector_width = W
  order = choose_order(GEMMLOOPSET)
  (N ≢ nothing) && (GEMMLOOPSET.loops[1] = nloop)
  (M ≢ nothing) && (GEMMLOOPSET.loops[2] = mloop)
  (K ≢ nothing) && (GEMMLOOPSET.loops[3] = kloop)
  order[5], last(order)
end
@generated function matmul_params(
  ::StaticInt{RS},
  ::StaticInt{RC},
  ::StaticInt{CLS},
) where {RS,RC,CLS}
  mᵣ, nᵣ = matmul_params(RS, RC, CLS)
  Expr(
    :tuple,
    Expr(:call, Expr(:curly, :StaticInt, mᵣ)),
    Expr(:call, Expr(:curly, :StaticInt, nᵣ)),
  )
end
matmul_params() = matmul_params(register_size(), register_count(), cache_linesize())

# function dotturbo(x,y)
#   s = zero(promote_type(eltype(x),eltype(y)))
#   @turbo for i ∈ eachindex(x,y)
#     s += x[i]*y[i]
#   end
#   s
# end
# dotturbo(Float64[1.0], Float64[2.0])


text:
push!(LOAD_PATH,"../")
using Documenter, GAlgebra

root = normpath(joinpath(dirname(pathof(GAlgebra)), ".."))
src = joinpath(root, "README.md")
dst = joinpath(root, "docs", "src",  "index.md")
cp(src, dst; force = true)

makedocs(
    # root = joinpath(dirname(pathof(GAlgebra)), ".."),
    # source = "./",
    # build = "./docs/build",
    modules = [GAlgebra],
    clean = false,
    sitename = "GAlgebra.jl",
    pages = Any[
        "Home" => "index.md",
        "API" => "api.md"
        # "Manual" => Any[
        #     "Plotting" => "man/plotting.md",
        #     "Compositing" => "man/compositing.md",
        #     "Backends" => "man/backends.md",
        #     "Themes" => "man/themes.md",
        # ]
    ]
)

deploydocs(
    repo   = "github.com/pygae/GAlgebra.jl.git"
)

ext: jl


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

Create a dictionary of results with string keys for each Scenario structure modeled.
"""
function reopt_results(m::JuMP.AbstractModel, p::REoptInputs; _n="")
	tstart = time()
    d = Dict{String, Any}()
    # TODO determine whether other results specific to electrical or thermal storage
    # systems warrant separate functions
    for b in p.s.storage.types.elec
        if p.s.storage.attr[b].max_kwh > 0
            add_electric_storage_results(m, p, d, b; _n)
        end
    end

    for b in p.s.storage.types.hot
        if p.s.storage.attr[b].max_kwh > 0
            add_hot_storage_results(m, p, d, b; _n)
        end
    end

    for b in p.s.storage.types.cold
        if p.s.storage.attr[b].max_kwh > 0
            add_cold_storage_results(m, p, d, b; _n)
        end
    end

    add_electric_tariff_results(m, p, d; _n)
    add_electric_utility_results(m, p, d; _n)
    add_financial_results(m, p, d; _n)
    add_electric_load_results(m, p, d; _n)

	if !isempty(p.techs.pv)
        add_pv_results(m, p, d; _n)
	end

    if "Wind" in p.techs.all
        add_wind_results(m, p, d; _n)
    end

    if "CHP" in p.techs.all
        add_chp_results(m, p, d; _n)
    end

	time_elapsed = time() - tstart
	@debug "Base results processing took $(round(time_elapsed, digits=3)) seconds."

	if !isempty(p.techs.gen) && isempty(_n)  # generators not included in multinode model
        tstart = time()
		add_generator_results(m, p, d)
        time_elapsed = time() - tstart
        @debug "Generator results processing took $(round(time_elapsed, digits=3)) seconds."
	end

	if !isempty(p.s.electric_utility.outage_durations) && isempty(_n)  # outages not included in multinode model
        tstart = time()
		add_outage_results(m, p, d)
        time_elapsed = time() - tstart
        @debug "Outage results processing took $(round(time_elapsed, digits=3)) seconds."
	end

    if !isempty(p.techs.boiler)
        add_existing_boiler_results(m, p, d)
    end

    if !isnothing(p.s.existing_chiller)
        add_existing_chiller_results(m, p, d)
    end

    if !isempty(p.techs.absorption_chiller)
        add_absorption_chiller_results(m, p, d)
    end

    if !isnothing(p.s.flexible_hvac)
        add_flexible_hvac_results(m, p, d)
    end

	return d
end

"""
    combine_results(bau::Dict, opt::Dict)

Combine two results dictionaries into one using BAU and optimal scenario results.
"""
function combine_results(p::REoptInputs, bau::Dict, opt::Dict, bau_scenario::BAUScenario)
    bau_outputs = (
        ("Financial", "lcc"),
        ("ElectricTariff", "year_one_energy_cost_before_tax"),
        ("ElectricTariff", "year_one_demand_cost_before_tax"),
        ("ElectricTariff", "year_one_fixed_cost_before_tax"),
        ("ElectricTariff", "year_one_min_charge_adder_before_tax"),
        ("ElectricTariff", "lifecycle_energy_cost_after_tax"),
        ("ElectricTariff", "lifecycle_demand_cost_after_tax"),
        ("ElectricTariff", "lifecycle_fixed_cost_after_tax"),
        ("ElectricTariff", "lifecycle_min_charge_adder_after_tax"),
        ("ElectricTariff", "lifecycle_export_benefit_after_tax"),
        ("ElectricTariff", "year_one_bill_before_tax"),
        ("ElectricTariff", "year_one_export_benefit_before_tax"),
        ("ElectricTariff", "year_one_coincident_peak_cost_before_tax"),
        ("ElectricTariff", "lifecycle_coincident_peak_cost_after_tax"),
        ("ElectricUtility", "year_one_to_load_series_kw"),
        ("ElectricUtility", "year_one_energy_supplied_kwh"),
        ("PV", "average_annual_energy_produced_kwh"),
        ("PV", "year_one_energy_produced_kwh"),
        ("PV", "lifecycle_om_cost_after_tax"),
        ("Generator", "average_annual_fuel_used_gal"),
        ("Generator", "lifecycle_fixed_om_cost_after_tax"),
        ("Generator", "lifecycle_variable_om_cost_after_tax"),
        ("Generator", "lifecycle_fuel_cost_after_tax"),
        ("Generator", "year_one_fuel_cost_before_tax"),
        ("Generator", "year_one_variable_om_cost_before_tax"),
        ("Generator", "year_one_fixed_om_cost_before_tax"),
        ("FlexibleHVAC", "temperatures_degC_node_by_time"),
        ("ExistingBoiler", "lifecycle_fuel_cost_after_tax" )
    )

    for t in bau_outputs
        if t[1] in keys(opt) && t[1] in keys(bau)
            if t[2] in keys(bau[t[1]])
                opt[t[1]][t[2] * "_bau"] = bau[t[1]][t[2]]
            end
        elseif t[1] == "PV" && !isempty(p.techs.pv)
            for pvname in p.techs.pv
                if pvname in keys(opt) && pvname in keys(bau)
                    if t[2] in keys(bau[pvname])
                        opt[pvname][t[2] * "_bau"] = bau[pvname][t[2]]
                    end
                end
            end
        end
    end
    opt["Financial"]["lifecycle_om_costs_before_tax_bau"] = bau["Financial"]["lifecycle_om_costs_after_tax"]
    opt["Financial"]["npv"] = round(opt["Financial"]["lcc_bau"] - opt["Financial"]["lcc"], digits=2)

    opt["ElectricLoad"]["bau_critical_load_met"] = bau_scenario.outage_outputs.bau_critical_load_met
    opt["ElectricLoad"]["bau_critical_load_met_time_steps"] = bau_scenario.outage_outputs.bau_critical_load_met_time_steps

    # TODO add FlexibleHVAC opex savings

    return opt
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
#=
    Target and port specifiction
=#
"""
    ports(ports...) ::Option

Define ports the scanner should scan for each host
"""
function ports(ports...) :: Option
    portlist = join(string.(collect(ports)), ",")
    return Option(
        function(scan::Scanner)
            append!(scan.args, ["-p", portlist])
        end)
end

"""
    target(targets...) ::Option

Define scanner targets
"""
function targets(targets...) :: Option
    return Option(
        function(s::Scanner)
            push!(s.args, string.(targets)...)
        end)
end

"""
    random_targets(::Int) ::Option

Amount of targets to randomly choose from the available targets
"""
function random_targets(amount::Int64) ::Option
    return Option("-iR ", string(amount))
end

"""
    exclude(target::String, others...) ::Option

Set targets (hosts, networks, etc) to exclude from the scan
At least one target must be provided
"""
function exclude(target::String, others...) ::Option
    excluded = string(target, join(string.(collect(others)), ","))
    return Option(
        function(s::Scanner)
            push!(s.args, "--exclude")
            push!(s.args, excluded)
        end
    )
end

"""
    exclude_file(filepath::String)

Set input files that define target exclusions

See also: `exclude`
"""
function exclude_file(filepath::String) ::Option
    return Option(
        function(s::Scanner)
            push!(s.args, "--excludefile")
            push!(s.args, filepath)
        end
    )
end

"""
    input_file(filepath::String) ::Option

Set the input file that define scan targets.
"""
function input_file(filepath::String) ::Option
    return Option("--iL", filepath)
end

#= end Target specification =#

function exclude_ports(ports...) ::Option

end

"""
    fastmode()

Faster scan that will scan fewer ports than the default scan
"""
function fastmode()::Option return Option("-F") end

function consecutive_port_scan()::Option  return Option("-r") end

function top_ports(number::Int64)::Option return Option("--top-ports", string(number)) end

function port_ratio(ratio::Float64) ::Option
    @assert ratio >= 1 and ratio >= 0 "port_ratio should be a value between 0 and 1"
    return Option("--port-ratio", @sprintf("%.1f", ratio))
end

#= end Port specification =#

#=
    Host Discovery
=#
"""
    listscan

-sL: List Scan - simply list targets to scan
"""
function listscan()::Option return Option("-sL") end
"""
    pingscan

-sn: Ping Scan - disable port scan
"""
function pingscan()::Option return Option("-sn") end
"""
    skip_discovery()

-Pn: Treat all hosts as online -- skip host discovery
"""
function skip_discovery()::Option return Option("-Pn") end

"""
    syn_discovery(ports...) :: Option

Set discovery mode to use SYN packets

If no ports are provided, SYN discovery will be enabled for all ports

Alternatively, use `syn(ports...)`
"""
function syn_discovery(ports...) ::Option
    return Option(string("-PS", join(string.(collect(ports)), ",")))
end
syn(ports...)::Option = syn_discovery()

"""
    ack_discovery(ports...)

Set discovery mode to use ACK packets

If no ports are provided, ACK discovery will be enabled for all ports

Alternatively, use `ack(ports...)`
"""
function ack_discovery(ports...) ::Option
    return Option(
        function(s::Scanner)
            push!(s.args, "-PA")
            push!(s.args, join(string.(collect(ports)), ","))
        end
    )
end
ack(ports...)::Option = ack_discovery(ports...)

function udp_discovery(ports...) ::Option
    return Option(
        function(s::Scanner)
            push!(s.args, "-PU")
            push!(s.args, join(string.(collect(ports)), ","))
        end
    )
end
udp(ports...) = udp_discovery(ports...)

function sctp_discovery(ports...) ::Option
    return Option(
        function(s::Scanner)
            push!(s.args, "-PS")
            push!(s.args, join(string.(collect(ports)), ","))
        end
    )
end
sctp(ports...)::Option = sctp_discovery(ports...)

"""
    discovery(modes...) ::Option

Utility option for grouping discovery mode definitions

Examples
```julia

scanner = Scanner(
    discovery(
        syn("22", "80"), # use syn packets for ports 22 and 80
        udp("9000"),     # use udp packets for port 9000
        ack()            # use ack packets in all ports
    )
)
```
"""
function discovery(modes...)
    @assert isempty(modes) || all(x->x isa Option, modes) "Provided arguments must be of type `Option`"
    return Option(
        function(scan::Scanner)
            [mode(scan) for mode in modes]
        end
    )
end

function icmp_echo_discovery()

end

function icmp_timestamp_discovery()

end

function icmp_netmask_discovery()

end

function ip_ping_discovery(protocols...)

end

function disabled_dns_resolution()
end

function forced_dns_resolution()
end

function dns_servers()
end

function system_dns()

end

"""
    traceroute()

Enable hop path tracing for each host
"""
function traceroute()::Option return Option("--traceroute") end

#= end Host Discovery =#

#=
    Scan Techniques
=#

function syn_scan()::Option     return Option("-sS") end

function connect_scan()::Option return Option("-sT") end

function ack_scan()::Option     return Option("-sA") end

function window_scan()::Option  return Option("-sW") end

function maimon_scan()::Option  return Option("-sM") end
"""
    updscan() ::Option

Enable scan technique to use UDP packets
"""
function udp_scan()::Option return Option("-sU") end

"""
    tcp_null_scan() ::Option

Set scan technique to use TCP null packets - TCP flag header = 0

If a `RST` packet is received, port is considered closed.

If no response is received, the port is considered open|filtered
"""
function tcp_null_scan()::Option return Option("-sN") end

function tcp_fin_scan()::Option  return Option("-sF") end

"""
    xmas_scan() ::Option

Set the scan technique to use TCP packets with FIN, PSH and URF flags set.

If a `RST` packet is received, port is considered closed.

If no response is received, the port is considered open|filtered
"""
function xmas_scan()::Option return Option("-sX") end

function ip_protocol_scan() ::Option
    return Option("-sO")
end

"""
    ftp_bounce_scan(ftp_relay::String) ::Option

Enables FTP relay host scan technique
"""
function ftp_bounce_scan(ftp_relay::String) return Option("-b", ftp_relay) end

@enum TCPFlag begin
    NULL = 0
    FIN  = 1
    SYN  = 2
    RST  = 4
    PSH  = 8
    ACK  = 16
    URG  = 32
    ECE  = 64
    CWR  = 128
    NS   = 256
end

function tcp_scan_flags(flags...) ::Option
    @assert all(x->x isa Int64 || x isa TCPFlag) "TCP Flags must be of type Int or TCPFlag"
    value = sum([Int(flag) for flag in flags])
    # format as lowercase hexadecimal
    return Option("--scanflags", @sprintf("%x", value))
end

function idle_scan(zombie::String, port::Int64=0) ::Option
    return Option(
        function(s::Scanner)
            push!(s.args, "-sI")
            if port > 0
                # format as <string>:<int>
                push!(s.args, @sprintf("%s:%d", zombie, port))
            else
                push!(s.args, @sprintf("%s", zombie))
            end
        end
    )
end

function sctp_init_scan()::Option return Option("-sY") end

function sctp_cookie_echo_scan()::Option return Option("-sZ") end

#= end Scan Techniques =#

#= Service/Version detection =#

"""
    service_info() ::Option

Probe open ports to determine service/version info
"""
function service_info(options...) ::Option
    return Option(
        function(scanner::Scanner)
            push!(scanner.args, "-sV")
            [option(scanner) for option in options]
        end
    )
end

"""
    version_intensity(level) ::Option

Set the level of intensity with which nmap should probe the open ports to get version information.
Intensity should be a value between 0 (light) and 9 (try all probes).

Default value is 7
"""
function version_intensity(level::Int) ::Option
    @assert level >= 0 && level <= 9 "--version-intensity accepts values from 0 (light) to 9 (try all probes)"
    return Option("--version-intensity", string(level))
end
version_intensity()::Option = version_intensity(7)
version_light()::Option     = version_intensity(2)
version_all()::Option       = version_intensity(9)

function version_trace() ::Option
    return Option("--version-trace")
end

#= end Service/Version detection =#

#= Script scan =#

function default_script() ::Option
    return Option("-sC")
end

"""
    script_args(args...; kwargs...) ::Option

For each argument, if a value is not provided or is empty, the key is used as argument

Arguments

 * args     :: Union{Pair, String}

Example
```julia
opt = NMAP.script_args("script-key.property"=>"localhost", "script.showall", user="username", pass=",{}=bar", whois="{whodb=nofollow+ripe}")

scanner(opt)
# > "--script-args=script-key.property=localhost,script.showall,user=username,pass=,{}=bar,whois={whodb=nofollow+ripe}"
```
"""
function script_args(args...; kwargs...) :: Option
    @assert all(x->x isa Pair || x isa String, args) """Provide positional arguments as Pairs or Strings: `script_args("some-script.property"=>"some value", "someKey.property"; user=foo)`"""

    args_list = []
    for pair in args
        k,v = pair isa Pair ? pair : (pair, "")
        push!(args_list, (k,v))
    end
    for pair in kwargs
        k,v = pair
        push!(args_list, (string(k), v))
    end
    scriptargs = []
    for tuple in args_list
        k,v = tuple
        k = strip(k)
        str = isnothing(v) || isempty(v) ? k : @sprintf("%s=%s", k, v)
        push!(scriptargs, str)
    end
    scriptargs = join(string.(strip.(scriptargs)), ",")
    return Option(
        function(scanner::Scanner)
            args = strip(scriptargs)
            push!(scanner.args, @sprintf("--script-args=%s", args))
        end
    )
end

function script(scripts...) ::Option
    @assert all(x->x isa String, scripts) """Provide script names as strings: `script("default")`"""
    list = join(string.(collect(scripts)), ",")

    return Option(
        function(scanner::Scanner)
            push!(scanner.args, @sprintf("--script=%s", list))
        end
    )
end

#= end Script scan =#

"""
    scandelay(delay::Int)

Sets the minimum time to wait between each probe sent to a host
"""
function scandelay(delay::Int) ::Option
end

#= IDS evasion, spoofind =#

"""
    decoy(decoys... String) ::Option

Cloak a scan with decoys

Performs a decoy scan which makes it seem to the remote host that your specified decoys
are also scanning the target network.

Argument `ME` can be used as one of the given `decoys` to represent the position of your real IP Address.
"""
function decoy(decoys...) ::Option
    @assert !isempty(decoys) && all(x->x isa String, decoys) "Provide decoy addresses or ME"
    return Option("-D", join(string.(collect(decoys)), ","))
end

"""
    spoof_source(ip::String)

Spoof the IP address of the machine running nmap with given `address`
"""
function spoof_source(ip::String) ::Option
    return Option("-S", ip)
end

"""
    spoof_mac(spoof::String)

Spoof MAC Address
Uses the given MAC address for all of the raw ethernet frames the scanner sends
"""
function spoof_mac(spoof::String) ::Option
    return Option("--spoof-mac", spoof)
end

"""
    interface(interface::String)

Set network interface used for scanning
"""
function interface(interface::String) ::Option
    return Option("-e", interface)
end

function badsum() ::Option
   return Option("--badsum")
end
#= end IDS =#

#= OS Detection =#

"""
    os_detection(options...) :: Option

Example

```julia
    Scan(
        os_detection(
            osscan_guess(),
            osscan_limit()
        )
    )
```
"""
function os_detection(options...) ::Option
    @assert isempty(options) || all(x->x isa Option, options) "Provided arguments must be of type `Option`"
    return Option(
        function(scan::Scanner)
            push!(scan.args, "-O")
            for option in options
                option(scan)
            end
        end
    )
end
function os(options...)::Option return os_detection(options...) end

function osscan_guess() ::Option
    return Option((x)->push!(x.args, "--osscan-guess"))
end

function osscan_limit() ::Option
    return Option((x)->push!(x.args, "--osscan-limit"))
end

#= end OS Detection =#

#= Output options =#

abstract type OutputFormat end
abstract type GreppableFormat   <:OutputFormat end
abstract type kIddi3Format      <:OutputFormat end
abstract type NormalFormat      <:OutputFormat end
abstract type XMLFormat         <:OutputFormat end

Base.string(::Type{XMLFormat})          = "-oX"
Base.string(::Type{NormalFormat})       = "-oN"
Base.string(::Type{kIddi3Format})       = "-oS"
Base.string(::Type{GreppableFormat})    = "-oG"

function output(::Type{O}, file::String="-")::Option where O<:OutputFormat
    return output(string(O), file)
end
function output(::Type{X}) where X<:XMLFormat
end
function output(format::String, file::String="-") ::Option
    return Option(
        function(s::Scanner)
            if !(format == string(XMLFormat))
                s.xml = false
            end
            push!(s.args, format, file)
        end
    )
end
function Base.in(::Type{O}, scanner::Scanner)::Option where O<:OutputFormat
    return string(O) in scanner.args
end

function only_open()::Option    return Option("--open") end
function packet_trace()::Option return Option("--packet-trace") end
function iflist() ::Option
    return Option(
        function(scanner::Scanner)
            scanner.debug = true
            push!(scanner.args, "--iflist")
        end
    )
end

function resume(logfile::String) ::Option
    return Option("--resume", logfile)
end

function reason() ::Option Option("--reason") end

function verbosity(level::Int=1) ::Option
    @assert level in range(1, stop=10) "Verbosity level should be an integer between 1 and 10"
    return Option(string("-", join(["v" for _ in 1:level])))
end
function verbose(level::Int=1)::Option verbosity(level)  end
function vv()::Option verbosity(2) end

function debugging(level::Int) ::Option
    @assert level in range(1, stop=10) "Debugging level should be an integer between 1 and 10"
    return Option(string("-", join(["d" for _ in 1:level])))
end
function debug(level::Int=1)::Option debugging(level)  end
function dd()::Option debugging(2) end

#= end Output options =#

#= Miscellaneous options =#

function ipv6_scanning()    return Option("-6") end
function aggressive_scan()  return Option("-A") end
function privileged()       return Option("--privileged") end
function unprivileged()     return Option("--unprivileged") end

#= end Miscellaneous options =#

#= Timing templates     =#
include("timing.jl")
#= end Timing templates =#


text:
using DataFrames
using CategoricalArrays
using Random

import Base: hash, size, getindex, length, show

struct ModalInstance <: AbstractVector{SubDataFrame}
    # worlds::Set{<:AbstractOntology} # TODO add set of worlds, for the semantics (DT)
    rows::Vector{DataFrameRow}
end

ModalInstance() = ModalInstance(DataFrameRow[])
DataFrames.nrow(::Array{DataFrameRow,1}) = 1

function show(io::IO, ::MIME"text/plain", mi::ModalInstance)
    println(io, "Modal instance with $(length(mi.rows)) frames")
    for i in 1:length(mi.rows)
        println(io, "Frame #$i")
        println(io, mi.rows[i])
    end
end

size(mi::ModalInstance) = (nrow(mi.rows),)
getindex(mi::ModalInstance, i::Int) = mi.rows[i]
length(mi::ModalInstance, fid::Int) = max([length(mi[fid][i]) for i in 1:length(mi[fid])]...)

mutable struct ModalFrame <: AbstractVector{SubDataFrame}
    # name::String # TODO add name?
    dimension::Int  # (0=static, 1=timeseries, 2=images, ecc.)
    data::DataFrame
    views::Vector{DataFrameRow} # intraframe instances (i.e., first, second, .., series)
end

dimension(mf::ModalFrame) = mf.dimension

@inline function ModalFrame(data::DataFrame)
    colnames = names(data)

    @assert length(colnames) > 0 && nrow(data) > 0 "There must be at least one column and one instance in data to compute its dimension."

    dimension = ndims(data[1, Symbol(colnames[1])])

    for colname in colnames
        curr_dimension = ndims(data[1, Symbol(colname)])
        @assert dimension == curr_dimension "Each attribute must have the same dimension."
    end

    # if dimension > 0
    #     typ = eltype(data[1, Symbol(colnams[1])])
    #     for colname in colnames
    #         curr_typ = eltype(data[1, Symbol(colname)])
    #         @assert typ == cur_typ "Each modal frame with higher-than-zero dimension must have the same element types."
    #     end
    # end

    views = Vector{DataFrameRow}(undef, nrow(data))
    for i in 1:nrow(data)
        views[i] = view(data, i, :)
    end
    ModalFrame(dimension, data, views)
end

size(mf::ModalFrame) = (nrow(mf.data),)
getindex(mf::ModalFrame, i::Int) = mf.views[i]

struct ClassificationDataset
    ldim::Int
    hdim::Int
    frames::Vector{ModalFrame}
    instances::Vector{ModalInstance}
    domains::Dict{Int, Dict{Symbol, Vector{Number}}}
    classes::CategoricalArray
    unique_classes::Vector
end

function ClassificationDataset(frames::Vector{ModalFrame}, classes::CategoricalArray)
    instances = ModalInstance[]
    domains = Dict{Int, Dict{Symbol, Vector{Number}}}()

    n_instances = size(frames[1])

    for i in 1:length(frames)
        curr_n_instances = size(frames[i])
        @assert n_instances == curr_n_instances "Each frame must have the same number of (intraframe) instances!"
    end

    @assert n_instances[1] == length(classes) "Each instance must have one and only one class!"

    # Computing the instances.
    # for each instance
    for i in 1:n_instances[1]
        instance = ModalInstance()
        # for each frame
        for f in 1:length(frames)
            # push the instance
            push!(instance.rows, frames[f][i])
        end
        push!(instances, instance)
    end

    # Computing the domains.
    # for each frame
    for f in 1:length(frames)
        attrs = names(frames[f][1])
        # for each attribute
        attr_domain = Dict{Symbol, Vector{Number}}()
        for attr in attrs
            # convert to symbol
            A = Symbol(attr)
            # for each instance
            # for i in 1:n_instances[1]

            # end
            attr_domain[A] = collect(Set([v for i in 1:n_instances[1] for v in frames[f][i][A]]))
            domains[f] = attr_domain
        end
    end

    ldim, hdim = extrema([dimension(i) for i in frames])

    ClassificationDataset(ldim, hdim, frames, instances, domains, classes, unique(classes))
end

instance(cs::ClassificationDataset, i::Int) = cs.instances[i] # TODO cs -> ds (either sup or unsupervised)
instance(cs::ClassificationDataset, i::Int, f::Int) = cs.instances[i][f]
function attributes(ds::ClassificationDataset)
    d = Dict{Int,Array{String,1}}()
    for (fid,frame) in enumerate(ds.frames)
        d[fid] = names(frame.data)
    end
    return d
end
attributes(ds::ClassificationDataset, fid::Int) = names(ds.frames[fid].data)

length(ds::ClassificationDataset) = length(ds.instances)

function show(io::IO, ::MIME"text/plain", ds::ClassificationDataset)
    println(io, "Classification dataset with $(length(ds.instances)) instances")

    for (i, inst) in enumerate(ds.instances)
        println(io, "Instance #$(i):")
        println(io, instance(ds, i))
        println(io)
    end
end

# TODO
# possible test: all([auslan.instances[i].rows[1][attr] === auslan.frames[1].data[i,attr] for i in 1:length(auslan.instances) for attr in attributes(auslan, 1)])
function transform!(ds::ClassificationDataset, f::Function, fid::Int; kwargs...)
    for attr in attributes(ds, fid)
        for i in 1:nrow(ds.frames[fid].data)
            ds.frames[fid].data[i,attr] = f(ds.frames[fid].data[i,attr]; kwargs...)
        end
    end
end

# function Base.minimum(ds::ClassificationDataset)
#     d = Dict{Int,Array{Float64,1}}()

#     for (fid,frame) in enumerate(ds.frames)
#         frame = ds.frames[fid]
#         if frame.dimension == 0
#             d[fid] = missing
#         else
#             d[fid] = [Base.minimum(vcat(frame.data[!,attr]...)) for attr in attributes(ds, fid)]
#         end
#     end
#     return d
# end

# function Base.maximum(ds::ClassificationDataset)
#     d = Dict{Int,Array{Float64,1}}()

#     for (fid,frame) in enumerate(ds.frames)
#         frame = ds.frames[fid]
#         if frame.dimension == 0
#             d[fid] = missing
#         else
#             d[fid] = [Base.maximum(vcat(frame.data[!,attr]...)) for attr in attributes(ds, fid)]
#         end
#     end
#     return d
# end

# TODO fix problem with select_attributes!
# TODO make it with bang!
function select_attributes(ds::ClassificationDataset, frame_attributes::Dict{Int, Array{Int, 1}})
    frames = collect(keys(frame_attributes))

    # Frame level.
    for i in 1:length(frames)
        ds.frames[frames[i]].data = ds.frames[frames[i]].data[:, frame_attributes[frames[i]]]
    end

    modal_frames = ModalFrame[]
    for i in 1:length(frames)
        push!(modal_frames, ModalFrame(ds.frames[frames[i]].data))
    end

    return ClassificationDataset(modal_frames, ds.classes)
end

# function transform!(ds::ClassificationDataset, f::Function; kwargs...)

# using MultivariateTimeSeries
# X, y = read_data_labeled(joinpath(dirname(pathof(GBDTs)), "..", "data", "auslan_youtube8"));
#
# # adesso voglio passare da MTS a ClassificationDataset
#
# my_frame = DataFrame()
#
# for attr in names(X)
#     insertcols!(my_frame, attr => Array{Float64,1}[])
# end
#
# # for each instance
# for i in 1:length(X)
#     array_data = Array{Float64, 2}(X[i])
#     push!(my_frame, [array_data[:, j] for j in 1:length(names(X))])
# end
#
# auslan = ClassificationDataset([ModalFrame(my_frame)], CategoricalArray(y))

function sample(ds::ClassificationDataset,
        range::Union{UnitRange{Int},Nothing};
        invert::Bool=false)

    n_instances = length(ds.instances)

    @assert range ⊆ 1:n_instances "The range must be a subset of 1:$(n_instances)."

    if ismissing(range)
        return ds
    else

        output_classes = Int[]
        timeseries_df = DataFrame()

        attrs = attributes(ds, 1) # TODO up to now, only timeseries

        for attr in attrs
            insertcols!(timeseries_df, attr => Array{Float64, 1}[]) # TODO also here Array{Float64, 1}[], i.e., timeseries
        end

        if invert == false
            for idx in range
                push!(timeseries_df, ds.instances[idx][1]) # ds.instances[idx][1] <-- [1] means the first frame TODO in future
                push!(output_classes, ds.classes[idx])
            end
        else
            for idx in setdiff(1:n_instances,range)
                push!(timeseries_df, ds.instances[idx][1]) # ds.instances[idx][1] <-- [1] means the first frame TODO in future
                push!(output_classes, ds.classes[idx])
            end
        end
    end

    return ClassificationDataset([ModalFrame(timeseries_df)], CategoricalArray(output_classes))
end

function sample(ds::ClassificationDataset,
    percentage::Float64;
    seed::Int=1,
    bias::Float64=0.0, # bias to uniform
    replace::Bool=true,
    invert::Bool=false)

    @assert 0 ≤ percentage && percentage ≤ 1 "Percentage p must be 0 ≤ p ≤ 1."

    n_instances = length(ds.instances)
    n_instances_per_class = countmap(ds.classes)
    classes = keys(n_instances_per_class)
    n_classes = length(n_instances_per_class)

    instances_per_class = Dict{CategoricalArrays.CategoricalValue{Int64, UInt32}, Array{Int, 1}}()

    for class in classes
        instances_per_class[class] = ModalInstance[]
    end

    instances_to_sample = Dict{CategoricalArrays.CategoricalValue{Int64, UInt32}, Int64}()

    for class in classes
        sample_size = Int(round(percentage * ((1-bias)*n_instances_per_class[class] + bias*n_instances / n_classes)))
        instances_to_sample[class] = sample_size
    end

    for i in 1:n_instances
        push!(instances_per_class[ds.classes[i]], i)
    end

    r = Random.seed!(seed)

    output_classes = Int[]
    timeseries_df = DataFrame()

    attrs = attributes(ds, 1) # TODO up to now, only timeseries

    for attr in attrs
        insertcols!(timeseries_df, attr => Array{Float64, 1}[]) # TODO also here Array{Float64, 1}[], i.e., timeseries
    end

    if replace
        for class in classes
            for i in 1:instances_to_sample[class]
                idx = instances_per_class[class][rand(r, 1:n_instances_per_class[class])]
                push!(timeseries_df, ds.instances[idx][1]) # ds.instances[idx][1] <-- [1] means the first frame TODO in future
                push!(output_classes, ds.classes[idx])
                # @show idx
                # @show ds.classes[idx]
            end
        end
    else
        for class in classes
            n_eligible = length(instances_per_class[class])
            selected = instances_per_class[class]

            # the selected elements are those at the end of then selected array
            for i in 1:instances_to_sample[class]
                idx = rand(r, 1:n_eligible)
                selected[idx], selected[n_eligible] = selected[n_eligible], selected[idx]
                n_eligible -= 1
            end

            # take the last instances_to_sample[class] instances
            if invert == false
                for i in n_instances_per_class[class]-instances_to_sample[class]+1:n_instances_per_class[class]
                    idx = selected[i]
                    push!(timeseries_df, ds.instances[idx][1]) # ds.instances[idx][1] <-- [1] means the first frame TODO in future
                    push!(output_classes, ds.classes[idx])
                end
            # take the first n_instances_per_class[class]-instances_to_sample[class] instances
            else
                for i in 1:n_instances_per_class[class]-instances_to_sample[class]
                    idx = selected[i]
                    push!(timeseries_df, ds.instances[idx][1]) # ds.instances[idx][1] <-- [1] means the first frame TODO in future
                    push!(output_classes, ds.classes[idx])
                end
            end
        end
    end

    return ClassificationDataset([ModalFrame(timeseries_df)], CategoricalArray(output_classes))
end

ds = SymbolicML.trip_no_trip("/home/edu/.julia/dev/SymbolicML/data/Siemens/Features");
# ds = SymbolicML.trip_no_trip("/home/edu/.julia/dev/SymbolicML/data/Siemens"; n_machine=42740);
# d = Dict{Int,Array{Int,1}}()
# d[1] = [1,2,4,5,9,19,22,28,29,30]
# SymbolicML.select_attributes!(ds, d)
SymbolicML.transform!(ds, SymbolicML.window_slice, 1; from=1,to=120)

grammar = @grammar begin
    b = G(bvec) | F(bvec) | G(implies(bvec,bvec))
    bvec = and(bvec, bvec)
    bvec = or(bvec, bvec)
    bvec = not(bvec)
    bvec = lt(rvec, rvec)
    bvec = lte(rvec, rvec)
    bvec = gt(rvec, rvec)
    bvec = gte(rvec, rvec)
    bvec = f_lt(x, xid, v, vid)
    bvec = f_lte(x, xid, v, vid)
    bvec = f_gt(x, xid, v, vid)
    bvec = f_gte(x, xid, v, vid)
    rvec = x[xid]
    # xid = |(["x_1","y_1","z_1","roll_1","pitch_1","yaw_1","thumbbend_1","forebend_1","middlebend_1","ringbend_1","littlebend_1","x_2","y_2","z_2","roll_2","pitch_2","yaw_2","thumbbend_2","forebend_2","middlebend_2","ringbend_2","littlebend_2"])
    xid = |("Fce", "Fsfc", "G2", "MAD", "T3-2", "T13-12", "T16-15", "T13-5", "T14-6", "T15-7")
    # xid = |("Ambient_air_humidity","Compr_IGV_position","Gas_fuel_valve_position","Ambient_air_temperature","Compressor_outlet_temperature","Exhaust_temperature_1","Exhaust_temperature_2","Exhaust_temperature_3","Exhaust_temperature_4","Exhaust_temperature_5","Exhaust_temperature_6","Exhaust_temperature_7","Exhaust_temperature_8","Exhaust_temperature_9","Exhaust_temperature_10","Exhaust_temperature_11","Exhaust_temperature_12","Exhaust_temperature_13","Exhaust_temperature_14","Exhaust_temperature_15","Exhaust_temperature_16","Compressor_outlet_pressure","Gas_fuel_mass_flow_rate","Compressor_inlet_air_mass_flow_rate","Rotational_speed","Power_output")
    # xid = |("Ambient_air_humidity","Compr_IGV_position","Gas_fuel_valve_position","Ambient_air_temperature","Compressor_outlet_temperature","Compressor_outlet_pressure","Gas_fuel_mass_flow_rate","Compressor_inlet_air_mass_flow_rate","Rotational_speed","Power_output")
    vid = |(1:10)
end

G(v) = all(v)                                                #globally
F(v) = any(v)                                                #eventually
f_lt(x, xid, v, vid) = lt(x[xid], v[xid][vid])               #feature is less than a constant
f_lte(x, xid, v, vid) = lte(x[xid], v[xid][vid])             #feature is less than or equal to a constant
f_gt(x, xid, v, vid) = gt(x[xid], v[xid][vid])               #feature is greater than a constant
f_gte(x, xid, v, vid) = gte(x[xid], v[xid][vid])             #feature is greater than or equal to a constant

#workarounds for slow dot operators:
implies(v1, v2) = (a = similar(v1); a .= v2 .| .!v1)         #implies
not(v) = (a = similar(v); a .= .!v)                          #not
and(v1, v2) = (a = similar(v1); a .= v1 .& v2)               #and
or(v1, v2) = (a = similar(v1); a .= v1 .| v2)                #or
lt(x1, x2) = (a = Vector{Bool}(undef,length(x1)); a .= x1 .< x2)   #less than
lte(x1, x2) = (a = Vector{Bool}(undef,length(x1)); a .= x1 .≤ x2)  #less than or equal to
gt(x1, x2) = (a = Vector{Bool}(undef,length(x1)); a .= x1 .> x2)   #greater than
gte(x1, x2) = (a = Vector{Bool}(undef,length(x1)); a .= x1 .≥ x2)  #greater than or equal to

balanced_ds = SymbolicML.sample(ds, 0.727; seed=1, bias=1.0, replace=false);
train_ds = SymbolicML.sample(balanced_ds, 0.8, seed=1, bias=1.0, replace=false);
test_ds = SymbolicML.sample(balanced_ds, 0.8, seed=1, bias=1.0, replace=false, invert=true);

# all_instances = union(train_ds.instances, test_ds.instances)

# A = for (i, inst) in enumerate(all_instances)
#     A = Bool[]
#     @show i
#     if inst in balanced_ds.instances
#         push!(A, true)
#     else
#         push!(A, false)
#     end
#     A
# end

const v = Dict{String,Vector{Float64}}()
mins, maxes = minimum(train_ds)[1], maximum(train_ds)[1]
for (i,xid) in enumerate(attributes(train_ds,1))#enumerate(names(X))
    v[xid] = collect(range(mins[i],stop=maxes[i],length=10))
end;

p = MonteCarlo(1000, 5)
# genetic = GeneticProgram(100,20,5,0.5,0.4,0.1)
using Random; Random.seed!(1)
model = induce_tree(grammar, :b, p, train_ds.frames[1], train_ds.classes, 5, min_members_per_branch=1, prevent_same_label=true);

train_ind = collect(1:length(train_ds.instances));
train_y_pred = classify(model, train_ds.frames[1], train_ind);
train_accuracy = count(train_y_pred .== train_ds.classes[train_ind]) / length(train_ind)

test_ind = collect(1:length(test_ds.instances));
test_y_pred = classify(model, test_ds.frames[1], test_ind);
test_accuracy = count(test_y_pred .== test_ds.classes[test_ind]) / length(test_ind)

# TODO now it works for the specific problem
function confusion_matrix(ŷ, y)
    @assert length(ŷ) == length(y) "The number of predictions must be equal to the number of true classes."
    # @show y
    # @show length(unique(y))
    # mat = zeros(Int, length(unique(y)), length(unique(y)))
    mat = zeros(Int, 2, 2)

    for i in 1:length(y)
        mat[Int(ŷ[i].ref), Int(y.refs[i])] += 1
    end

    return mat
end

function performance_metrics(M)
    return tn, fn, fp, tp = M[1,1], M[1,2], M[2,1], M[2,2]
end

function sensitivity(M)
    tn, fn, fp, tp = performance_metrics(M)
    return tp/(tp+fn)
end

function specificity(M)
    tn, fn, fp, tp = performance_metrics(M)
    return tn/(tn+fp)
end

@show train_accuracy
@show test_accuracy

ds = SymbolicML.trip_no_trip("/home/edu/.julia/dev/SymbolicML/data/Siemens"; n_machine=42740);

const v = Dict{String,Vector{Float64}}()

for seed in 1:10
    for window_sz in 1:2
        for hour in 24:-1:19
            start_t = 1440 - hour*60 + 1
            end_t = 1440 - hour*60 + (window_sz*60)

            balanced_ds = SymbolicML.sample(ds; seed=seed, percentage=0.727, bias=1.0, replace=false);
            SymbolicML.transform!(balanced_ds, SymbolicML.window_slice, 1; from=start_t,to=end_t)
            # @show countmap(balanced_ds.classes)
            train_ds = SymbolicML.sample(balanced_ds, seed=seed, percentage=0.8, bias=1.0, replace=false);
            # @show countmap(train_ds.classes)
            test_ds = SymbolicML.sample(balanced_ds, seed=seed, percentage=0.8, bias=1.0, replace=false, invert=true);
            # @show countmap(test_ds.classes)

            # local v = Dict{String,Vector{Float64}}()
            mins, maxes = minimum(train_ds)[1], maximum(train_ds)[1]
            for (i,xid) in enumerate(attributes(train_ds,1))#enumerate(names(X))
                v[xid] = collect(range(mins[i],stop=maxes[i],length=10))
            end

            p = MonteCarlo(500, 5)

            Random.seed!(seed)

            model = induce_tree(grammar, :b, p, train_ds.frames[1], train_ds.classes, 7, min_members_per_branch=1);

            # train_ind = collect(1:length(train_ds.instances));
            # train_y_pred = classify(model, train_ds.frames[1], train_ind);
            # train_accuracy = count(train_y_pred .== train_ds.classes[train_ind]) / length(train_ind)

            test_ind = collect(1:length(test_ds.instances));
            test_y_pred = classify(model, test_ds.frames[1], test_ind);
            test_accuracy = round(count(test_y_pred .== test_ds.classes[test_ind]) / length(test_ind), digits=4)

            test_M = confusion_matrix(test_y_pred, test_ds.classes[test_ind])
            test_sens = round(sensitivity(test_M),digits=4)
            test_spec = round(specificity(test_M),digits=4)

            @show (seed,window_sz,hour,test_accuracy,test_sens,test_spec)
        end
    end
end

for seed in 1:10
    # @show seed
    for window_sz in 1:2
        # @show window_sz
        for hour in 24:-1:19
            # @show hour
            start_t = 1440 - hour*60 + 1
            end_t = 1440 - hour*60 + (window_sz*60)

            balanced_ds = SymbolicML.sample(ds, 0.727; seed=seed, bias=1.0, replace=false);
            SymbolicML.transform!(balanced_ds, SymbolicML.window_slice, 1; from=start_t,to=end_t)
            test_M = zeros(Int, 2, 2)
            for r in 1:length(balanced_ds.instances)
                # @show r
                train_ds = SymbolicML.sample(balanced_ds, r:r, invert=true);
                test_ds = SymbolicML.sample(balanced_ds, r:r, invert=false);

                mins, maxes = minimum(train_ds)[1], maximum(train_ds)[1]
                for (i,xid) in enumerate(attributes(train_ds,1))#enumerate(names(X))
                    v[xid] = collect(range(mins[i],stop=maxes[i],length=10))
                end

                p = MonteCarlo(500, 5)

                Random.seed!(seed)

                model = induce_tree(grammar, :b, p, train_ds.frames[1], train_ds.classes, 7, min_members_per_branch=1);

                test_ind = collect(1:length(test_ds.instances));
                test_y_pred = classify(model, test_ds.frames[1], test_ind);
                # test_accuracy = round(count(test_y_pred .== test_ds.classes[test_ind]) / length(test_ind), digits=4)

                # @show test_y_pred
                # @show convert(Int, test_y_pred[1])
                # # @show test_ds.classes[r]
                # @show test_ds.classes[test_ind][1]
                # @show Int(test_y_pred[1].ref)
                # @show test_ds.classes[test_ind]
                # @show test_ds.classes[r]
                # tn
                if convert(Int, test_y_pred[1]) == 0 && convert(Int, test_ds.classes[test_ind][1]) == 0
                    test_M[1,1] += 1
                # tp
                elseif convert(Int, test_y_pred[1]) == 1 && convert(Int, test_ds.classes[test_ind][1]) == 1
                    test_M[2,2] +=1
                # fn
                elseif convert(Int, test_y_pred[1]) == 0 && convert(Int, test_ds.classes[test_ind][1]) == 1
                    test_M[1,2] += 1
                # fp
                elseif convert(Int, test_y_pred[1]) == 1 && convert(Int, test_ds.classes[test_ind][1]) == 0
                    test_M[2,1] += 1
                end

                # @show test_M
            end

            # @show test_M

            test_accuracy = (test_M[1,1] + test_M[2,2])/length(balanced_ds.instances)
            test_sens = round(sensitivity(test_M),digits=4)
            test_spec = round(specificity(test_M),digits=4)
            @show (seed,window_sz,hour,test_accuracy,test_sens,test_spec)
        end
    end
end


text:
@test_throws DimensionMismatch shuffleobs((X, rand(149)))

@testset "typestability" begin
    for var in vars
        @test typeof(@inferred(shuffleobs(var))) <: SubArray
    end
    for tup in tuples
        @test typeof(@inferred(shuffleobs(tup))) <: Tuple
    end
end

@testset "Array and SubArray" begin
    for var in vars
        @test size(shuffleobs(var)) == size(var)
    end
    # tests if all obs are still present and none duplicated
    @test sum(shuffleobs(Y1)) == 120
end

@testset "Tuple of Array and SubArray" begin
    for var in ((X,yv), (Xv,y), tuples...)
        @test_throws MethodError shuffleobs(var...)
        @test typeof(shuffleobs(var)) <: Tuple
        @test all(map(x->(typeof(x)<:SubArray), shuffleobs(var)))
        @test all(map(x->(numobs(x)===15), shuffleobs(var)))
    end
    # tests if all obs are still present and none duplicated
    # also tests that both paramter are shuffled identically
    x1, y1, z1 = shuffleobs((X1,Y1,X1))
    @test vec(sum(x1,dims=2)) == fill(120,10)
    @test vec(sum(z1,dims=2)) == fill(120,10)
    @test sum(y1) == 120
    @test all(x1' .== y1)
    @test all(z1' .== y1)
end

@testset "SparseArray" begin
    for var in (Xs, ys)
        @test typeof(shuffleobs(var)) <: SubArray
        @test numobs(shuffleobs(var)) == numobs(var)
    end
    # tests if all obs are still present and none duplicated
    @test vec(sum(getobs(shuffleobs(sparse(X1))),dims=2)) == fill(120,10)
    @test sum(getobs(shuffleobs(sparse(Y1)))) == 120
end

@testset "Tuple of SparseArray" begin
    for var in ((Xs,ys), (X,ys), (Xs,y), (Xs,Xs), (XX,X,ys))
        @test_throws MethodError shuffleobs(var...)
        @test typeof(shuffleobs(var)) <: Tuple
        @test numobs(shuffleobs(var)) == numobs(var)
    end
    # tests if all obs are still present and none duplicated
    # also tests that both paramter are shuffled identically
    x1, y1 = getobs(shuffleobs((sparse(X1),sparse(Y1))))
    @test vec(sum(x1,dims=2)) == fill(120,10)
    @test sum(y1) == 120
    @test all(x1' .== y1)
end

@testset "RNG" begin
    # tests reproducibility
    explicit_shuffle = shuffleobs(MersenneTwister(42), (X, y))
    @test explicit_shuffle == shuffleobs(MersenneTwister(42), (X, y))
end


text:
module TwoSourceSound
	using QuadGK
	using Cubature
	using PyCall
	const SPEED_OF_SOUND = 343

	struct TwoSourceOpts
		f0::Float64
		f1::Float64
		d::Float64
		center::Tuple{Float64, Float64}
	end

	export TwoSourceOpts

	function calcAvgAmplitudeIntegral(a::Float64, b::Float64, f::Float64, df::Float64)
		k = 2 * pi / SPEED_OF_SOUND

		return -(b*cos(a*f*k) + a*cos(b*f*k))/(a*b*k*df)
	end

	function calcAvgAmplitude(opts::TwoSourceOpts, p::Tuple{Float64, Float64})
		df = opts.f1 - opts.f0
		spk0x = opts.center[1] - opts.d / 2
		spk1x = opts.center[1] + opts.d / 2
		a = sqrt((p[1] - spk0x)^2 + (p[2] - opts.center[2])^2)
		b = sqrt((p[1] - spk1x)^2 + (p[2] - opts.center[2])^2)

		return calcAvgAmplitudeIntegral(a, b, opts.f1, df) - calcAvgAmplitudeIntegral(a, b, opts.f0, df)
	end

	function calcAvgAmplitudeDevIntegral(a::Float64, b::Float64, avgAmp::Float64, f::Float64, df::Float64)
		k = 2 * pi / SPEED_OF_SOUND

		m = a*avgAmp*b*f*k + a*cos(b*f*k) + b*cos(a*f*k)
		n = sign(sin(a*f*k) + sin(b*f*k) - avgAmp)
		l = a * b * k * df

		return (-m)/l
	end

	function calcAvgAmplitudeDev(opts::TwoSourceOpts, p::Tuple{Float64, Float64})
		df = opts.f1 - opts.f0
		spk0x = opts.center[1] - opts.d / 2
		spk1x = opts.center[1] + opts.d / 2
		a = sqrt((p[1] - spk0x)^2 + (p[2] - opts.center[2])^2)
		b = sqrt((p[1] - spk1x)^2 + (p[2] - opts.center[2])^2)
		avgAmp = calcAvgAmplitude(opts, p)

		f(x) = abs(sin(2*pi*a*x/SPEED_OF_SOUND) + sin(2*pi*b*x/SPEED_OF_SOUND) - avgAmp)

		return (1/df) * quadgk(f, opts.f0, opts.f1, rtol=0.001)[1]
	end

	function calcAvgPlaneAmplitudeDev(d::Float64, fRange::Tuple{Float64, Float64}, xRange::Tuple{Float64, Float64}, yRange::Tuple{Float64, Float64})
		df = fRange[2] - fRange[1]
		dx = xRange[2] - xRange[1]
		dy = yRange[2] - yRange[1]

		A_0(f, x, y) = sin(2*pi*sqrt(x^2 + y^2)*f/SPEED_OF_SOUND) + sin(2*pi*sqrt((x - d)^2 + y^2)*f/SPEED_OF_SOUND)
		A_bar_indef(f, x, y) = -(SPEED_OF_SOUND/(2*pi*df)) * ((cos(2*pi*f*sqrt((x-d)^2 + y^2)))/(sqrt((x-d)^2 + y^2)) + (cos(2*pi*f*sqrt(x^2 + y^2)))/(sqrt(x^2 + y^2)))
		A_bar(x, y) = A_bar_indef(fRange[2], x, y) - A_bar_indef(fRange[1], x, y)

		func(y, x, f) = abs(A_0(f, x, y) - A_bar(x, y))
		func_cubature(x) = func(x[1], x[2], x[3])

		return 1/(dy * dx * df) * hcubature(func_cubature, (yRange[1], xRange[1], fRange[1]), (yRange[2], xRange[2], fRange[2]), reltol=0.01)[1]
	end

	export calcAvgAmplitudeDev
end


text:
typealias FirstOrderSolver Union{AcceleratedGradientDescent, ConjugateGradient, GradientDescent,
                                 MomentumGradientDescent, BFGS, LBFGS}
typealias SecondOrderSolver Union{Newton, NewtonTrustRegion}
# Multivariate optimization
function check_kwargs(kwargs, fallback_method)
    kws = Dict{Symbol, Any}()
    method = nothing
    for kwarg in kwargs
        if kwarg[1] != :method
            kws[kwarg[1]] = kwarg[2]
        else
            method = kwarg[2]
        end
    end

    if method == nothing
        method = fallback_method
    end
    kws, method
end

function optimize{F<:Function}(f::F, initial_x::Array; kwargs...)
    checked_kwargs, method = check_kwargs(kwargs, NelderMead())
    optimize(f, initial_x, method, Options(; checked_kwargs...))
end

function optimize{F<:Function, G<:Function}(f::F, g!::G, initial_x::Array; kwargs...)
    checked_kwargs, method = check_kwargs(kwargs, BFGS())
    optimize(f, g!, initial_x, method, Options(;checked_kwargs...))
end
function optimize(d::DifferentiableFunction, initial_x::Array; kwargs...)
    checked_kwargs, method = check_kwargs(kwargs, BFGS())
    optimize(d, initial_x, method, Options(checked_kwargs...))
end

function optimize{F<:Function, G<:Function, H<:Function}(f::F,
                  g!::G,
                  h!::H,
                  initial_x::Array; kwargs...)
    checked_kwargs, method = check_kwargs(kwargs, Newton())
    optimize(f, g!, h!, initial_x, method, Options(checked_kwargs...))
end

function optimize(d::TwiceDifferentiableFunction, initial_x::Array; kwargs...)
    checked_kwargs, method = check_kwargs(kwargs, Newton())
    optimize(d, initial_x, method, Options(;kwargs...))
end

function optimize(d, initial_x::Array,
                  method::Optimizer, options::Options = Options())
    optimize(d, initial_x, method, options)
end
optimize(d::Function, initial_x, options::Options) = optimize(d, initial_x, NelderMead(), options)
optimize(d::DifferentiableFunction, initial_x, options::Options) = optimize(d, initial_x, BFGS(), options)
optimize(d::TwiceDifferentiableFunction, initial_x, options::Options) = optimize(d, initial_x, Newton(), options)

function optimize{F<:Function, G<:Function}(f::F,
                  g!::G,
                  initial_x::Array,
                  method::Optimizer,
                  options::Options = Options())
    d = DifferentiableFunction(f, g!)
    optimize(d, initial_x, method, options)
end
function optimize{F<:Function, G<:Function}(f::F,
                  g!::G,
                  initial_x::Array,
                  options::Options)
    d = DifferentiableFunction(f, g!)
    optimize(d, initial_x, BFGS(), options)
end

function optimize{F<:Function, G<:Function, H<:Function}(f::F,
                  g!::G,
                  h!::H,
                  initial_x::Array,
                  method::Optimizer,
                  options::Options = Options())
    d = TwiceDifferentiableFunction(f, g!, h!)
    optimize(d, initial_x, method, options)
end
function optimize{F<:Function, G<:Function, H<:Function}(f::F,
                  g!::G,
                  h!::H,
                  initial_x::Array,
                  options)
    d = TwiceDifferentiableFunction(f, g!, h!)
    optimize(d, initial_x, Newton(), options)
end

function optimize{F<:Function, T, M <: Union{FirstOrderSolver, SecondOrderSolver}}(f::F,
                  initial_x::Array{T},
                  method::M,
                  options::Options)
    if !options.autodiff
        if M <: FirstOrderSolver
            d = DifferentiableFunction(f)
        else
            error("No gradient or Hessian was provided. Either provide a gradient and Hessian, set autodiff = true in the Options if applicable, or choose a solver that doesn't require a Hessian.")
        end
    else
        gcfg = ForwardDiff.GradientConfig(initial_x)
        g! = (x, out) -> ForwardDiff.gradient!(out, f, x, gcfg)

        fg! = (x, out) -> begin
            gr_res = DiffBase.DiffResult(zero(T), out)
            ForwardDiff.gradient!(gr_res, f, x, gcfg)
            DiffBase.value(gr_res)
        end

        if M <: FirstOrderSolver
            d = DifferentiableFunction(f, g!, fg!)
        else
            hcfg = ForwardDiff.HessianConfig(initial_x)
            h! = (x, out) -> ForwardDiff.hessian!(out, f, x, hcfg)
            d = TwiceDifferentiableFunction(f, g!, fg!, h!)
        end
    end

    optimize(d, initial_x, method, options)
end

function optimize(d::DifferentiableFunction,
                  initial_x::Array,
                  method::Newton,
                  options::Options)
    if !options.autodiff
        error("No Hessian was provided. Either provide a Hessian, set autodiff = true in the Options if applicable, or choose a solver that doesn't require a Hessian.")
    else
        hcfg = ForwardDiff.HessianConfig(initial_x)
        h! = (x, out) -> ForwardDiff.hessian!(out, d.f, x, hcfg)
    end
    optimize(TwiceDifferentiableFunction(d.f, d.g!, d.fg!, h!), initial_x, method, options)
end

function optimize(d::DifferentiableFunction,
                  initial_x::Array,
                  method::NewtonTrustRegion,
                  options::Options)
    if !options.autodiff
        error("No Hessian was provided. Either provide a Hessian, set autodiff = true in the Options if applicable, or choose a solver that doesn't require a Hessian.")
    else
        hcfg = ForwardDiff.HessianConfig(initial_x)
        h! = (x, out) -> ForwardDiff.hessian!(out, d.f, x, hcfg)
    end
    optimize(TwiceDifferentiableFunction(d.f, d.g!, d.fg!, h!), initial_x, method, options)
end

update_g!(d, state, method) = nothing

function update_g!{M<:Union{FirstOrderSolver, Newton}}(d, state, method::M)
    # Update the function value and gradient
    state.f_x_previous, state.f_x = state.f_x, d.fg!(state.x, state.g)
    state.f_calls, state.g_calls = state.f_calls + 1, state.g_calls + 1
end

update_h!(d, state, method) = nothing

# Update the Hessian
function update_h!(d, state, method::SecondOrderSolver)
    d.h!(state.x, state.H)
    state.h_calls += 1
end

after_while!(d, state, method, options) = nothing

function optimize{T, M<:Optimizer}(d, initial_x::Array{T}, method::M, options::Options)
    t0 = time() # Initial time stamp used to control early stopping by options.time_limit

    if length(initial_x) == 1 && typeof(method) <: NelderMead
        error("Use optimize(f, scalar, scalar) for 1D problems")
    end

    state = initial_state(method, options, d, initial_x)

    tr = OptimizationTrace{typeof(method)}()
    tracing = options.store_trace || options.show_trace || options.extended_trace || options.callback != nothing
    stopped, stopped_by_callback, stopped_by_time_limit = false, false, false

    x_converged, f_converged = false, false
    g_converged = if typeof(method) <: NelderMead
        nmobjective(state.f_simplex, state.m, state.n) < options.g_tol
    elseif  typeof(method) <: ParticleSwarm || typeof(method) <: SimulatedAnnealing
        g_converged = false
    else
        vecnorm(state.g, Inf) < options.g_tol
    end

    converged = g_converged
    iteration = 0

    options.show_trace && print_header(method)
    trace!(tr, state, iteration, method, options)

    while !converged && !stopped && iteration < options.iterations
        iteration += 1

        update_state!(d, state, method) && break # it returns true if it's forced by something in update! to stop (eg dx_dg == 0.0 in BFGS)
        update_g!(d, state, method)
        x_converged, f_converged,
        g_converged, converged = assess_convergence(state, options)
        # We don't use the Hessian for anything if we have declared convergence,
        # so we might as well not make the (expensive) update if converged == true
        !converged && update_h!(d, state, method)

        # If tracing, update trace with trace!. If a callback is provided, it
        # should have boolean return value that controls the variable stopped_by_callback.
        # This allows for early stopping controlled by the callback.
        if tracing
            stopped_by_callback = trace!(tr, state, iteration, method, options)
        end

        # Check time_limit; if none is provided it is NaN and the comparison
        # will always return false.
        stopped_by_time_limit = time()-t0 > options.time_limit ? true : false

        # Combine the two, so see if the stopped flag should be changed to true
        # and stop the while loop
        stopped = stopped_by_callback || stopped_by_time_limit ? true : false
    end # while

    after_while!(d, state, method, options)

    return MultivariateOptimizationResults(state.method_string,
                                            initial_x,
                                            state.x,
                                            Float64(state.f_x),
                                            iteration,
                                            iteration == options.iterations,
                                            x_converged,
                                            options.x_tol,
                                            f_converged,
                                            options.f_tol,
                                            g_converged,
                                            options.g_tol,
                                            tr,
                                            state.f_calls,
                                            state.g_calls,
                                            state.h_calls)
end

# Univariate Options
function optimize{F<:Function, T <: AbstractFloat}(f::F,
                                      lower::T,
                                      upper::T;
                                      method = Brent(),
                                      rel_tol::Real = sqrt(eps(T)),
                                      abs_tol::Real = eps(T),
                                      iterations::Integer = 1_000,
                                      store_trace::Bool = false,
                                      show_trace::Bool = false,
                                      callback = nothing,
                                      show_every = 1,
                                      extended_trace::Bool = false)
    show_every = show_every > 0 ? show_every: 1
    if extended_trace && callback == nothing
        show_trace = true
    end
    if show_trace
        @printf "Iter     Function value   Gradient norm \n"
    end
    optimize(f, Float64(lower), Float64(upper), method;
             rel_tol = rel_tol,
             abs_tol = abs_tol,
             iterations = iterations,
             store_trace = store_trace,
             show_trace = show_trace,
             show_every = show_every,
             callback = callback,
             extended_trace = extended_trace)
end

function optimize{F<:Function}(f::F,
                  lower::Real,
                  upper::Real;
                  kwargs...)
    optimize(f,
             Float64(lower),
             Float64(upper);
             kwargs...)
end

function optimize{F<:Function}(f::F,
                  lower::Real,
                  upper::Real,
                  mo::Union{Brent, GoldenSection};
                  kwargs...)
    optimize(f,
             Float64(lower),
             Float64(upper),
             mo;
             kwargs...)
end


text:
"""
    build_fνmax_fast(f::AbstractArray{T}, n::Int)::Array{T, 1} where T <: Number

Description
-------------
Constructs array of partial sums of one or two-dimensional array `f` with
`n` partial sum terms.
This assumes, that the array is symmetric around the mid index.

Examples
-------------
```
julia> LadderDGA.build_fνmax_fast([3,2,1,2,3],2)
[5, 11]
julia> arr = [3 3 3 3 3; 3 2 2 2 3; 3 2 1 2 3; 3 2 2 2 3; 3 3 3 3 3
julia> LadderDGA.build_fνmax_fast(arr,2)
[17, 65]
```
"""
function build_fνmax_fast(f::AbstractArray{T}, n::Int)::Array{T, 1} where T <: Number
    n_iν = minimum(size(f))
    lo = n
    up = n_iν - n + 1
    f_νmax  = Array{T, 1}(undef, lo)
    build_fνmax_fast!(f_νmax, f, lo, up)
    return f_νmax
end

"""
    build_fνmax_fast!(f::AbstractArray{T,1}, lo::Int, up::Int)::Array{T, 1} where T <: Number
    build_fνmax_fast!(f::AbstractArray{T,2}, lo::Int, up::Int)::Array{T, 1} where T <: Number

Description
-------------
Constructs array of partial sums of one or two-dimensional array `f` starting
at with `lo` and ending wie `up` summands.
This assumes, that the array is symmetric around the mid index.
For a more convenient, but possible slower version, see also [`build_fνmax_fast`](@ref build_fνmax_fast!)
"""
function build_fνmax_fast!(f_νmax::AbstractArray{T,1}, f::AbstractArray{T,1}, lo::Int, up::Int) where T <: Number
    f_νmax[1] = sum(f[lo:up])
    for i in 2:length(f_νmax)
        lo = lo - 1
        up = up + 1
        @inbounds f_νmax[i] = f_νmax[i-1] + f[lo] + f[up]
    end
end

function build_fνmax_fast!(f_νmax::AbstractArray{T,1}, f::AbstractArray{T,2}, lo::Int, up::Int) where T <: Number
    f_νmax[1] = sum(f[lo:up, lo:up])
    for i in 2:length(f_νmax)
        lo = lo - 1
        up = up + 1
        @inbounds f_νmax[i] = f_νmax[i-1] + sum(f[lo, lo:up]) + sum(f[up, lo:up]) +
                    sum(f[(lo+1):(up-1),lo]) + sum(f[(lo+1):(up-1),up])
    end
end

#TODO: write function to tune for best range/exps. use usable_ω for that purpose
default_fit_range(arr::AbstractArray) = default_fit_range(length(arr))
default_fit_range(s::Int) = ceil(Int,s/3):floor(Int, s/2)

"""
    get_sum_helper(range, sP::SimulationParameters)

Description
-------------
Construct helper for (improved) sums from setting in
`sP` over a given fit range `range`.
"""
function get_sum_helper(range, sP::SimulationParameters, type)
    tc = if type == :f sP.tc_type_f else sP.tc_type_b end
    fitRange = (type == :f) ? (default_fit_range(range)) : (1:length(sP.bosonic_tail_coeffs))

    sumHelper = if tc == :nothing
        Naive()
    elseif tc == :richardson
        (type == :f) ? Richardson(fitRange, sP.fermionic_tail_coeffs, method=:bender) : Richardson(fitRange, sP.bosonic_tail_coeffs, method=:bender)
    elseif tc == :coeffs
        #TODO: move coeffs to SeriesAccelleration
        #Naive()
        (type == :f) ? Richardson(fitRange, sP.fermionic_tail_coeffs, method=:bender) : Richardson(fitRange, sP.bosonic_tail_coeffs, method=:bender)
    else
        @error("Unrecognized tail correction, falling back to naive sums!")
        Naive()
    end
    return sumHelper
end

"""
    sum_freq(f, dims::AbstractArray{Int,1}, type::T, β::Float64;
                  corr::Float64=0.0) where T <: SumHelper

Description
-------------
Wrapper for estimation of initite Matsubara sums from a finite set
of samples. Computes ``\\frac{1}{\\beta}^D \\sum_{n} f(i\\omega_n)``.
For convenience, a correction term (being the analytic sum over a
previously subtracted tail) can also be provided throught the `corr` argument.
For a faster version of 1D arrays, see also [`sum_freq_1D`](@ref sum_freq_1D).

Arguments
-------------
- **`f`**    : array containing function of Matsubara frequencies
- **`dims`** : dimensions over which to sum
- **`type`** : SumHelper object, used for estimation of limit
- **`β`**    : Inverse Temperature for normalization
- **`corr`** : When tail coefficients are know, the subtraction and
               subsequent addition of the analytic infinite sum over
               the tail can be added back here.

Examples
-------------
```
julia> arr = ones(3,5,5)
julia> sum_freq(arr, [2,3], Naive(), 1.0)
3×1×1 Array{Float64, 3}:
    [:, :, 1] =
     25.0
     25.0
     25.0
```
"""
function sum_freq(f::AbstractArray{T1}, dims::Array{Int,1}, type::T2, β::Float64;
        corr::Float64=0.0) where {T1 <: Real, T2 <: SumHelper}
    length(dims) == ndims(f) && return  sum_freq_full(f,type,β,corr=corr)
    res = mapslices(x -> esum_c(build_fνmax_fast(x, npartial_sums(type)) .+ corr, type), f, dims=dims)
    return res/(β^length(dims))
end

#TODO: provide interface for tail subtraction in sum_freq

function sum_freq(f::AbstractArray{T1}, dims::Array{Int,1}, type::T2, β::Float64;
        corr::Float64=0.0) where {T1 <: Complex, T2 <: SumHelper}
    length(dims) == ndims(f) && return sum_freq_full(f,type,β,corr=corr)
    res_re = mapslices(x -> esum_c(build_fνmax_fast(x,npartial_sums(type)) .+ corr, type), real.(f), dims=dims)
    res_im = mapslices(x -> esum_c(build_fνmax_fast(x,npartial_sums(type)) .+ corr, type), imag.(f), dims=dims)
    res = res_re + res_im*im
    return res/(β^length(dims))
end

"""
    sum_freq_full(f, type::T, β::Float64; corr::Float64=0.0) where T <: SumHelper

Description
-------------
Fast wrapper of [`sum_freq`](@ref sum_freq), lacking the `dims` argument.
"""
sum_freq_full(f::AbstractArray{Float64}, type::T2, β::Float64;corr::Float64=0.0) where {T2 <: SumHelper} =
    esum_c(build_fνmax_fast(f, npartial_sums(type)) .+ corr, type)/(β^ndims(f))

function sum_freq_full(f::AbstractArray{Complex{Float64}}, type::T2, β::Float64;corr::Float64=0.0) where {T2 <: SumHelper}
    tmp = build_fνmax_fast(f, npartial_sums(type)) .+ corr
    return (esum_c(real.(tmp), type) .+ esum_c(imag.(tmp), type).*im)/(β^ndims(f))
end

function find_usable_interval(arr::Array{Float64,1}; sum_type::Union{Symbol,Tuple{Int,Int}}=:common, reduce_range_prct::Float64 = 0.0)
    mid_index = Int(ceil(length(arr)/2))
    if sum_type == :full
        return 1:length(arr)
    elseif typeof(sum_type) == Tuple{Int,Int}
        return default_sum_range(mid_index, sum_type)
    end

    darr = diff(arr; dims=1)
    if arr[mid_index] < 0.0
        res = [mid_index]
        return res
    end
    # interval for condition 1 (positive values)
    cond1_intervall_range = 1
    # find range for positive values
    @inbounds while (cond1_intervall_range < mid_index) &&
        (arr[(mid_index-cond1_intervall_range)] > 0) &&
        (arr[(mid_index+cond1_intervall_range)] > 0)
        cond1_intervall_range = cond1_intervall_range + 1
    end

    # interval for condition 2 (monotonicity)
    cond2_intervall_range = 1
    # find range for first turning point
   @inbounds while (cond2_intervall_range < mid_index-1) &&
        (darr[(mid_index-cond2_intervall_range)] > 0) &&
        (darr[(mid_index+cond2_intervall_range)] < 0)
        cond2_intervall_range = cond2_intervall_range + 1
    end
    intervall_range = minimum([cond1_intervall_range, cond2_intervall_range])
    range = ceil(Int64, intervall_range*(1-reduce_range_prct))
    if length(arr)%2 == 1
        res = ((mid_index-range+1):(mid_index+range-2) .+ 1)
    else
        res = ((mid_index-range+1):(mid_index+range-2) .+ 2)
    end

    if length(res) < 1
        @warn "   ---> WARNING: could not determine usable range. Defaulting to single frequency!"
        res = [mid_index]
    end
    return res
end

"""
    find_usable_γ(arr; threshold=50, prct_red=0.05)

Usable νₙ range for γ. γ(ωₙ) is usable, for a range in which `γ[n]/γ[n-1]` does not
exceed some threshold value. See also [`extend_γ!`](@ref extend_γ!)

Arguments
-------------
- **`arr`**    : 1D γ(νₙ) slice for fixed ω and q
- **`threshold`** : Optional, default `50`, unusable if `arr[i]/arr[i±1] > threshold`
- **`prct_red`**  : Optional, default `0.05`, reduce unusble range by this amount
"""
function find_usable_γ(arr; threshold=50, prct_red=0.05)
    indh = ceil(Int, length(arr)/2)
    i = indh + floor(Int,indh/5)
    red = ceil(Int,length(arr)*prct_red)
    lo = 1
    up = length(arr)
    while i <= length(arr)
        change = abs(real(arr[i]/arr[i-1]))
        (change > threshold) && (up = i-1; break)
        i += 1
    end
    i = indh - floor(Int,indh/5)
    while i > 0
        change = abs(real(arr[i]/arr[i+1]))
        (change > threshold) && (lo = i+1; break)
        i -= 1
    end
    (lo > 1) && (lo += red)
    (up < length(arr)) && (up -= red)
    return lo,up
end

"""
    extend_γ!(arr::AbstractArray{Complex{Float64},1})

Extends γ for fermionic frequencies, for which the sum extrapolation failed.
See also (`extend_γ(::AbstractArray, ::AbstractArray)`) for a version, that uses
the naive sum as reference.

Description
-------------
Assuming that ∀ ωₙ,q  γ[ωₙ,q,νₙ] → 1 + 0i for ν → ∞, this function cuts off
typical overfitting spikes for large νₙ. It is furthermore assumed, that
γ[ωₙ,q,νₙ] /  γ[ωₙ,q,νₙ ± ν_{n+1}] never exceeds 50.
If such a large change is detected all γ[ωₙ,q,νₙ] after that point are
set to df⋅γ[ωₙ,q,νₙ ∓ ν{n}], df being the difference between the last
accepted elements.
The tail behavior is approximated by the 0th to 2nd Taylor terms.
This will most likely lead to a discontinuiety at the first altered frequency
in the second derivative.
Since it is know, that the real part tends towards 1 and the imaginary part
towards 0, a frequency depndended dampening factor is multiplied with the
first and second derivative. I.e. `arr[i] = (1 + df/h + ddf/h^2) arr[i-1]`
and `df = weight * df`, `ddf = weight^2 ddf` at each iteration.

Arguments
-------------
- **`arr`**    : 1D array containing sum-extrapolated γ[νₙ] (for fixed ωₙ and q)
- **`h`**      : 2 pi/ beta
- **`weight`** : Optional, default `0.9`, see description above for effect of weight.

Examples
-------------
```
extend_γ!(view(γ,1,1,:), 0.5235987755982988)
 """
function extend_γ!(arr::AbstractArray{Complex{Float64},1}, h::Float64; weight=0.9)
    indh = ceil(Int, length(arr)/2)
    lo,up = find_usable_γ(arr)
    # left
    i = lo
    df = -conj(arr[i] - arr[i+1])/h
    ddf = -conj(arr[i+2] - 2*arr[i+1] + arr[i])/(2*h^2)
    wi = weight
    while i > 1
        i -= 1
        arr[i] = (ddf*wi + df*wi + 1) * arr[i+1]
        #df = df * wi
        #ddf = ddf * wi
        wi = 10*(arr[i] - 1.0 + 0im)
    end
    # right
    i = up
    df = -conj(arr[i] - arr[i-1])/h
    ddf = -conj(arr[i-2] - 2*arr[i-1] + arr[i])/(2*h^2)
    wi = weight
    while i < length(arr)
        i += 1
        arr[i] = (ddf*wi + df*wi + 1) * arr[i-1]
        #df = df * wi
        #ddf = ddf * wi
        wi = 10*(arr[i] - 1.0 + 0im)
    end
end

"""
    extend_γ!(arr::AbstractArray{Complex{Float64},1}, ref::AbstractArray{Complex{Float64},1})

Extends γ for fermionic frequencies, for which the sum extrapolation failed.
See also (`extend_γ(::AbstractArray)`) for a version, that does not require
the naive sum as input.

Description
-------------
Assuming that ∀ ωₙ,q  γ[ωₙ,q,νₙ] → 0 for ν → ∞, this function cuts off
typical overfitting spikes for large νₙ. It is furthermore assumed, that
γ[ωₙ,q,νₙ] /  γ[ωₙ,q,νₙ ± ν_{n+1}] never exceeds 50.
If such a large change is detected all γ[ωₙ,q,νₙ] after that point are
set to df⋅γ[ωₙ,q,νₙ ∓ ν{n}], df being the difference between the last
accepted elements.
This will most likely lead to a discontinuiety at the first altered frequency.
See also (`extend_γ(::AbstractArray, ::AbstractArray)`) for a version, that uses
the naive sum as reference.

Arguments
-------------
- **`arr`**    : 1D array containing sum-extrapolated γ[νₙ] (for fixed ωₙ and q)
- **`ref`**    : 1D array containing naively summed γ[νₙ] (for fixed ωₙ and q)

Examples
-------------
```
extend_γ!(view(γ,1,1,:), view(γ_ref,1,1,:))
 """
function extend_γ!(arr::AbstractArray{Complex{Float64},1}, ref::AbstractArray{Complex{Float64},1})
    indh = ceil(Int, length(arr)/2)
    i = floor(Int,indh/5)
    override= false
    # right
    while i <= length(arr)
        change = abs(arr[i]/arr[i-1])
        (!override && change > 50) && (arr[i-1]=ref[i-1];override = true)
        override && (arr[i] = ref[i])
        i += 1
    end
    # left
    i = floor(Int,indh/5)
    override= false
    while i > 0
        change = abs(arr[i]/arr[i+1])
        (!override && change > 50) && (arr[i+1]=ref[i+1];override = true)
        override && (arr[i] = ref[i])
        i -= 1
    end
end


text:
module Junctions

using ..Blocks
using ..Nodes
using ..Graphs: Counter, add!

export Junction
export junctions, left, right

struct Junction
    left::Node{Block}
    right::Node{Block}
end

blocks(j::Junction) = (left=(j.left.block,j.left.strand), right=(j.right.block,j.right.strand))

# --------------------------------
# extension of base operators

Base.hash(j::Junction) = Base.hash(blocks(j)) ⊻ Base.hash(blocks(reverse(j)))
Base.isequal(j₁::Junction, j₂::Junction) = Base.isequal(blocks(j₁), blocks(j₂)) || Base.isequal(blocks(j₁), blocks(reverse(j₂)))

# --------------------------------
# custom operators

function left(j::Junction)
    return (
        block  = j.left.block,
        strand = j.left.strand,
    )
end

function right(j::Junction)
    return (
        block  = j.right.block,
        strand = j.right.strand,
    )
end

function reverse(j::Junction)
    return Junction(
        Node{Block}(j.right.block, !j.right.strand),
        Node{Block}(j.left.block,  !j.left.strand),
    )
end

function junctions(paths)
    js = Dict{Junction, Counter}()
    function push(key, name)
        if key ∈ keys(js)
            add!(js[key], name)
        else
            js[key] = Counter([name => 1])
        end
    end

    for path in paths
        for (i, node) ∈ enumerate(path.node[2:end])
            push(Junction(path.node[i], node), path.name)
        end

        if path.circular && length(path) > 1
            push(Junction(path.node[end],path.node[1]), path.name)
        end
    end

    return js
end

end

max_stars_repo_name: neherlab/pangraph

max_issues_repo_name: neherlab/pangraph

max_forks_repo_name: neherlab/pangraph


text:
const LOG_2π = log(2π)
const newaxis = [CartesianIndex()]

heappush!(h, v) = insert!(h, searchsortedfirst(h, v), v)

mutable struct UnionFind{T <: Integer}
    parent:: Vector{T}  # parent[root] is the negative of the size
    label::Dict{Int, Int}
    cnt::Int64

    function UnionFind{T}(nodes::T) where T<:Integer
        if nodes <= 0
            throw(ArgumentError("invalid argument for nodes: $nodes"))
        end

        parent = -ones(T, nodes)
        label = Dict([(i, i) for i in 1 : nodes])
        new{T}(parent, label, nodes)
    end
end

UnionFind(nodes::Integer) = UnionFind{typeof(nodes)}(nodes)
group(uf::UnionFind, x)::Int = uf.label[root(uf, x)]
members(uf::UnionFind, x::Int) = collect(keys(filter(n->n.second == x, uf.label)))

function root(uf::UnionFind{T}, x::T)::T where T<:Integer
    if uf.parent[x] < 0
        return x
    else
        # uf.parent[x] = root{T}(uf, uf.parent[x])
        # return uf.parent[x]
        return uf.parent[x] = root(uf, uf.parent[x])
    end
end

function issame(uf::UnionFind{T}, x::T, y::T)::Bool where T<:Integer
    return root(uf, x) == root(uf, y)
end

function Base.size(uf::UnionFind{T}, x::T)::T where T<:Integer
    return -uf.parent[root(uf, x)]
end

function unite!(uf::UnionFind{T}, x::T, y::T)::Bool where T<:Integer
    x = root(uf, x)
    y = root(uf, y)
    if x == y
        return false
    end
    if uf.parent[x] > uf.parent[y]
        x, y = y, x
    end
    # unite smaller tree(y) to bigger one(x)
    uf.parent[x] += uf.parent[y]
    uf.parent[y] = x
    uf.cnt += 1
    uf.label[y] = uf.cnt
    for i in members(uf, group(uf, x))
        uf.label[i] = uf.cnt
    end
    return true
end

#error function
#if absolute of `x` is smaller than 2.4, we use Taylor expansion.
#other wise, we use Continued fraction expansion
function erf(x)
    absx = abs(x)
    if absx<2.4
        c=1
        a=0
        for i in 1 : 40
            a+=(x/(2i-1)*c)
            c=-c*x^2/i
        end
        return a*2/sqrt(π)
    else
        if absx>1e50
            a = 1
        else
            y = absx*sqrt(2)
            a = 0
            for i in 40:-1:1
                a=i/(y+a)
            end
            a=1-exp(-x^2)/(y+a)*sqrt(2/π)
        end
        if x<0
            return -a
        else
            return a
        end
    end
end
erfc(x) = 1-erf(x)
norm_cdf(x) = 1/2*erfc(-x/sqrt(2))

function process_parameters(dim, mean, cov)
    if dim === nothing
        if mean === nothing
            if cov === nothing
                dim = 1
            else
                cov = convert(Array{Float64}, cov)
                if ndims(cov) < 2
                    dim = 1
                else
                    dim = size(cov, 1)
                end
            end
        else
            mean = convert(Array{Float64}, mean)
            dim = length(mean)
        end
    else
        !isa(dim, Number) && throw(DimensionMismatch("dimension of random variable must be a scalar"))
    end

    if mean === nothing
        mean = zeros(dim)
    end
    mean = convert(Array{Float64}, mean)

    if cov === nothing
        cov = [1.0]
    end
    cov = convert(Array{Float64}, cov)

    if dim == 1
        mean = reshape(mean, 1)
        cov = reshape(cov, 1, 1)
    end

    if ndims(mean) != 1 || size(mean, 1) != dim
        throw(ArgumentError("array `mean` must be vector of length $dim"))
    end
    if ndims(cov) == 0
        cov = cov * Matrix{Float64}(I, dim, dim)
    elseif ndims(cov) == 1
        cov = diag(cov)
    else
        size(cov) != (dim, dim) && throw(DimensionMismatch("array `cov` must be at most two-dimensional, but ndims(cov) = $(ndims(cov))"))
    end
    return dim, mean, cov
end

function process_quantiles(x, dim)
    x = convert(Array{Float64}, x)

    if ndims(x) == 0
        x = [x]
    elseif ndims(x) == 1
        if dim == 1
            x = x[:, :]
        else
            x = x[newaxis, :]
        end
    end
    return x
end

function pinv_1d(v; _eps=1e-5)
    return [(abs(x)<_eps) ? 0 : 1/x for x in v]
end

function psd_pinv_decomposed_log_pdet(mat; cond=nothing, rcond=nothing)
    u, s = eigvecs(mat), eigvals(mat)

    if rcond !== nothing
        cond = rcond
    end
    if cond === nothing || cond == -1
        cond = 1e6 * Base.eps()
    end
    _eps = cond * maximum(abs.(s))

    if minimum(s) < -_eps
        throw(ArgumentError("the covariance matrix must be positive semidefinite"))
    end
    s_pinv = pinv_1d(s, _eps=_eps)
    U = u .* sqrt.(s_pinv')
    log_pdet = sum(log.(s[findall(s.>_eps)]))

    return U, log_pdet
end

function squeeze_output(out)
    if length(out) == 1
        out = Float64(out...)
    else
        out = vec(out)
    end
    return out
end

function _logpdf(x, mean, prec_U, log_det_cov)
    dim = size(x, ndims(x))
    dev = x - mean'
    tmp = (dev * prec_U).^2
    maha = sum(tmp, dims=ndims(tmp))
    maha = dropdims(maha, dims=ndims(tmp))
    return -0.5 * ((dim*LOG_2π+log_det_cov).+maha)
end

function pdf(x, mean, cov)
    dim, mean, cov = process_parameters(nothing, mean, cov)
    x = process_quantiles(x, dim)
    prec_U, log_det_cov = psd_pinv_decomposed_log_pdet(cov)
    out = exp.(_logpdf(x, mean, prec_U, log_det_cov))
    return squeeze_output(out)
end

logpdf(x, mean, cov) = log.(pdf(x, mean, cov))


text:
@testset "Test ext" begin
    internal = IS.InfrastructureSystemsInternal()
    @test isnothing(internal.ext)
    ext = IS.get_ext(internal)
    ext["my_value"] = 1
    @test IS.get_ext(internal)["my_value"] == 1

    internal2 = JSON2.read(JSON2.write(internal), IS.InfrastructureSystemsInternal)
    @test internal.uuid == internal2.uuid
    @test internal.ext == internal2.ext

    IS.clear_ext(internal)
    @test isnothing(internal.ext)
end


text:
module Functions

#=
two-variable function
x = 5.90133
y = 0.5
fmin = -43.3159
-30 <= x, y <= 30
=#
function chichinadze(pars)
    @assert(length(pars) == 2)
    x = pars[1]
    y = pars[2]
    return x^2 - 12.0x + 11.0 + 10.0 * cos(pi * x / 2.0) + 8.0 * sin(5 * pi * x) -
           (1 / sqrt(5)) * exp(-0.5 * (y - 0.5)^2)
end

#=
n-variable function
x = [1, 1, ..., 1]
fmin = 0
-10 <= x_i <= 10
=#
function levy(pars)
    d = length(pars)

    w = zeros(Float64, d)

    for ii = 1:d
        w[ii] = 1 + (pars[ii] - 1.0) / 4.0
    end

    term1 = (sin(pi * w[1]))^2
    term3 = (w[d] - 1)^2.0 * (1 + (sin(2.0 * pi * w[d]))^2.0)

    sum = 0
    for ii = 1:(d-1)
        wi = w[ii]
        new = (wi - 1)^2 * (1 + 10 * (sin(pi * wi + 1))^2)
        sum = sum + new
    end

    return term1 + sum + term3
end

#=
n-variable function
x = [0, 0, ..., 0]
fmin = 0
-32.768 <= x_i <= 32.768
=#
function ackley(xx)
    d = length(xx)
    a = 20
    b = 0.2
    c = 2 * pi

    sum1 = 0
    for i in 1:d
        sum1 = sum1 + xx[i]^2
    end

    sum2 = 0
    for i in 1:d
        sum2 = sum2 + cos(c * xx[i])
    end

    term1 = -a * exp(-b*sqrt(sum1/d))
    term2 = -exp(sum2/d)

    y = term1 + term2 + a + exp(1)

    return(y)
end

end

max_stars_repo_name: jbytecode/MCCGA

max_issues_repo_name: jbytecode/MCCGA

max_forks_repo_name: jbytecode/MCCGA


text:
"""
## glob_to_axial

This subroutine transforms the global end reactions
into an axial force for rod fin_els (2- or 3-d).

### Function
```julia
glob_to_axial(glob, coord)
```

### Arguments
```julia
* glob::Vector{Float64}      : Globale forces and moments
* coord::Matrix{Float64}     : Nodal coordinates
```

### REturn value
```julia
* ::Float64                  : Axial force
```
"""
function glob_to_axial(glob::Vector{Float64}, coord::LinearAlgebra.Adjoint{Float64,Array{Float64,2}})
  #
  # This subroutine transforms the global end reactions
  # into an axial force for rod fin_els (2- or 3-d).
  #
  local axial = 0.0
  local ndim = size(coord,2)
  local add = 0.0
  for i in 1:ndim
    add = add+(coord[2,i]-coord[1,i])^2
  end
  ell = sqrt(add)
  for i in 1:ndim
    axial = axial+(coord[2,i]-coord[1,i])/ell*glob[ndim+i]
  end
  axial
end


text:
const data_path = pwd() * "/data/"
const inflated_cost = 1e6
const number_of_runs_to_generate_heuristic = 1000

mutable struct MountainCarRTAAPlanner <: Planner
    mountaincar::MountainCar
    num_expansions::Int64
    actions::Vector{MountainCarAction}
    residuals::Matrix{Float64}
    heuristic::Matrix{Float64}
    params::MountainCarParameters
end

function MountainCarRTAAPlanner(
    mountaincar::MountainCar,
    num_expansions::Int64,
    params::MountainCarParameters,
)
    heuristic =
        fill(0.0, (mountaincar.position_discretization, mountaincar.speed_discretization))
    residuals =
        fill(0.0, (mountaincar.position_discretization, mountaincar.speed_discretization))
    MountainCarRTAAPlanner(
        mountaincar,
        num_expansions,
        getActions(mountaincar),
        residuals,
        heuristic,
        params,
    )
end

function generateHeuristic!(planner::MountainCarRTAAPlanner; max_steps = 1e4, num_runs = number_of_runs_to_generate_heuristic, cache = true)
    heuristic_path = getHeuristicFilePath(planner)
    if isfile(heuristic_path) && cache
        loadHeuristic!(planner)
    else
        println("Generating Heuristic")
        clearHeuristic!(planner)
        rng = MersenneTwister(0)
        for i = 1:num_runs
            state = init(planner.mountaincar; rng = rng)
            num_steps = 0
            while !checkGoal(planner.mountaincar, state) && num_steps < max_steps
                num_steps += 1
                best_action, info = act(planner, state)
                updateResiduals!(planner, info)
                state, cost = step(planner.mountaincar, state, best_action, planner.params)
            end
        end
        planner.heuristic = deepcopy(planner.residuals)
        clearResiduals!(planner)
        if cache
            saveHeuristic(planner)
        end
        println("Generated Heuristic")
    end
end

function generateHeuristic!(planners::Vector{MountainCarRTAAPlanner}; max_steps = 1e4)
    for planner in planners
        generateHeuristic!(planner, max_steps = max_steps)
    end
end

function clearHeuristic!(planner::MountainCarRTAAPlanner)
    planner.heuristic = fill(
        0.0,
        (
            planner.mountaincar.position_discretization,
            planner.mountaincar.speed_discretization,
        ),
    )
end

function clearResiduals!(planner::MountainCarRTAAPlanner)
    planner.residuals = fill(
        0.0,
        (
            planner.mountaincar.position_discretization,
            planner.mountaincar.speed_discretization,
        ),
    )
end

function clearResiduals!(planners::Vector{MountainCarRTAAPlanner})
    for planner in planners
        clearResiduals!(planner)
    end
end

function getHeuristic(planner::MountainCarRTAAPlanner, state::MountainCarDiscState)
    heuristic = planner.heuristic[state.disc_position+1, state.disc_speed+1]
    # adding 1 for julia indexing
    residual = planner.residuals[state.disc_position+1, state.disc_speed+1]
    heuristic + residual
end

function getLookaheadValue(planner::MountainCarRTAAPlanner, state::MountainCarDiscState)
    _, info = act(planner, state)
    info["best_node_f"]
end

function getSuccessors(
    planner::MountainCarRTAAPlanner,
    state::MountainCarDiscState,
    action::MountainCarAction,
)
    next_disc_state, cost = step(planner.mountaincar, state, action, planner.params)
    next_disc_state, cost
end

function checkGoal(planner::MountainCarRTAAPlanner, state::MountainCarDiscState)
    checkGoal(planner.mountaincar, state)
end

function getActions(planner::MountainCarRTAAPlanner, state::MountainCarDiscState)
    planner.actions
end

function act(planner::MountainCarRTAAPlanner, state::MountainCarDiscState)
    best_action, info = astar(planner, planner.num_expansions, state)
    return best_action, info
end

function act(
    planner::MountainCarRTAAPlanner,
    state::MountainCarDiscState,
    num_expansions::Int64,
)
    best_action, info = astar(planner, num_expansions, state)
    return best_action, info
end

function updateResiduals!(planner::MountainCarRTAAPlanner, info::Dict)
    best_node_f = info["best_node_f"]
    for node in info["closed"]
        state = node.state
        g = node.g
        h = planner.heuristic[state.disc_position+1, state.disc_speed+1]
        # adding 1 for julia indexing
        planner.residuals[state.disc_position+1, state.disc_speed+1] = best_node_f - g - h
    end
end

function saveHeuristic(planner::MountainCarRTAAPlanner)
    heuristic_path = getHeuristicFilePath(planner)
    save(heuristic_path, "heuristic", planner.heuristic + planner.residuals)
end

function getHeuristicFilePath(planner::MountainCarRTAAPlanner)
    heuristic_file_name =
        "heuristic" *
        "_$(planner.params.theta1)_$(planner.params.theta2)_$(planner.mountaincar.rock_c)" *
        ".jld"
    data_path * heuristic_file_name
end

function loadHeuristic!(planner::MountainCarRTAAPlanner)
    heuristic_path = getHeuristicFilePath(planner)
    if isfile(heuristic_path)
        planner.heuristic = load(heuristic_path)["heuristic"]
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Random, Distributions, Plots
	Random.seed!(0)
end

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing5.04"

# ╔═╡ fa68607e-22e5-11eb-0558-c9a4d9f77426
n, N, alpha = 3, 10^7, 0.1

# ╔═╡ 0a6fcc80-3a99-11eb-3de9-13429e2152ee
begin
	myT(nObs) = rand(Normal())/sqrt(rand(Chisq(nObs-1))/(nObs-1))
	mcQuantile = quantile([myT(n) for _ in 1:N],alpha)
	analyticQuantile = quantile(TDist(n-1),alpha)
end

# ╔═╡ 0a990730-3a99-11eb-0e9d-439c6e34d87c
Text("Quantile from Monte Carlo: $(mcQuantile)")

# ╔═╡ 0a998098-3a99-11eb-209b-2f467a34d644
Text("Analytic qunatile: $(analyticQuantile)")

# ╔═╡ 0a9fd40c-3a99-11eb-387b-05b58e298e6e
begin
	xGrid = -5:0.1:5
	plot(xGrid, pdf.(Normal(), xGrid), c=:black, label="Normal Distribution")
	scatter!(xGrid, pdf.(TDist(1) ,xGrid),
		c=:blue, msw=0, label="DOF = 1")
	scatter!(xGrid, pdf.(TDist(3), xGrid),
		c=:red, msw=0, label="DOF = 3")
	scatter!(xGrid, pdf.(TDist(100),xGrid),
		c=:green, msw=0, label="DOF = 100",
		xlims=(-4,4), ylims=(0,0.5), xlabel="X", ylabel="Density")
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing5.04"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═fa68607e-22e5-11eb-0558-c9a4d9f77426
# ╠═0a6fcc80-3a99-11eb-3de9-13429e2152ee
# ╠═0a990730-3a99-11eb-0e9d-439c6e34d87c
# ╠═0a998098-3a99-11eb-209b-2f467a34d644
# ╠═0a9fd40c-3a99-11eb-387b-05b58e298e6e
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12