text:
"""
	gmtgmtregress(cmd0::String="", arg1=[], kwargs...)

Linear regression of 1-D data sets.

Full option list at [`gmtregress`](http://gmt.soest.hawaii.edu/doc/latest/gmtregress.html)

Parameters
----------

- **A** : **all_slopes** : -- Str or List --        Flags = min/max/inc

    Instead of determining a best-fit regression we explore the full range of regressions.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/gmtregress.html#a)
- **C** : **confidence_level** : -- Int --      Flags = level

    Set the confidence level (in %) to use for the optional calculation of confidence
    bands on the regression [95].
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/gmtregress.html#c)
- **E** : **type_of_regression** : -- Str --    Flags = x|y|o|r

    Type of linear regression, i.e., select the type of misfit we should calculate.
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/gmtregress.html#e)
- **F** : **column_combination** : -- Str --   Flags = x|y|m|l|c

    Append a combination of the columns you wish returned;
    [`-F`](http://gmt.soest.hawaii.edu/doc/latest/gmtregress.html#f)
- **N** : **norm** : -- Int or Str --           Flags = 1|2|r|w

    Selects the norm to use for the misfit calculation.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/gmtregress.html#n)
- **S** : **restrict** : -- Str or [] --        Flags = [r]

    Restricts which records will be output.
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/gmtregress.html#s)
- **T** : **equi_space** : -- List or Str --     Flags = [min/max/]inc[+a|n]] or file|list

    Evaluate the best-fit regression model at the equidistant points implied by the arguments.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/gmtregress.html#t)
- **W** : **weighted** : -- Str or [] --     Flags = [w][x][y][r]

    Specifies weighted regression and which weights will be provided.
    [`-W`](http://gmt.soest.hawaii.edu/doc/latest/gmtregress.html#w)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_swap_xy)
"""
function gmtregress(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("gmtregress", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd = parse_V_params("", d)
	cmd, = parse_b(cmd, d)
	cmd, = parse_d(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_g(cmd, d)
	cmd, = parse_h(cmd, d)
	cmd, = parse_i(cmd, d)
	cmd, = parse_o(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'A', d, [:A :all_slopes])
	cmd = add_opt(cmd, 'C', d, [:C :confidence_level])
	cmd = add_opt(cmd, 'E', d, [:E :type_of_regression])
	cmd = add_opt(cmd, 'F', d, [:F :column_combination])
	cmd = add_opt(cmd, 'N', d, [:N :norm])
	cmd = add_opt(cmd, 'S', d, [:S :restrict])
	cmd = add_opt(cmd, 'T', d, [:T :equi_space])
	cmd = add_opt(cmd, 'W', d, [:W :weighted])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "gmtregress", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmtregress(arg1=[], cmd0::String=""; kw...) = gmtregress(cmd0, arg1; kw...)


text:
function rhs_to_args(ex::Expr)
    indsall = Union{Symbol,Expr}[]
    symorig, symgend, isconjs = Union{Symbol,Expr}[], Symbol[], Bool[]

    function exparse(x, isconj)
        if x isa Number
            x
        elseif @capture(x, -(rhs__))
            y = exparse.(rhs, isconj)
            :(-($(y...)))
        elseif @capture(x, +(rhs__))
            y = exparse.(rhs, isconj)
            :(+($(y...)))
        elseif @capture(x, *(rhs__))
            y = exparse.(rhs, isconj)
            :(*($(y...)))
        elseif @capture(x, conj(rhs_))
            y = exparse(rhs, !isconj)
            :(conj($y))
        elseif @capture(x, sym_[ind__])
            any(x -> isa(x, Integer), ind) && error("NCON style is unsupported")
            @gensym new
            push!(symorig, sym)
            push!(symgend, new)
            push!(isconjs, isconj)
            append!(indsall, ind)
            :($new[$(ind...)])
        else
            @gensym new
            push!(symorig, x)
            push!(symgend, new)
            push!(isconjs, isconj)
            new
        end
    end

    exreplaced = exparse(ex, false)
    return exreplaced, symorig, symgend, isconjs, indsall
end

function make_only_product(ex::Expr, sym::Symbol)
    hassym(x) =
        if @capture(x, -(y__) | +(y__) | *(y__))
            any(hassym.(y))
        elseif @capture(x, conj(y_))
            hassym(y)
        elseif @capture(x, $sym) || @capture(x, $sym[__])
            true
        else
            false
        end

    return MacroTools.postwalk(ex) do x
        if @capture(x, -(y__))
            @assert 1 ≤ length(y) ≤ 2
            if length(y) == 1
                x
            elseif hassym(first(y))
                first(y)
            elseif hassym(last(y))
                :(-$(last(y)))
            else
                x
            end
        elseif @capture(x, +(y__))
            @assert 1 ≤ length(y)
            y = filter(hassym, y)
            @assert length(y) ≤ 1
            length(y) == 1 ? first(y) : x
        else
            x
        end
    end
end

function make_scalar_first(ex::Expr)
    return MacroTools.postwalk(ex) do x
        if @capture(x, *(y__))
            tensors = Expr[]
            notensors = Union{Expr,Symbol,Number}[]

            for z in y
                if @capture(z, _[__]) || @capture(z, conj(_[__]))
                    push!(tensors, z)
                else
                    push!(notensors, z)
                end
            end

            isempty(tensors) && return x

            if isempty(notensors)
                return x
            elseif length(notensors) == 1
                return :(*($(notensors[1]), $(tensors...)))
            else
                return :(*(*($(notensors...)), $(tensors...)))
            end
        else
            return x
        end
    end
end

ext: jl


text:
using Printf

# Abstract types can't be instantiated like Structs
# but they can have subtypes
abstract type Animal end

# Dog inherits from Animal
struct Dog <: Animal
    name::String
    bark::String
end

# Cat inherits from Animal
struct Cat <: Animal
    name::String
    meow::String
end

function makeSound(animal::Dog)
    println("$(animal.name) says $(animal.bark)")
end

function makeSound(animal::Cat)
    println("$(animal.name) says $(animal.meow)")
end

function abstract_type_examples()
    println("--------------------------------------")
    println("Abstract type examples")

    println()
    bowser = Dog("Bowser", "Ruff")
    muffin = Cat("Muffin", "Meow")

    makeSound(bowser)
    makeSound(muffin)
end

max_stars_repo_name: James-P-D/JuliaDump

max_issues_repo_name: James-P-D/JuliaDump

max_forks_repo_name: James-P-D/JuliaDump


text:
module TestCat
    using Base.Test
    using DataFrames

    #
    # hcat
    #

    dvint = @data([1, 2, NA, 4])
    dvstr = @data(["one", "two", NA, "four"])

    df2 = DataFrame(Any[dvint, dvstr])
    df3 = DataFrame(Any[dvint])
    df4 = convert(DataFrame, [1:4 1:4])
    df5 = DataFrame(Any[@data([1,2,3,4]), dvstr])

    dfh = hcat(df3, df4)
    @test size(dfh, 2) == 3
    @test names(dfh) == [:x1, :x1_1, :x2]
    @test isequal(dfh[:x1], df3[:x1])
    @test isequal(dfh, [df3 df4])
    @test isequal(dfh, DataFrames.hcat!(DataFrame(), df3, df4))

    dfh3 = hcat(df3, df4, df5)
    @test names(dfh3) == [:x1, :x1_1, :x2, :x1_2, :x2_1]
    @test isequal(dfh3, hcat(dfh, df5))
    @test isequal(dfh3, DataFrames.hcat!(DataFrame(), df3, df4, df5))

    @test isequal(df2, DataFrames.hcat!(df2))

    #
    # vcat
    #

    null_df = DataFrame(Int, 0, 0)
    df = DataFrame(Int, 4, 3)

    # Assignment of rows
    df[1, :] = df[1, :]
    df[1:2, :] = df[1:2, :]

    # Broadcasting assignment of rows
    df[1, :] = 1

    # Assignment of columns
    df[1] = zeros(4)

    # Broadcasting assignment of columns
    df[:, 1] = 1
    df[1] = 3
    df[:x3] = 2

    vcat([])
    vcat(null_df)
    vcat(null_df, null_df)
    vcat(null_df, df)
    vcat(df, null_df)
    vcat(df, df)
    vcat(df, df, df)
    @test vcat(DataFrame[]) == DataFrame()

    alt_df = deepcopy(df)
    vcat(df, alt_df)

    # Don't fail on non-matching types
    df[1] = zeros(Int, nrow(df))
    vcat(df, alt_df)

    # Don't fail on non-matching names
    names!(alt_df, [:A, :B, :C])
    vcat(df, alt_df)

    dfr = vcat(df4, df4)
    @test size(dfr, 1) == 8
    @test names(df4) == names(dfr)
    @test isequal(dfr, [df4; df4])

    dfr = vcat(df2, df3)
    @test size(dfr) == (8,2)
    @test names(df2) == names(dfr)
    @test isna(dfr[8,:x2])

    # Eltype promotion
    @test eltypes(vcat(DataFrame(a = [1]), DataFrame(a = [2.1]))) == [Float64]
    @test eltypes(vcat(DataFrame(a = [NA]), DataFrame(a = [2.1]))) == [Float64]

    # Minimal container type promotion
    dfa = DataFrame(a = @pdata([1, 2, 2]))
    dfb = DataFrame(a = @pdata([2, 3, 4]))
    dfc = DataFrame(a = @data([2, 3, 4]))
    dfd = DataFrame(Any[2:4], [:a])
    @test vcat(dfa, dfb)[:a] == @pdata([1, 2, 2, 2, 3, 4])
    @test vcat(dfa, dfc)[:a] == @pdata([1, 2, 2, 2, 3, 4])
    # ^^ container may flip if container promotion happens in Base/DataArrays
    dc = vcat(dfd, dfc)
    @test vcat(dfc, dfd) == dc

    # Zero-row DataFrames
    dfc0 = similar(dfc, 0)
    @test vcat(dfd, dfc0, dfc) == dc
    @test eltypes(vcat(dfd, dfc0)) == eltypes(dc)

    # Missing columns
    rename!(dfd, :a, :b)
    dfda = DataFrame(b = @data([2, 3, 4, NA, NA, NA]),
                     a = @pdata([NA, NA, NA, 1, 2, 2]))
    @test isequal(vcat(dfd, dfa), dfda)

    # Alignment
    @test isequal(vcat(dfda, dfd, dfa), vcat(dfda, dfda))
end


text:
using Statistics
using LinearAlgebra
using Random

using Distributions

include("da.jl")
import .DA

include("ens_forecast.jl")
import .ens_forecast

include("models.jl")
import .Models

include("integrators.jl")
import .Integrators

Random.seed!(1)

D = 40
models = [Models.lorenz96_err, Models.lorenz96_err2, Models.lorenz96_err3,
          Models.lorenz96_err4]

model_errs_prescribed = [nothing, nothing, nothing, nothing]
model_true = Models.lorenz96_true
n_models = length(models)
obs_ops = [I(D), I(D), I(D), I(D)]
H = I(D)
ens_sizes = [20, 20, 20, 20]
model_sizes = [D, D, D, D]
integrator = Integrators.rk4
da_method = DA.ensrf
localization_radius = 4
localization = DA.gaspari_cohn_localization(localization_radius, D, cyclic=true)
x0 = randn(D)
t0 = 0.0
Δt = 0.05
outfreq = 1
transient = 2000
x = integrator(models[1], x0, t0, transient*outfreq*Δt, Δt, inplace=false)
R = Symmetric(diagm(0.25*ones(D)))
ens_errs = [Symmetric(diagm(0.25*ones(D))), Symmetric(diagm(0.25*ones(D))),
            Symmetric(diagm(0.25*ones(D))), Symmetric(diagm(0.25*ones(D)))]
gen_ensembles = false
all_orders = true
assimilate_obs = true
save_Q_hist = false
save_P_hist = false

leads = 1
x0 = x[end, :]

n_cycles = 3000*leads
ρ = 1e-3

window = 4

infos = Vector(undef, n_models)
for model=1:n_models
    model_errs = [0.1*diagm(ones(D))]

    ens_size = sum(ens_sizes)
    ensembles = [x0 .+ rand(MvNormal(R), ens_size)]

    info = ens_forecast.da_cycles(x0=x0, ensembles=ensembles, models=[models[model]],
                                  model_true=model_true, obs_ops=[obs_ops[model]], H_true=H,
                                  model_errs=model_errs,
                                  model_errs_prescribed=[model_errs_prescribed[model]],
                                  integrator=integrator, da_method=da_method,
                                  localization=localization, ens_sizes=[ens_size], Δt=Δt,
                                  window=window, n_cycles=n_cycles, outfreq=outfreq,
                                  model_sizes=model_sizes, R=R, ens_errs=ens_errs, ρ=ρ,
                                  gen_ensembles=gen_ensembles,
                                  assimilate_obs=assimilate_obs, save_analyses=false,
                                  leads=leads, save_Q_hist=save_Q_hist,
                                  save_P_hist=save_P_hist)
    infos[model] = info
end

ensembles = [x0 .+ rand(MvNormal(R), ens_sizes[model]) for model=1:n_models]

model_errs = [0.1*diagm(ones(D)) for model=1:n_models]

info_mm = ens_forecast.da_cycles(x0=x0, ensembles=ensembles, models=models,
                                 model_true=model_true, obs_ops=obs_ops, H_true=H,
                                 model_errs=model_errs,
                                 model_errs_prescribed=model_errs_prescribed,
                                 integrator=integrator, da_method=da_method,
                                 localization=localization, ens_sizes=ens_sizes, Δt=Δt,
                                 window=window, n_cycles=n_cycles, outfreq=outfreq,
                                 model_sizes=model_sizes, R=R, ens_errs=ens_errs, ρ=ρ,
                                 all_orders=all_orders, combine_forecasts=true,
                                 gen_ensembles=gen_ensembles, assimilate_obs=assimilate_obs,
                                 leads=leads, save_Q_hist=save_Q_hist)

ensembles = [x0 .+ rand(MvNormal(R), ens_sizes[model]) for model=1:n_models]

model_errs = [0.1*diagm(ones(D)) for model=1:n_models]

info_mm2 = ens_forecast.da_cycles(x0=x0, ensembles=ensembles, models=models,
                                  model_true=model_true, obs_ops=obs_ops, H_true=H,
                                  model_errs=model_errs,
                                  model_errs_prescribed=model_errs_prescribed,
                                  integrator=integrator, da_method=da_method,
                                  localization=localization, ens_sizes=ens_sizes, Δt=Δt,
                                  window=window, n_cycles=n_cycles, outfreq=outfreq,
                                  model_sizes=model_sizes, R=R, ens_errs=ens_errs, ρ=ρ,
                                  combine_forecasts=false, gen_ensembles=gen_ensembles,
                                  assimilate_obs=assimilate_obs, leads=leads,
                                  save_Q_hist=save_Q_hist)

max_stars_repo_name: eviatarbach/mmda

max_issues_repo_name: eviatarbach/mmda

max_forks_repo_name: eviatarbach/mmda


text:
"""

# Available functions
- `enable_output()`
- `disable_output()`
- `set_voltage(voltage)`
- `get_voltage()`
- `set_current_limit(current)`
- `get_current_limit()`
- `set_channel(channel)`
- `get_channel()`
"""
struct AgilentE36312A <: PowerSupply end

"""
    set_channel(obj::Instr{AgilentE36312A}, chan)

This will set the global channel on a device.

Any commands like set_voltage that affect the
device with respect to a specific channel will be impacted
by this command.

By setting the channel on a device, all subsequent commands
will operate on that channel unless they allow for an optional channel argument.

Parameters:
  - obj
    - must be a Power Supply Instrument
  - chan
    - This can be a string or int: 1, 2, 3 .. to n
    where n is the total number of channels

Supported Instruments:
   - Power supply

Returns:
  Nothing
"""
set_channel(obj::Instr{AgilentE36312A}, chan) = write(obj, "INST:NSEL $chan")

"""
    get_channel(obj::Instr{AgilentE36312A}; v=false)

This will return the global or default channel of a device.

Allows you to see what the global channel is set to at the
moment

Parameters:
  - obj
    - must be a Power Supply Instrument
  - chan
    - This can be a string or int: 1, 2, 3 .. to n
    where n is the total number of channels
  - v
    - optional boolean flag argument that is set to false
    - setting to `true` will print the verbose channel name

Supported Instruments:
   - Power supply

Returns:
  String
    - {"1"|"2"|...}
    - when v == true: "P6V", .. etc
"""
get_channel(obj::Instr{AgilentE36312A}; v=false) = v ? query(obj, "INST:SEL?") : strip(query(obj, "INST:NSEL?"), '+')

function psu_chan(obj, num, cmd; float=false)
    original_chan = get_channel(obj)
    num != 0 && set_channel(obj, num)
    ans = ""
    if occursin("?", cmd)
        if float
            ans = f_query(obj, cmd)
        else
            ans = query(obj, cmd)
        end
    else
        write(obj, cmd)
    end
    num != 0 && set_channel(obj, original_chan)
    ans
end

"""
    enable_output(obj::Instr{AgilentE36312A})

This will enable an output on a device.

If the device has multiple channels it will enable the
device for the currently selected channel. To see the
channel that will effected use the `get_channel` function.

If you want to enable a different channel, first use
`set_channel` to choose the channel. Running this function
subsequently will enable that channel

Arguments:
  - obj
    - must be a Power Supply Instrument
Supported Instruments:
   - Power supply
"""
enable_output(obj::Instr{AgilentE36312A}) = write(obj, ":OUTPUT:STATe ON")

"""
    disable_output(obj::Instr{AgilentE36312A})

This will disable an output on a device.

If the device has multiple channels it will disable the
device for the currently selected channel. To see the
channel that will effected use the `get_channel` function.

If you want to disable a different channel, first use
`set_channel` to choose the channel. Running this function
subsequently will disable that channel

Arguments:
  - obj
    - must be a Power Supply Instrument
Supported Instruments:
   - Power supply
"""
disable_output(obj::Instr{AgilentE36312A}) = write(obj, ":OUTPUT:STATE OFF")

"""
    get_output(obj::Instr{AgilentE36312A})

This will return the state of an output on a device.

If the device has multiple channels is will display the
state of the currently selected channel. To see the
channel that will read use the `get_channel` function.

If you want to see the state of a different channel, first use
`set_channel` to choose the channel. Running this function
subsequently will disable that channel

Arguments:
  - obj
    - must be a Power Supply Instrument
Supported Instruments:
   - Power supply

Returns:
  String: {"0"|"1"}
"""
function get_output(obj::Instr{AgilentE36312A})
    output = query(obj, ":OUTPUT:STATE?")
    if output == "0"
        return false
    elseif output == "1"
        return true
    else
        error("Undefined output: $output")
    end
end

"""
    set_voltage(obj::Instr{AgilentE36312A}, num::Voltage; chan=0)

This will change the voltage output voltage of a device.

Parameters:
  - obj
    - must be a Power Supply Instrument
  - num
    - integer or decimal of the desired voltage
  - chan
    - This is an optional parameter
    - If not provided it will use the default channel (see `set_channel`)
    - Otherwise this can be a string or int: 1, 2, 3 .. to n
    where n is the total number of channels

Supported Instruments:
   - Power supply

Returns:
  Nothing
"""
set_voltage(obj::Instr{AgilentE36312A}, num::Voltage; chan=0) = psu_chan(obj, chan, "SOURCE:VOLTAGE $(raw(num))")

"""
    get_voltage(obj::Instr{AgilentE36312A}; chan=0)

This will return the voltage of a device's channel.

Parameters:
  - obj
    - must be a Power Supply Instrument
  - chan
    - This is an optional parameter
    - If not provided it will use the default channel (see `set_channel`)
    - Otherwise this can be a string or int: 1, 2, 3 .. to n
    where n is the total number of channels

Supported Instruments:
   - Power supply

Returns:
  Voltage
"""
get_voltage(obj::Instr{AgilentE36312A}; chan=0) = psu_chan(obj, chan, "SOURCE:VOLTAGE?"; float=true) * V

"""
    set_current_limit(obj::Instr{AgilentE36312A}, num::Current; chan=0)

This will change the current limit of a device on a given
channel

Parameters:
  - obj
    - must be a Power Supply Instrument
  - num
    - Desired current limit of type Unitful Amps: 1.0u"A"
  - chan
    - This is an optional parameter
    - If not provided it will use the default channel (see `set_channel`)
    - Otherwise this can be a string or int: 1, 2, 3 .. to n
    where n is the total number of channels

Supported Instruments:
   - Power supply

Returns:
  Nothing
"""
set_current_limit(obj::Instr{AgilentE36312A}, num::Current; chan=0) = psu_chan(obj, chan, "SOURCE:CURRENT $(raw(num))")

"""
    get_current_limit(obj::Instr{AgilentE36312A}; chan=0)

This will return the current limit of a device.

Parameters:
  - obj
    - must be a Power Supply Instrument
  - chan
    - This is an optional parameter
    - If not provided it will use the default channel (see `set_channel`)
    - Otherwise this can be a string or int: 1, 2, 3 .. to n
    where n is the total number of channels

Supported Instruments:
   - Power supply

Returns:
  Current Limit
"""
get_current_limit(obj::Instr{AgilentE36312A}; chan=0) = psu_chan(obj, chan, "SOURCE:CURRENT?"; float=false) * A


text:
function countingsort!(v::Vector{<:Integer})
    mini, maxi = extrema(v)

    keys = zeros(Int, maxi-mini+1)

    for key in v
        keys[key-mini+1] += 1
    end

    i = 1
    for (key, amount) in enumerate(keys)
        v[i:i+amount-1] .= key+mini-1
        i += amount
    end
    return v
end

v = [144, 89, 4, 9, 95, 12, 86, 25]
println("unsorted array = ",v)

countingsort!(v)

println("counting sorted array = ", v)


text:
module FilterCircuit

using TinyModia
using DifferentialEquations
@usingModiaPlot
using Unitful

setLogMerge(true)

include("../models/ElectricTest.jl")

# res = Resistor
# @showModel(res)

@showModel(Capacitor)

Filter = Model(
    R = Resistor | Map(R=0.5u"Ω"),
    C = Capacitor | Map(C=2.0u"F", v=Map(init=0.1u"V")),
    V = ConstantVoltage | Map(V=10.0u"V"),
    ground = Ground,
    connect = :[
      (V.p, R.p)
      (R.n, C.p)
      (C.n, V.n, ground.p)
    ]
)

@showModel(Filter)

#=
Filter2 = Model(
    r = 1.0u"Ω",
    c = 1.0u"F",
    v = 1.0u"V",
    R = Resistor | Map(R=:(up.r)),
    C = Capacitor | Map(C=:(up.c)),
    V = ConstantVoltage | Map(V=:(up.v)),
    ground = Ground,
    connect = :[
      (V.p, R.p)
      (R.n, C.p)
      (C.n, V.n, ground.p)
    ]
)

Cpar = Map(C = 5.0u"F")

TwoFilters = Model( f1 = Filter | Map( R = Map(R = 10.0u"Ω"), C = Cpar), f2 = Filter)

VoltageDividerAndFilter = TwoFilters | Map(f1 = Map(C = Redeclare | Resistor | (R = 20.0u"Ω", start = Map(v = 0u"V"))))

println("Build")
@time Filters = Model(
    filters = [Filter for i in 1:10]
)

setLogMerge(false)

model = @instantiateModel(Filters, logDetails=false, aliasReduction=false)

println("Simulate")
@time simulate!(model, Tsit5(), stopTime = 50, requiredFinalStates =
    [9.999999294887072, 9.999999294887072, 9.999999294887072, 9.999999294887072, 9.999999294887072, 9.999999294887072, 9.999999294887072, 9.999999294887072, 9.999999294887072, 9.999999294887072])
plot(model, [("filters_1.R.v", "filters_1.C.v"), ("filters_2.R.v", "filters_2.C.v")])
=#
end


text:
"""
    Fermi.ConfigurationInteraction.DetOperations.Determinant

structure holding alpha and beta strings for a determinant as integers.

# Fields

    α integer representing an alpha string
    β integer representing an beta string
"""
struct Determinant
    α::Integer
    β::Integer
end

struct CASCI{T} <: AbstractCIWavefunction
    ref::Fermi.HartreeFock.RHF
    energy::T
    dets::Array{Determinant,1}
    coef::Array{T,1}
end

function select_precision(A::String)
    implemented = Dict{Any,Any}("single" => Float32,
                                "double" => Float64)
    try
        return implemented[A]
    catch KeyError
        throw(Fermi.InvalidFermiOptions("Invalid precision: $(A)"))
    end
end

function select_algorithm(A::String)
    implemented = Dict{Any,Any}("sparse" => SparseHamiltonian(),
                                "aci" => ACI())
    try
        return implemented[A]
    catch KeyError
        throw(Fermi.InvalidFermiOptions("Invalid CI algorithm: $(A)"))
    end
end

# Implementations
include("DetOperations.jl")
include("MatrixElement.jl")
include("SparseHamiltonian.jl")
include("ACI.jl")

# Most general function. It defines the precision and call a precision-specific function
"""
    Fermi.CoupledCluster.RCCSD()

Compute a RCCSD wave function using Fermi.CurrentOptions data.
"""
function CASCI()
    @output "Selecting CAS precision...\n"
    prec = select_precision(Fermi.CurrentOptions["precision"])
    CASCI{prec}()
end

"""
    Fermi.CoupledCluster.RCCSD{T}()

Compute a RCCSD wave function for a given precision T (Float64 or Float32)
"""
function CASCI{T}() where T <: AbstractFloat
    @output "Selecting CAS Algorithm...\n"
    alg = select_algorithm(Fermi.CurrentOptions["ci_alg"])
    CASCI{T}(alg)
end


text:
#' ## Model complexity
#'
#' In the last notebook, we saw that we could customize a model by adding a
#' parameter. Doing so, we were able to fit that model to a data point.
#' This fit was perfect, insofar as numerics would allow.
#'
#' In the next notebook, we'll see that as we add more data to our data set,
#' fitting a model to our data usually becomes more challenging and the result
#' will be less perfect.
#'
#' For one thing, we will find that we can add complexity to our model to
#' capture added complexity in the data. We can do this by adding more
#' parameters to our model. We'll see that for a data set with two data points,
#' we can again get a "perfect" fit to our model by adding a second parameter
#' to our model.
#'
#' However, we can't simply add a parameter to our model every time we add a
#' data point to our data set, since this will lead to a phenomenon called
#' **overfitting**.
#'
#' In the image below, we depict a data set that is close to linear, and
#' models that exhibit underfitting, fitting well, and overfitting, from left
#' to right:
#'
#' ![model_fitting](data/model_fitting.png)
#'
#'
#' In the first image, the model accounts for the slope along which the data
#' falls, but not the offset.
#'
#' In the second image, the model accounts for both the slope and offset of
#' the data. Adding this second parameter (the offset) to the model creates
#' a much better fit.
#'
#' However, we can imagine that a model can have too many parameters, where
#' we begin to fit  not only the high level features of the data, but also
#' the noise. This overfitting is depicted in the third image.

#' Our aim will be to fit the data well, but avoiding *over* fitting the data!


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#

# Dummy feature
struct EmptyObject3DFeature <: Modia3D.AbstractObject3DFeature
end
const emptyObject3DFeature = EmptyObject3DFeature()

struct EmptyTwoObject3DObject <: Modia3D.AbstractTwoObject3DObject
end
const emptyTwoObject3DObject = EmptyTwoObject3DObject()

mutable struct InteractionManner
    gripper::Bool
    movable::Bool
    lockable::Bool
    movablePos::Union{Int64,Nothing}
    originPos::Int64
    actualPos::Int64
    InteractionManner() = new(false, false, false, nothing, 0, 0)

    function InteractionManner(interactionBehavior::InteractionBehavior)
        new(false, false, false, nothing, 0, 0)
    end
end

"""
    Object3D(;
        parent      = nothing,
        translation = Modia3D.ZeroVector3D,
        rotation    = Modia3D.ZeroVector3D,
        feature     = nothing)

Generate a new Object3D object, that is a coordinate system (= frame, object3D) with associated feature that is described relatively to its (optional) parent Object3D.

# Arguments

All arguments have default values. Thus, not all arguments must be defined.

- `parent`: Parent Object3D. If `parent` is present, the Object3D is defined relatively to the parent. If `parent` is not present, the Object3D is either a reference object (such as the world-Object3D), or the object is connected later with a joint to another Object3D.

- `translation`: Vector from the origin of the parent to the origin of the Object3D, resolved in the parent coordinate system in [m].
  - Example: `translation = [0.0, 0.5, 0.0]` is a relative translation of 0.5 m in y-direction of the parent. If nothing is defined, the default value (= no translation) is taken.

- `rotation`: Vector `[angleX, angleY, angleZ]` to rotate the parent coordinate system along the x-axis with `angleX`, the y-axis with `angleY` and the z-axis with `angleZ` to arrive at the Object3D coordinate system. Angles can be provided in radians or with unit `u"°"` (degree).
  - Example: `rotation = [0.0, pi/2, 0.0]` or `rotation = [0.0, 90u"°", 0.0]` defines that a rotation around the y-axis of the parent coordinate system with 90 degrees arrives at the Object3D. If nothing is defined, the default value (= no rotation) is taken.

- `feature`: Defines the (optional) property associated with the Object3D by a constructor call. Supported constructors:
    - `Scene`: A [Scene](@ref) feature marks the root (world, origin) Object3D. It has no parent Object3D and allows to define global properties, such as the gravity field.
    - `Visual`: A [Visual](@ref) feature defines a shape used for visualization.
    - `Solid`: A [Solid](@ref) feature defines the solid properties of an Object3D, like mass, inertia tensor, collision behavior.
    - `nothing`: No feature is associated with the Object3D. This might be useful for helper Object3Ds, for example to mark a point on a shape and connecting it later via a revolute joint.

# Example

```julia
model = Model(
    world      = Object3D(feature = Scene()),
    worldFrame = Object3D(parent  = :world,
                          feature = Visual(shape=CoordinateSystem())),
    support    = Object3D(parent  = :world, translation = [0.0, 0.5, 0.0],
                          feature = Visual(shape = Sphere(diameter=0.1)))
)
```
"""
mutable struct Object3D <: Modia3D.AbstractObject3D
    # Tree connection structure of Object3D (for efficient processing of Object3Ds)
    parent::Object3D                 # Parent Object3D (if parent===Object3D, no parent is yet defined)
    children::Vector{Object3D}       # All Object3Ds, where Object3D is the parent
    isRootObj::Bool                  # = true, if it is a root obj of a super obj
    interactionManner::InteractionManner

    # Joint properties, defining the relative motion from parent to Object3D
    joint::Modia3D.AbstractJoint     # ::Fix, ::Revolute, ...

    # Efficient access of joint properties
    jointKind::JointKind             # Kind of joint
    jointIndex::Int                  # scene.<jointType>[jointIndex] = joint
    ndof::Int                        # Number of degrees of freedom
    canCollide::Bool                 # = false, if no collision parent/Object3D
    r_rel::SVector{3,Float64}        # Relative position vector from frame of parent Object3D to origin of Object3D frame, resolved in parent frame in [m]
                                     # (if parent===Object3D, r_rel=Modia3D.ZeroVector3D)
    R_rel::SMatrix{3,3,Float64,9}    # Rotation matrix from frame of parent Object3D to Object3D frame.
    r_abs::SVector{3,Float64}        # Absolute position vector from origin of world frame to origin of Object3D frame, resolved in world frame in [m]
    R_abs::SMatrix{3,3,Float64,9}    # Absolute rotation matrix from world frame to Object3D frame
    v0::SVector{3,Float64}           # Absolute velocity of Object3D origin, resolved in world frame (= der(r_abs)) in [m/s]
    w::SVector{3,Float64}            # Absolute angular velocity of Object3D, resolved in Object3D in [rad/s]
    a0::SVector{3,Float64}           # Absolute acceleration of Object3D origin, resolved in world frame (= der(v0)) in [m/s^2]
    z::SVector{3,Float64}            # Absolute angular acceleration of Object3D, resolved in Object3D in [rad/s^2]
    f::SVector{3,Float64}            # Cut-force resolved in Object3D in [N]
    t::SVector{3,Float64}            # Cut-torque resolved in Object3D in [N*m]

    # Mass properties.
    #   The root of each super object has potentially hasMass=true. All other Object3Ds have hasMass=false.
    #   The initial (fixed) mass properties defined in the ModiaLang model are stored in feature.
    hasMass::Bool                 # = false, if m and I_CM are zero. = true, otherwise.
    m::Float64                    # Mass in [kg]
    r_CM::SVector{3,Float64}      # Position vector from Object3D to Center of Mass resolved in Object3D in [m]
    I_CM::SMatrix{3,3,Float64,9}  # Inertia matrix at Center of Mass resolved in a frame in the center of mass that is parallel to Object3D in [kg*m^2]

    # Additional information associated with Object3D
    feature::Union{Modia3D.AbstractObject3DFeature, Modia3D.AbstractSceneOptions}  # Optional feature associated with Object3D
    twoObject3Dobject::Vector{Modia3D.AbstractTwoObject3DObject}  # Optional AbstractTwoObject3DObject object associated with Object3D
    hasCutJoint::Bool          # = true if it has a cut joint
    hasForceElement::Bool      # = true if it has a force element
    hasChildJoint::Bool        # = true if its child has a joint
    computeAcceleration::Bool  # = true if acceleration needs to be computed

    # internal shortcuts to avoid costly dynamic multiple dispatch
    shapeKind::Shapes.ShapeKind
    shape::Modia3D.AbstractShape
    visualMaterial::Shapes.VisualMaterial

    # = True: Coordinate system of Object3D is always visualized
    # = False: Coordinate system of Object3D is never visualized
    # = Inherited: Coordinate system of Object3D is visualized, if SceneOptions(visualizeFrames=true)
    visualizeFrame::Modia3D.Ternary
    visualizationFrame::Vector{Object3D} # If to be visualized, the Object3D holding the coordinate system.

    # additional 3D Shapes
    fileMeshConvexShapes::Vector{Object3D} # a graphical decomposition of a 3D mesh must be stored additionally

    # if enabled, all AABBs are visualized and stored in world Object3D
    AABBVisu::Vector{Object3D}
    # if enabled, contact and support points are set to world Object3D
    # for visualizing contact points (each contact has two points)
    contactVisuObj1::Vector{Object3D}
    contactVisuObj2::Vector{Object3D}
    # for visualizing support points (each contact has 6 support points, 3 on each collision obj)
    supportVisuObj1A::Vector{Object3D}
    supportVisuObj2A::Vector{Object3D}
    supportVisuObj3A::Vector{Object3D}
    supportVisuObj1B::Vector{Object3D}
    supportVisuObj2B::Vector{Object3D}
    supportVisuObj3B::Vector{Object3D}

    path::String

    ###--------------------- Object3D constructor ------------------------------
    ## Constructor 1
    # constructor for Modia interface (only keyword arguments)
    # calls constructor 2 or 3 or 5
    function Object3D(;
        parent::Union{Object3D,Nothing}=nothing,
        path::String="",
        translation::AbstractVector = Modia3D.ZeroVector3D,
        rotation::Union{AbstractVector,Frames.RotationMatrix,Nothing} = nothing,
        feature::Any=nothing,
        kwargs...)
        if length(kwargs) > 0
            @warn "Object3D: path=$path, kwargs = $(kwargs...)"
        end

        # create Object3D with created feature
        fixed = true
        v_start = Modia3D.ZeroVector3D
        w_start = Modia3D.ZeroVector3D
        w_startVariables = WCartesian

        visualizeFrame = Modia3D.Inherited

        if isnothing(feature)
            feature = emptyObject3DFeature
        else
            feature = feature
        end

        if !isnothing(parent)
            # with parent -> call constructor 3 or 5
            if !isnothing(rotation) && (typeof(rotation) != Frames.RotationMatrix)
                rotation = Frames.rot123(rotation[1], rotation[2], rotation[3])
            end
            obj = Object3D(parent, feature, fixed=fixed, r=translation, R=rotation, v_start=v_start, w_start=w_start, w_startVariables=w_startVariables, visualizeFrame=visualizeFrame, path=path)
        else
            # no parent -> call constructor 2
            obj = Object3D(feature, visualizeFrame=visualizeFrame, path=path)
        end

        if typeof(feature) == Modia3D.Visual && typeof(feature.shape) == Modia3D.FileMesh && feature.shape.convexPartition
            createConvexPartition(obj, feature, feature.shape)
        end

        if typeof(feature) == Modia3D.Solid && typeof(feature.shape) == Modia3D.FileMesh && feature.shape.convexPartition
            createConvexPartition(obj, feature, feature.shape)
        end

        return obj
    end

    ###--------------------- Object3D constructor ------------------------------
    ## Constructor 2
    # Object3D constructor: with feature and without parent (less keywords)
    function Object3D(feature;
                    visualizeFrame::Union{Modia3D.Ternary,Bool} = Modia3D.Inherited, interactionBehavior::InteractionBehavior = Modia3D.NoInteraction,
                    path::String="")::Object3D
        obj = Object3DWithoutParent( new(), visualizeFrame = visualizeFrame, interactionBehavior = interactionBehavior, path=path)

        obj.feature = feature
        (obj.shapeKind, obj.shape, obj.visualMaterial) = setShapeKind(feature)
        return obj
    end

    ###--------------------- Object3D constructor ------------------------------
    ## Constructor 3
    # Object3D constructor: with feature and with parent
    function Object3D(parent::Object3D,
                      feature::Modia3D.AbstractObject3DFeature = emptyObject3DFeature;
                      fixed::Bool = true,
                      interactionBehavior::InteractionBehavior = Modia3D.NoInteraction,
                      r::AbstractVector = Modia3D.ZeroVector3D,
                      R::Union{Frames.RotationMatrix,Nothing} = nothing,
                      q::Union{Frames.Quaternion,Nothing} = nothing,
                      v_start::AbstractVector = Modia3D.ZeroVector3D,
                      w_start::AbstractVector = Modia3D.ZeroVector3D,
                      w_startVariables::WStartVariables = WCartesian,
                      visualizeFrame::Union{Modia3D.Ternary,Bool} = Modia3D.Inherited,
                      path::String="")::Object3D

        if !isnothing(R) && !isnothing(q)
            error("Modia3D.Object3D: either R or q must be nothing but both have a value.")
        end
        if !isnothing(R)
            Modia3D.assertRotationMatrix(R)
        elseif !isnothing(q)
            Modia3D.assertQuaternion(q)
        end

        r_rel = SVector{3,Float64}(r)
        R_rel = !isnothing(R) ? R : (!isnothing(q) ? Modia3D.from_q(q) : Modia3D.NullRotation)
        r_abs = parent.r_abs + r_rel
        R_abs = R_rel*parent.R_abs

        visualizeFrame2 = typeof(visualizeFrame) == Modia3D.Ternary ? visualizeFrame : (visualizeFrame ? Modia3D.True : Modia3D.False)

        (shapeKind, shape, visualMaterial) = setShapeKind(feature)
        obj = new(parent, Vector{Object3D}[],
              false, InteractionManner(interactionBehavior), FixedJoint(), FixKind, 0, 0, true,
              r_rel, R_rel, r_abs, R_abs,
              Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D,
              false, 0.0, Modia3D.ZeroVector3D, SMatrix{3,3,Float64,9}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
              feature, Modia3D.AbstractTwoObject3DObject[],
              false, false, false, false,
              shapeKind, shape, visualMaterial,
              visualizeFrame2,
              Vector{Object3D}[],
              Vector{Object3D}[], Vector{Object3D}[],
              Vector{Object3D}[], Vector{Object3D}[], Vector{Object3D}[], Vector{Object3D}[],
              Vector{Object3D}[], Vector{Object3D}[], Vector{Object3D}[], Vector{Object3D}[],
              path)

        if fixed
            push!(parent.children, obj)

        else
            # FreeMotionCardan pushes obj on parent.children
            rot_start = !isnothing(R) ? rot123fromR(R) : !isnothing(q) ? rot123fromR(Modia3D.from_q(q)) : Modia3D.ZeroVector3D
            if w_startVariables == WCardan123
                # convert Cardan123 -> Cartesian
                om_start = wfromrot123(rot_start, w_start)
            else
                om_start = w_start
            end
            FreeMotionCardan(; r   = r_rel,
                               rot = rot_start,
                               v   = SVector{3,Float64}(v_start),
                               w   = SVector{3,Float64}(om_start) )
        end

        return obj
    end

    ###--------------------- Object3D constructor -----------------------------
    ## Constructor 4
    # Constructor used only for internal purposes (not to be directly used by the user)
    # used by e.g. copyObject3D
    function Object3D(
            parent::Object3D, r_rel::SVector{3,Float64},
            R_rel::SMatrix{3,3,Float64,9}, r_abs::SVector{3,Float64},
            R_abs::SMatrix{3,3,Float64,9}, feature::Modia3D.AbstractObject3DFeature,
            visualizeFrame::Modia3D.Ternary, path::String="")

        (shapeKind, shape, visualMaterial) = setShapeKind(feature)
        new(parent, Vector{Object3D}[], false,
        InteractionManner(Modia3D.NoInteraction), FixedJoint(), FixKind, 0, 0, true,
        r_rel, R_rel, r_abs, R_abs,
        Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D, Modia3D.ZeroVector3D,
        false, 0.0, Modia3D.ZeroVector3D, SMatrix{3,3,Float64,9}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0),
        feature, Modia3D.AbstractTwoObject3DObject[],
        false, false, false, false, shapeKind, shape, visualMaterial, visualizeFrame,
        Vector{Object3D}[],
        Vector{Object3D}[], Vector{Object3D}[],
        Vector{Object3D}[], Vector{Object3D}[], Vector{Object3D}[], Vector{Object3D}[],
        Vector{Object3D}[], Vector{Object3D}[], Vector{Object3D}[], Vector{Object3D}[],
        path)
    end
end

function setShapeKind(feature)
    if typeof(feature) == Modia3D.Solid || typeof(feature) == Modia3D.Visual
        shapeKind = Modia3D.getShapeKind(feature.shape)
        shape = feature.shape
        visualMaterial = feature.visualMaterial
        if shapeKind == Modia3D.UndefinedShapeKind
            shape = Modia3D.Sphere()
        end
        if isnothing(visualMaterial)
            visualMaterial = Modia3D.VisualMaterial()
            @warn("No visualMaterial defined for ", shape)
        end
        return shapeKind, shape, visualMaterial
    else
        return Modia3D.UndefinedShapeKind, Modia3D.Sphere(), Modia3D.VisualMaterial()
    end
end

# Object3DWithoutParent is called from constructor 2 (for objs without a parent)
function Object3DWithoutParent(obj::Object3D;
                               visualizeFrame::Union{Modia3D.Ternary,Bool} = Modia3D.Inherited, interactionBehavior::InteractionBehavior = Modia3D.NoInteraction,
                               path::String="")
    obj.parent             = obj
    obj.children           = Vector{Object3D}[]
    obj.isRootObj          = false
    obj.interactionManner  = InteractionManner(interactionBehavior)
    obj.joint              = FixedJoint()
    obj.jointKind          = FixKind
    obj.jointIndex         = 0
    obj.ndof               = 0
    obj.canCollide         = true
    obj.r_rel              = Modia3D.ZeroVector3D
    obj.R_rel              = Modia3D.NullRotation
    obj.r_abs              = Modia3D.ZeroVector3D
    obj.R_abs              = Modia3D.NullRotation
    obj.v0                 = Modia3D.ZeroVector3D
    obj.w                  = Modia3D.ZeroVector3D
    obj.a0                 = Modia3D.ZeroVector3D
    obj.z                  = Modia3D.ZeroVector3D
    obj.f                  = Modia3D.ZeroVector3D
    obj.t                  = Modia3D.ZeroVector3D
    obj.hasMass            = false
    obj.m                  = 0.0
    obj.r_CM               = Modia3D.ZeroVector3D
    obj.I_CM               = SMatrix{3,3,Float64,9}(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
    obj.twoObject3Dobject  = Modia3D.AbstractTwoObject3DObject[]
    obj.hasCutJoint        = false
    obj.hasForceElement    = false
    obj.hasChildJoint      = false
    obj.computeAcceleration = false
    obj.visualizeFrame     = typeof(visualizeFrame) == Modia3D.Ternary ? visualizeFrame : (visualizeFrame ? Modia3D.True : Modia3D.False)
    obj.visualizationFrame = Vector{Object3D}[]
    obj.fileMeshConvexShapes = Vector{Object3D}[]
    obj.AABBVisu           = Vector{Object3D}[]
    obj.contactVisuObj1    = Vector{Object3D}[]
    obj.contactVisuObj2    = Vector{Object3D}[]
    obj.supportVisuObj1A   = Vector{Object3D}[]
    obj.supportVisuObj2A   = Vector{Object3D}[]
    obj.supportVisuObj3A   = Vector{Object3D}[]
    obj.supportVisuObj1B   = Vector{Object3D}[]
    obj.supportVisuObj2B   = Vector{Object3D}[]
    obj.supportVisuObj3B   = Vector{Object3D}[]
    obj.path               = path
    return obj
end

# addObject3DVisualizationFrame! calls constructor 4
#  - internally used for visualizing a frame by a coordinate system
#    its called: "name of obj" + ".visualizationFrame"
function addObject3DVisualizationFrame!(obj::Object3D,
                                        feature::Modia3D.AbstractObject3DFeature, name)
    push!(obj.visualizationFrame, Object3D(obj.parent,
                                    obj.r_rel, obj.R_rel, obj.r_abs,
                                    obj.R_abs, feature, obj.visualizeFrame) )
end

#    its called: "name of obj" + ".mesh[i]"
# createConvexPartition(obj::Object3D, feature, mesh) = nothing
function createConvexPartition(obj::Object3D, feature, mesh) # feature Visual
    (head,ext) = splitext(mesh.filename)
    if ext == ".obj" && mesh.convexPartition
        convexDecompositionDirectory = joinpath(dirname(mesh.filename),"convexSolids_" * basename(mesh.filename))
        contentDir = readdir(convexDecompositionDirectory)
        # names = Vector{AbstractString}()
        i = 1
        for name in contentDir
            path = String(Symbol(obj.path, ".", "mesh[", i , "]"))
            (head,ext) = splitext(name)
            if ext == ".obj"
                fileMesh = Modia3D.FileMesh(filename = joinpath(convexDecompositionDirectory, name),
                scale = mesh.scaleFactor, useMaterialColor = mesh.useMaterialColor, smoothNormals = mesh.smoothNormals, convexPartition = false)

                feature = createFileFeature(feature, fileMesh)

                push!(obj.fileMeshConvexShapes, Object3D(obj, feature, path = path) )
                i = i + 1
            end
        end
    end
    return nothing
end

function createFileFeature(feature::Shapes.Visual, fileMesh)
    return Modia3D.Visual(shape = fileMesh, visualMaterial = feature.visualMaterial)
end

function createFileFeature(feature::Shapes.Solid, fileMesh)
    return Modia3D.Solid(shape=fileMesh, massProperties=nothing, solidMaterial=feature.solidMaterial, collision=feature.collision, contactMaterial=feature.contactMaterial, collisionSmoothingRadius=feature.collisionSmoothingRadius, visualMaterial=feature.visualMaterial)
end

function addAABBVisuToWorld!(world::Object3D, AABB::Vector{Vector{Basics.BoundingBox}})
    k = 0
    @inbounds for i = 1:length(AABB)
        for j = 1:length(AABB[i])
            k = k + 1
            name = String(Symbol(world.path, ".", "AABBVisu", "[", i, "][",j,"]"))
            aabb = AABB[i][j]
            feature = Modia3D.Visual(shape = Modia3D.Box(
                    lengthX = abs(aabb.x_max - aabb.x_min), lengthY = abs(aabb.y_max - aabb.y_min), lengthZ = abs(aabb.z_max - aabb.z_min)),
                visualMaterial = Modia3D.VisualMaterial(color="grey96", transparency=0.8))
            push!(world.AABBVisu,  Object3D(world, feature, path = name) )
        end
    end

end

function addContactVisuObjToWorld!(world::Object3D, nVisualContSupPoints, defaultContactSphereDiameter)
    world.contactVisuObj1 = fill(Object3D(), nVisualContSupPoints)
    world.contactVisuObj2 = fill(Object3D(), nVisualContSupPoints)
    @inbounds for i = 1:length(world.contactVisuObj1)
        name1 = String(Symbol(world.path, ".", "contactVisuObj1", "[", i, "]"))
        name2 = String(Symbol(world.path, ".", "contactVisuObj2", "[", i, "]"))

        feature1 = Modia3D.Visual(shape = Modia3D.Sphere(diameter = defaultContactSphereDiameter), visualMaterial = Modia3D.VisualMaterial(color="Red",   transparency=1.0))
        feature2 = Modia3D.Visual(shape = Modia3D.Sphere(diameter = defaultContactSphereDiameter), visualMaterial = Modia3D.VisualMaterial(color="Black",   transparency=1.0))

        world.contactVisuObj1[i] =  Object3D(world, feature1, path = name1)
        world.contactVisuObj2[i] =  Object3D(world, feature2, path = name2)
    end
end

function addSupportVisuObjToWorld!(world::Object3D, nVisualContSupPoints, defaultContactSphereDiameter)
    world.supportVisuObj1A = fill(Object3D(), nVisualContSupPoints)
    world.supportVisuObj2A = fill(Object3D(), nVisualContSupPoints)
    world.supportVisuObj3A = fill(Object3D(), nVisualContSupPoints)
    world.supportVisuObj1B = fill(Object3D(), nVisualContSupPoints)
    world.supportVisuObj2B = fill(Object3D(), nVisualContSupPoints)
    world.supportVisuObj3B = fill(Object3D(), nVisualContSupPoints)
    @inbounds for i = 1:length(world.supportVisuObj1A)
        name1 = String(Symbol(world.path, ".", "supportVisuObj1A", "[", i, "]"))
        name2 = String(Symbol(world.path, ".", "supportVisuObj2A", "[", i, "]"))
        name3 = String(Symbol(world.path, ".", "supportVisuObj3A", "[", i, "]"))
        name4 = String(Symbol(world.path, ".", "supportVisuObj1B", "[", i, "]"))
        name5 = String(Symbol(world.path, ".", "supportVisuObj2B", "[", i, "]"))
        name6 = String(Symbol(world.path, ".", "supportVisuObj3B", "[", i, "]"))

        featureA = Modia3D.Visual(shape = Modia3D.Sphere(diameter = defaultContactSphereDiameter), visualMaterial = Modia3D.VisualMaterial(color="Red", transparency=1.0))

        featureB = Modia3D.Visual(shape = Modia3D.Sphere(diameter = defaultContactSphereDiameter), visualMaterial = Modia3D.VisualMaterial(color="Black", transparency=1.0))

        world.supportVisuObj1A[i] = Object3D(world, featureA, path = name1)
        world.supportVisuObj2A[i] = Object3D(world, featureA, path = name2)
        world.supportVisuObj3A[i] = Object3D(world, featureA, path = name3)
        world.supportVisuObj1B[i] = Object3D(world, featureB, path = name4)
        world.supportVisuObj2B[i] = Object3D(world, featureB, path = name5)
        world.supportVisuObj3B[i] = Object3D(world, featureB, path = name6)
    end
end

# Inquire properties of a Object3D
fullName(             obj::Object3D) = obj.path # ModiaLang.SimulationModel.modelName + Object3D.path (wird ca. 98x verwendet)
instanceName(         obj::Object3D) = obj.path
fullName(obj) = error("fullName not implemented for ", typeof(obj))
instanceName(obj) = error("instanceName not implemented for ", typeof(obj))
hasParent(            obj::Object3D) = !(obj.parent === obj)
hasNoParent(          obj::Object3D) =   obj.parent === obj
hasChildren(          obj::Object3D) = length(obj.children) > 0
hasNoChildren(        obj::Object3D) = length(obj.children) == 0
isWorld(              obj::Object3D) = hasNoParent(obj) && typeof(obj.feature) == Modia3D.SceneOptions
isNotWorld(           obj::Object3D) = !(isWorld(obj))
isMovable(            obj::Object3D) = obj.interactionManner.movable
isLockable(           obj::Object3D) = obj.interactionManner.lockable
isFixed(              obj::Object3D) = obj.ndof == 0
isNotFixed(           obj::Object3D) = obj.ndof > 0
isFree(               obj::Object3D) = obj.ndof == 6
isNotFree(            obj::Object3D) = obj.ndof < 6
hasJoint(             obj::Object3D) = obj.ndof > 0 && obj.ndof < 6
hasNoJoint(           obj::Object3D) = isFixed(obj) || isFree(obj)
isCoordinateSystem(   obj::Object3D) = typeof(obj.feature) == Shapes.CoordinateSystem
isNotCoordinateSystem(obj::Object3D) = typeof(obj.feature) != Shapes.CoordinateSystem
hasForceElement(      obj::Object3D) = obj.hasForceElement
hasChildJoint(        obj::Object3D) = obj.hasChildJoint
needsAcceleration(    obj::Object3D) = obj.computeAcceleration
objectHasMass(        obj::Object3D) = obj.hasMass
isRootObject(         obj::Object3D) = obj.isRootObj
objectHasMovablePos(  obj::Object3D) = !isnothing(obj.interactionManner.movablePos)

featureHasMass(obj::Object3D) = featureHasMass(obj.feature)
featureHasMass(feature::Modia3D.AbstractObject3DFeature) = false
featureHasMass(feature::Modia3D.AbstractSceneOptions) = false
featureHasMass(feature::Shapes.Solid)                 = !isnothing(feature.massProperties)

isVisible(obj::Object3D, renderer::Modia3D.AbstractRenderer) = isVisible(obj.feature, renderer)
isVisible(feature::Modia3D.AbstractObject3DFeature, renderer::Modia3D.AbstractRenderer) = false

isVisible(feature::Modia3D.AbstractSceneOptions, renderer::Modia3D.AbstractRenderer) = false

canCollide(feature::Modia3D.AbstractObject3DFeature) = false

function canCollide(obj::Object3D)::Bool
    if typeof(obj.feature) == Modia3D.Shapes.Solid
        if !isnothing(obj.feature.shape) && !isnothing(obj.feature.contactMaterial)
            if typeof(obj.feature.contactMaterial) == String
                return obj.feature.contactMaterial != ""
            else
                return true
    end; end; end
    return false
end

function driveJoint!(obj::Object3D)::Nothing
    obj.isDriven = true
    return nothing
end

"""    rootObject3D(obj) - returns the root Object3D of all parents of obj"""
function rootObject3D(obj::Object3D)::Object3D
    obj1 = obj
    while hasParent(obj1)
        obj1 = obj1.parent
    end
    return obj1
end

"""    removeChild!(obj, child) - Remove child from obj.children"""
function removeChild!(obj::Object3D, child::Object3D)::Nothing
    children = obj.children
    for i in eachindex(children)
        if children[i] === child
            deleteat!(children,i)
            return nothing
    end; end
    error("\nError from Modia3D.Composition.removeChild!(", obj, ", ", child, ")\n",
            child, " is not a child of object ", obj)
end

# Print Object3D
function Base.show(io::IO, obj::Object3D)
    print(io,"Object3D(path=\"", obj.path, "\"")
    commaNeeded = true

    if hasParent(obj)
        if commaNeeded
            print(io,", ")
        end
        print(io, "parent=", obj.parent.path)
    end
    if hasJoint(obj)
        if commaNeeded
            print(io,", ")
        end
        print(io, "joint=", typeof(obj.joint), "(path = \"", obj.joint.path, "\", ...)")
    end
    if commaNeeded
        print(io,", ")
    end
    print(io, "feature=", typeof(obj.feature), "(...))")
end
#=
function Base.show(io::IO, obj::Object3D)
    commaNeeded = true
    if obj.parent === obj || hasJoint(obj)
        commaNeeded = false
        print(io,"Object3D(")
    else
        if obj.parent.path == ""
            print(io,"Object3D(", Modia3D.fullName(obj.parent))
        else
            print(io,"Object3D(", obj.parent.path)
        end
    end

    if typeof(obj.feature) != EmptyObject3DFeature
        if commaNeeded
            print(io,", ")
        end
        typeName = string( typeof(obj.feature) )
        s = findfirst("Modia3D", typeName) # search(typeName, "Modia3D")
        if s.start == 1
            typeName = "Modia3D." * Basics.trailingPartOfName( typeName )
        end
        #print(io, typeName, "(...)")
        print(io, typeName, "(")
        Base.show(io, obj.feature)
        print(io, ")")
    end

    if obj.r_rel != Modia3D.ZeroVector3D
        print(io, "; r_rel = ", obj.r_rel)
        if obj.path == ""
            print(io, ")")
        else
            print(io,", path = \"", obj.path, "\")")
        end
    else
        if obj.path == ""
            print(io, ")")
        else
            print(io, "; path = \"", obj.path, "\")")
        end
    end
    if obj.path == ""
        print(io, " # ", Modia3D.instanceName(obj), "\n")
    end
end
=#


text:
"""
    articulation(g)

Compute the [articulation points](https://en.wikipedia.org/wiki/Biconnected_component)
of a connected graph `g` and return an array containing all cut vertices.

# Examples
```jldoctest
julia> using LightGraphs

julia> articulation(star_graph(5))
1-element Array{Int64,1}:
 1

julia> articulation(path_graph(5))
3-element Array{Int64,1}:
 2
 3
 4
```
"""
function articulation end
@traitfn function articulation(g::AG::(!IsDirected)) where {T, AG<:AbstractGraph{T}}
    s = Vector{Tuple{T, T, T}}()
    is_articulation_pt = falses(nv(g))
    low = zeros(T, nv(g))
    pre = zeros(T, nv(g))

    @inbounds for u in vertices(g)
        pre[u] != 0 && continue
        v = u
        children = 0
        wi::T = zero(T)
        w::T = zero(T)
        cnt::T = one(T)
        first_time = true

        while !isempty(s) || first_time
            first_time = false
            if  wi < 1
                pre[v] = cnt
                cnt += 1
                low[v] = pre[v]
                v_neighbors = outneighbors(g, v)
                wi = 1
            else
                wi, u, v = pop!(s)
                v_neighbors = outneighbors(g, v)
                w = v_neighbors[wi]
                low[v] = min(low[v], low[w])
                if low[w] >= pre[v] && u != v
                    is_articulation_pt[v] = true
                end
                wi += 1
            end
            while wi <= length(v_neighbors)
                w = v_neighbors[wi]
                if pre[w] == 0
                    if u == v
                        children += 1
                    end
                    push!(s, (wi, u, v))
                    wi = 0
                    u = v
                    v = w
                    break
                elseif w != u
                    low[v] = min(low[v], pre[w])
                end
                wi += 1
            end
            wi < 1 && continue
        end

        if children > 1
            is_articulation_pt[u] = true
        end
    end

    articulation_points = Vector{T}()

    for u in findall(is_articulation_pt)
        push!(articulation_points, T(u))
    end

    return articulation_points
end


text:
const binary_operators = Set()

"""
    BinaryOperation{X, Y, Z, O, A, B, IA, IB, IΩ, G} <: AbstractOperation{X, Y, Z, G}

An abstract representation of a binary operation on `AbstractField`s.
"""
struct BinaryOperation{X, Y, Z, O, A, B, IA, IB, IΩ, G} <: AbstractOperation{X, Y, Z, G}
      op :: O
       a :: A
       b :: B
      ▶a :: IA
      ▶b :: IB
     ▶op :: IΩ
    grid :: G

    """
        BinaryOperation{X, Y, Z}(op, a, b, ▶a, ▶b, ▶op, grid)

    Returns an abstract representation of the binary operation `op(▶a(a), ▶b(b))`,
    followed by interpolation by `▶op` to `(X, Y, Z)`, where `▶a` and `▶b` interpolate
    `a` and `b` to a common location.
    """
    function BinaryOperation{X, Y, Z}(op, a, b, ▶a, ▶b, ▶op, grid) where {X, Y, Z}

        any((X, Y, Z) .=== Nothing) && throw(ArgumentError("Nothing locations are invalid! " *
                                                           "Cannot construct BinaryOperation at ($X, $Y, $Z)."))

        return new{X, Y, Z, typeof(op), typeof(a), typeof(b), typeof(▶a), typeof(▶b),
                   typeof(▶op), typeof(grid)}(op, a, b, ▶a, ▶b, ▶op, grid)
    end
end

@inline Base.getindex(β::BinaryOperation, i, j, k) = β.▶op(i, j, k, β.grid, β.op, β.▶a, β.▶b, β.a, β.b)

#####
##### BinaryOperation construction
#####

"""Create a binary operation for `op` acting on `a` and `b` with locations `La` and `Lb`.
The operator acts at `Lab` and the result is interpolated to `Lc`."""
function _binary_operation(Lc, op, a, b, La, Lb, Lab, grid)
     ▶a = interpolation_operator(La, Lab)
     ▶b = interpolation_operator(Lb, Lab)
    ▶op = interpolation_operator(Lab, Lc)
    return BinaryOperation{Lc[1], Lc[2], Lc[3]}(op, a, b, ▶a, ▶b, ▶op, grid)
end

"""Return an expression that defines an abstract `BinaryOperator` named `op` for `AbstractField`."""
function define_binary_operator(op)
    return quote
        import Oceananigans.Grids: AbstractGrid
        import Oceananigans.Fields: AbstractField

        local location = Oceananigans.Fields.location
        local FunctionField = Oceananigans.Fields.FunctionField
        local AF = AbstractField

        @inline $op(i, j, k, grid::AbstractGrid, ▶a, ▶b, a, b) =
            @inbounds $op(▶a(i, j, k, grid, a), ▶b(i, j, k, grid, b))

        """
            $($op)(Lc, Lab, a, b)

        Returns an abstract representation of the operator `$($op)` acting on `a` and `b` at
        location `Lab`, and subsequently interpolated to location `Lc`.
        """
        function $op(Lc::Tuple, Lop::Tuple, a, b)
            La = location(a)
            Lb = location(b)
            grid = Oceananigans.AbstractOperations.validate_grid(a, b)
            return Oceananigans.AbstractOperations._binary_operation(Lc, $op, a, b, La, Lb, Lop, grid)
        end

        $op(Lc::Tuple, a, b) = $op(Lc, Lc, a, b)
        $op(Lc::Tuple, a::Number, b) = $op(Lc, location(b), a, b)
        $op(Lc::Tuple, a, b::Number) = $op(Lc, location(a), a, b)
        $op(Lc::Tuple, a::AF{X, Y, Z}, b::AF{X, Y, Z}) where {X, Y, Z} = $op(Lc, location(a), a, b)

        # Sugar for mixing in functions of (x, y, z)
        $op(Lc::Tuple, a::Function, b::AbstractField) = $op(Lc, FunctionField(Lc, a, b.grid), b)
        $op(Lc::Tuple, a::AbstractField, b::Function) = $op(Lc, a, FunctionField(Lc, b, a.grid))

        # Sugary versions with default locations
        $op(a::AF, b::AF) = $op(location(a), a, b)
        $op(a::AF, b::Number) = $op(location(a), a, b)
        $op(a::Number, b::AF) = $op(location(b), a, b)

        $op(a::AF, b::Function) = $op(location(a), a, FunctionField(location(a), b, a.grid))
        $op(a::Function, b::AF) = $op(location(b), FunctionField(location(b), a, b.grid), b)
    end
end

"""
    @binary op1 op2 op3...

Turn each binary function in the list `(op1, op2, op3...)`
into a binary operator on `Oceananigans.Fields` for use in `AbstractOperations`.

Note: a binary function is a function with two arguments: for example, `+(x, y)` is a binary function.

Also note: a binary function in `Base` must be imported to be extended: use `import Base: op; @binary op`.

Example
=======

```jldoctest
julia> using Oceananigans, Oceananigans.AbstractOperations, Oceananigans.Grids

julia> plus_or_times(x, y) = x < 0 ? x + y : x * y
plus_or_times (generic function with 1 method)

julia> @binary plus_or_times
6-element Array{Any,1}:
 :+
 :-
 :/
 :^
 :*
 :plus_or_times

julia> c, d = (Field(Center, Center, Center, CPU(), RegularRectilinearGrid(size=(1, 1, 1), extent=(1, 1, 1))) for i = 1:2);

julia> plus_or_times(c, d)
BinaryOperation at (Center, Center, Center)
├── grid: RegularRectilinearGrid{Float64, Periodic, Periodic, Bounded}(Nx=1, Ny=1, Nz=1)
│   └── domain: x ∈ [0.0, 1.0], y ∈ [0.0, 1.0], z ∈ [-1.0, 0.0]
└── tree:
    plus_or_times at (Center, Center, Center) via identity
    ├── Field located at (Center, Center, Center)
    └── Field located at (Center, Center, Center)
"""
macro binary(ops...)
    expr = Expr(:block)

    for op in ops
        defexpr = define_binary_operator(op)
        push!(expr.args, :($(esc(defexpr))))

        add_to_operator_lists = quote
            push!(Oceananigans.AbstractOperations.operators, Symbol($op))
            push!(Oceananigans.AbstractOperations.binary_operators, Symbol($op))
        end

        push!(expr.args, :($(esc(add_to_operator_lists))))
    end

    return expr
end

#####
##### Architecture inference for BinaryOperation
#####

architecture(β::BinaryOperation) = architecture(β.a, β.b)

function architecture(a, b)
    arch_a = architecture(a)
    arch_b = architecture(b)

    arch_a === arch_b && return arch_a
    isnothing(arch_a) && return arch_b
    isnothing(arch_b) && return arch_a

    throw(ArgumentError("Operation involves fields on architectures $arch_a and $arch_b"))

    return nothing
end

#####
##### Nested computations
#####

function compute_at!(β::BinaryOperation, time)
    compute_at!(β.a, time)
    compute_at!(β.b, time)
    return nothing
end

#####
##### GPU capabilities
#####

"Adapt `BinaryOperation` to work on the GPU via CUDAnative and CUDAdrv."
Adapt.adapt_structure(to, binary::BinaryOperation{X, Y, Z}) where {X, Y, Z} =
    BinaryOperation{X, Y, Z}(Adapt.adapt(to, binary.op), Adapt.adapt(to, binary.a),  Adapt.adapt(to, binary.b),
                             Adapt.adapt(to, binary.▶a), Adapt.adapt(to, binary.▶b), Adapt.adapt(to, binary.▶op),
                             binary.grid)


text:
"""
Stuctured Mesh Generation package.
Provides tools for structured mesh definition
generation and export to abaqus.

Copyright 2019 Jonathan Hahn

Licensed under MIT License, see LICENSE.md
"""
module StructuredMeshing

using LinearAlgebra

import Base.angle
import Base.show

# types

export MeshDef, Mesh

# functions

## mesh definition
export emptyMeshDef, addVertex!, extrude!, transitionextrude!, defineCartesian, connect!, rotate!

## mesh generation
export mesh

## abaqus export
export writeabq

## visualisation
export show

## debugging
export blockdef2fun2, blockdef2fun

# include

include("definition.jl")
include("transfinite.jl")
include("generate.jl")
include("show.jl")
include("helper.jl")
include("abaqus.jl")

end


text:
using ExcelReaders
using MultivariateStats
using Plots
using StatsBase

#Leitura dos dados
dados = readxlsheet("WDI01.xlsx","Data", "D1:BK5")

# Criação das variáveis
gbr = dados[:,1];
gmx = dados[:,2];
cbr = gbr.*dados[:,3];
cmx = gmx.*dados[:,4];

# Transformação logarítmica
gbr = log.(gbr);
gmx = log.(gmx);
cbr = log.(cbr);
cmx = log.(cmx);

# Séries para estimação

T=length(gbr);
ano = LinRange(1960,2017,T);
time = LinRange(1,58,58);
time2 = time.^2;
X = [ones(58) time time2];   # tendência quadrática
# X = [ones(58) time];       # tendência linear

# Regressões

y = gbr;
a = llsq(X, y; bias=false);
gbrs = X * a;
gbrc = gbr - gbrs;

y = gmx;
a = llsq(X, y; bias=false);
gmxs = X * a;
gmxc = gmx - gmxs;

y = cbr;
a = llsq(X, y; bias=false);
cbrs = X * a;
cbrc = cbr - cbrs;

y = cmx;
a = llsq(X, y; bias=false);
cmxs = X * a;
cmxc = cmx - cmxs;

# Gráficos

p1 = plot(ano, [gbr gbrs], title="BRA-PIB")

p2 = plot(ano, gbrc, title="BRA-PIB-ciclo")

p3 = plot(ano, [cbr cbrs], title="BRA-Con")

p4 = plot(ano, cbrc, title="BRA-Con-ciclo")

p5 = plot(ano, [gmx gmxs], title="MEX-PIB")

p6 = plot(ano, gmxc,  title="MEX-PIB-ciclo")

p7 = plot(ano, [cmx cmxs], title="MEX-Con")

p8 = plot(ano, cmxc, title="MEX-Con-ciclo")

plot(p1,p2,p3,p4,p5,p6,p7,p8, layout=(4,2), legend=false)


text:
include("lib/component_types.jl")
include("lib/cartesian_tuple.jl")
include("lib/aliases.jl")
include("lib/tuple_behavior.jl")
include("lib/point_behavior.jl")
include("lib/vector_behavior.jl")
include("lib/conversion.jl")


text:
num_machines() = static(1)
num_sockets() = static(1)

function _get_num_threads()::Int
  (Sys.CPU_THREADS)::Int >> (Sys.ARCH !== :aarch64)
end

const _get_num_cores = _get_num_threads

let nc = static(_get_num_threads())
  global num_l1cache() = nc
  global num_cores() = nc
end
let syst = static((Sys.CPU_THREADS)::Int)
  global sys_threads() = syst
  global num_threads() = syst
end
@static if Sys.ARCH === :aarch64
  num_l2cache() = static(1)
  num_l3cache() = static(0)
else
  num_l2cache() = num_l1cache()
  num_l3cache() = static(1)
end
num_l4cache() = static(0)

if Sys.ARCH === :aarch64 && Sys.isapple()
  cache_size(::Union{Val{1},StaticInt{1}}) = StaticInt{131072}()
else
  cache_size(::Union{Val{1},StaticInt{1}}) = StaticInt{32768}()
end
cache_associativity(::Union{Val{1},StaticInt{1}}) = StaticInt{0}()
cache_type(::Union{Val{1},StaticInt{1}}) = Val{:Data}()
cache_inclusive(::Union{Val{1},StaticInt{1}}) = False()

if Sys.ARCH === :aarch64 && Sys.isapple()
  cache_size(::Union{Val{2},StaticInt{2}}) = StaticInt{12582912}()
else
  cache_size(::Union{Val{2},StaticInt{2}}) = StaticInt{262144}()
end
cache_associativity(::Union{Val{2},StaticInt{2}}) = StaticInt{0}()
cache_type(::Union{Val{2},StaticInt{2}}) = Val{:Unified}()
cache_inclusive(_) = False()
@static if Sys.isapple() && Sys.ARCH === :aarch64
  cache_linesize(_) = StaticInt{128}() # assume...
else
  cache_linesize(_) = StaticInt{64}() # assume...
end
cache_size(_) = StaticInt{0}()

cache_type(::Union{Val{3},StaticInt{3}}) = Val{:Unified}()
cache_size(::Union{Val{3},StaticInt{3}}) = num_cores() * StaticInt{1441792}()

_extra_init() = nothing


text:
using Documenter, RAMS

makedocs(
    sitename="RAMS.jl",
    pages = [
        "Index" => "index.md",
        "API Reference" => "apiref.md"
    ]
)

deploydocs(
    repo="github.com/lsterzinger/RAMS.jl.git",
    devbranch="main"
)


text:
using Images, TestImages, StaticArrays, ImageTracking, Random, OffsetArrays, CoordinateTransformations, BenchmarkTools

img1 = Gray{Float64}.(testimage("mandrill"))
trans = Translation(0, -3.0)
img2 = warp(img1, trans, axes(img1))

corners = imcorner(img1, method=shi_tomasi)
I = findall(!iszero, corners)
r, c = (getindex.(I, 1), getindex.(I, 2))
points = map((ri, ci) -> SVector{2}(Float64(ri), Float64(ci)), r, c)

flow, indicator  = optical_flow(img1, img2, points, LucasKanade( 11, 4, 20, 0.000001))

@btime flow, indicator  =  optical_flow($img1, $img2, $points, LucasKanade( 11, 4, 20, 0.000001))

(sum(indicator)/ length(indicator)) * 100

@code_warntype  optical_flow(img1, img2, points, LucasKanade( 11, 4, 20, 0.000001))


text:
using Test
using ScottishTaxBenefitModel
using .ModelHousehold:
    Household,
    Person,
    People_Dict,
    default_bu_allocation,
    get_benefit_units,
    get_head,
    get_spouse,
    is_single,
    num_children,
    pers_is_carer,
    pers_is_disabled,
    search

using .GeneralTaxComponents:
    WEEKS_PER_MONTH,
    WEEKS_PER_YEAR

using .ExampleHouseholdGetter
using .Definitions

using .TimeSeriesUtils: fy_from_bits
using .HistoricBenefits:
    RATIO_BENS,
    make_benefit_ratios!

using .STBIncomes
using .Intermediate:
    MTIntermediate,
    apply_2_child_policy,
    make_intermediate

using .NonMeansTestedBenefits:
    calc_widows_benefits,
    calc_state_pension,
    calc_child_benefit!,
    calc_pip,
    calc_attendance_allowance,
    calc_dla,
    calc_esa,
    calc_maternity_allowance,
    calc_carers_allowance,
    calc_jsa,
    calc_pre_tax_non_means_tested!,
    calc_post_tax_non_means_tested!

using .IncomeTaxCalculations:
    calc_income_tax!

using .STBParameters:
    AttendanceAllowance,
    ChildBenefit,
    DisabilityLivingAllowance,
    CarersAllowance,
    PersonalIndependencePayment,
    ContributoryESA,
    WidowsPensions,
    BereavementSupport,
    RetirementPension,
    JobSeekersAllowance,
    NonMeansTestedSys

using .Results:
    LMTResults,
    BenefitUnitResult,
    init_household_result

using .ExampleHelpers

using Dates

## FIXME don't need both
sys = get_system( scotland=true )
ruksys = get_system( scotland=false )
itsys_scot = get_default_it_system( year=2019, scotland=true )

@testset "CB" begin
    cb = sys.nmt_bens.child_benefit
    sph = get_example( single_parent_hh )
    hhres = init_household_result( sph )
    bu = get_benefit_units( sph )[1]
    head = get_head( bu )
    hp = head.pid
    bures = hhres.bus[1]
    @test num_children( bu ) == 2
    calc_child_benefit!( bures, bu, cb )
    @test bures.pers[hp].income[GUARDIANS_ALLOWANCE] == 0
    @test bures.pers[hp].income[CHILD_BENEFIT] ≈ 34.40
    # one of two is now non-related
    head.relationships[320190010302] = Other_relative
    calc_child_benefit!( bures, bu, cb )
    @test bures.pers[hp].income[GUARDIANS_ALLOWANCE] == 17.20
    @test bures.pers[hp].income[CHILD_BENEFIT] ≈ 34.40
    # cap - 2k over so £20 pw reduction
    head.income = Incomes_Dict{Float64}() # clear other incomes
    head.income[wages] = 52_000/WEEKS_PER_YEAR

    hhres = init_household_result( sph )
    bures = hhres.bus[1]
    # for total income calculation
    calc_income_tax!(
        bures,
        head,
        nothing,
        itsys_scot )

    calc_child_benefit!( bures, bu, cb )
    @test bures.pers[hp].income[GUARDIANS_ALLOWANCE] == 17.20
    @test bures.pers[hp].income[CHILD_BENEFIT] ≈ (34.40 - 20)
end

@testset "PIP" begin
    pip = sys.nmt_bens.pip
    sph = get_example( single_parent_hh )
    head = get_head( sph )
    disable_seriously!( head )
    head.pip_daily_living_type = standard_pip
    head.pip_mobility_type = enhanced_pip
    pl,pm = calc_pip( head, pip )
    @test (pl+pm) ≈ 119.90
end

@testset "AA" begin
    aa = sys.nmt_bens.attendance_allowance
    sph = get_example( single_parent_hh )
    head = get_head( sph )
    disable_seriously!( head )
    head.attendance_allowance_type = high
    a = calc_attendance_allowance( head, aa )
    @test a ≈ 87.65
end

@testset "DLA" begin
    dla = sys.nmt_bens.dla
    sph = get_example( single_parent_hh )
    head = get_head( sph )
    disable_seriously!( head )
    head.dla_self_care_type = mid
    dd,dm = calc_dla( head, dla )
    @test dd ≈ 58.70
    @test dm == 0
    head.dla_mobility_type = high
    dd,dm = calc_dla( head, dla )
    @test dd ≈ 58.70
    @test dm ≈ 61.20
end

@testset "ESA" begin
    esa = sys.nmt_bens.esa
    sph = get_example( single_parent_hh )
    head = get_head( sph )
    head.esa_type = contributory_jsa
    e = calc_esa( head, esa )
    @test e ≈ 73.10
    head.age = 22
    e = calc_esa( head, esa )
    @test e ≈ 57.90
    # should put them in `has_limited_capactity_for_work_activity`
    disable_seriously!( head )
    head.age = 50
    e = calc_esa( head, esa )
    @test e ≈ 73.10 + 38.55
end

@testset "Maternity" begin
    mat = sys.nmt_bens.maternity
    sph = get_example( single_parent_hh )
    head = get_head( sph )
    head.income[maternity_allowance] = 123
    m = calc_maternity_allowance( head, mat )
    @test m ≈ 148.68

end

@testset "Carers" begin
    care = sys.nmt_bens.carers
    sph = get_example( single_parent_hh )
    head = get_head( sph )
    hp = head.pid
    unemploy!( head )
    head.hours_of_care_given = 60
    hhres = init_household_result( sph )
    bures = hhres.bus[1]
    # for total income calculation
    calc_income_tax!(
        bures,
        head,
        nothing,
        itsys_scot )
    pres = bures.pers[hp]
    c = calc_carers_allowance( head, pres, care )
    @test c ≈ 66.15
    employ!( head )
    hhres = init_household_result( sph )
    bures = hhres.bus[1]
    calc_income_tax!(
        bures,
        head,
        nothing,
        itsys_scot )
    pres = bures.pers[hp]
    c = calc_carers_allowance( head, pres, care )
    @test c == 0.0
end

@testset "JSA" begin
    sph = get_example( single_parent_hh )
    head = get_head( sph )
    hp = head.pid
    unemploy!( head )
    head.jsa_type = contributory_jsa
    head.age = 23
    j = calc_jsa( head, sys.nmt_bens.jsa, sys.hours_limits )
    @test j ≈ 57.90
    head.age = 50
    j = calc_jsa( head, sys.nmt_bens.jsa, sys.hours_limits  )
    @test j ≈ 73.10
    employ!( head )
    j = calc_jsa( head, sys.nmt_bens.jsa, sys.hours_limits  )
    @test j ≈ 0.0
end

@testset "Widows" begin
    bp = sys.nmt_bens.bereavement
    wp = sys.nmt_bens.widows_pension
    sph = get_example( single_parent_hh )
    fy = fy_from_bits( sph.interview_year, sph.interview_month )
    @test fy == 2019

    bu = get_benefit_units( sph )[1]
    head = get_head( bu )
    # old style
    head.income[bereavement_allowance_or_widowed_parents_allowance_or_bereavement] = 100.0
    head.bereavement_type = widowed_parents
    make_benefit_ratios!( head, sph.interview_year, sph.interview_month )

    p = calc_widows_benefits( head, true, bp, wp )
    @test p ≈ wp.standard_rate*100.0/119.9 # sys.widows_pension.standard_rate
    head.bereavement_type = bereavement_allowance
    p = calc_widows_benefits( head, true, bp, wp )
    @test p ≈ bp.lump_sum_higher*2/3 + bp.higher

end

@testset "State Pension" begin
    rp = sys.nmt_bens.pensions
    hh = get_example( single_hh )
    fy = fy_from_bits( hh.interview_year, hh.interview_month )
    bu = get_benefit_units( hh )[1]
    head = get_head( bu )
    head.income[state_pension] = 100.0
    make_benefit_ratios!( head, hh.interview_year, hh.interview_month )
    head.age = age_now(65)
    # female over pension age but too young to have been on old pension
    p = calc_state_pension(
        head,
        rp,
        sys.age_limits )
    @test p == rp.new_state_pension
    head.age = age_now(60)
    # under pension age
    p = calc_state_pension(
        head,
        rp,
        sys.age_limits )
    @test p == 0 # under pension age
    head.age = age_now(70)
    p = calc_state_pension(
        head,
        rp,
        sys.age_limits )
    @test p ≈ rp.cat_a*100.0/129.2
end

@testset "Pre-Tax NMT" begin
    for exn in instances( SS_Examples )
        hh = get_example( exn )
        println( "on $exn")
        hhres = init_household_result( hh )
        calc_pre_tax_non_means_tested!(
            hhres, # :: HouseholdResult,
            hh,    #    :: Household,
            sys.nmt_bens, #   :: NonMeansTestedSys,
            sys.hours_limits, #  :: HoursLimits,
            sys.age_limits ) # :: AgeLimits )
    end
end

@testset "Post-Tax NMT" begin
    for exn in instances( SS_Examples )
        hh = get_example( exn )
        println( "on $exn")
        hhres = init_household_result( hh )
        calc_post_tax_non_means_tested!(
            hhres, # :: HouseholdResult,
            hh,    #    :: Household,
            sys.nmt_bens, #   :: NonMeansTestedSys,
            sys.age_limits ) # :: AgeLimits )
    end
end


text:
module TestStructArrays

using Test
using BangBang
using StructArrays

@testset "push!!" begin
    xs = StructVector{ComplexF64}(([1.0, 2.0], [3.0, 4.0]))
    ys = push!!(xs, 5 + 6im)
    @test ys == [1 + 3im, 2 + 4im, 5 + 6im]
    @test_broken xs === ys
    xs = StructVector{Complex{Int}}(([1, 2], [3, 4]))
    ys = push!!(xs, 5.0 + 6.0im)
    @test ys == [1 + 3im, 2 + 4im, 5 + 6im]
    @test_broken eltype(ys) == ComplexF64
end

@testset "append!!" begin
    xs = StructVector{ComplexF64}(([1.0, 2.0], [3.0, 4.0]))
    ys = append!!(xs, [5 + 6im])
    @test ys == [1 + 3im, 2 + 4im, 5 + 6im]
    @test_broken xs === ys
    xs = StructVector{Complex{Int}}(([1, 2], [3, 4]))
    ys = append!!(xs, 5.0 + 6.0im)
    @test ys == [1 + 3im, 2 + 4im, 5 + 6im]
    @test_broken eltype(ys) == ComplexF64
end

@testset "type inference" begin
    @test (@inferred singletonof(StructVector, (a=1,))) == StructVector((a=[1]))
end

end  # module


text:
# TODO: check trees are on the same leave set

"""
    RF(tree1::T, tree2::T)::Int64 where T <:GeneralNode

Calculate the Robinson-Foulds distance between the two trees.
In its current form the function assumes the trees have identical leave sets.

Returns result of algorithm as integer.

* `tree1` : tree used to determine RF distance.

* `tree2` : tree used to determine RF distance.
"""
function RF(tree1::T, tree2::T)::Int64 where T <: GeneralNode
    bt3 = get_bipartitions(tree1)
    bt4 = get_bipartitions(tree2)
    length(bt3)+length(bt4) - 2* length(intersect(bt3, bt4))
end

"""
    get_bipartitions(tree::T)::Vector{Tuple} where T <:GeneralNode

Get a vector of all bipartions of `tree`.

Returns a vector containing Tuples of sets representing the bipartions.
"""
function get_bipartitions(tree::T)::Vector{Tuple} where T <:GeneralNode
    po_vect= post_order(tree)[1:end-1]
    bt = Vector{Tuple}(undef, length(po_vect))
    Base.Threads.@threads for ind in eachindex(po_vect)
        elem = po_vect[ind]::T
        inset = Set{Int64}()
        outset = Set{Int64}()
        @simd for elem2 in po_vect
            startswith(elem2.binary, elem.binary) ? push!(inset, elem2.num) : push!(outset, elem2.num)
        end # for
        @inbounds bt[ind] = (inset, outset)
    end # for
    bt
end

"""
    BHV_bounds(tree1::T, tree2::T)::Tuple{Float64, Float64} where T <:GeneralNode

This function calculates the lower and upper bounds of the geodesic in the
Billera-Holmes-Vogtman space.

Returns tuple of floats.
"""
function BHV_bounds(tree1::T, tree2::T)::Tuple{Float64, Float64} where T <:GeneralNode
    res_upper_1 = Threads.Atomic{Float64}(0.0)
    res_upper_2 = Threads.Atomic{Float64}(0.0)
    res_upper_3 = Threads.Atomic{Float64}(0.0)
    po::Vector{T} = post_order(tree1)
    Base.Threads.@threads for node in po[1:end-1]
        nom = find_num(tree2, node.num)
        if node.mother.num == nom.mother.num
            Threads.atomic_add!(res_upper_3, (node.inc_length-nom.inc_length)^2)
        else
            Threads.atomic_add!(res_upper_2 , nom.inc_length^2)
            Threads.atomic_add!(res_upper_1, node.inc_length^2)
        end # if

    end # for
    res_low::Float64 = res_upper_1[]+res_upper_2[]+res_upper_3[]
    res_high::Float64 = (sqrt(res_upper_2[])+sqrt(res_upper_1[]))^2+res_upper_3[]

    sqrt(res_low), sqrt(res_high)
end


text:
module TestMisc

include("preamble.jl")

# The only name crash of Python's double-underscore methods and
# Julia's methods is module-level `__init__`.  Let's make sure to
# never define it since it would be confusing if someone finds that
# `PyBase.__init__` exists and it does not have any effect.  I should
# define it in a dummy submodule if I need it.
@test !isdefined(PyBase, :__init__)

end  # module


text:
using LOLTools.LeagueV4

api_key = get(ENV, "RIOT_TOKEN", "")
platform = "kr"
league = LeagueV4.challengerleagues(api_key, platform, "RANKED_SOLO_5x5")
@info sort(league.entries, by = x -> x.leaguePoints)[1:2]


text:
module MetidaReports

    using StatsModels, DataFrames
    import Base: ht_keyindex

    include("bioeq.jl")

end # module


text:
using JuMP
using GLPKMathProgInterface
using Ipopt

n = 3;
lowRes = -20;
hiRes  = 20;
m = Model(solver=IpoptSolver())
@defVar(m, lowRes <= α[i=1:n]  <=hiRes)
@defVar(m, lowRes <= β1[i=1:n] <=hiRes)
@defVar(m, lowRes <= β2[i=1:n] <=hiRes)
@defVar(m, lowRes <= β3[i=1:n] <=hiRes)
@defVar(m, lowRes <= β4[i=1:n] <=hiRes)

@defVar(m, lowRes <= A011 <= hiRes)
@defVar(m, lowRes <= A012 <= hiRes)
@defVar(m, lowRes <= A013 <= hiRes)
@defVar(m, lowRes <= A021 <= hiRes)
@defVar(m, lowRes <= A022 <= hiRes)
@defVar(m, lowRes <= A023 <= hiRes)
@defVar(m, lowRes <= A031 <= hiRes)
@defVar(m, lowRes <= A032 <= hiRes)
@defVar(m, lowRes <= A033 <= hiRes)

@defVar(m, lowRes <= A111 <= hiRes)
@defVar(m, lowRes <= A112 <= hiRes)
@defVar(m, lowRes <= A113 <= hiRes)
@defVar(m, lowRes <= A121 <= hiRes)
@defVar(m, lowRes <= A122 <= hiRes)
@defVar(m, lowRes <= A123 <= hiRes)
@defVar(m, lowRes <= A131 <= hiRes)
@defVar(m, lowRes <= A132 <= hiRes)
@defVar(m, lowRes <= A133 <= hiRes)

@defVar(m, lowRes <= B011 <= hiRes)
@defVar(m, lowRes <= B012 <= hiRes)
@defVar(m, lowRes <= B013 <= hiRes)
@defVar(m, lowRes <= B021 <= hiRes)
@defVar(m, lowRes <= B022 <= hiRes)
@defVar(m, lowRes <= B023 <= hiRes)
@defVar(m, lowRes <= B031 <= hiRes)
@defVar(m, lowRes <= B032 <= hiRes)
@defVar(m, lowRes <= B033 <= hiRes)

@defVar(m, lowRes <= B111 <= hiRes)
@defVar(m, lowRes <= B112 <= hiRes)
@defVar(m, lowRes <= B113 <= hiRes)
@defVar(m, lowRes <= B121 <= hiRes)
@defVar(m, lowRes <= B122 <= hiRes)
@defVar(m, lowRes <= B123 <= hiRes)
@defVar(m, lowRes <= B131 <= hiRes)
@defVar(m, lowRes <= B132 <= hiRes)
@defVar(m, lowRes <= B133 <= hiRes)

A0 = [A011 A012
     A021 A022]
A1 = [A111 A112
     A121 A122]
B0 = [B011 B012
      B021 B022]
B1 = [B111 B112
      B121 B122]
e = [1 1]'

#@setObjective(m, Max, 0)
@addConstraints(m,begin
                  sum{α[i], i=1:n}==1
                  sum{β1[i], i=1:n}==1
                  sum{β1[i], i=1:n}==1
                  sum{β2[i], i=1:n}==0
                  sum{β3[i], i=1:n}==0
                  sum{β4[i], i=1:n}==0
                  sum{β4[i], i=1:n}==0
                end)
@addNLConstraints(m,begin
                  B111*β1[1]+B112*β1[1]+B113*β1[1] + B121*β1[2]+B122*β1[2]+B123*β1[2] + B131*β1[3]+B132*β1[3]+B133*β1[3] ==0
                  B111*β2[1]+B112*β2[1]+B113*β2[1] + B121*β2[2]+B122*β2[2]+B123*β2[2] + B131*β2[3]+B132*β2[3]+B133*β2[3] ==1
                  B111*β3[1]+B112*β3[1]+B113*β3[1] + B121*β3[2]+B122*β3[2]+B123*β3[2] + B131*β3[3]+B132*β3[3]+B133*β3[3] ==0
                  B111*β4[1]+B112*β4[1]+B113*β4[1] + B121*β4[2]+B122*β4[2]+B123*β4[2] + B131*β4[3]+B132*β4[3]+B133*β4[3] ==0
                  A011*α[1]+A012*α[1]+A013*α[1]    + A021*α[2]+A022*α[2]+A023*α[2]    + A031*α[3]+A032*α[3]+A033*α[3] ==1//2
                  B011*α[1]+B012*α[1]+B013*α[1]    + B021*α[2]+B022*α[2]+B023*α[2]    + B031*α[3]+B032*α[3]+B033*α[3] ==1
                  A111*β1[1]+A112*β1[1]+A113*β1[1] + A121*β1[2]+A122*β1[2]+A123*β1[2] + A131*β1[3]+A132*β1[3]+A133*β1[3]==1
                  A111*β2[1]+A112*β2[1]+A113*β2[1] + A121*β2[2]+A122*β2[2]+A123*β2[2] + A131*β2[3]+A132*β2[3]+A133*β2[3]==0
                  A111*β3[1]+A112*β3[1]+A113*β3[1] + A121*β3[2]+A122*β3[2]+A123*β3[2] + A131*β3[3]+A132*β3[3]+A133*β3[3]==-1
                  A111*β4[1]+A112*β4[1]+A113*β4[1] + A121*β4[2]+A122*β4[2]+A123*β4[2] + A131*β4[3]+A132*β4[3]+A133*β4[3]==0
                end)

@addNLConstraint(m,α[1]  * (B011+B012+B013)*(B011+B012+B013)  + α[2]  * (B021+B022+B123)*(B021+B022+B123) + α[3]  * (B031+B032+B133)*(B031+B032+B133)  == 3//2)
@addNLConstraint(m,β1[1] * (B111+B112+B013)*(B111+B112+B013)  + β1[2] * (B121+B122+B123)*(B121+B122+B123) + β1[3] * (B131+B132+B133)*(B131+B132+B133)  == 1)
@addNLConstraint(m,β2[1] * (B111+B112+B013)*(B111+B112+B013)  + β2[2] * (B121+B122+B123)*(B121+B122+B123) + β2[3] * (B131+B132+B133)*(B131+B132+B133)  == 0)
@addNLConstraint(m,β3[1] * (B111+B112+B013)*(B111+B112+B013)  + β3[2] * (B121+B122+B123)*(B121+B122+B123) + β3[3] * (B131+B132+B133)*(B131+B132+B133)  == -1)
@addNLConstraint(m,β4[1] * (B111+B112+B013)*(B111+B112+B013)  + β4[2] * (B121+B122+B123)*(B121+B122+B123) + β4[3] * (B131+B132+B133)*(B131+B132+B133)  == 2)

@addNLConstraint(m,β1[1] * (B111*(B111+B112+B113) + B112*(B121+B122+B123) + B113*(B131+B132+B133))  + β1[2] * (B121*(B111+B112+B113) + B122*(B121+B122+B123) + B123*(B131+B132+B133))  + β1[3] * (B131*(B111+B112+B113) + B132*(B121+B122+B123) + B133*(B131+B132+B133))  == 0)
@addNLConstraint(m,β2[1] * (B111*(B111+B112+B113) + B112*(B121+B122+B123) + B113*(B131+B132+B133))  + β2[2] * (B121*(B111+B112+B113) + B122*(B121+B122+B123) + B123*(B131+B132+B133))  + β2[3] * (B131*(B111+B112+B113) + B132*(B121+B122+B123) + B133*(B131+B132+B133))  == 0)
@addNLConstraint(m,β3[1] * (B111*(B111+B112+B113) + B112*(B121+B122+B123) + B113*(B131+B132+B133))  + β3[2] * (B121*(B111+B112+B113) + B122*(B121+B122+B123) + B123*(B131+B132+B133))  + β3[3] * (B131*(B111+B112+B113) + B132*(B121+B122+B123) + B133*(B131+B132+B133))  == 0)
@addNLConstraint(m,β4[1] * (B111*(B111+B112+B113) + B112*(B121+B122+B123) + B113*(B131+B132+B133))  + β4[2] * (B121*(B111+B112+B113) + B122*(B121+B122+B123) + B123*(B131+B132+B133))  + β4[3] * (B131*(B111+B112+B113) + B132*(B121+B122+B123) + B133*(B131+B132+B133))  == 1)

@addNLConstraint(m,(1/2)*((A111*(B011+B012+B013) + A112*(B021+B022+B023) + A113*(B031+B032+B033))*β1[1] +
      (A121*(B011+B012+B013) + A122*(B021+B022+B023) + A123*(B031+B032+B033))β1[2] +
      (A131*(B011+B012+B013) + A132*(B021+B022+B023) + A133*(B031+B032+B033))β1[3]) +
      (1/3)*((A111*(B011+B012+B013) + A112*(B021+B022+B023) + A113*(B031+B032+B033))*β3[1] +
            (A121*(B011+B012+B013) + A122*(B021+B022+B023) + A123*(B031+B032+B033))β3[2] +
            (A131*(B011+B012+B013) + A132*(B021+B022+B023) + A133*(B031+B032+B033))β3[3]) == 0)

print(m)

status = solve(m)

println("Objective value: ", getObjectiveValue(m))
println("β1 = ", getValue(β1))
println("β2 = ", getValue(β2))
println("β3 = ", getValue(β3))
println("β4 = ", getValue(β4))


text:
require("ZMQ")
using ZMQ

@assert length(ZMQ.version()) == 3

major, minor, patch = ZMQ.version()

ctx=ZMQContext(1)

@assert typeof(ctx) == ZMQContext

try
	ZMQContext(-1)
	@assert false
catch ex
	@assert typeof(ex) == ZMQStateError
end

ZMQ.close(ctx)

#try to create socket with expired context
try
	ZMQSocket(ctx, ZMQ_PUB)
	@assert false
catch ex
	@assert typeof(ex) == ZMQStateError
end

ctx2=ZMQContext(1)
s=ZMQSocket(ctx2, ZMQ_PUB)
@assert typeof(s) == ZMQSocket
ZMQ.close(s)

#trying to close already closed socket
try
	ZMQ.close(s)
catch ex
	@assert typeof(ex) == ZMQStateError
end

s1=ZMQSocket(ctx2, ZMQ_REP)
if major == 2
	ZMQ.set_hwm(s1, 1000)
else
	ZMQ.set_sndhwm(s1, 1000)
end
ZMQ.set_linger(s1, 1)
ZMQ.set_identity(s1, "abcd")

@assert ZMQ.get_identity(s1)::String == "abcd"
if major == 2
	@assert ZMQ.get_hwm(s1)::Integer == 1000
else
	@assert ZMQ.get_sndhwm(s1)::Integer == 1000
end
@assert ZMQ.get_linger(s1)::Integer == 1
@assert ZMQ.ismore(s1) == false

s2=ZMQSocket(ctx2, ZMQ_REQ)
@assert ZMQ.get_type(s1) == ZMQ_REP
@assert ZMQ.get_type(s2) == ZMQ_REQ

ZMQ.bind(s1, "tcp://*:5555")
ZMQ.connect(s2, "tcp://localhost:5555")

ZMQ.send(s2, ZMQMessage("test request"))
@assert (bytestring(ZMQ.recv(s1)) == "test request")
ZMQ.send(s1, ZMQMessage("test response"))
@assert (bytestring(ZMQ.recv(s2)) == "test response")

ZMQ.send(s2, ZMQMessage("another test request"))
msg = ZMQ.recv(s1)
o=convert(IOStream, msg)
seek(o, 0)
@assert (takebuf_string(o)=="another test request")

ZMQ.close(s1)
ZMQ.close(s2)
ZMQ.close(ctx2)


text:
#!/usr/bin/env julia

template(text) =
    """
     $(border(text, '_'))
    < $text >
     $(border(text, '-'))
            \\   ^__^
             \\  (oo)\\_______
                (__)\\       )\\/\\
                    ||----w |
                    ||     ||
    """

border(text, char) = repeat(char, length(text) + 2)

if length(ARGS) == 1
    text = ARGS[1]
    println(template(text))
else
    println("Usage: cowsay [message]")
end

max_stars_repo_name: gaborbata/rosetta-cow

max_issues_repo_name: gaborbata/rosetta-cow

max_forks_repo_name: gaborbata/rosetta-cow


text:
using LearnBase: getobs, nobs, default_obsdim

@test typeof(LearnBase.getobs) <: Function
@test typeof(LearnBase.getobs!) <: Function
@test typeof(LearnBase.gettargets) <: Function
@test typeof(LearnBase.datasubset) <: Function

@testset "getobs and nobs" begin

    @testset "array" begin
        a = rand(2,3)
        @test nobs(a) == 3
        @test @inferred getobs(a, 1) == a[:,1]
        @test @inferred getobs(a, 2) == a[:,2]
        @test @inferred getobs(a, 1:2) == a[:,1:2]
        @test @inferred getobs(a, 1, obsdim=1) == a[1,:]
        @test @inferred getobs(a, 2, obsdim=1) == a[2,:]
        @test @inferred getobs(a, 2, obsdim=nothing) ≈ a[:,2]
    end

    @testset "tuple" begin
        # A dataset with 3 observations, each with 2 input features
        X, Y = rand(2, 3), rand(3)
        dataset = (X, Y)
        @test nobs(dataset) == 3
        if VERSION >= v"1.6"
            o = @inferred getobs(dataset, 2)
        else
            o = getobs(dataset, 2)
        end
        @test o[1] == X[:,2]
        @test o[2] == Y[2]

        if VERSION >= v"1.6"
            o = @inferred getobs(dataset, 1:2)
        else
            o = getobs(dataset, 1:2)
        end

        @test o[1] == X[:,1:2]
        @test o[2] == Y[1:2]
    end

    @testset "named tuple" begin
        X, Y = rand(2, 3), rand(3)
        dataset = (x=X, y=Y)
        @test nobs(dataset) == 3
        if VERSION >= v"1.6"
            o = @inferred getobs(dataset, 2)
        else
            o = getobs(dataset, 2)
        end
        @test o.x == X[:,2]
        @test o.y == Y[2]

        if VERSION >= v"1.6"
            o = @inferred getobs(dataset, 1:2)
        else
            o = getobs(dataset, 1:2)
        end
        @test o.x == X[:,1:2]
        @test o.y == Y[1:2]
    end

    @testset "dict" begin
        X, Y = rand(2, 3), rand(3)
        dataset = Dict("X" => X, "Y" => Y)
        @test nobs(dataset) == 3

        # o = @inferred getobs(dataset, 2) # not inferred
        o = getobs(dataset, 2)
        @test o["X"] == X[:,2]
        @test o["Y"] == Y[2]

        o = getobs(dataset, 1:2)
        @test o["X"] == X[:,1:2]
        @test o["Y"] == Y[1:2]
    end
end


text:
function mask_vals(img,mask=ones(Float16, (size(img)[1],size(img)[2])))
    vals = []
    for i in 1:size(img)[1]
        for j in 1:size(img)[2]
            if mask[i,j] == 1
                push!(vals,img[i,j])
            end
        end
    end
    return vals
end

"""
A mask can be applied to a datacube cube so that each image of the datacube has 0s wherever the mask is dark and the original value wherever the mask is lit.
# Example:
```
datacube = rand(1:10.0, 10,10,10)
mask = rand(0:1, 10, 10)
apply_mask(datacube, mask)
```
"""
function apply_mask(data, mask = ones((size(data)[1], size(data)[2])))
    #use multiple dispatch here
    if ndims(data) == 2
        return data .* mask
    end
    if typeof(data) == VectorOfArray{Any,3,Array{Any,1}}
        for i in 1:length(data)
            data[i] = data[i] .* mask
        end
        return data
    end
    masked_data = copy(data)
    for i in 1:size(data)[3]
        masked_data[:, :, i] = data[:, :, i] .* mask
    end
    return masked_data
end

add_masks = function(mask1, mask2)
    mask_sum =  BitArray(zeros(size(mask1)[1], size(mask2[2])))
    mask_sum           = mask1 .| BitArray(mask2)
    return sparse(Array{Int8}(mask_sum))
end


text:
"""
defines power from generators at each bus, varies in contingencies
"""
function expression_bus_generation(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    if !haskey(var(pm, nw), :bus_pg)
        var(pm, nw)[:bus_pg] = Dict{Int,Any}()
    end
    if !haskey(var(pm, nw), :bus_qg)
        var(pm, nw)[:bus_qg] = Dict{Int,Any}()
    end

    bus = ref(pm, nw, :bus, i)
    bus_gens = ref(pm, nw, :bus_gens, i)

    expression_bus_generation(pm, nw, i, bus_gens)
end

"""
defines power from non-generator components at each bus, static in contingencies
"""
function expression_bus_withdrawal(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    if !haskey(var(pm, nw), :bus_wdp)
        var(pm, nw)[:bus_wdp] = Dict{Int,Any}()
    end
    if !haskey(var(pm, nw), :bus_wdq)
        var(pm, nw)[:bus_wdq] = Dict{Int,Any}()
    end

    bus = ref(pm, nw, :bus, i)
    bus_loads = ref(pm, nw, :bus_loads, i)
    bus_shunts = ref(pm, nw, :bus_shunts, i)
    bus_storage = ref(pm, nw, :bus_storage, i)

    bus_pd = Dict(k => ref(pm, nw, :load, k, "pd") for k in bus_loads)
    bus_qd = Dict(k => ref(pm, nw, :load, k, "qd") for k in bus_loads)

    bus_gs = Dict(k => ref(pm, nw, :shunt, k, "gs") for k in bus_shunts)
    bus_bs = Dict(k => ref(pm, nw, :shunt, k, "bs") for k in bus_shunts)

    expression_bus_withdrawal(pm, nw, i, bus_storage, bus_pd, bus_qd, bus_gs, bus_bs)
end

""
function expression_branch_power_yt_from_goc(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    if !haskey(var(pm, nw), :p)
        var(pm, nw)[:p] = Dict{Tuple{Int,Int,Int},Any}()
    end
    if !haskey(var(pm, nw), :q)
        var(pm, nw)[:q] = Dict{Tuple{Int,Int,Int},Any}()
    end

    branch = ref(pm, nw, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    g, b = calc_branch_y(branch)
    tr, ti = calc_branch_t(branch)
    g_fr = branch["g_fr"]
    b_fr = branch["b_fr"]
    tm = branch["tap"]

    if branch["transformer"]
        expression_branch_power_yt_from_goc(pm, nw, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    else
        expression_branch_power_yt_from(pm, nw, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    end
end

""
function expression_branch_power_yt_to(pm::AbstractPowerModel, i::Int; nw::Int=pm.cnw)
    if !haskey(var(pm, nw), :p)
        var(pm, nw)[:p] = Dict{Tuple{Int,Int,Int},Any}()
    end
    if !haskey(var(pm, nw), :q)
        var(pm, nw)[:q] = Dict{Tuple{Int,Int,Int},Any}()
    end

    branch = ref(pm, nw, :branch, i)
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    g, b = calc_branch_y(branch)
    tr, ti = calc_branch_t(branch)
    g_to = branch["g_to"]
    b_to = branch["b_to"]
    tm = branch["tap"]

    expression_branch_power_yt_to(pm, nw, f_bus, t_bus, f_idx, t_idx, g, b, g_to, b_to, tr, ti, tm)
end