text:
export BrickletIndustrialDigitalIn4EdgeCountConfig
struct BrickletIndustrialDigitalIn4EdgeCountConfig
    edge_type::Integer
    debounce::Integer
end

export BrickletIndustrialDigitalIn4Identity
struct BrickletIndustrialDigitalIn4Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletIndustrialDigitalIn4
"""
4 galvanically isolated digital inputs
"""
mutable struct BrickletIndustrialDigitalIn4 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletIndustrialDigitalIn4(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_industrial_digital_in_4")
        deviceInternal = package.BrickletIndustrialDigitalIn4(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_value
"""
    $(SIGNATURES)

Returns the input value with a bitmask. The bitmask is 16bit long, *true*
refers to high and *false* refers to low.

For example: The value 3 or 0b0011 means that pins 0-1 are high and the other
pins are low.

If no groups are used (see :func:`Set Group`), the pins correspond to the
markings on the IndustrialDigital In 4 Bricklet.

If groups are used, the pins correspond to the element in the group.
Element 1 in the group will get pins 0-3, element 2 pins 4-7, element 3
pins 8-11 and element 4 pins 12-15.
"""
function get_value(device::BrickletIndustrialDigitalIn4)
    return device.deviceInternal.get_value()
end

export set_group
"""
    $(SIGNATURES)

Sets a group of Digital In 4 Bricklets that should work together. You can
find Bricklets that can be grouped together with :func:`Get Available For Group`.

The group consists of 4 elements. Element 1 in the group will get pins 0-3,
element 2 pins 4-7, element 3 pins 8-11 and element 4 pins 12-15.

Each element can either be one of the ports ('a' to 'd') or 'n' if it should
not be used.

For example: If you have two Digital In 4 Bricklets connected to port A and
port B respectively, you could call with ``['a', 'b', 'n', 'n']``.

Now the pins on the Digital In 4 on port A are assigned to 0-3 and the
pins on the Digital In 4 on port B are assigned to 4-7. It is now possible
to call :func:`Get Value` and read out two Bricklets at the same time.

Changing the group configuration resets all edge counter configurations
and values.
"""
function set_group(device::BrickletIndustrialDigitalIn4, group)
    device.deviceInternal.set_group(group)
end

export get_group
"""
    $(SIGNATURES)

Returns the group as set by :func:`Set Group`
"""
function get_group(device::BrickletIndustrialDigitalIn4)
    return device.deviceInternal.get_group()
end

export get_available_for_group
"""
    $(SIGNATURES)

Returns a bitmask of ports that are available for grouping. For example the
value 5 or 0b0101 means: Port A and port C are connected to Bricklets that
can be grouped together.
"""
function get_available_for_group(device::BrickletIndustrialDigitalIn4)
    return device.deviceInternal.get_available_for_group()
end

export set_debounce_period
"""
    $(SIGNATURES)

Sets the debounce period of the :cb:`Interrupt` callback.

For example: If you set this value to 100, you will get the interrupt
maximal every 100ms. This is necessary if something that bounces is
connected to the Digital In 4 Bricklet, such as a button.
"""
function set_debounce_period(device::BrickletIndustrialDigitalIn4, debounce)
    device.deviceInternal.set_debounce_period(debounce)
end

export get_debounce_period
"""
    $(SIGNATURES)

Returns the debounce period as set by :func:`Set Debounce Period`.
"""
function get_debounce_period(device::BrickletIndustrialDigitalIn4)
    return device.deviceInternal.get_debounce_period()
end

export set_interrupt
"""
    $(SIGNATURES)

Sets the pins on which an interrupt is activated with a bitmask.
Interrupts are triggered on changes of the voltage level of the pin,
i.e. changes from high to low and low to high.

For example: An interrupt bitmask of 9 or 0b1001 will enable the interrupt for
pins 0 and 3.

The interrupts use the grouping as set by :func:`Set Group`.

The interrupt is delivered with the :cb:`Interrupt` callback.
"""
function set_interrupt(device::BrickletIndustrialDigitalIn4, interrupt_mask)
    device.deviceInternal.set_interrupt(interrupt_mask)
end

export get_interrupt
"""
    $(SIGNATURES)

Returns the interrupt bitmask as set by :func:`Set Interrupt`.
"""
function get_interrupt(device::BrickletIndustrialDigitalIn4)
    return device.deviceInternal.get_interrupt()
end

export get_edge_count
"""
    $(SIGNATURES)

Returns the current value of the edge counter for the selected pin. You can
configure the edges that are counted with :func:`Set Edge Count Config`.

If you set the reset counter to *true*, the count is set back to 0
directly after it is read.

The edge counters use the grouping as set by :func:`Set Group`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function get_edge_count(device::BrickletIndustrialDigitalIn4, pin, reset_counter)
    return device.deviceInternal.get_edge_count(pin, reset_counter)
end

export set_edge_count_config
"""
    $(SIGNATURES)

Configures the edge counter for the selected pins. A bitmask of 9 or 0b1001 will
enable the edge counter for pins 0 and 3.

The edge type parameter configures if rising edges, falling edges or
both are counted if the pin is configured for input. Possible edge types are:

* 0 = rising
* 1 = falling
* 2 = both

Configuring an edge counter resets its value to 0.

If you don't know what any of this means, just leave it at default. The
default configuration is very likely OK for you.

The edge counters use the grouping as set by :func:`Set Group`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function set_edge_count_config(device::BrickletIndustrialDigitalIn4, selection_mask, edge_type, debounce)
    device.deviceInternal.set_edge_count_config(selection_mask, edge_type, debounce)
end

export get_edge_count_config
"""
    $(SIGNATURES)

Returns the edge type and debounce time for the selected pin as set by
:func:`Set Edge Count Config`.

.. versionadded:: 2.0.1\$nbsp;(Plugin)
"""
function get_edge_count_config(device::BrickletIndustrialDigitalIn4, pin)
    return device.deviceInternal.get_edge_count_config(pin)
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletIndustrialDigitalIn4)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletIndustrialDigitalIn4, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
"""
    f∂f(f,x)

returns f and ∂f/∂x evaluated in `x`, using `ForwardDiff.jl`, `DiffResults.jl` and `StaticArrays.jl` to calculate everything in one pass.
"""
function f∂f(f,x::T) where T
    _f(z) = f(only(z))
    x_vec =   SVector(x)
    ∂result = DiffResults.GradientResult(x_vec)
    _∂f =  ForwardDiff.gradient!(∂result, _f,x_vec)
    fx =  DiffResults.value(_∂f)
    ∂f∂x = only(DiffResults.gradient(_∂f))
    return fx,∂f∂x
end

"""
    f∂f∂2f(f,x)

returns f,∂f/∂x,and ∂²f/∂²x and evaluated in `x`, using `ForwardDiff.jl`, `DiffResults.jl` and `StaticArrays.jl` to calculate everything in one pass.
"""
function f∂f∂2f(f,x::T) where T
    _f(z) = f(only(z))
    x_vec =   SVector(x)
    ∂result = DiffResults.HessianResult(x_vec)
    _∂f =  ForwardDiff.hessian!(∂result, _f,x_vec)
    fx =  DiffResults.value(_∂f)
    ∂f∂x = only(DiffResults.gradient(_∂f))
    ∂²f∂²x =  only(DiffResults.hessian(_∂f))
    return fx,∂f∂x,∂²f∂²x
end

"""
    ad_newton(f,x0;kwargs...)

    Performs newton root finding (via Roots.jl), but the derivatives are calculated with ForwardDiff.jl
    kwargs are passed to `Roots.find_zero`
"""
function ad_newton(f,x0;kwargs...)
    function f0(x)
        fx,∂f∂x = f∂f(f,x)
        return fx,fx/∂f∂x
    end
    return Roots.find_zero(f0,x0,Roots.Newton();kwargs...)
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ 9820cc3a-fe50-11ea-3fb6-991358edb7ff
using Pkg, DrWatson

# ╔═╡ 98210cc2-fe50-11ea-0fec-cfe81a89d0cb
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample, GLM
	using StatisticalRethinking
end

# ╔═╡ c5f141c2-fe4f-11ea-03d2-e5eb5d2349f6
md"## Clip-06-02-06s.jl"

# ╔═╡ 9831d8b8-fe50-11ea-38fd-c99a2b5fd0bc
md"### Snippet 6.1"

# ╔═╡ 98328740-fe50-11ea-3008-27d3dda98cd2
begin
	N = 100
	df = DataFrame(
		height = rand(Normal(10, 2), N),
		leg_prop = rand(Uniform(0.4, 0.5), N),
	)
	df.leg_left = df.leg_prop .* df.height + rand(Normal(0, 0.02), N)
	df.leg_right = df.leg_prop .* df.height + rand(Normal(0, 0.02), N)
end;

# ╔═╡ 9842396a-fe50-11ea-217d-2bba5fa44fb5
md"### Snippet 6.2"

# ╔═╡ 9842dbb8-fe50-11ea-158e-cbc2873c64bd
stan6_1 = "
data {
  int <lower=1> N;
  vector[N] H;
  vector[N] LL;
  vector[N] LR;
}
parameters {
  real a;
  real bL;
  real bR;
  real <lower=0> sigma;
}
model {
  vector[N] mu;
  mu = a + bL * LL + bR * LR;
  a ~ normal(10, 100);
  bL ~ normal(2, 10);
  bR ~ normal(2, 10);
  sigma ~ exponential(1);
  H ~ normal(mu, sigma);
}
";

# ╔═╡ 9851b4a8-fe50-11ea-083b-17e3c182a55a
begin
	m6_1s = SampleModel("m6.1s", stan6_1,
		method=StanSample.Sample(num_samples=1000))
	m_6_1_data = Dict(
	  :H => df[:, :height],
	  :LL => df[:, :leg_left],
	  :LR => df[:, :leg_right],
	  :N => size(df, 1)
	)
	rc6_1s = stan_sample(m6_1s, data=m_6_1_data)
	success(rc6_1s) && (part6_1s = read_samples(m6_1s, :particles))
end

# ╔═╡ eb921ad4-fe50-11ea-2450-7dfe19847755
if success(rc6_1s)
	(s0, p0) = plot_model_coef([m6_1s], [:a, :bL, :bR, :sigma];
		title="Multicollinearity between bL and bR")
	p0
end

# ╔═╡ b825edf6-00eb-11eb-039c-cfcb859dc43d
s0

# ╔═╡ 985265b0-fe50-11ea-0594-05ac311d2e87
if success(rc6_1s)
	post6_1s_df = read_samples(m6_1s, :dataframe)

	# Fit a linear regression

	m = lm(@formula(bL ~ bR), post6_1s_df)

	# estimated coefficients from the model

	coefs = coef(m)

	fig1 = plot(xlabel="bR", ylabel="bL", lab="bL ~ bR")
	plot!(post6_1s_df[:, :bR], post6_1s_df[:, :bL])
	fig2 = density(part6_1s.bR.particles + part6_1s.bL.particles, xlabel="sum of bL and bR",
		ylabel="Density", lab="bL + bR")
	plot(fig1, fig2, layout=(1, 2))
end

# ╔═╡ 9862bdb4-fe50-11ea-1402-bbac3257c25d
md"## End of clip-06-02-06s.jl"

# ╔═╡ Cell order:
# ╟─c5f141c2-fe4f-11ea-03d2-e5eb5d2349f6
# ╠═9820cc3a-fe50-11ea-3fb6-991358edb7ff
# ╠═98210cc2-fe50-11ea-0fec-cfe81a89d0cb
# ╟─9831d8b8-fe50-11ea-38fd-c99a2b5fd0bc
# ╠═98328740-fe50-11ea-3008-27d3dda98cd2
# ╟─9842396a-fe50-11ea-217d-2bba5fa44fb5
# ╠═9842dbb8-fe50-11ea-158e-cbc2873c64bd
# ╠═9851b4a8-fe50-11ea-083b-17e3c182a55a
# ╠═eb921ad4-fe50-11ea-2450-7dfe19847755
# ╠═b825edf6-00eb-11eb-039c-cfcb859dc43d
# ╠═985265b0-fe50-11ea-0594-05ac311d2e87
# ╟─9862bdb4-fe50-11ea-1402-bbac3257c25d


text:
const DO_NOT_TRY_URL_LIST = String[
    ]


text:
function compute_one_form(t, x, p, Θ)
    @assert size(x) == size(p)
    @assert length(t) == size(x,2) == size(p,2)

    tp = zeros(eltype(x), size(x,1))

    for k in 1:size(x,3)
        for j in 1:size(x,2)
            Θ(t[j], x[:,j,k], tp)
            p[:,j,k] .= tp
        end
    end
end

function compute_velocity(γ, γ̇)
    @assert size(γ) == size(γ̇)

    n = size(γ,2)
    k = collect(0:div(n-1,2)+1)
    κ = 2π .* 1im .* k ./ n

    for l in 1:size(γ,1)
        γ̂ = rfft(γ[l,:])
        γ̇[l,:] .= irfft(γ̂ .* κ, n)
    end
end

function compute_correction(t, x, λ, σ, g)
    @assert size(x) == size(λ) == size(σ)
    @assert length(t) == size(x,2) == size(σ,2)

    tσ = zeros(eltype(x), size(x,1))

    for k in 1:size(x,3)
        for j in 1:size(x,2)
            g(t[j], x[:,j,k], λ[:,j,k], tσ)
            σ[:,j,k] .= tσ
        end
    end
end

function compute_loop_integral(p::Array{T,2}, v::Array{T,2}) where {T}
    local result = zero(T)
    local error  = zero(T)

    for i in 1:size(v,2)
        for l in 1:size(v,1)
            result, error = compensated_summation(p[l,i] * v[l,i], result, error)
        end
    end

    return result
end


text:
using Plots,Mux,Interact,CSV,DataFrames

function addplot(varsym,varnames)
    options = Observable{Any}(Dict{String,Symbol}()) # create an Observable var 'options' as Dict{String,Symbol}
    map((varname,varsym) -> options.val[varname] = varsym, varnames, varsym) # map the 'varnames' and 'vaesym' to 'options'
    plotvar = dropdown(options.val;label = nothing,multiple = true) # create a widget of type 'dropdown' to define the var's to plot
    # plt = Interact.@map plot(collect(-50:50),&mOutput ,label = "Model output")
end

function updateVarToPlot(datain)
    # println("updateVarToPlot function input type is : $(typeof(datain))")
    # println("updateVarToPlot function input is : $datain")
    if isa(datain,DataFrame)
        varsym = names(datain) # parse the 'DataFrame' fields names from 'data.val' as Symbol
        varnames = string.(names(datain)) # parse the 'DataFrame' fields names from 'data.val' as String
        addplot(varsym,varnames)
        # println(varsym,varnames)
    else
        varsym = nothing
        varnames = nothing
        # println("nothing")
    end
end

function createLayout()
    data = Observable{Any}(DataFrame) # create an Observible var 'data' of type 'DataFrame'
    plt = Observable{Any}(plot()) # create an Observible var 'data' of type 'DataFrame'

    loadbutton = filepicker(label = "load CSV file"; multiple = false, accept = "*.csv") #call the load file widget
    map!(CSV.read, data, loadbutton) # map the path from the 'loadbutton' to the CSV.read and save the output to 'data'
    map!(updateVarToPlot,,plt,data)
    # Interact.@map updateVarToPlot(&data)
    # Observables.onany(updateVarToPlot,data)
    # (varsym,varnames) = Observables.onany(updateVarToPlot,data)

    # addPlotButton = button(label = "Add a Plot")
    # on(_->addplot(data),addPlotButton)

    # mOutput = Interact.@map myModel(&p1s,&p2s)
    # plt = Interact.@map plot(collect(-50:50),&mOutput ,label = "Model output")
    # wdg = Widget(["p1" => p1s, "p2" => p2s, "addPlotButton" => addPlotButton, "loadFileButton" => loadFileButton], output = mOutput)
    # @layout! wdg hbox(plt,vbox(:p1, :p2, :loadFileButton, :addPlotButton))
end

function serverLauout(destinationPort)
    try
        WebIO.webio_serve(page("/",req -> createLayout()),destinationPort)
    catch e
        if isa(e,IOError)
            # sleep and then try later
            sleep(0.1)
            serverLauout(destinationPort)
        else
            throw(e)
        end
    end
end

serverLauout(8000)

max_stars_repo_name: kobied/interactTests

max_issues_repo_name: kobied/interactTests

max_forks_repo_name: kobied/interactTests


text:
"""
    $SIGNATURES

Return `(::Bool, ::OKVarlist, ::NotOkayVarList)` on whether all other variables (besides `loovar::Symbol`)
attached to factor `fct::Symbol` are all initialized -- i.e. `fct` is usable.

Notes:
- Special carve out for multihypo cases, see issue 427, where at least one hypothesis should be available, but not all required at first.

Development Notes
* TODO get faster version of isInitialized for database version

Related

doautoinit!, initManual!, isInitialized, isMultihypo
"""
function factorCanInitFromOtherVars(dfg::AbstractDFG,
                                    fct::Symbol,
                                    loovar::Symbol;
                                    solveKey::Symbol=:default)
  #
  # all variables attached to this factor
  varsyms = DFG.getNeighbors(dfg, fct)

  # which element is being solved for
  sfidx = (1:length(varsyms))[varsyms .== loovar][1]
  # list of factors to use in init operation
  fctlist = Symbol[]
  # list fo variables that cannot be used
  faillist = Symbol[]
  isinit = Bool[]
  for vsym in varsyms
    # check each variable one by one
    xi = DFG.getVariable(dfg, vsym)
    isi = isInitialized(xi, solveKey)
    push!(isinit, isi)
    if !isi
      push!(faillist, vsym)
    end
  end

  ## determine if this factor can be used
  # priors and general n-ary cases
  canuse = length(varsyms)==1 || (length(faillist)==1 && loovar in faillist)
  ## special multihypo case (at least one hypothesis is available or initializing first hypo)
  fctnode = getFactor(dfg, fct)
  # @show canuse, isMultihypo(fctnode), isinit
  if !canuse && isMultihypo(fctnode)
    # multihypo=[1;0.5;0.5] : sfidx=1, isinit=[0,1,0] -- true
    # multihypo=[1;0.5;0.5] : sfidx=1, isinit=[0,0,1] -- true
    # multihypo=[1;0.5;0.5] : sfidx=2|3, isinit=[1,0,0] -- true
    mhp = getMultihypoDistribution(fctnode).p
    allmhp,certainidx,uncertnidx = getHypothesesVectors(mhp)
    if isLeastOneHypoAvailable(sfidx, certainidx, uncertnidx, isinit)
      # special case works
      @info "allowing init from incomplete set of previously initialized hypotheses, fct=$fct"
      canuse = true
    end
  end

  # should add the factor for use?
  if canuse
    push!(fctlist, fct)
  end

  # return if can use, the factor in an array, and the non-initialized variables attached to the factor
  return (canuse, fctlist, faillist)::Tuple{Bool, Vector{Symbol}, Vector{Symbol}}
end

"""
    $(SIGNATURES)

EXPERIMENTAL: initialize target variable `xi` based on connected factors in the
factor graph `fgl`.  Possibly called from `addFactor!`, or `doCliqAutoInitUp!` (?).

Notes:
- Special carve out for multihypo cases, see issue 427.

Development Notes:
> Target factor is first (singletons) or second (dim 2 pairwise) variable vertex in `xi`.
* TODO use DFG properly with local operations and DB update at end.
* TODO get faster version of `isInitialized` for database version.
* TODO: Persist this back if we want to here.
"""
function doautoinit!( dfg::AbstractDFG,
                      xi::DFGVariable;
                      solveKey::Symbol=:default,
                      singles::Bool=true,
                      N::Int=getSolverParams(dfg).N, #maximum([length(getPoints(getBelief(xi, solveKey))); getSolverParams(dfg).N]),
                      logger=ConsoleLogger() )
  #
  didinit = false
  # don't initialize a variable more than once
  if !isInitialized(xi, solveKey)
    with_logger(logger) do
      @info "try doautoinit! of $(xi.label)"
    end
    # get factors attached to this variable xi
    vsym = xi.label
    neinodes = DFG.getNeighbors(dfg, vsym)
    # proceed if has more than one neighbor OR even if single factor
    if (singles || length(neinodes) > 1)
      # Which of the factors can be used for initialization
      useinitfct = Symbol[]
      # Consider factors connected to $vsym...
      for xifct in neinodes
        canuse, usefct, notusevars = factorCanInitFromOtherVars(dfg, xifct, vsym, solveKey=solveKey)
        if canuse
          union!(useinitfct, usefct)
        end
      end
      with_logger(logger) do
        @info "init with useinitfct $useinitfct"
      end
      # println("Consider all singleton (unary) factors to $vsym...")
      # calculate the predicted belief over $vsym
      if length(useinitfct) > 0
        with_logger(logger) do
          @info "do init of $vsym"
        end
        # FIXME ensure a product of only partial densities and returned pts are put to proper dimensions
        fcts = map(fx->getFactor(dfg,fx), useinitfct)
        bel,ipc = propagateBelief(dfg, getVariable(dfg,vsym), fcts, solveKey=solveKey, logger=logger, N=N)
        # while the propagate step might allow large point counts, the graph should stay restricted to N
        bel_ = Npts(bel) == getSolverParams(dfg).N ? bel : resample(bel, getSolverParams(dfg).N)
        # @info "MANIFOLD IS" bel.manifold isPartial(bel) string(bel._partial) string(getPoints(bel, false)[1])
        setValKDE!(xi, bel_, true, ipc, solveKey=solveKey) # getPoints(bel, false)
        # Update the estimates (longer DFG function used so cloud is also updated)
        setVariablePosteriorEstimates!(dfg, xi.label, solveKey)
        # Update the data in the event that it's not local
        # TODO perhaps usecopy=false
        updateVariableSolverData!(dfg, xi, solveKey, true; warn_if_absent=false)
        # deepcopy graphinit value, see IIF #612
        updateVariableSolverData!(dfg, xi.label, getSolverData(xi, solveKey), :graphinit, true, Symbol[]; warn_if_absent=false)
        didinit = true
      end
    end
  end
  return didinit
end

function doautoinit!( dfg::AbstractDFG,
                      Xi::Vector{<:DFGVariable};
                      solveKey::Symbol=:default,
                      singles::Bool=true,
                      N::Int=getSolverParams(dfg).N,
                      logger=ConsoleLogger()  )
  #
  #
  # Mighty inefficient function, since we only need very select fields nearby from a few neighboring nodes
  # do double depth search for variable nodes

  didinit = true

  # loop over all requested variables that must be initialized
  for xi in Xi
    didinit &= doautoinit!(dfg, xi, solveKey=solveKey, singles=singles, N=N, logger=logger)
  end
  return didinit
end

function doautoinit!( dfg::AbstractDFG,
                      xsyms::Vector{Symbol};
                      solveKey::Symbol=:default,
                      singles::Bool=true,
                      N::Int=getSolverParams(dfg).N,
                      logger=ConsoleLogger()  )
  #
  verts = getVariable.(dfg, xsyms)
  return doautoinit!(dfg, verts, solveKey=solveKey, singles=singles, N=N, logger=logger)
end
function doautoinit!( dfg::AbstractDFG,
                      xsym::Symbol;
                      solveKey::Symbol=:default,
                      singles::Bool=true,
                      N::Int=getSolverParams(dfg).N,
                      logger=ConsoleLogger()  )
  #
  return doautoinit!(dfg, [getVariable(dfg, xsym);], solveKey=solveKey, singles=singles, N=N, logger=logger)
end

"""
    $(TYPEDSIGNATURES)

Method to manually initialize a variable using a set of points.

Notes
- Disable automated graphinit on `addFactor!(fg, ...; graphinit=false)
  - any un-initialized variables will automatically be initialized by `solveTree!`

Example:

```julia
# some variable is added to fg
addVariable!(fg, :somepoint3, ContinuousEuclid{2})

# data is organized as (row,col) == (dimension, samples)
pts = randn(2,100)
initManual!(fg, :somepoint3, pts)

# manifold management should be done automatically.
# note upgrades are coming to consolidate with Manifolds.jl, see RoME #244

## it is also possible to initManual! by using existing factors, e.g.
initManual!(fg, :x3, [:x2x3f1])
```

DevNotes
- TODO better document graphinit and treeinit.
"""
function initManual!( variable::DFGVariable,
                      ptsArr::ManifoldKernelDensity,
                      solveKey::Symbol=:default;
                      dontmargin::Bool=false,
                      N::Int=length(getPoints(ptsArr)) )
  #
  @debug "initManual! $(getLabel(variable))"
  if !(solveKey in listSolveKeys(variable))
    @debug "$(getLabel(variable)) needs new VND solveKey=$(solveKey)"
    varType = getVariableType(variable)
    setDefaultNodeData!(variable, 0, N, getDimension(varType), solveKey=solveKey,
    initialized=false, varType=varType, dontmargin=dontmargin)
  end
  setValKDE!(variable, ptsArr, true, solveKey=solveKey)
  return nothing
end
function initManual!( dfg::AbstractDFG,
                      label::Symbol,
                      belief::ManifoldKernelDensity,
                      solveKey::Symbol=:default;
                      dontmargin::Bool=false,
                      N::Int=getSolverParams(dfg).N  )
  #
  variable = getVariable(dfg, label)
  initManual!(variable, belief, solveKey, dontmargin=dontmargin, N=N)
  return nothing
end
function initManual!( dfg::AbstractDFG,
                      label::Symbol,
                      usefcts::Vector{Symbol},
                      solveKey::Symbol=:default;
                      dontmargin::Bool=false,
                      N::Int=getSolverParams(dfg).N )
  #
  pts = predictbelief(dfg, label, usefcts, solveKey=solveKey)[1]
  vert = getVariable(dfg, label)
  Xpre = manikde!(getManifold(getVariableType(vert)), pts )
  initManual!(vert, Xpre, solveKey, dontmargin=dontmargin, N=N )
  # setValKDE!(vert, Xpre, true, solveKey=solveKey)
  # return nothing
end

function initManual!( dfg::AbstractDFG,
                      sym::Symbol,
                      pts::AbstractVector{P},
                      solveKey::Symbol=:default;
                      bw=nothing ) where {P}
  #
  var = getVariable(dfg, sym)
  M = getManifold(var)
  pp = manikde!(M, pts, bw=bw)
  initManual!(var,pp, solveKey)
end

const initVariableManual! = initManual!

"""
    $SIGNATURES

Set solveKey values of `dest::AbstractDFG` according to `initKey::Symbol=:graphinit` values.

Notes
- Some flexibility for using two DFGs and different key values, see Examples and code for details.
- Can also be specific with `varList::Vector{Symbol}`.
- Returns `dest` graph.
- Uses the supersolve mechanism.

Examples
```julia
resetInitialValues!(fg)
resetInitialValues!(fg1,fg2)  # into 1 from 2
resetInitialValues!(fg1,fg1,:myotherinit)  # use different init value into solveKey :default
resetInitialValues!(fg1,fg1,:graphinit, :mysolver) # not into solveKey=:default but :mysolver
resetInitialValues!(fg1, varList=[:x1;:l3])  # Specific variables only

# Into `fgNew` object, leaving `fg` untouched
fgNew = deepcopy(fg)
resetInitialValues!(fgNew,fg)
```

Related

initManual!, graphinit (keyword)
"""
function resetInitialValues!(dest::AbstractDFG,
                            src::AbstractDFG=dest,
                            initKey::Symbol=:graphinit,
                            solveKey::Symbol=:default;
                            varList::AbstractVector{Symbol}=ls(dest))
  #
  for vs in varList
    vnd = getSolverData(getVariable(src, vs), initKey)
    # guess we definitely want to use copy to preserve the initKey memory
    updateVariableSolverData!(dest,vs,vnd,solveKey,true; warn_if_absent=false)
  end
  return dest
end
const resetInitValues! = resetInitialValues!

"""
    $SIGNATURES

Ensure that no variables set as `solvable=1` are floating free without any connected `solvable=1` factors.  If any found, then set those 'free' variable's `solvable=solvableFallback` (default `0`).

Related

[`initAll!`](@ref)
"""
function ensureSolvable!(dfg::AbstractDFG; solvableTarget::Int=1, solvableFallback::Int=0)
  # workaround in case isolated variables occur
  solvVars = ls(dfg, solvable=solvableTarget)
  varHasFact = (x->length(ls(dfg,x, solvable=solvableTarget))==0).(solvVars)
  blankVars = solvVars[findall(varHasFact)]
  if 0 < length(blankVars)
    @warn("solveTree! dissallows solvable variables without any connected solvable factors -- forcing solvable=0 on $(blankVars)")
    (x->setSolvable!(dfg, x, solvableFallback)).(blankVars)
  end
  return blankVars
end

"""
    $SIGNATURES

Perform `graphinit` over all variables with `solvable=1` (default).

Related

ensureSolvable!, (EXPERIMENTAL 'treeinit')
"""
function initAll!(dfg::AbstractDFG,
                  solveKey::Symbol=:default;
                  solvable::Int=1,
                  N::Int=getSolverParams(dfg).N )
  #
  # allvarnodes = getVariables(dfg)
  syms = intersect(getAddHistory(dfg), ls(dfg, solvable=solvable) )
  # syms = ls(dfg, solvable=solvable) # |> sortDFG

  # May have to first add the solveKey VNDs if they are not yet available
  for sym in syms
    var = getVariable(dfg, sym)
    # does SolverData exist for this solveKey?
    if !( solveKey in listSolveKeys(var) )
      varType = getVariableType(var)
      # accept complete defaults for a novel solveKey
      setDefaultNodeData!(var, 0, N, getDimension(varType), solveKey=solveKey,
                          initialized=false, varType=varType, dontmargin=false)
    end
  end

  # do the init
  repeatCount = 0
  repeatFlag = true
  while repeatFlag
    repeatFlag = false
    repeatCount += 1
    if 10 < repeatCount
      @info "not able to initialize all variables via the factor graph, abort autoinit."
      break;
    end
    for sym in syms
      var = getVariable(dfg, sym)
      # is this SolverData initialized?
      if !isInitialized(var, solveKey)
        @info "$(var.label) is not initialized, and will do so now..."
        doautoinit!(dfg, [var;], solveKey=solveKey, singles=true)
        !isInitialized(var, solveKey) ? (repeatFlag = true) : nothing
      end
    end
  end
  nothing
end


text:
##############################################################################
##
## FixedEffect
##
##############################################################################

struct FixedEffect{R <: AbstractVector{<:Integer}, I <: AbstractVector{<:Real}}
	refs::R                 # refs must be between 0 and n
	interaction::I          # the continuous interaction
	n::Int                  # Number of potential values (= maximum(refs))
	function FixedEffect{R, I}(refs, interaction, n) where {R <: AbstractVector{<:Integer}, I <: AbstractVector{<: Real}}
		length(refs) == length(interaction) || throw(DimensionMismatch(
			"cannot match refs of length $(length(refs)) with interaction of length $(length(interaction))"))
		return new(refs, interaction, n)
	end
end

function FixedEffect(args...; interaction::AbstractVector = uweights(length(args[1])))
	g = GroupedArray(args..., sort = nothing)
	FixedEffect{typeof(g.groups), typeof(interaction)}(g.groups, interaction, g.ngroups)
end

Base.show(io::IO, ::FixedEffect) = print(io, "Fixed Effects")

function Base.show(io::IO, ::MIME"text/plain", fe::FixedEffect)
	print(io, fe, ':')
	print(io, "\n  refs (", length(fe.refs), "-element ", typeof(fe.refs), "):")
	print(io, "\n    [", string.(Int.(fe.refs[1:min(5, length(fe.refs))])).*", "..., "... ]")
	if fe.interaction isa UnitWeights
		print(io, "\n  interaction (UnitWeights):")
		print(io, "\n    none")
	else
		print(io, "\n  interaction (", length(fe.interaction), "-element ", typeof(fe.interaction), "):")
		print(io, "\n    [", (sprint(show, x; context=:compact=>true)*", " for x in fe.interaction[1:min(5, length(fe.interaction))])..., "... ]")
	end
end

Base.size(fe::FixedEffect) = size(fe.refs)
Base.length(fe::FixedEffect) = length(fe.refs)
Base.eltype(::FixedEffect{R,I}) where {R,I} = eltype(I)

Base.getindex(fe::FixedEffect, ::Colon) = fe

@propagate_inbounds function Base.getindex(fe::FixedEffect, esample)
	@boundscheck checkbounds(fe.refs, esample)
	@boundscheck checkbounds(fe.interaction, esample)
	@inbounds refs = fe.refs[esample]
	@inbounds interaction = fe.interaction[esample]
	return FixedEffect{typeof(fe.refs), typeof(fe.interaction)}(refs, interaction, fe.n)
end
##############################################################################
##
## Find connected components
##
##############################################################################
# Return a vector of sets that contains the indices of each unique value
function refsrev(fe::FixedEffect)
	out = Vector{Int}[Int[] for _ in 1:fe.n]
	for i in eachindex(fe.refs)
		push!(out[fe.refs[i]], i)
	end
	return out
end

# Returns a vector of all components
# A component is a vector that, for each fixed effect,
# contains all the refs that are included in the component.
function components(fes::AbstractVector{<:FixedEffect})
	refs_vec = Vector{Int}[fe.refs for fe in fes]
	refsrev_vec = Vector{Vector{Int}}[refsrev(fe) for fe in fes]
	visited = falses(length(refs_vec[1]))
	out = Vector{Set{Int}}[]
	for i in eachindex(visited)
		if !visited[i]
			# obs not visited yet, so create new component
			component_vec = Set{Int}[Set{Int}() for _ in 1:length(refsrev_vec)]
			# visit all obs in the same components
			tovisit = Set{Int}(i)
			while !isempty(tovisit)
				for (component, refs, refsrev) in zip(component_vec, refs_vec, refsrev_vec)
					ref = refs[i]
					# if group is not in component yet
					if ref ∉ component
						# add group to the component
						push!(component, ref)
						# visit other observations in same group
						union!(tovisit, refsrev[ref])
					end
				end
				# mark obs as visited
				i = pop!(tovisit)
				visited[i] = true
			end
			push!(out, component_vec)
		end
	end
	return out
end

##############################################################################
##
## normalize! a vector of fixedeffect coefficients using connected components
##
##############################################################################

function normalize!(fecoefs::AbstractVector{<: Vector{<: Real}}, fes::AbstractVector{<:FixedEffect}; kwargs...)
	# The solution is generally not unique. Find connected components and scale accordingly
	idx = findall(fe -> isa(fe.interaction, UnitWeights), fes)
	length(idx) >= 2 && rescale!(view(fecoefs, idx), view(fes, idx))
	return fecoefs
end

function rescale!(fecoefs::AbstractVector{<: Vector{<: Real}}, fes::AbstractVector{<:FixedEffect})
	for component_vec in components(fes)
		m = 0.0
		# demean all fixed effects except the first
		for j in length(fecoefs):(-1):2
			fecoef, component = fecoefs[j], component_vec[j]
			mj = 0.0
			for k in component
				mj += fecoef[k]
			end
			mj = mj / length(component)
			for k in component
				fecoef[k] -= mj
			end
			m += mj
		end
		# rescale the first fixed effects
		fecoef, component = fecoefs[1], component_vec[1]
		for k in component
			fecoef[k] += m
		end
	end
end

function full(fecoefs::AbstractVector{<: Vector{<: Real}}, fes::AbstractVector{<:FixedEffect})
	[fecoef[fe.refs] for (fecoef, fe) in zip(fecoefs, fes)]
end

max_stars_repo_name: matthieugomez/FixedEffects

max_issues_repo_name: matthieugomez/FixedEffects

max_forks_repo_name: matthieugomez/FixedEffects


text:
using RDatasets, DataFrames, FixedEffectModels, Base.Test

x = dataset("plm", "Cigar")
x[:id1] = x[:State]
x[:id2] = x[:Year]
x[:pid1] = pool(x[:id1])
x[:ppid1] = pool(div(x[:id1], 10))
x[:pid2] = pool(x[:id2])
x[:y] = x[:Sales]
x[:x1] = x[:Price]
x[:z1] = x[:Pimin]
x[:x2] = x[:NDI]
x[:w] = x[:Pop]

##############################################################################
##
## coefficients
##
##############################################################################

# simple
@test_approx_eq_eps coef(reg(y ~ x1, x))  [139.73446,-0.22974] 1e-4
@test_approx_eq_eps coef(reg(y ~ x1, x, weight = :w))  [137.72495428982756,-0.23738] 1e-4

x[:SalesInt] = round(Int64, x[:Sales])
@test_approx_eq_eps coef(reg(SalesInt ~ Price, x)) [139.72674,-0.2296683205] 1e-4

# absorb
@test_approx_eq_eps coef(reg(y ~ x1 |> pid1, x))  [-0.20984] 1e-4
@test reg(y ~ x1 |> pid1, x).iterations == 1

@test_approx_eq_eps coef(reg(y ~ x1 |> pid1 + pid2, x))   [-1.08471] 1e-4
@test_approx_eq_eps coef(reg(y ~ x1 |> pid1 + pid1&id2 , x))    [-0.53470] 1e-4
@test_approx_eq_eps coef(reg(y ~ x1 |> pid1*id2 , x))    [-0.53470] 1e-4
#@test isempty(coef(reg(y ~ 0 |> pid1*x1, x)))
@test_approx_eq_eps coef(reg(y ~ x1 |> ppid1&pid2 , x))    [-1.44255] 1e-4

@test_approx_eq_eps  coef(reg(y ~ x1 |> pid1&id2, x))  [13.993028174622104,-0.5804357763515606] 1e-4
@test_approx_eq_eps  coef(reg(y ~ x1 |> id2&pid1 , x))  [13.993028174622104,-0.5804357763515606] 1e-4
@test_approx_eq_eps  coef(reg(y ~ 1 |> id2&pid1 , x))  [174.4084407796102] 1e-4
@test_approx_eq_eps  coef(reg(y ~ x1 |> pid1&id2 + pid2&id1, x))  [51.2358,-0.5797] 1e-4
@test_approx_eq_eps  coef(reg(y ~ x1 + x2 |> pid1&id2 + pid2&id1, x))  [-46.4464,-0.2546, -0.005563] 1e-4
@test_approx_eq_eps  coef(reg(y ~ 0 + x1 + x2 |> pid1&id2 + pid2&id1, x))   [-0.21226562244177932,-0.004775616634862829] 1e-4

# recheck these two below
@test_approx_eq_eps coef(reg(y ~ z1 |> x1&x2&pid1 , x))    [122.98713,0.30933] 1e-4
@test_approx_eq_eps coef(reg(y ~ z1 |> (x1&x2)*pid1 , x))    [0.421406] 1e-4

# absorb + weights
@test_approx_eq_eps coef(reg(y ~ x1 |> pid1, x, weight = :w)) [-0.21741] 1e-4
@test_approx_eq_eps coef(reg(y ~ x1 |> pid1 + pid2, x,  weight = :w))  [-0.88794] 1e-3
@test_approx_eq_eps coef(reg(y ~ x1 |> pid1 + pid1&id2, x, weight = :w))  [-0.461085492] 1e-4

# iv
@test_approx_eq_eps coef(reg(y ~ (x1 = z1), x))  [138.19479,-0.20733] 1e-4
@test_approx_eq_eps coef(reg(y ~ x2 + (x1 = z1), x))   [137.45096,0.00516,-0.76276] 1e-4
@test_approx_eq_eps coef(reg(y ~ x2 + (x1 = z1 + w), x))  [137.57335,0.00534,-0.78365] 1e-4
## multiple endogeneous variables
@test_approx_eq_eps coef(reg(y ~ (x1 + x2 = z1 + w), x))  [139.544, .8001, -.00937] 1e-4
@test_approx_eq_eps coef(reg(y ~ 1 + (x1 + x2 = z1 + w), x))  [139.544, .8001, -.00937] 1e-4
result = [196.576, 0.00490989, -2.94019, -3.00686, -2.94903, -2.80183, -2.74789, -2.66682, -2.63855, -2.52394, -2.34751, -2.19241, -2.18707, -2.09244, -1.9691, -1.80463, -1.81865, -1.70428, -1.72925, -1.68501, -1.66007, -1.56102, -1.43582, -1.36812, -1.33677, -1.30426, -1.28094, -1.25175, -1.21438, -1.16668, -1.13033, -1.03782]
@test_approx_eq_eps coef(reg(y ~ x2 + (x1&pid2 = z1&pid2), x)) result 1e-4

# iv + weight
@test_approx_eq_eps coef(reg(y ~ (x1 = z1), x, weight = :w))  [137.03637,-0.22802] 1e-4

# iv + weight + absorb
@test_approx_eq_eps coef(reg(y ~ (x1 = z1) |> pid1, x)) [-0.20284] 1e-4
@test_approx_eq_eps coef(reg(y ~ (x1 = z1) |> pid1, x, weight = :w)) [-0.20995] 1e-4
@test_approx_eq_eps coef(reg(y ~ x2 + (x1 = z1), x))   [137.45096580480387,0.005169677634275297,-0.7627670265757879] 1e-4
@test_approx_eq_eps coef(reg(y ~ x2 + (x1 = z1) |> pid1, x))    [0.0011021722526916768,-0.3216374943695231] 1e-4

# non high dimensional factors
@test_approx_eq_eps coef(reg(y ~ x1 + pid2 |> pid1, x))[1]   -1.08471 1e-4
@test_approx_eq_eps coef(reg(y ~ x1 + pid2 |> pid1, x,  weight = :w))[1]   -0.88794 1e-4
@test_approx_eq_eps coef(reg(y ~ x2 + (x1 = z1) + pid2 |> pid1, x))[1]   -0.00525 1e-4

##############################################################################
##
## collinearity
## add more tests
##
##############################################################################
# ols
x[:x12] = x[:x1]
model = reg(y ~ x1 + x12, x)
@test_approx_eq_eps coef(model) [139.7344639806166,-0.22974688593485126,0.0] 1e-4

# iv
x[:x22] = x[:x2]
model = reg(y ~  x22 + x2 + (x1 = z1), x)
@test  coef(model)[2] == 0 || coef(model)[3] == 0

x[:zz1] = x[:z1]
model = reg(y ~  zz1 + (x1 = x2 + z1), x)
@test coef(model)[2] != 0.0

# catch when IV underidentified : re-try when 0.5
#@test_throws ErrorException reg(y ~ x1 + (x2 + w = x2), x)

# catch continuous variables in fixed effects
@test_throws ErrorException reg(y ~ x1 |> x2, x)
@test_throws ErrorException reg(y ~ x1 |> x2 + pid1, x)

##############################################################################
##
## std errors
##
##############################################################################

# Simple
@test_approx_eq_eps stderr(reg(y ~ x1 , x)) [1.52126,0.01889] 1e-4
# Stata ivreg
@test_approx_eq_eps stderr(reg(y ~ (x1 = z1) , x)) [1.53661,0.01915] 1e-4
# Stata areg
@test_approx_eq_eps stderr(reg(y ~ x1 |> pid1, x))  [0.00980] 1e-4

# White
# Stata reg
@test_approx_eq_eps stderr(reg(y ~ x1 , x, VcovWhite())) [1.68679,0.01670] 1e-4
# Stata ivreg
@test_approx_eq_eps stderr(reg(y ~ (x1 = z1) , x, VcovWhite()))  [1.63305,0.01674] 1e-4
# Stata areg
@test_approx_eq_eps stderr(reg(y ~ x1 |> pid1, x, VcovWhite())) [0.01100] 1e-4

# cluster
@test_approx_eq_eps stderr(reg(y ~ x1 , x, VcovCluster(:pid1)))[2]  0.03792 1e-4
@test_approx_eq_eps stderr(reg(y ~ x1 |> pid1, x, VcovCluster(:pid2)))  [0.02205] 1e-4
# stata reghxe
@test_approx_eq_eps stderr(reg(y ~ x1 |> pid1, x, VcovCluster(:pid1)))  [0.03573] 1e-4
# no reference
@test_approx_eq_eps stderr(reg(y ~ x1, x, VcovCluster([:pid1, :pid2])))[1]  6.17025 1e-4
# no reference
@test_approx_eq_eps stderr(reg(y ~ x1 |> pid1, x, VcovCluster([:pid1, :pid2])))[1]   0.04037 1e-4

@test_throws ErrorException reg(y ~ x1 , x, VcovCluster(:State))

##############################################################################
##
## subset
##
##############################################################################

result = reg(y ~ x1 + pid1, x, subset = x[:State] .<= 30)
@test length(result.esample) == size(x, 1)
@test_approx_eq_eps coef(result)  coef(reg(y ~ x1 + pid1, x[x[:State] .<= 30, :])) 1e-4
@test_approx_eq_eps vcov(result)  vcov(reg(y ~ x1 + pid1, x[x[:State] .<= 30, :])) 1e-4
@test_throws ErrorException reg(y ~ x1, x, subset = x[:pid2] .<= 30)

##############################################################################
##
## F Stat
##
##############################################################################

@test_approx_eq_eps reg(y ~ x1, x).F  147.82425 1e-4
@test_approx_eq_eps reg(y ~ x1 |> pid1, x).F  458.45825 1e-4
@test_approx_eq_eps reg(y ~ (x1 = z1), x).F  117.17329 1e-4

@test_approx_eq_eps reg(y ~ x1, x, VcovCluster(:pid1)).F    36.70275 1e-4
@test_approx_eq_eps reg(y ~ (x1 = z1), x, VcovCluster(:pid1)).F  39.67227 1e-4

#  xtivreg2
@test_approx_eq_eps reg(y ~ (x1 = z1) |> pid1, x).F  422.46444 1e-4

##############################################################################
##
## F_kp r_kp statistics for IV. Difference degrees of freedom.
##
##############################################################################

@test_approx_eq_eps reg(y ~ (x1 = z1), x).F_kp    52248.79247 1e-4
@test_approx_eq_eps reg(y ~ x2 + (x1 = z1), x).F_kp    5159.812208193612 1e-4
@test_approx_eq_eps reg(y ~ (x1 = z1 + CPI), x).F_kp   27011.44080 1e-4
@test_approx_eq_eps reg(y ~ (x1 = z1) |> pid1, x).F_kp 100927.75710 1e-4

# exactly same with ivreg2
@test_approx_eq_eps reg(y ~ (x1 = z1), x, VcovWhite()).F_kp    23160.06350 1e-4
@test_approx_eq_eps reg(y ~ (x1 = z1) |> pid1, x, VcovWhite()).F_kp  37662.82808 1e-4
@test_approx_eq_eps reg(y ~ CPI + (x1 = z1), x, VcovWhite()).F_kp    2093.46609 1e-4
@test_approx_eq_eps reg(y ~ (x1 = z1 + CPI), x, VcovWhite()).F_kp  16418.21196 1e-4

# like in ivreg2 but += 5 difference. there is a degrees of freedom difference. not sure where.
@test_approx_eq_eps reg(y ~ (x1 = z1), x, VcovCluster(:pid1)).F_kp     7249.88606 1e-4
@test_approx_eq_eps reg(y ~ CPI + (x1 = z1), x, VcovCluster(:pid1)).F_kp   538.40393 1e-4
@test_approx_eq_eps reg(y ~ CPI + (x1 = z1), x, VcovCluster([:pid1, :pid2])).F_kp   423.00342 1e-4
@test_approx_eq_eps reg(y ~ (x1 = z1 + CPI), x, VcovCluster(:pid1)).F_kp      4080.66081 1e-4
@test_approx_eq_eps reg(y ~  (x1 = z1 + CPI), x, VcovCluster([:pid1, :pid2])).F_kp     2877.94477 1e-4

##############################################################################
##
## Test unbalanced panel
##
## corresponds to abdata in Stata, for instance reghxe wage emp [w=indoutpt], a(id year)
##
##############################################################################

x = dataset("plm", "EmplUK")
x[:id1] = x[:Firm]
x[:id2] = x[:Year]
x[:pid1] = pool(x[:id1])
x[:pid2] = pool(x[:id2])
x[:y] = x[:Wage]
x[:x1] = x[:Emp]
x[:w] = x[:Output]

for method in [:cholesky, :qr, :lsmr]
	# absorb
	@test_approx_eq_eps coef(reg(y ~ x1 |> pid1, x, method = method))  [-0.11981270017206136] 1e-4
	@test_approx_eq_eps  coef(reg(y ~ x1 |> pid1&id2, x, method = method))   [-315.0000747500431,-0.07633636891202833] 1e-4
	#
	@test_approx_eq_eps  coef(reg(y ~ x1 |> id2&pid1 , x, method = method))   [-315.0000747500431,-0.07633636891202833] 1e-4

	@test_approx_eq_eps  coef(reg(y ~ 1 |> id2&pid1 , x, method = method))   [-356.40430526316396] 1e-4
	@test_approx_eq_eps coef(reg(y ~ x1 |> pid1, x, method = method, weight = :w))  [-0.11514363590574725] 1e-4

	# absorb + weights
	@test_approx_eq_eps coef(reg(y ~ x1 |> pid1 + pid2, x, method = method))    [-0.04683333721137311] 1e-4
	@test_approx_eq_eps coef(reg(y ~ x1 |> pid1 + pid2, x, method = method,  weight = :w))   [-0.043475472188120416] 1e-3

	## the last two ones test an ill conditioned model matrix
	@test_approx_eq_eps coef(reg(y ~ x1 |> pid1 + pid1&id2 , x, method = method))     [-0.122354] 1e-4
	@test reg(y ~ x1 |> pid1 + pid1&id2 , x, method = method).iterations <= 30

	@test_approx_eq_eps coef(reg(y ~ x1 |> pid1 + pid1&id2, x, method = method, weight = :w))  [-0.11752306001586807] 1e-4
	@test reg(y ~ x1 |> pid1 + pid1&id2, x, method = method, weight = :w).iterations <= 50
end


text:
# on-device array type

export AbstractDeviceArray

## device array

"""
    AbstractDeviceArray{T, N} <: DenseArray{T, N}

Supertype for `N`-dimensional GPU arrays (or array-like types) with elements of type `T`.
Instances of this type are expected to live on the device, see [`AbstractGPUArray`](@ref)
for device-side objects.
"""
abstract type AbstractDeviceArray{T, N} <: DenseArray{T, N} end

Base.IndexStyle(::AbstractDeviceArray) = IndexLinear()

@inline function Base.iterate(A::AbstractDeviceArray, i=1)
    if (i % UInt) - 1 < length(A)
        (@inbounds A[i], i + 1)
    else
        nothing
    end
end

function Base.sum(A::AbstractDeviceArray{T}) where T
    acc = zero(T)
    for elem in A
        acc += elem
    end
    acc
end


text:
# look at average days per year which are snow covered during the timeseries (look at different elevations)
# look at amount that snow contributes to runoff (Gesamtvolumen pro Jahr)

using DelimitedFiles
using Plots
using Statistics
using StatsPlots
using Plots.PlotMeasures
using CSV
using Dates
using DocStringExtensions

path_45 = "/home/sarah/Master/Thesis/Data/Projektionen/new_station_data_rcp45/rcp45/"
Name_Projections_45 = readdir(path_45)
path_85 = "/home/sarah/Master/Thesis/Data/Projektionen/new_station_data_rcp85/rcp85/"
Name_Projections_85 = readdir(path_85)
"""
Computes the mean yearly number of snow covered days at each elevation of the catchment for all best parameter sets.

$(SIGNATURES)

The function returns the mean yearly number of snow covered days as an array (12xnr_runs) It takes as input the snow cover data, a timeseries and the elevation zones of the catchment.
"""
function snow_covered_days(Snow_Cover::Array{Float64,2}, Timeseries, Elevations, nr_runs)
    # Snow_Cover is [elevations*runs, timeseries]
    Years = collect(Dates.year(Timeseries[1]): Dates.year(Timeseries[end]))
    Days_Snow_Cover = zeros(nr_runs*length(Elevations))
    #Date_max_Annual_Discharge = Float64[]
    for (i, Current_Year) in enumerate(Years)
        Dates_Current_Year = filter(Timeseries) do x
                                          Dates.Year(x) == Dates.Year(Current_Year) &&
                                          (Dates.Month(x) == Dates.Month(9) ||
                                          Dates.Month(x) == Dates.Month(10) ||
                                          Dates.Month(x) == Dates.Month(11) ||
                                          Dates.Month(x) == Dates.Month(12))
                                    end
        Dates_Next_Year = filter(Timeseries) do x
                                          Dates.Year(x) == Dates.Year(Current_Year+1) &&
                                          (Dates.Month(x) == Dates.Month(1) ||
                                          Dates.Month(x) == Dates.Month(2) ||
                                          Dates.Month(x) == Dates.Month(3) ||
                                          Dates.Month(x) == Dates.Month(4) ||
                                          Dates.Month(x) == Dates.Month(5) ||
                                          Dates.Month(x) == Dates.Month(6) ||
                                          Dates.Month(x) == Dates.Month(7) ||
                                          Dates.Month(x) == Dates.Month(8))
                                    end

        append!(Dates_Current_Year, Dates_Next_Year)
        Snow_Cover_Current_Year = Snow_Cover[:, indexin(Dates_Current_Year, Timeseries)]
        Nr_Days_Snow = Float64[]
        for run in 1:nr_runs*length(Elevations)
            append!(Nr_Days_Snow, length(findall(x->x > 0.5, Snow_Cover_Current_Year[run,:])))
        end

        # put all the years into one array
        Days_Snow_Cover = hcat(Days_Snow_Cover, Nr_Days_Snow)
    end
    # the first column is only zeros, the last column does not contain a whole year, so both have to be deleted
    Days_Snow_Cover = Days_Snow_Cover[:,2:end-1]
    # then get the mean value of snow covered days during the timeperiod, dims=2 takes mean over columns (takes mean over one row)
    Mean_Nr_Days_Snow = mean(Days_Snow_Cover, dims = 2)
    # now reshape the array so that results of each elevation are in one row
    Mean_Nr_Days_Snow = reshape(Mean_Nr_Days_Snow,length(Elevations), nr_runs)

    return Mean_Nr_Days_Snow
end

"""
Computes the mean yearly number of snow covered days at each elevation of the catchment for all projections.

$(SIGNATURES)

The function returns the mean yearly number of snow covered days as an array (12x(14*nr_runs) It takes as input the path to the projection, the catchment name and the elvations of the catchment.
"""
function snow_cover_days_projections(path_to_projections, Catchment_Name, Elevations_Catchment, nr_runs)
    Name_Projections = readdir(path_to_projections)
    Timeseries_Past = collect(Date(1981,1,1):Day(1):Date(2010,12,31))
    Timeseries_End = readdlm("/home/sarah/Master/Thesis/Data/Projektionen/End_Timeseries_45_85.txt",',')
    Mean_Nr_Days_Snow_Future_All_Proj = zeros(length(Elevations_Catchment))

    if path_to_projections[end-2:end-1] == "45"
        index = 1
        rcp = "45"
        print(rcp, " ", path_to_projections)
    elseif path_to_projections[end-2:end-1] == "85"
        index = 2
        rcp="85"
        print(rcp, " ", path_to_projections)
    end

    for i in 13:14
        @time begin
        name = Name_Projections[i]
        print("run ", i)
        println("name ", name, "\n")
        Timeseries_Future = collect(Date(Timeseries_End[i,index]-29,1,1):Day(1):Date(Timeseries_End[i,index],12,31))
        #Snow_Cover_Past = readdlm(path_to_projections*name*"/"*Catchment_Name*"/300_model_results_snow_cover_past_2010.csv", ',')
        Snow_Cover_Future = readdlm(path_to_projections*name*"/"*Catchment_Name*"/300_model_results_snow_cover_future_2100.csv", ',')
        #Mean_Nr_Days_Snow_Past = snow_covered_days(Snow_Cover_Past, Timeseries_Past, Elevations_Catchment, nr_runs)
        Mean_Nr_Days_Snow_Future = snow_covered_days(Snow_Cover_Future, Timeseries_Future, Elevations_Catchment, nr_runs)
        #Mean_Nr_Days_Snow_Past = transpose(Mean_Nr_Days_Snow_Past)
        Mean_Nr_Days_Snow_Future = transpose(Mean_Nr_Days_Snow_Future)

        # open("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/Mean_Nr_Days_Snow_Past_"*rcp*".csv", "a") do io
        #         writedlm(io, Mean_Nr_Days_Snow_Past,",")
        # end
        open("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/Mean_Nr_Days_Snow_Future_"*rcp*".csv", "a") do io
                writedlm(io, Mean_Nr_Days_Snow_Future,",")
        end
        end
    end
end

# Elevation_Gailtal = collect(500:200:2700)
# Elevation_Pitten = collect(500:200:1500)
# Elevation_Palten = collect(700:200:2500)
Elevation_Silbertal = collect(700.0:200:2700.0)
Elevation_Defreggental = collect(1100:200:3500)
Elevation_Catchment = Elevation_Silbertal
#Timeseries_Past = collect(Date(1981,1,1):Day(1):Date(2010,12,31))
#Timeseries_End = readdlm("/home/sarah/Master/Thesis/Data/Projektionen/End_Timeseries_45_85.txt",',')
#snow_cover_days_projections(path_45, "IllSugadin", Elevation_Silbertal, 300)

# Catchment_Name = "IllSugadin"
# Mean_Nr_Days_Snow_Past_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/Mean_Nr_Days_Snow_Past_45.csv", ',')
# Mean_Nr_Days_Snow_Future_45 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/Mean_Nr_Days_Snow_Future_45.csv", ',')
# Mean_Nr_Days_Snow_Past_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/Mean_Nr_Days_Snow_Past_85.csv", ',')
# Mean_Nr_Days_Snow_Future_85 = readdlm("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/Mean_Nr_Days_Snow_Future_85.csv", ',')
# Mean_Nr_Days_Snow_Past_45 = transpose(Mean_Nr_Days_Snow_Past_45)
# Mean_Nr_Days_Snow_Past_85 = transpose(Mean_Nr_Days_Snow_Past_85)
# Mean_Nr_Days_Snow_Future_45 = transpose(Mean_Nr_Days_Snow_Future_45)
# Mean_Nr_Days_Snow_Future_85 = transpose(Mean_Nr_Days_Snow_Future_85)

"""
Plots the mean yearly number of snow covered days at each elevation of the catchment.

$(SIGNATURES)

The function plots the relative and absolute changes in snow covered days as well as a comparison of absolute values of past adn future.
"""
function plot_snow_cover(Mean_Nr_Days_Snow_Past_45, Mean_Nr_Days_Snow_Future_45, Mean_Nr_Days_Snow_Past_85, Mean_Nr_Days_Snow_Future_85, Catchment_Name, Elevations_Catchment)
    #print(Elevations_Catchment, "\n")
    Farben45=palette(:blues)
    Farben85=palette(:reds)
    plot()
    for (i, elevation) in enumerate(Elevations_Catchment)
        boxplot!(relative_error(Mean_Nr_Days_Snow_Future_45[i,:], Mean_Nr_Days_Snow_Past_45[i,:])*100, size=(2000,800), leg=false, color="blue")
        boxplot!(relative_error(Mean_Nr_Days_Snow_Future_85[i,:], Mean_Nr_Days_Snow_Past_85[i,:])*100,size=(2000,800), left_margin = [5mm 0mm], leg=false, color="red")
    end
    ylabel!("Relative Change [%]")
    title!("Relative Change in Annual Number of Snow Days per Elevation in "*Catchment_Name)
    xticks!([1.5:2:length(Elevations_Catchment)*2-0.5;], string.(Elevations_Catchment))
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/"*Catchment_Name*"_rel_change.png")
    plot()
    for (i, elevation) in enumerate(Elevations_Catchment)
        boxplot!(Mean_Nr_Days_Snow_Future_45[i,:] - Mean_Nr_Days_Snow_Past_45[i,:], size=(2000,800), leg=false, color="blue")
        boxplot!(Mean_Nr_Days_Snow_Future_85[i,:] - Mean_Nr_Days_Snow_Past_85[i,:],size=(2000,800), left_margin = [5mm 0mm], leg=false, color="red")
    end
    ylabel!("Absolute Change [Days]")
    title!("Absolute Change in Annual Number of Snow Days per Elevation in "*Catchment_Name)
    xticks!([1.5:2:length(Elevations_Catchment)*2-0.5;], string.(Elevations_Catchment))
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/"*Catchment_Name*"_abs_change.png")
    # plot true data
    plot()
    for (i, elevation) in enumerate(Elevations_Catchment)
        boxplot!(Mean_Nr_Days_Snow_Past_45[i,:], size=(2000,800), leg=false, color=[Farben45[1]])
        boxplot!(Mean_Nr_Days_Snow_Future_45[i,:],size=(2000,800), left_margin = [5mm 0mm], leg=false, color=[Farben45[2]])
    end
    ylabel!("Annual Number of Days with Snow")
    title!("Annual Number of Snow Days per Elevation in "*Catchment_Name*" Past vs. Future RCP 4.5")
    xticks!([1.5:2:length(Elevations_Catchment)*2-0.5;], string.(Elevations_Catchment))
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/"*Catchment_Name*"_comparison_past_futureRCP45.png")

    plot()
    for (i, elevation) in enumerate(Elevations_Catchment)
        boxplot!(Mean_Nr_Days_Snow_Past_85[i,:], size=(2000,800), leg=false, color=[Farben85[1]])
        boxplot!(Mean_Nr_Days_Snow_Future_85[i,:],size=(2000,800), left_margin = [5mm 0mm], leg=false, color=[Farben85[2]])
    end
    ylabel!("Annual Number of Days with Snow")
    title!("Annual Number of Snow Days per Elevation in "*Catchment_Name* " Past vs. Future RCP 8.5")
    xticks!([1.5:2:length(Elevations_Catchment)*2-0.5;], string.(Elevations_Catchment))
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/"*Catchment_Name*"_comparison_past_futureRCP85.png")

    # plot differently
    plot()
    mean_change_45 = Float64[]
    max_change_45 = Float64[]
    min_change_45 = Float64[]
    mean_change_85 = Float64[]
    max_change_85 = Float64[]
    min_change_85 = Float64[]
    for (i, elevation) in enumerate(Elevations_Catchment)
        change_45 = relative_error(Mean_Nr_Days_Snow_Future_45[i,:], Mean_Nr_Days_Snow_Past_45[i,:])*100
        change_85 = relative_error(Mean_Nr_Days_Snow_Future_85[i,:], Mean_Nr_Days_Snow_Past_85[i,:])*100
        append!(mean_change_45, mean(change_45))
        append!(max_change_45, maximum(change_45))
        append!(min_change_45, minimum(change_45))
        append!(mean_change_85, mean(change_85))
        append!(max_change_85, maximum(change_85))
        append!(min_change_85, minimum(change_85))
    end

    plot(Elevations_Catchment, mean_change_45,  olor=[Farben45[2]], label="RCP 4.5", ribbon = (mean_change_45 - min_change_45, max_change_45 - mean_change_45))
    plot!(Elevations_Catchment, mean_change_85, color=[Farben85[2]], label="RCP 8.5", ribbon = (mean_change_85 - min_change_85, max_change_85 - mean_change_85), size=(1600,800))
    ylabel!("Relative Change [%]")
    xlabel!("Elevation Zones")
    title!("Relative Change in Annual Number of Snow Days per Elevation in "*Catchment_Name)
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/"*Catchment_Name*"_rel_change_2.png")
end

#plot_snow_cover(Mean_Nr_Days_Snow_Past_45, Mean_Nr_Days_Snow_Future_45,  Mean_Nr_Days_Snow_Past_85, Mean_Nr_Days_Snow_Future_85, Catchment_Name, Elevation_Catchment)

"""
Computes the mean yearly amount of snow melt as contribution to hydrological system in the timeseries.

$(SIGNATURES)

The function returns the mean yearly amount of snow melt [mm]. It takes as input, snowstorage and corresponding timeseries.
"""
function snow_contribution_yearly(Snow_Melt, Timeseries)
    Years = collect(Dates.year(Timeseries[1]): Dates.year(Timeseries[end]))
    Snow_Contribution_All_Years = Float64[]
    #Date_max_Annual_Discharge = Float64[]
    for (i, Current_Year) in enumerate(Years)
        Dates_Current_Year = filter(Timeseries) do x
                                          Dates.Year(x) == Dates.Year(Current_Year) &&
                                          (Dates.Month(x) == Dates.Month(10) ||
                                          Dates.Month(x) == Dates.Month(11) ||
                                          Dates.Month(x) == Dates.Month(12))
                                    end
        Dates_Next_Year = filter(Timeseries) do x
                                          Dates.Year(x) == Dates.Year(Current_Year+1) &&
                                          (Dates.Month(x) == Dates.Month(1) ||
                                          Dates.Month(x) == Dates.Month(2) ||
                                          Dates.Month(x) == Dates.Month(3) ||
                                          Dates.Month(x) == Dates.Month(4) ||
                                          Dates.Month(x) == Dates.Month(5) ||
                                          Dates.Month(x) == Dates.Month(6) ||
                                          Dates.Month(x) == Dates.Month(7) ||
                                          Dates.Month(x) == Dates.Month(8) ||
                                          Dates.Month(x) == Dates.Month(9))
                                    end

        append!(Dates_Current_Year, Dates_Next_Year)
        Snow_Contribution_Yearly = sum(Snow_Melt[indexin(Dates_Current_Year, Timeseries)])
        append!(Snow_Contribution_All_Years, Snow_Contribution_Yearly)
    end
    #because last years data only includes have the winter
    Snow_Contribution_All_Years = Snow_Contribution_All_Years[1:end-1]
    #returns the mean volume of precipitation stored as snow
    return mean(Snow_Contribution_All_Years)::Float64
end

"""
Computes the monthly snow melt of the past and future and the relative changes, of the projections of the path and the different parameter sets

$(SIGNATURES)

The function returns relative change in monthly sum of snow melt, the monthly sum of snow melt of the past, and the monthly mean snow melt of the future
"""
function change_monthly_snowmelt(path_to_projections, Catchment_Name)
    Name_Projections_45 = readdir(path_to_projections)
    Timeseries_Past = collect(Date(1981,1,1):Day(1):Date(2010,12,31))
    Timeseries_End = readdlm("/home/sarah/Master/Thesis/Data/Projektionen/End_Timeseries_45_85.txt",',')
    if path_45[end-2:end-1] == "45"
        index = 1
        rcp = "45"
        print(rcp, " ", path_to_projections)
    elseif path_45[end-2:end-1] == "85"
        index = 2
        rcp="85"
        print(rcp, " ", path_to_projections)
    end
    average_monthly_Discharge_past = Float64[]
    average_monthly_Discharge_future = Float64[]
    error_average_monthly_Discharge_all_runs = Float64[]
    all_months_all_runs = Float64[]
    for (i, name) in enumerate(Name_Projections_45)
        Timeseries_Future_45 = collect(Date(Timeseries_End[i,index]-29,1,1):Day(1):Date(Timeseries_End[i,index],12,31))
        Past_Discharge_45 = readdlm(path_to_projections*name*"/"*Catchment_Name*"/300_model_results_snow_melt_snow_redistr_past_2010.csv", ',')
        Future_Discharge_45 = readdlm(path_to_projections*name*"/"*Catchment_Name*"/300_model_results_snow_melt_snow_redistr_future_2100.csv", ',')
        println(size(Past_Discharge_45)[1])
        for run in 1:size(Past_Discharge_45)[1]
            # computes sum of monthly snow melt of each month
            Monthly_Discharge_past, Month = monthly_snowmelt(Past_Discharge_45[run,:], Timeseries_Past)
            Monthly_Discharge_future, Month_future = monthly_snowmelt(Future_Discharge_45[run,:], Timeseries_Future_45)
            for month in 1:12
                # computes the average mean monthly discharge over all years for each month
                current_Month_Discharge = Monthly_Discharge_past[findall(x->x == month, Month)]
                current_Month_Discharge_future = Monthly_Discharge_future[findall(x->x == month, Month_future)]
                current_Month_Discharge = mean(current_Month_Discharge)
                current_Month_Discharge_future = mean(current_Month_Discharge_future)
                error = relative_error(current_Month_Discharge_future, current_Month_Discharge)
                append!(average_monthly_Discharge_past, current_Month_Discharge)
                append!(average_monthly_Discharge_future, current_Month_Discharge_future)
                append!(all_months_all_runs, month)
                append!(error_average_monthly_Discharge_all_runs, error)
            end
        end
    end
    return error_average_monthly_Discharge_all_runs, average_monthly_Discharge_past, average_monthly_Discharge_future, all_months_all_runs
end

"""
Computes the monthly snow storage of the past and future and the relative changes, of the projections of the path and the different parameter sets

$(SIGNATURES)

The function returns relative change in monthly sum of snow melt, the monthly sum of snow melt of the past, and the monthly mean snow melt of the future
"""
function change_monthly_snowstorage(path_to_projections, Catchment_Name)
    Name_Projections_45 = readdir(path_to_projections)
    Timeseries_Past = collect(Date(1981,1,1):Day(1):Date(2010,12,31))
    Timeseries_End = readdlm("/home/sarah/Master/Thesis/Data/Projektionen/End_Timeseries_45_85.txt",',')
    if path_45[end-2:end-1] == "45"
        index = 1
        rcp = "45"
        print(rcp, " ", path_to_projections)
    elseif path_45[end-2:end-1] == "85"
        index = 2
        rcp="85"
        print(rcp, " ", path_to_projections)
    end
    average_monthly_Discharge_past = Float64[]
    average_monthly_Discharge_future = Float64[]
    error_average_monthly_Discharge_all_runs = Float64[]
    all_months_all_runs = Float64[]
    for (i, name) in enumerate(Name_Projections_45)
        Timeseries_Future_45 = collect(Date(Timeseries_End[i,index]-29,1,1):Day(1):Date(Timeseries_End[i,index],12,31))
        Past_Discharge_45 = readdlm(path_to_projections*name*"/"*Catchment_Name*"/300_model_results_snow_storage_past_2010_distribution_600_2700_new.csv", ',')
        Future_Discharge_45 = readdlm(path_to_projections*name*"/"*Catchment_Name*"/300_model_results_snow_storage_future_2100_distribution_600_2700_new.csv", ',')
        println("future ",size(Future_Discharge_45))
        println("past ",size(Past_Discharge_45))
        for run in 1:size(Past_Discharge_45)[1]
            # computes sum of monthly snow melt of each month
            Monthly_Discharge_past, Month = monthly_discharge(Past_Discharge_45[run,:], Timeseries_Past)
            Monthly_Discharge_future, Month_future = monthly_discharge(Future_Discharge_45[run,:], Timeseries_Future_45)
            for month in 1:12
                # computes the average mean monthly discharge over all years for each month
                current_Month_Discharge = Monthly_Discharge_past[findall(x->x == month, Month)]
                current_Month_Discharge_future = Monthly_Discharge_future[findall(x->x == month, Month_future)]
                current_Month_Discharge = mean(current_Month_Discharge)
                current_Month_Discharge_future = mean(current_Month_Discharge_future)
                error = relative_error(current_Month_Discharge_future, current_Month_Discharge)
                append!(average_monthly_Discharge_past, current_Month_Discharge)
                append!(average_monthly_Discharge_future, current_Month_Discharge_future)
                append!(all_months_all_runs, month)
                append!(error_average_monthly_Discharge_all_runs, error)
            end
        end
    end
    return error_average_monthly_Discharge_all_runs, average_monthly_Discharge_past, average_monthly_Discharge_future, all_months_all_runs
end

"""
Computes the mean yearly amount of snow melt as contribution to hydrological system for past and future of all 14 simulations.

$(SIGNATURES)

The function returns the mean yearly amount of snow melt [mm] in past and future of all 14 simulations with 100 runs. It takes as input the path to projection and the catchment name.
"""
function snow_contribution_yearly_projections(path_to_projections, Catchment_Name)
    Name_Projections = readdir(path_to_projections)
    Timeseries_Past = collect(Date(1981,1,1):Day(1):Date(2010,12,31))
    Timeseries_End = readdlm("/home/sarah/Master/Thesis/Data/Projektionen/End_Timeseries_45_85.txt",',')

    Mean_Annual_Snow_Storage_Past_All_Proj = Float64[]
    Mean_Annual_Snow_Storage_Future_All_Proj = Float64[]

    if path_to_projections[end-2:end-1] == "45"
        index = 1
        rcp = "45"
        print(rcp, " ", path_to_projections)
    elseif path_to_projections[end-2:end-1] == "85"
        index = 2
        rcp="85"
        print(rcp, " ", path_to_projections)
    end

    for (i, name) in enumerate(Name_Projections)
        Timeseries_Future = collect(Date(Timeseries_End[i,index]-29,1,1):Day(1):Date(Timeseries_End[i,index],12,31))
        Snow_Storage_Past = readdlm(path_to_projections*name*"/"*Catchment_Name*"/300_model_results_snow_melt_snow_redistr_past_2010.csv", ',')
        Snow_Storage_Future = readdlm(path_to_projections*name*"/"*Catchment_Name*"/300_model_results_snow_melt_snow_redistr_future_2100.csv", ',')
        for run in 1:size(Snow_Storage_Past)[1]
            Mean_Annual_Snow_Storage_Past = snow_contribution_yearly(Snow_Storage_Past[run,:], Timeseries_Past)
            Mean_Annual_Snow_Storage_Future = snow_contribution_yearly(Snow_Storage_Future[run,:], Timeseries_Future)
            append!(Mean_Annual_Snow_Storage_Past_All_Proj, Mean_Annual_Snow_Storage_Past)
            append!(Mean_Annual_Snow_Storage_Future_All_Proj, Mean_Annual_Snow_Storage_Future)
        end
    end
    return Mean_Annual_Snow_Storage_Past_All_Proj::Array{Float64,1}, Mean_Annual_Snow_Storage_Future_All_Proj::Array{Float64,1}
end

#

"""
Plots the mean yearly amount of snow melt as contribution to hydrological system for past and future.

$(SIGNATURES)

It plots the relative and absolute change in snow melt and a figure comparing past to future snow melt.
"""
function plot_snow_storage_contribution(Snow_Storage_Past45, Snow_Storage_Future45, Snow_Storage_Past85, Snow_Storage_Future85, Catchment_Name)
    Farben45=palette(:blues)
    Farben85=palette(:reds)
    rcps = ["RCP 4.5", "RCP 8.5"]
    plot()
    boxplot!([rcps[1]], relative_error(Snow_Storage_Future45, Snow_Storage_Past45)*100, size=(2000,800), leg=false, color="blue")
    boxplot!([rcps[2]], relative_error(Snow_Storage_Future85, Snow_Storage_Past85)*100,size=(2000,800), left_margin = [5mm 0mm], leg=false, color="red")
    violin!([rcps[1]], relative_error(Snow_Storage_Future45, Snow_Storage_Past45)*100, size=(2000,800), leg=false, color="blue", alpha=0.6)
    violin!([rcps[2]], relative_error(Snow_Storage_Future85, Snow_Storage_Past85)*100,size=(2000,800), left_margin = [5mm 0mm], leg=false, color="red", alpha=0.6)
    ylabel!("Relative Change [%]")
    title!("Relative Change in Annual (Oct-Sep) Snow Melt in "*Catchment_Name)
    relative_change = boxplot!()
    #savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/rel_change_snow_storage.png")

    plot()
    boxplot!([rcps[1]], Snow_Storage_Future45 - Snow_Storage_Past45, size=(2000,800), leg=false, color="blue")
    boxplot!([rcps[2]], Snow_Storage_Future85 - Snow_Storage_Past85,size=(2000,800), left_margin = [5mm 0mm], leg=false, color="red")
    violin!([rcps[1]], Snow_Storage_Future45 - Snow_Storage_Past45, size=(2000,800), leg=false, color="blue", alpha=0.6)
    violin!([rcps[2]], Snow_Storage_Future85 - Snow_Storage_Past85,size=(2000,800), left_margin = [5mm 0mm], leg=false, color="red", alpha=0.6)
    ylabel!("Absolute Change [mm]")
    title!("Absolute Change in Annual (Oct-Sep) Snow Melt in "*Catchment_Name)
    absolute_change = boxplot!()
    plot(relative_change, absolute_change, size=(2000,800), left_margin = [5mm 0mm])
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/"*Catchment_Name*"_change_snow_melt_new.png")

    # plot real values
    plot()
    boxplot!([rcps[1]*" Past"], Snow_Storage_Past45, size=(2000,800), leg=false, color=[Farben45[1]])
    boxplot!([rcps[1]*" Future"], Snow_Storage_Future45, size=(2000,800), left_margin = [5mm 0mm], leg=false, color=[Farben45[2]])
    violin!([rcps[1]*" Past"], Snow_Storage_Past45, size=(2000,800), leg=false, color=[Farben45[1]], alpha=0.6)
    violin!([rcps[1]*" Future"], Snow_Storage_Future45,size=(2000,800), left_margin = [5mm 0mm], leg=false, color=[Farben45[2]], alpha=0.6)
    ylabel!("Snow Melt Contribution [mm]")
    title!("Annual Snow Melt (Oct-Sep) RCP 4.5 in "*Catchment_Name)
    relative_change = boxplot!()
    #savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/rel_change_snow_storage.png")

    plot()
    boxplot!([rcps[2]*" Past"], Snow_Storage_Past85, size=(2000,800), leg=false, color=[Farben85[1]])
    boxplot!([rcps[2]*" Future"], Snow_Storage_Future85 ,size=(2000,800), left_margin = [5mm 0mm], leg=false, color=[Farben85[2]])
    violin!([rcps[2]*" Past"], Snow_Storage_Past85, size=(2000,800), leg=false, color=[Farben85[1]], alpha=0.6)
    violin!([rcps[2]*" Future"], Snow_Storage_Future85,size=(2000,800), left_margin = [5mm 0mm], leg=false, color=[Farben85[2]], alpha=0.6)
    ylabel!("Snow Melt Contribution [mm]")
    title!(" Annual (Oct-Sep) Snow Melt RCP 8.5 in "*Catchment_Name)
    absolute_change = boxplot!()
    plot(relative_change, absolute_change, size=(2000,800), left_margin = [5mm 0mm])
    savefig("/home/sarah/Master/Thesis/Results/Projektionen/"*Catchment_Name*"/PastvsFuture/Snow_Cover/"*Catchment_Name*"_snow_melt_past_future_new.png")
end

max_stars_repo_name: sarah-hanus/hbv-mountain

max_issues_repo_name: sarah-hanus/hbv-mountain

max_forks_repo_name: sarah-hanus/hbv-mountain


text:
import FuzzyCompletions: complete_path, completion_text, score
import Distributed
import .PkgCompat: package_completions
using Markdown
import REPL

###
# RESPONSES FOR AUTOCOMPLETE & DOCS
###

function format_path_completion(completion)
    replace(replace(completion_text(completion), "\\ " => " "), "\\\\" => "\\")
end

responses[:completepath] = function response_completepath(🙋::ClientRequest)
    path = 🙋.body["query"]
    pos = lastindex(path)

    results, loc, found = complete_path(path, pos)
    isenough(x) = x ≥ -0.1
    filter!(isenough ∘ score, results) # too many candiates otherwise. -0.1 instead of 0 to enable autocompletions for paths: `/` or `/asdf/`

    start_utf8 = let
        # REPLCompletions takes into account that spaces need to be prefixed with `\` in the shell, so it subtracts the number of spaces in the filename from `start`:
        # https://github.com/JuliaLang/julia/blob/c54f80c785a3107ae411267427bbca05f5362b0b/stdlib/REPL/src/REPLCompletions.jl#L270

        # we don't use prefixes, so we need to reverse this.

        # this is from the Julia source code:
        # https://github.com/JuliaLang/julia/blob/c54f80c785a3107ae411267427bbca05f5362b0b/stdlib/REPL/src/REPLCompletions.jl#L195-L204
        if Base.Sys.isunix() && occursin(r"^~(?:/|$)", path)
            # if the path is just "~", don't consider the expanded username as a prefix
            if path == "~"
                dir, prefix = homedir(), ""
            else
                dir, prefix = splitdir(homedir() * path[2:end])
            end
        else
            dir, prefix = splitdir(path)
        end

        loc.start + count(isequal(' '), prefix)
    end
    stop_utf8 = nextind(path, pos) # advance one unicode char, js uses exclusive upper bound

    scores = [max(0.0, score(r)) for r in results]
    formatted = format_path_completion.(results)

    # sort on score. If a tie (e.g. both score 0.0), sort on dir/file. If a tie, sort alphabetically.
    perm = sortperm(collect(zip(.-scores, (!isdirpath).(formatted), formatted)))

    msg = UpdateMessage(:completion_result,
        Dict(
            :start => start_utf8 - 1, # 1-based index (julia) to 0-based index (js)
            :stop => stop_utf8 - 1, # idem
            :results => formatted[perm]
            ), 🙋.notebook, nothing, 🙋.initiator)

    putclientupdates!(🙋.session, 🙋.initiator, msg)
end

function package_name_to_complete(str)
	matches = match(r"(import|using) ([a-zA-Z0-9]+)$", str)
	matches === nothing ? nothing : matches[2]
end

responses[:complete] = function response_complete(🙋::ClientRequest)
    try require_notebook(🙋) catch; return; end
    query = 🙋.body["query"]
    pos = lastindex(query) # the query is cut at the cursor position by the front-end, so the cursor position is just the last legal index

    workspace = WorkspaceManager.get_workspace((🙋.session, 🙋.notebook))

    results, loc, found = if package_name_to_complete(query) !== nothing
        p = package_name_to_complete(query)
        cs = package_completions(p) |> sort
        [(c,"package",true) for c in cs], (nextind(query, pos-length(p)):pos), true
    else
        if will_run_code(🙋.notebook) && isready(workspace.dowork_token)
            # we don't use eval_format_fetch_in_workspace because we don't want the output to be string-formatted.
            # This works in this particular case, because the return object, a `Completion`, exists in this scope too.
            Distributed.remotecall_eval(Main, workspace.pid, :(PlutoRunner.completion_fetcher(
                $query, $pos,
                getfield(Main, $(QuoteNode(workspace.module_name))),
                )))
        else
            # We can at least autocomplete general julia things:
            PlutoRunner.completion_fetcher(query, pos, Main)
        end
    end

    start_utf8 = loc.start
    stop_utf8 = nextind(query, pos) # advance one unicode char, js uses exclusive upper bound

    msg = UpdateMessage(:completion_result,
        Dict(
            :start => start_utf8 - 1, # 1-based index (julia) to 0-based index (js)
            :stop => stop_utf8 - 1, # idem
            :results => results
            ), 🙋.notebook, nothing, 🙋.initiator)

    putclientupdates!(🙋.session, 🙋.initiator, msg)
end

responses[:docs] = function response_docs(🙋::ClientRequest)
    require_notebook(🙋)
    query = 🙋.body["query"]

    # Expand string macro calls to their macro form:
    # `html"` should yield `@html_str` and
    # `Markdown.md"` should yield `@Markdown.md_str`. (Ideally `Markdown.@md_str` but the former is easier)
    if endswith(query, "\"") && query != "\""
        query = "@$(query[begin:end-1])_str"
    end

    doc_html, status = if REPL.lookup_doc(Symbol(query)) isa Markdown.MD
        # available in Base, no need to ask worker
        doc_md = REPL.lookup_doc(Symbol(query))
        (repr(MIME("text/html"), doc_md), :👍)
    else
        workspace = WorkspaceManager.get_workspace((🙋.session, 🙋.notebook))

        if will_run_code(🙋.notebook) && isready(workspace.dowork_token)
            Distributed.remotecall_eval(Main, workspace.pid, :(PlutoRunner.doc_fetcher(
                $query,
                getfield(Main, $(QuoteNode(workspace.module_name))),
            )))
        else
            (nothing, :⌛)
        end
    end

    msg = UpdateMessage(:doc_result,
        Dict(
            :status => status,
            :doc => doc_html,
            ), 🙋.notebook, nothing, 🙋.initiator)

    putclientupdates!(🙋.session, 🙋.initiator, msg)
end


text:
using GaussianProcesses
using HDF5

# Load the possible parameter combinations

# Points at which to create images
parameters = readdlm("parameters.dat")

drop_out = parameters[end, :]
parameters = parameters[1:end-1, :]

# Column headers
# nu, mass, r_c, T_10, q, logM_gas, ksi, incl, PA = row

d = 2 # dimensions
n = size(parameters, 1) # number of observations

# Load the sampled visibilities from HDF5
fid = h5open("grid_VV_old.hdf5", "r")
reals = read(fid, "real")
imags = read(fid, "imag")
VV = reals + imags .* im # Complex visibility
close(fid)

# Choose our visibility
vis = squeeze(VV[1, :], 1)

# Subtract the mean value
vis_norm = vis .- mean(vis)

vis_drop_out = vis_norm[end]
vis_norm = vis_norm[1:end-1]

# # Plot a histogram of these
# import PyPlot.plt
# using LaTeXStrings
#
# fig = plt[:figure]()
# ax = fig[:gca](projection="3d")
#
# ax[:scatter](parameters[:,1], parameters[:,2], real(vis_norm), color="k")
# ax[:scatter](drop_out[1], drop_out[2], real(vis_drop_out), color="r")
#
# plt[:show]()
# fig[:savefig]("3d_plot.png")
#
# quit()

#Choose mean and covariance function
mZero = MeanZero()                             # Zero mean function

# Log of length scale
# ll = Float64[-3, -1, -0.5, 1, -1, 1, 2, 0, -1]
ll = Float64[-4.0, -0.5]
logSigma = 0.1

kern = SE(ll, logSigma)

# log standard deviation of observation noise (this is optional)
# logObsNoise = -5.
gp = GP(parameters', real(vis_norm), mZero, kern) #, logObsNoise)   # Fit the GP

optimize!(gp; ftol=1e-20)

# nus = linspace(345.7956, 345.7975, 5)
# masses = linspace(0.8, 1.3, 5)
# r_cs = linspace(300, 600, 5)
# T_10s = linspace(60, 130, 5)
# qs = linspace(0.3, 0.7, 5)
# logM_gass = linspace(-3, -1.5, 5)
# ksis = linspace(0.2, 0.5, 5)
# incls = linspace(35, 55, 5)
# PAs = linspace(145, 160, 5)

println()
println("GP mean: ", gp.m)
println("GP kernel: ", gp.k)

x_predict = drop_out

mu, sigma2 = predict(gp, x_predict')
sigma = sqrt(sigma2)

println("mu: $mu, sigma: $sigma")
println("dropout: ", real(vis_drop_out))

max_stars_repo_name: iancze/DiskEmulator

max_issues_repo_name: iancze/DiskEmulator

max_forks_repo_name: iancze/DiskEmulator


text:
using Test

@testset "./src/LocalPolytope/adjacency_decomposition.jl" begin

using BellScenario

@testset "rotate_facet()" begin
    @testset "simplex rotation" begin
        F = [-1,0,0,0]
        G1 = [0,0,-1,0]
        G2 = [0,0,1,1]

        xbar = [1,0,0]

        # neighboring facets
        @test LocalPolytope.rotate_facet(F,G1,xbar) == [0,0,-1,0]
        @test LocalPolytope.rotate_facet(F,G2,xbar) == [1,0,1,1]
    end
end

end


text:
import PredictMDAPI

PredictMDAPI.get_history(::AbstractFittable; kwargs...)::Nothing = nothing

PredictMDAPI.get_underlying(::AbstractFittable; kwargs...)::Nothing = nothing

PredictMDAPI.parse_functions!(::AbstractFittable)::Nothing = nothing

PredictMDAPI.set_feature_contrasts!(::AbstractFittable,
                                    ::AbstractFeatureContrasts)::Nothing = nothing

PredictMDAPI.set_max_epochs!(::AbstractFittable,::Integer)::Nothing = nothing


text:
# Get pattern

function _barcode_pattern_code128(code::Vector{<:AbstractString})

    m = match(r"^START (A|B|C)$", first(code))
    m !== nothing || throw(
        ArgumentError(
            "First element of `code` should be either `START A`, `START B` or `START C`",
        ),
    )

    last(code) == "STOP" ||
        throw(ArgumentError("Last element of `code` should be `STOP`"))

    count(code .== "STOP") == 1 ||
        throw(ArgumentError("""There should be only one "STOP" in the encoding sequence"""))

    count(match.(r"^START [A|B|C]$", code) .!== nothing) == 1 ||
        throw(
            ArgumentError(
                "There should be only one r\"START [A|B|C]\" in the encoding sequence"
            )
        )

    subtype = Symbol("code128$(lowercase(m.captures[1]))")
    nextsubtype = subtype

    # initialize pattern with the quiet zone, which is at least 10x, where x is
    # the width of each module, assumed here to be one bit. We use 11x just to have the
    # same length as the other symbols (except the double bar termination symbol).
    quiet_zone = "0"^11
    pattern = [quiet_zone]

    nrow = findfirst(==(first(code)), CODE128[:, subtype])
    push!(pattern, CODE128.pattern[nrow])
    chk_sum = CODE128.value[nrow]
    multiplier = 0

    for c in code[2:end]

        multiplier += 1

        if c == "CHECKSUM"
            push!(pattern, CODE128.pattern[rem(chk_sum, 103)+1])
            chk_sum += multiplier * chk_sum
        else
            nrow = findfirst(==(c), CODE128[:, nextsubtype])
            nrow === nothing && throw(
                ArgumentError(
                    "$c is not part of subtype $(titlecase(string(nextsubtype))) of Code128",
                ),
            )
            push!(pattern, CODE128.pattern[nrow])
            chk_sum += multiplier * CODE128.value[nrow]

            nextsubtype = subtype
            m = match(r"^CODE (A|B|C)$", c)
            if m !== nothing
                subtype = nextsubtype = Symbol("code128$(lowercase(m.captures[1]))")
            end

            m = match(r"^SHIFT (A|B)$", c)
            if m !== nothing
                nextsubtype = Symbol("code128$(lowercase(m.captures[1]))")
            end

        end
    end
    # "END" bars
    push!(pattern, "11")

    # Quiet zone
    push!(pattern, quiet_zone)

    pattern = prod(pattern)
    return pattern
end

"""
    barcode_pattern(code::Vector{<:AbstractString}, encoding_type::Symbol)

Return the binary pattern for a given vector of codes following the specifications
determined by the `encoding_type`.

Currently, only Code128 specification is available.

If `encoding_type` is either `:code128a`, `:code128b`, or `:code128c`, it returns the
encoding following the corresponding subtype. If `encoding_type` is `:code128`, it returns
an optimized encoding, possibily mixing different subtypes. This strategy follows
the specifications in "GS1 General Specifications, Version 13, Issue 1, Jan-2013,
Section 5.4.7.7. Use of Start, Code Set, and Shift symbols to Minimize Symbol Length
(Informative), pages 268 to 269."

It is assumed that `code` starts with either `START A`, `START B`, or `START C` and that
it ends with `STOP`.

The checksum must either be already computed or one can add an element `CHECKSUM`
for the check sum to be computed and included where this directive appears.

# Examples

```jldoctest
julia> pattern = barcode_pattern(
           ["START C", "00", "01", "32", "CHECKSUM", "STOP"],
           :code128
       )
"000000000001101001110011011001100110011011001100011011010111101110110001110101100000000000"

julia> pattern = barcode_pattern(
           ["START A", "A", "B", "C", "CHECKSUM", "STOP"],
           :code128
       )
"000000000001101000010010100011000100010110001000100011011011001100110001110101100000000000"
```
"""
function barcode_pattern(code::Vector{<:AbstractString}, encoding_type::Symbol)
    if encoding_type in (:code128, :code128a, :code128b, :code128c)
        pattern = _barcode_pattern_code128(code)
    else
        throw(
            ArgumentError(
                "Pattern encoding for encoding type `$encoding_type` not implemented"
            )
        )
    end
    return pattern
end

"""
    barcode_pattern(msg::AbstractString, encoding_type::Symbol)

Generate the barcode pattern for the given `msg`, according to the encoding Code128.

If `encoding_type` is set to either `:code128a`, `:code128b`, or `:code128c`, it returns
the encoding following the corresponding subtype. It throws an `ArgumentError` if the given
`msg` is not suitable for encoding with the given subtype.

If `encoding_type` is set to `:code128`, it attempts to use either of these subtypes or
a combination of them. It throws an `ArgumentError` if not possible to do so.

The encoding is returned as a vector of string patterns, with each element corresponding
to the encoding of each symbol in `code`, and appended with the proper quiet zones and the
`START`, `STOP`, checksum, and ending bar patterns.
"""
function barcode_pattern(msg::AbstractString, encoding_type::Symbol)
    code = encode(msg, encoding_type)
    return barcode_pattern(code, encoding_type)
end

function _split_pattern_code128(pattern::AbstractString)
    all(d -> d in ('0', '1'), pattern) || throw(
        ArgumentError(
            "A barcode pattern should be made of 0's and 1's only"
        )
    )

    firstbar = findfirst(==('1'), pattern)
    lastbar = findlast(==('1'), pattern)
    lastbar - firstbar ≥ 34 || throw(
        ArgumentoError(
            "Barcode pattern is too short to be a Code128 barcode pattern"
        )
    )
    pattern[lastbar-1:lastbar] == "11" || throw(
        ArgumentError(
            "Barcode pattern does not end properly with a double bar"
        )
    )
    rem(lastbar - firstbar - 1, 11) == 0 || throw(
        ArgumentError(
            "Barcode pattern contents should have a length multiple of 11"
        )
    )

    core_patterns = [pattern[i:i+10] for i in firstbar:11:lastbar-2]

    first(core_patterns) in CODE128[104:106, :pattern] || throw(
        ArgumentError(
            "Barcode pattern should start with a pattern for either `START A`, " *
            "`START B`, or `START C`"
        )
    )

    last(core_patterns) == CODE128[107, :pattern] || throw(
        ArgumentError(
            "Barcode pattern should end with the pattern for `STOP`"
        )
    )

    count(core_patterns .== CODE128[107, :pattern]) == 1 ||
        throw(ArgumentError("There should be only one pattern for \"STOP\""))

    count([p in CODE128[104:106, :pattern] for p in core_patterns]) == 1 ||
        throw(ArgumentError("There should be only one pattern for r\"START [A|B|C]\""))

    splited_pattern = [
        pattern[begin:firstbar-1]
        core_patterns
        "11"
        pattern[lastbar+1:end]
    ]

    return splited_pattern
end

function _barcode_depattern_code128(pattern::AbstractString)

    splited_pattern = _split_pattern_code128(pattern)

    # initialize code with the subtype START code
    code = CODE128.code128a[CODE128[:, :pattern] .== splited_pattern[2]]
    nextsubtype = subtype = Symbol("code128$(lowercase(first(code)[end]))")

    nrow = findfirst(==(first(code)), CODE128[:, subtype])
    chk_sum = CODE128.value[nrow]
    multiplier = 0

    # skip quiet zone and START [A|B|C]
    # stop before CHECKSUM, STOP, END, and last quiet zone
    for p in splited_pattern[3:end-4]

        multiplier += 1

        nrow = findfirst(==(p), CODE128[:, :pattern])
        nrow === nothing && throw(
            ArgumentError(
                "$p is not a valid CODE128 pattern",
            ),
        )
        push!(code, CODE128[nrow, nextsubtype])
        chk_sum += multiplier * CODE128.value[nrow]

        nextsubtype = subtype

        p == CODE128.pattern[99] && (
            (subtype == :code128a && (nextsubtype = :code128b)) ||
            (subtype == :code128b && (nextsubtype = :code128a))
        )
        p == CODE128.pattern[100] && subtype in (:code128a, :code128b) &&
            (subtype = nextsubtype = :code128c)
        p == CODE128.pattern[101] && subtype in (:code128a, :code128c) &&
            (subtype = nextsubtype = :code128b)
        p == CODE128.pattern[102] && subtype in (:code128b, :code128c) &&
            (subtype = nextsubtype = :code128a)
    end

    chk_sum = chk_sum % 103
    CODE128.value[CODE128.pattern .== splited_pattern[end-3]][1] == chk_sum ||
        @warn "Barcode checksum does not match computed checksum"
    push!(code, "CHECKSUM $chk_sum")
    push!(code, "STOP")

    return code
end

function barcode_depattern(pattern::AbstractString, encoding_type::Symbol)
    encoding_type == :code128 || throw(
        ArgumentError(
            "Decoding type `$(Meta.quot(mode))` not implemented"
        )
    )
    encoding_type == :code128 && (code = _barcode_depattern_code128(pattern))

    return code
end

function barcode_decode(pattern::AbstractString, encoding_type::Symbol)
    encoding_type == :code128 || throw(
        ArgumentError(
            "Decoding type `$(Meta.quot(mode))` not implemented"
        )
    )
    encoding_type == :code128 && (code = _barcode_depattern_code128(pattern))

    msg = decode(code, encoding_type)

    return msg
end


text:
using BenchmarkTools
using BenchmarkCIExample

SUITE = BenchmarkGroup()

setup_sum() = randn(10_000)

let suite = SUITE["sum(xs)"] = BenchmarkGroup()
    suite["base"] = @benchmarkable sum(xs) setup = (xs = setup_sum())
    suite["naive_sum"] = @benchmarkable naive_sum(xs) setup = (xs = setup_sum())
end


text:
function quaternion_mult(	p0::Float64,
						p1::Float64,
						p2::Float64,
						p3::Float64,
						q0::Float64,
						q1::Float64,
						q2::Float64,
						q3::Float64)

	r0::Float64 = p0 * q0 - p1 * q1 - p2 * q2 - p3 * q3
	r1::Float64 = p0 * q1 + p1 * q0 - p2 * q3 + p3 * q2
	r2::Float64 = p0 * q2 + p2 * q0 + p1 * q3 - p3 * q1
	r3::Float64 = p0 * q3 - p1 * q2 + p2 * q1 + p3 * q0

	return (r0, r1, r2, r3)
end