text:
"""
This function finds the hourly availability data for projects
to be passed to expected and actual market clearing modules.
"""
function find_project_availability_data(project::P,
                                     availability_df::DataFrames.DataFrame,
                                     num_invperiods::Int64,
                                     num_hours::Int64) where P <: Project{<:BuildPhase}
    type = get_type(get_tech(project))
    zone = get_zone(get_tech(project))
        if in(get_name(project), names(availability_df))
            availability_raw = availability_df[:, Symbol(get_name(project))]
        elseif in("$(type)_$(zone)", names(availability_df))
            availability_raw = availability_df[:, Symbol("$(type)_$(zone)")]
        else
            error("project availiability data not found!")
        end

        data_years = unique(availability_df.Year)

        data_hours = count(i -> i == data_years[1], availability_df.Year)

        @assert data_hours == num_hours

        availability_input = zeros(num_invperiods, data_hours)

        # Fill raw availability data into yearly and hourly values
        for y in 1:num_invperiods
            availability_input[y, 1:data_hours] = availability_raw[1:data_hours]
        end

        return availability_input
end

"""
This function extracts the technical details for generators
to be passed to expected and actual market clearing modules.
"""
function get_technical_details(project::P) where P <: GeneratorEMIS{<: BuildPhase}

    project_type = "generator"
    min_input = 0.0
    max_input = 0.0
    efficiency_in = 0.0
    efficiency_out = 0.0
    min_storage = 0.0
    max_storage = 0.0
    init_storage = 0.0

    return project_type, min_input, max_input, efficiency_in, efficiency_out, min_storage, max_storage, init_storage
end

"""
This function extracts the technical details for storage units
to be passed to expected and actual market clearing modules.
"""
function get_technical_details(project::P) where P <: StorageEMIS{<: BuildPhase}

    tech = get_tech(project)

    project_type = "storage"

    input_power_limits = get_input_active_power_limits(tech)
    min_input = input_power_limits[:min]
    max_input = input_power_limits[:max]

    efficiency = get_efficiency(tech)
    efficiency_in = efficiency[:in]
    efficiency_out = efficiency[:out]

    storage_capacity = get_storage_capacity(tech)
    min_storage = storage_capacity[:min]
    max_storage = storage_capacity[:max]

    init_storage = get_soc(tech)

    return project_type, min_input, max_input, efficiency_in, efficiency_out, min_storage, max_storage, init_storage
end

function get_marginal_cost_energy(product::T) where T <: Product
    return
end

function get_marginal_cost_energy(product::Energy)
    marginal_cost_energy = get_marginal_cost(product)
    return marginal_cost_energy
end

"""
This function returns the project's marginal cost of energy product to be passed to economic dispatch and CEM.
"""
function get_project_energy_cost(project::P) where P <: Project{<: BuildPhase}
    marginal_cost_energy = 0.0
    for product in find_operating_products(get_products(project))
        marginal_cost_temp = get_marginal_cost_energy(product)
        if !isnothing(marginal_cost_temp)
            marginal_cost_energy = marginal_cost_temp
        end
    end
    return marginal_cost_energy
end

function get_marginal_cost_reserveup(product::T) where T <: Product
    return
end

function get_marginal_cost_reserveup(product::OperatingReserve{ReserveUpEMIS})
    marginal_cost_reserveup = get_marginal_cost(product)
    return marginal_cost_reserveup
end

"""
This function returns the project's marginal cost of reserve up product to be passed to economic dispatch and CEM.
"""
function get_project_reserve_up_cost(project::P, voll::Float64) where P <: Project{<: BuildPhase}
    marginal_cost_reserve_up = voll
    for product in find_operating_products(get_products(project))
        marginal_cost_temp = get_marginal_cost_reserveup(product)
        if !isnothing(marginal_cost_temp)
            marginal_cost_reserve_up = marginal_cost_temp
        end
    end
    return marginal_cost_reserve_up
end

function get_max_reserveup(product::T) where T <: Product
    return
end

function get_max_reserveup(product::OperatingReserve{ReserveUpEMIS})
    max_reserveup = get_max_limit(product)
    return max_reserveup
end

"""
This function returns the project's maximum reserve up participation limit to be passed to economic dispatch and CEM.
"""
function get_project_max_reserveup(project::P) where P <: Project{<: BuildPhase}
    max_reserveup = 0.
    for product in find_operating_products(get_products(project))
        max_reserveup_temp = get_max_reserveup(product)
        if !isnothing(max_reserveup_temp)
            max_reserveup = max_reserveup_temp * get_maxcap(project)
        end
    end
    return max_reserveup
end

function get_marginal_cost_reservedown(product::T) where T <: Product
    return
end

function get_marginal_cost_reservedown(product::OperatingReserve{ReserveDownEMIS})
    marginal_cost_reserve_down = get_marginal_cost(product)
    return marginal_cost_reserve_down
end

"""
This function returns the project's marginal cost of reserve down product to be passed to economic dispatch and CEM.
"""
function get_project_reserve_down_cost(project::P, voll::Float64) where P <: Project{<: BuildPhase}
    marginal_cost_reserve_down = voll
    for product in find_operating_products(get_products(project))
        marginal_cost_temp = get_marginal_cost_reservedown(product)
        if !isnothing(marginal_cost_temp)
            marginal_cost_reserve_down = marginal_cost_temp
        end
    end
    return marginal_cost_reserve_down
end

function get_max_reservedown(product::T) where T <: Product
    return
end

function get_max_reservedown(product::OperatingReserve{ReserveDownEMIS})
    max_reservedown = get_max_limit(product)
    return max_reservedown
end

"""
This function returns the project's maximum reserve down participation limit to be passed to economic dispatch and CEM.
"""
function get_project_max_reservedown(project::P) where P <: Project{<: BuildPhase}
    max_reservedown = 0.
    for product in find_operating_products(get_products(project))
        max_reservedown_temp = get_max_reservedown(product)
        if !isnothing(max_reservedown_temp)
            max_reservedown = max_reservedown_temp * get_maxcap(project)
        end
    end
    return max_reservedown
end

"""
This function returns the project's derating factor to be passed to CEM and capacity market clearing module.
Returns 0 if there is no capacity market participation.
"""
function get_project_derating(project::P) where P <: Project{<: BuildPhase}
    derating_factor = 0.
    for product in get_products(project)
        derating_temp = get_derating(product)
        if !isnothing(derating_temp)
            derating_factor = derating_temp
        end
    end
    return derating_factor
end

"""
This function returns the project's capacity market bid to be passed to the capacity market clearing module.
Returns 0 if there is no capacity market participation.
"""
function get_project_capacity_market_bid(project::P) where P <: Project{<: BuildPhase}
    capacity_bid = 0.
    for product in get_products(project)
        capacity_bid_temp = get_capacity_bid(product)
        if !isnothing(capacity_bid_temp)
            capacity_bid = capacity_bid_temp
        end
    end
    return capacity_bid
end

"""
This function returns the project's REC energy output to be passed the REC market market clearing module.
Returns 0 if there is no REC market participation.
"""
function get_project_rec_output(project::P) where P <: Project{<: BuildPhase}
    rec_output = 0.
    for product in get_products(project)
        rec_output_temp = get_rec_certificates(product)
        if !isnothing(rec_output_temp)
            rec_output = rec_output_temp
        end
    end
    return rec_output
end

"""
This function returns the project's REC market bid to be passed the REC market market clearing module.
Returns 0 if there is no REC market participation.
"""
function get_project_rec_market_bid(project::P) where P <: Project{<: BuildPhase}
    rec_bid = 0.
    for product in get_products(project)
        rec_bid_temp = get_rec_bid(product)
        if !isnothing(rec_bid_temp)
            rec_bid = rec_bid_temp
        end
    end
    return rec_bid
end

"""
This function creates the MarketProject struct to be passed to CEM price projection and endogeneous Economic Dispatch models.
"""
function populate_market_project(project::P,
                                reserve_up_cost::Float64,
                                reserve_down_cost::Float64,
                                project_type::String,
                                min_input::Float64,
                                max_input::Float64,
                                efficiency_in::Float64,
                                efficiency_out::Float64,
                                min_storage::Float64,
                                max_storage::Float64,
                                init_storage::Float64,
                                availability_input::Array{Float64, 2},
                                existing_units::Int64,
                                units_inqueue::Vector{Float64},
                                remaining_lag_time::Int64,
                                remaining_life_time::Int64,
                                iteration_year::Int64,
                                num_invperiods::Int64) where P <: Project{<: BuildPhase}

    finance_data = get_finance_data(project)

    market_project = MarketProject(
        get_name(project),                                                                       # name
        project_type,                                                                            # is project of storage type
        get_type(get_tech(project)),                                                             # technology type
        get_fixed_OM_cost(finance_data),                                                           # annualfixed O&M costs
        get_queue_cost(finance_data),                                                             # Queue cost
        get_project_energy_cost(project),                                                        # marginal cost of energy
        reserve_up_cost,                                                                         # marginal cost of reserve down
        reserve_down_cost,                                                                       # marginal cost of reserve up
        get_investment_cost(finance_data)[iteration_year:iteration_year + num_invperiods - 1],    # yearly investment cost
        get_discount_rate(finance_data),                                                          # discount rate
        get_mincap(project),                                                                     # minimum capacity
        get_maxcap(project),                                                                     # maximum capacity
        min_input,                                                                               # minimum input power
        max_input,                                                                               # maximum input power
        efficiency_in,                                                                           # input efficiency
        efficiency_out,                                                                          # output efficiency
        min_storage,                                                                             # minimum storage capacity
        max_storage,                                                                             # maximum storage capacity
        init_storage,                                                                            # initial storage level
        availability_input,                                                                      # Hourly availability
        get_project_derating(project),                                                           # de-rating factor
        get_ramp_limits(get_tech(project)),                                                       # ramp limits
        get_project_max_reserveup(project),                                                      # maximum reserve up limit
        get_project_max_reservedown(project),                                                    # maximum reserve down limit
        existing_units,                                                                          # existing units
        units_inqueue,                                                                           # units in queue
        get_lag_time(finance_data),                                                               # construction lead time
        remaining_lag_time,                                                                      # remaining construction time
        1,                                                                                # maximum units
        1,                                                                                # base cost units
        get_capex_years(finance_data),                                                            # capital cost recovery years
        get_life_time(finance_data),                                                              # total life_time
        remaining_life_time,                                                                      # remaining life_time
        in(:Capacity, get_name.(get_products(project))),                                         # eligible for capacity markets
        in(:REC, get_name.(get_products(project))),                                              # eligible for rps compliance
        get_zone(get_tech(project)),                                                             # project zone
        [get_ownedby(finance_data)])                                                             # owned by

    return market_project
end

"""
This function processes data of Existing projects to create the MarketProject struct
passed to CEM price projection and endogeneous Economic Dispatch models.
"""
function create_market_project(project::P,
                              pricecap_energy::Float64,
                              pricecap_reserveup::Float64,
                              pricecap_reservedown::Float64,
                              max_peak_loads::AxisArrays.AxisArray{Float64, 1},
                              iteration_year::Int64,
                              num_hours::Int64,
                              num_invperiods::Int64,
                              availability_df::DataFrames.DataFrame) where P <: Project{Existing}

    finance_data = get_finance_data(project)

    queue_time = length(get_queue_cost(finance_data))

    # Get technical characteristics based on project type
    project_type,
    min_input,
    max_input,
    efficiency_in,
    efficiency_out,
    min_storage,
    max_storage,
    init_storage = get_technical_details(project)

    availability_input = find_project_availability_data(project, availability_df, num_invperiods, num_hours)

    # Calculate remaining life_time
    remaining_life_time = min(get_life_time(finance_data), get_end_life_year(project) - iteration_year + 1)

    units_inqueue = zeros(queue_time)

    # Existing project have completed their queue time
    units_inqueue[queue_time] = 1.0

    existing_units = 1
    remaining_lag_time = 0

    market_project = populate_market_project(project,
                                             get_project_reserve_up_cost(project, pricecap_reserveup),
                                             get_project_reserve_down_cost(project, pricecap_reservedown),
                                             project_type,
                                             min_input,
                                             max_input,
                                             efficiency_in,
                                             efficiency_out,
                                             min_storage,
                                             max_storage,
                                             init_storage,
                                             availability_input,
                                             existing_units,
                                             units_inqueue,
                                             remaining_lag_time,
                                             remaining_life_time,
                                             iteration_year,
                                             num_invperiods)
    return market_project

end

"""
This function processes data of Option projects to create the MarketProject struct
passed to CEM price projection and endogeneous Economic Dispatch models.
"""
function create_market_project(project::P,
                              pricecap_energy::Float64,
                              pricecap_reserveup::Float64,
                              pricecap_reservedown::Float64,
                              max_peak_loads::AxisArrays.AxisArray{Float64, 1},
                              iteration_year::Int64,
                              num_hours::Int64,
                              num_invperiods::Int64,
                              availability_df::DataFrames.DataFrame) where P <: Project{Option}

    # Get technical characteristics based on project type
    project_type,
    min_input,
    max_input,
    efficiency_in,
    efficiency_out,
    min_storage,
    max_storage,
    init_storage = get_technical_details(project)

    availability_input = find_project_availability_data(project, availability_df, num_invperiods, num_hours)

    finance_data = get_finance_data(project)

    queue_time = length(get_queue_cost(finance_data))

    # End of life_time for option projects can be up to the end of horizon.
    remaining_life_time = num_invperiods
    units_inqueue = zeros(queue_time)
    existing_units = 0
    remaining_lag_time = get_lag_time(finance_data) + queue_time

    market_project = populate_market_project(project,
                                             get_project_reserve_up_cost(project, pricecap_reserveup),
                                             get_project_reserve_down_cost(project, pricecap_reservedown),
                                             project_type,
                                             min_input,
                                             max_input,
                                             efficiency_in,
                                             efficiency_out,
                                             min_storage,
                                             max_storage,
                                             init_storage,
                                             availability_input,
                                             existing_units,
                                             units_inqueue,
                                             remaining_lag_time,
                                             remaining_life_time,
                                             iteration_year,
                                             num_invperiods)

    market_project.name = "option_$(market_project.tech_type)_$(market_project.zone)"

    if market_project.derating_factor == 0.0
        derating = 0.5                              # Dummy derating for maximum capacity if capacity product doesn't exist
    else
        derating = market_project.derating_factor
    end

    modified_max_gen = max_peak_loads[market_project.zone] / derating
    market_project.max_new_options = round(modified_max_gen / market_project.max_gen)

    return market_project

end

"""
This function processes data of Planned projects to create the MarketProject struct
passed to CEM price projection and endogeneous Economic Dispatch models.
"""
function create_market_project(project::P,
                              pricecap_energy::Float64,
                              pricecap_reserveup::Float64,
                              pricecap_reservedown::Float64,
                              max_peak_loads::AxisArrays.AxisArray{Float64, 1},
                              iteration_year::Int64,
                              num_hours::Int64,
                              num_invperiods::Int64,
                              availability_df::DataFrames.DataFrame) where P <: Project{Planned}

    # Get technical characteristics based on project type
    project_type,
    min_input,
    max_input,
    efficiency_in,
    efficiency_out,
    min_storage,
    max_storage,
    init_storage = get_technical_details(project)

    availability_input = find_project_availability_data(project, availability_df, num_invperiods, num_hours)

    finance_data = get_finance_data(project)

    queue_time = length(get_queue_cost(finance_data))

    # Number of construction years left
    remaining_lag_time = get_construction_year(project) - iteration_year

    units_inqueue = zeros(queue_time)

    # Planned units have completed their queue time
    units_inqueue[queue_time] = 1.0

    # Calculate remaining life_time
    remaining_life_time = get_life_time(finance_data) + remaining_lag_time

    existing_units = 0

    market_project = populate_market_project(project,
                                             get_project_reserve_up_cost(project, pricecap_reserveup),
                                             get_project_reserve_down_cost(project, pricecap_reservedown),
                                             project_type,
                                             min_input,
                                             max_input,
                                             efficiency_in,
                                             efficiency_out,
                                             min_storage,
                                             max_storage,
                                             init_storage,
                                             availability_input,
                                             existing_units,
                                             units_inqueue,
                                             remaining_lag_time,
                                             remaining_life_time,
                                             iteration_year,
                                             num_invperiods)
    return market_project

end

"""
This function processes data of Queue projects to create the MarketProject struct
passed to CEM price projection and endogeneous Economic Dispatch models.
"""
function create_market_project(project::P,
                              pricecap_energy::Float64,
                              pricecap_reserveup::Float64,
                              pricecap_reservedown::Float64,
                              max_peak_loads::AxisArrays.AxisArray{Float64, 1},
                              iteration_year::Int64,
                              num_hours::Int64,
                              num_invperiods::Int64,
                              availability_df::DataFrames.DataFrame) where P <: Project{Queue}

    # Get technical characteristics based on project type
    project_type,
    min_input,
    max_input,
    efficiency_in,
    efficiency_out,
    min_storage,
    max_storage,
    init_storage = get_technical_details(project)

    availability_input = find_project_availability_data(project, availability_df, num_invperiods, num_hours)

    finance_data = get_finance_data(project)

    queue_time = length(get_queue_cost(finance_data))

    units_inqueue = zeros(queue_time)

    queue_year = iteration_year - get_decision_year(project)

    # Calculate number of years left in queue
    if queue_year <= queue_time
        units_inqueue[queue_year] = 1.0
    end

    remaining_lag_time = get_lag_time(finance_data) + queue_time - queue_year

    # Calculate remaining life_time
    remaining_life_time = get_life_time(finance_data) + get_lag_time(finance_data) + queue_time - queue_year

    existing_units = 0

     market_project = populate_market_project(project,
                                             get_project_reserve_up_cost(project, pricecap_reserveup),
                                             get_project_reserve_down_cost(project, pricecap_reservedown),
                                             project_type,
                                             min_input,
                                             max_input,
                                             efficiency_in,
                                             efficiency_out,
                                             min_storage,
                                             max_storage,
                                             init_storage,
                                             availability_input,
                                             existing_units,
                                             units_inqueue,
                                             remaining_lag_time,
                                             remaining_life_time,
                                             iteration_year,
                                             num_invperiods)
    return market_project

end


text:
##################################################
## Create Bulge
## 4  cases DoubleShift/SingleShift x Descending/Twisted Q Factorizations
## RFactorization type doesn't enter here

# ## The bulge is created by  (A-rho1) * (A - rho2) * e_1 where rho1 and rho2 are eigenvalue or random
# ## for real case, we take the real part of this result

const IT_COUNT = 15  # if in 15 steps deflation does not occur, this kicks in a random rotator

"""
   create_bulge(state)

Finds m=1 or 2 shifts (m=1 in the CSS case, m=2 in the RDS case) based on the eigenvalues of the lower 2x2 block (using `stop_index`). The
vector `x = alpha (A - rho_1 I) e_1` or `x = alpha (A-rho_1 I) (A-rho_2 I) e1` is found. From this, one or two core transforms
are found so that `U_1' x = gamma e_1` or `U_1' U_2' x = gamma e_1`. The values `U_1` or `U_1`, `U_2` are store in `state`.

"""
function create_bulge(QF::QFactorization{T,S,VV}, RF, storage,  ctr) where {T, S <: Real,VV}

    A = storage.A

    if iszero(ctr.it_count)
        ctr.it_count = IT_COUNT

        t = rand(T) * pi

        re1, ie1 = sincos(t)
        re2, ie2 = re1, -ie1

        i = ctr.start_index
        j = i + 1

        U =  Rotator(re1, ie1, i)
        V = Rotator(re2, ie2, j)

    else

        Δ = ctr.stop_index
        diagonal_block!(A, QF, RF, Δ, Δ) #+1)
        e1,  e2 = eigen_values(A)

        l1r, l1i  =  real(e1), imag(e1)
        l2r, l2i =  real(e2), imag(e2)

        delta = ctr.start_index

        diagonal_block!(A, QF, RF,  delta, delta)
        bk11, bk12 = A[1,1], A[1,2]
        bk21, bk22 = A[2,1], A[2,2]

        diagonal_block!(A, QF, RF, delta+1, delta+1)
        bk32 = A[2,1]

        # compute `x`
        c1 = -l1i * l2i + l1r*l2r -l1r*bk11 -l2r * bk11 + bk11^2 + bk12 * bk21
        c2 = -l1r * bk21 - l2r * bk21 + bk11* bk21 + bk21 * bk22
        c3 = bk21 * bk32

        c, s, nrm = givensrot(c2, c3)
        i = ctr.start_index #1
        j = i + 1
        V = Rotator(c, -s, j)

        c, s, tmp = givensrot(c1, nrm)
        U = Rotator(c, -s, i)
    end

    storage.VU[2] = U
    storage.VU[1] = V

    return nothing
end

## CSS case
function create_bulge(QF::QFactorization{T, S, VV}, RF, storage, ctr) where {T, S <: Complex, VV}
    ray = true # true seems to take  fewer steps than false

    A = storage.A

    if iszero(ctr.it_count)
        ctr.it_count = IT_COUNT

        t = rand(T) * pi
        if ray
            shift = complex(cos(t), sin(t))
        else
            shift = complex(cos(t), zero(T))
        end

    else

        Delta = ctr.stop_index
        diagonal_block!(A, QF, RF, Delta, Delta)

        if ray
            # Wilkinson
            e1, e2 = eigen_values(A[1,1], A[1,2], A[2,1], A[2,2])
            shift = norm(A[2,2] - e1) < norm(A[2,2] - e2) ? e1 : e2
        else
            shift = A[2,2]
        end

    end

    delta =  ctr.start_index
    diagonal_block!(A, QF, RF, delta,  delta)
    c,s,nrm = givensrot(A[1,1] - shift, A[2,1])
    R = Rotator(c, s, delta)
    storage.VU[1] = R'

    nothing

end

##################################################
##
##  Twisted case

function create_bulge(QF::QFactorizationTwisted{T, S}, RF, storage, ctr, m) where {T, S}

    if iszero(ctr.it_count)
        ctr.it_count = IT_COUNT

        for i in 1:m
            storage.VU[i] = random_rotator.(S, ctr.start_index + m - i)
        end
    else
        _create_bulge(QF, RF, storage, ctr, m)
    end

    return nothing

end

## Real
function create_bulge(QF::QFactorizationTwisted{T, S}, RF, storage, ctr, m) where   {T, S <: Real}

    @assert m >= 2  # for even we need atleast a quadratic poly

    if iszero(ctr.it_count)
        ctr.it_count = IT_COUNT

        t = rand(T) * pi

        re1, ie1 = sincos(t)
        re2, ie2 = re1, -ie1

        i = ctr.start_index
        j = i + 1

        U =  Rotator(re1, ie1, i)
        V = Rotator(re2, ie2, j)

        storage.VU[1], storage.VU[2] = V, U
        return nothing
    end

    delta, Delta = ctr.start_index, ctr.stop_index
    n = Delta - delta
    A = storage.A

    offset = delta - 1

    diagonal_block!(A, QF, RF, Delta-m+1, Delta)

    Am = view(A, 2:m+1, 2:m+1)

    if m == 2

        rhos = eigen_values(Am)

    else

        rhos = _eigvals(Am)  ## sorted by LinearAlgebra.sorteig!

    end

    x = zeros(S, m+1)
    x[1] = one(S)

    diagonal_block!(A, QF, RF, delta, delta + m - 1)
    Am = view(A, 1:m+1, 1:m+1)

    i = 1
    while i <= m
        rhoi = rhos[i]
        pi = (i==1) ? :left : QF.Q.pv[i-1]
        pj = QF.Q.pv[i]

        if i < m && isapprox(rhoi, conj(rhos[i+1]))
            ## XXX deal with left right...
            x = (Am^2 - 2real(rhos[i])*Am + norm(rhos[i])^2*I) * x
            if pi == :right
                x = Am \ x
            end
            if pj == :right
                x = Am \ x
            end
            i += 1
        else
            x = (Am - real(rhos[i])*I) * x
            if pi == :right
                x = Am \ x
            end
        end

        i += 1
    end

    # now roll up
    _rollup!(storage.VU, x, delta-1)
    return nothing

end

## complex
function _create_bulge(QF::QFactorizationTwisted{T, S}, RF, storage, ctr, m) where   {T, S <:  Complex}

    delta, Delta = ctr.start_index, ctr.stop_index
    n = Delta -  delta
    A = storage.A

    diagonal_block!(A,QF, RF, Delta-m+1, Delta)

    if m == 1

        Am = view(A, 1:2, 1:2)
        e1, e2 = eigen_values(Am)

        if norm(Am[2,2] - e1) < norm(Am[2,2] - e2)
            rhos = [e1]
        else
            rhos = [e2]
        end

    elseif m == 2

        Am = view(A, 2:3, 2:3)
        rhos = eigen_values(Am)

    else

        Am = view(A, 2:m+1, 2:m+1)
        rhos = _eigvals(Am)

    end

    x = zeros(S, m+1)
    x[1] = one(S)

    diagonal_block!(A, QF, RF, delta, delta + m - 1)
    Am = view(A, 1:m+1, 1:m+1)

    for i in 1:m
        x = (Am - rhos[i] * I) * x
        pi = (i==1) ? :left : QF.Q.pv[i-1]
        if  pi == :right
            x = Am \ x
        end
    end

    # now roll up
    _rollup!(storage.VU, x, delta - 1)

    return nothing
end

function _rollup!(VU, x, offset=0)
    r = x[end]
    m = length(x)
    for i in (m-1):-1:1
        c, s, r = givensrot(x[i], r)
        U = Rotator(c, s, i+offset)
        VU[m-i] = U'
    end
end

ext: jl


text:
function is_key_pressed(key::Int)
    return ccall((:sfKeyboard_isKeyPressed, libcsfml_window), Int32, (Int32,), key) == 1
end

baremodule KeyCode
    const UNKNOWN = -1
    const A = 0
    const B = 1
    const C = 2
    const D = 3
    const E = 4
    const F = 5
    const G = 6
    const H = 7
    const I = 8
    const J = 9
    const K = 10
    const L = 11
    const M = 12
    const N = 13
    const O = 14
    const P = 15
    const Q = 16
    const R = 17
    const S = 18
    const T = 19
    const U = 20
    const V = 21
    const W = 22
    const X = 23
    const Y = 24
    const Z = 25
    const NUM0 = 26
    const NUM1 = 27
    const NUM2 = 28
    const NUM3 = 29
    const NUM4 = 30
    const NUM5 = 31
    const NUM6 = 32
    const NUM7 = 33
    const NUM8 = 34
    const NUM9 = 35
    const ESCAPE = 36
    const LCONTROL = 37
    const LSHIFT = 38
    const LALT = 39
    const LSYSTEM = 40
    const RCONTROL = 41
    const RSHIFT = 42
    const RALT = 43
    const RSYSTEM = 44
    const MENU = 45
    const LBRACKET = 46
    const RBRACKET = 47
    const SEMICOLON = 48
    const COMMA = 49
    const PERIOD = 50
    const QUOTE = 51
    const SLASH = 52
    const BACKSLASH = 53
    const TILDE = 54
    const EQUAL = 55
    const DASH = 56
    const SPACE = 57
    const RETURN = 58
    const BACK = 59
    const TAB = 60
    const PAGEUP = 61
    const PAGEDOWN = 62
    const END = 63
    const HOME = 64
    const INSERT = 65
    const DELETE = 66
    const ADD = 67
    const SUBTRACT = 68
    const MULTIPLY = 69
    const DIVIDE = 70
    const LEFT = 71
    const RIGHT = 72
    const UP = 73
    const DOWN = 74
    const NUMPAD0 = 75
    const NUMPAD1 = 76
    const NUMPAD2 = 77
    const NUMPAD3 = 78
    const NUMPAD4 = 79
    const NUMPAD5 = 80
    const NUMPAD6 = 81
    const NUMPAD7 = 82
    const NUMPAD8 = 83
    const NUMPAD9 = 84
    const F1 = 85
    const F2 = 86
    const F3 = 87
    const F4 = 88
    const F5 = 89
    const F6 = 90
    const F7 = 91
    const F8 = 92
    const F9 = 93
    const F10 = 94
    const F11 = 95
    const F12 = 96
    const F13 = 97
    const F14 = 98
    const F15 = 99
    const PAUSE = 100
end


text:
using ParticleFilters

const N_PARTICLES = 1000

function initial_particle_belief(pomdp::OCPOMDP)
    d0 = initialstate_distribution(pomdp)
    p = Vector{OCState}(N_PARTICLES)
    for i=1:N_PARTICLES
        p[i] = rand(rng, d0)
    end
    return ParticleCollection{OCState}(p)
end

function init_updater(pomdp::OCPOMDP, N::Int64 = N_PARTICLES)
    filter = SIRParticleFilter(pomdp, N)
    return filter
end

function POMDPs.action(policy::QMDPPolicy, b::ParticleCollection{OCState})
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
"""
    LikelihoodRatioTest

Results of MixedModels.likelihoodratiotest

## Fields
* `formulas`: Vector of model formulae
* `models`: NamedTuple of the `dof` and `deviance` of the models
* `tests`: NamedTuple of the sequential `dofdiff`, `deviancediff`, and resulting `pvalues`

## Properties
* `deviance`
* `pvalues`

"""
struct LikelihoodRatioTest
    formulas::AbstractVector{String}
    models::NamedTuple{(:dof,:deviance)}
    tests::NamedTuple{(:dofdiff,:deviancediff,:pvalues)}
end

Base.propertynames(lrt::LikelihoodRatioTest, private = false) = (
    :deviance,
    :formulas,
    :models,
    :pvalues,
    :tests,
)

function Base.getproperty(lrt::LikelihoodRatioTest, s::Symbol)
    if s == :dof
        lrt.models.dof
    elseif s == :deviance
        lrt.models.deviance
    elseif s == :pvalues
        lrt.tests.pvalues
    elseif s == :formulae
        lrt.formulas
    else
        getfield(lrt, s)
    end
end

# backward syntactic but not type compatibility
Base.getindex(lrt::LikelihoodRatioTest, s::Symbol) = getfield(lrt,s)

"""
    likelihoodratiotest(m::LinearMixedModel...)

Likeihood ratio test applied to a set of nested models.

Note that nesting of the models is not checked.  It is incumbent on the user to check this.
"""
function likelihoodratiotest(m::LinearMixedModel...)
    allequal(getproperty.(getproperty.(m,:optsum),:REML)) ||
        throw(ArgumentError("Models must all be fit with the same objective (i.e. alll ML or all REML)"))
    if any(getproperty.(getproperty.(m,:optsum),:REML))
        allequal(coefnames.(m))  ||
                throw(ArgumentError("Likelihood-ratio tests for REML-fitted models are only valid when the fixed-effects specifications are identical"))
    end
    _likelihoodratiotest(m...)
end

function likelihoodratiotest(m::GeneralizedLinearMixedModel...)
    # TODO: test that all models are fit with same fast/nAGQ option?
    glms = getproperty.(m,:resp);
    allequal(Distribution.(glms)) ||
        throw(ArgumentError("Models must be fit to the same distribution"))
    allequal(string.(Link.(glms))) ||
        throw(ArgumentError("Models must have the same link function"))

    _likelihoodratiotest(m...)
end

function _likelihoodratiotest(m::Vararg{T}) where T <: MixedModel
    m = collect(m)   # change the tuple to an array
    dofs = dof.(m)
    formulas = String.(Symbol.(getproperty.(m,:formula)))
    ord = sortperm(dofs)
    dofs = dofs[ord]
    formulas = formulas[ord]
    devs = objective.(m)[ord]
    dofdiffs = diff(dofs)
    devdiffs = .-(diff(devs))
    pvals = ccdf.(Chisq.(dofdiffs), devdiffs)

    LikelihoodRatioTest(
        formulas,
        (dof = dofs, deviance = devs),
        (dofdiff = dofdiffs, deviancediff = devdiffs, pvalues = pvals)
    )
end

function _array_union_nothing(arr::Array{T}) where T
    Array{Union{T,Nothing}}(arr)
end

function _prepend_0(arr::Array{T}) where T
    pushfirst!(copy(arr), -zero(T))
end

function Base.show(io::IO, lrt::LikelihoodRatioTest; digits=2)
    println(io, "Model Formulae")

    for (i, f) in enumerate(lrt.formulas)
        println(io, "$i: $f")
    end

    # the following was adapted from StatsModels#162
    # from nalimilan
    Δdf = lrt.tests.dofdiff
    Δdev = lrt.tests.deviancediff

    nc = 6
    nr = length(lrt.formulas)
    outrows = Matrix{String}(undef, nr+1, nc)

    outrows[1, :] = ["",
                    "model-dof",
                    "deviance",
                    "χ²",
                    "χ²-dof",
                    "P(>χ²)"] # colnms

    outrows[2, :] = ["[1]",
                    string(lrt.dof[1]),
                    Ryu.writefixed(lrt.deviance[1], 4),
                    " "," ", " "]

    for i in 2:nr
        outrows[i+1, :] = ["[$i]",
                           string(lrt.dof[i]),
                           Ryu.writefixed(lrt.deviance[i], 4),
                           Ryu.writefixed(Δdev[i-1], 4),
                           string(Δdf[i-1]),
                           string(StatsBase.PValue(lrt.pvalues[i-1]))]
    end
    colwidths = length.(outrows)
    max_colwidths = [maximum(view(colwidths, :, i)) for i in 1:nc]
    totwidth = sum(max_colwidths) + 2*5

    println(io, '─'^totwidth)

    for r in 1:nr+1
        for c in 1:nc
            cur_cell = outrows[r, c]
            cur_cell_len = length(cur_cell)

            padding = " "^(max_colwidths[c]-cur_cell_len)
            if c > 1
                padding = "  "*padding
            end

            print(io, padding)
            print(io, cur_cell)
        end
        print(io, "\n")
        r == 1 && println(io, '─'^totwidth)
    end
    print(io, '─'^totwidth)

    nothing
end


text:
using MDSplus

function status_ok(stat)
    (stat & 1) == 1
end

function get_signal_length(sig)
    s = [0]
    idesc = descr(Int32, size(s), s)
    stat = MdsValue("SIZE($sig)", idesc)
    return s[1]
end

function run_example(shot)
    socket = MdsConnect("atlas.gat.com")

    stat = MdsOpen("EFIT01",shot)

    s = get_signal_length("\\AMINOR")

    data = zeros(Float64, s)
    sigdesc = descr(data)
    stat = MdsValue("\\AMINOR", sigdesc)

    time = zeros(Float64, s)
    timedesc = descr(time)
    stat = MdsValue("DIM_OF(\\AMINOR)", timedesc)

    for i=1:10
        println(@sprintf("%2i X: %6.2f Y: %6.2f",i,time[i],data[i]))
    end

    MdsDisconnect()
end

run_example(107000)


text:
#=
# Snippet used to generate reference derivatives
using AnalyticBSplines
function analytic_derivatives(kl, kr, a, b, N,
                              mll=kl, mlr=kl,
                              mrl=kr, mrr=kr)
    tl = LinearKnotSet(kl, a//1, b//1, N, mll, mlr)
    tr = LinearKnotSet(kr, a//1, b//1, N, mrl, mrr)

    L = gen_basis(tl, mll==kl ? 1//1 : 0, mlr==kl ? 1//1 : 0)
    R = gen_basis(tr, mrr==kr ? 1//1 : 0, mrr==kr ? 1//1 : 0)

    [derop(L, R, tl, o) for o=1:min(kl,kr)-1]
end
=#

@testset "Derivatives" begin
    @testset "Derivative operators" begin
        @testset "k = 3" begin
            @testset "One basis" begin
                k = 3
                t = LinearKnotSet(k, 0, 1, 2)
                B = BSpline(t)
                D = Derivative(axes(B,1))

                @test B'D*B ≈ [  -1//2    5//12   1//12   0
                                 -5//12   0//1    1//3   1//12
                                 -1//12  -1//3    0//1   5//12
                                 0     -1//12  -5//12  1//2]

                @test B'D'D*B ≈ [ -8//3   2//1   2//3    0
                                  2//1  -8//3   0//1   2//3
                                  2//3   0//1  -8//3   2//1
                                  0     2//3   2//1  -8//3]
            end
            @testset "kl = 3, kr = 4" begin
                tl = LinearKnotSet(3, 0, 1, 2)
                tr = LinearKnotSet(4, 0, 1, 2)

                N = max(CompactBases.num_quadrature_points(3,1),
                        CompactBases.num_quadrature_points(4,1))

                L = BSpline(tl, N)
                R = BSpline(tr, N)
                D = Derivative(axes(R,1))

                @test L'D*R ≈ [ -3//5   17//40   3//20    1//40  0//1
                                -7//20  -3//20   1//5     1//4   1//20
                                -1//20  -1//4   -1//5     3//20  7//20
                                0//1   -1//40  -3//20  -17//40  3//5 ]

                @test L'D'D*R ≈ [ -3//1   7//4   1//1   1//4   0//1
                                  5//2  -5//2  -1//1   1//2   1//2
                                  1//2   1//2  -1//1  -5//2   5//2
                                  0//1   1//4   1//1   7//4  -3//1]
            end
            @testset "kl = 4, kr = 3" begin
                tl = LinearKnotSet(4, 0, 1, 2)
                tr = LinearKnotSet(3, 0, 1, 2)

                N = max(CompactBases.num_quadrature_points(4,1),
                        CompactBases.num_quadrature_points(3,1))

                L = BSpline(tl, N)
                R = BSpline(tr, N)
                D = Derivative(axes(R,1))

                @test L'D*R ≈ [  -2//5    7//20   1//20   0//1
                                 -17//40   3//20   1//4    1//40
                                 -3//20  -1//5    1//5    3//20
                                 -1//40  -1//4   -3//20  17//40
                                 0//1   -1//20  -7//20   2//5]

                @test L'D'D*R ≈ [ -3//1   5//2   1//2   0//1
                                  7//4  -5//2   1//2   1//4
                                  1//1  -1//1  -1//1   1//1
                                  1//4   1//2  -5//2   7//4
                                  0//1   1//2   5//2  -3//1]
            end
        end
        @testset "k = 7" begin
            k = 7
            t = LinearKnotSet(k, 0, 1, 3)
            B = BSpline(t)
            D = Derivative(axes(B,1))

            @test B'D*B ≈ [      -1//2          10625//29568     283021//2395008    23827//1197504     479//199584        5//27216          1//149688         0                0
                                 -10625//29568          0//1          35867//177408    118891//1064448   26167//709632    11017//1419264     7141//8515584       1//946176         0
                                 -283021//2395008   -35867//177408         0//1           5129//39424     11645//101376   245743//4257792    20455//1216512    7141//8515584       1//149688
                                 -23827//1197504  -118891//1064448    -5129//39424          0//1         11035//118272   218689//2128896   245743//4257792   11017//1419264       5//27216
                                 -479//199584    -26167//709632    -11645//101376    -11035//118272        0//1         11035//118272     11645//101376    26167//709632      479//199584
                                 -5//27216     -11017//1419264  -245743//4257792  -218689//2128896  -11035//118272        0//1           5129//39424    118891//1064448   23827//1197504
                                 -1//149688     -7141//8515584   -20455//1216512  -245743//4257792  -11645//101376    -5129//39424          0//1         35867//177408   283021//2395008
                                 0               -1//946176     -7141//8515584   -11017//1419264  -26167//709632  -118891//1064448   -35867//177408        0//1         10625//29568
                                 0                 0               -1//149688        -5//27216      -479//199584   -23827//1197504  -283021//2395008  -10625//29568         1//2]

            @test B'D'D*B ≈ [   -108//11       16509//2464      169651//66528      2347//4752        377//5544         7//1188           1//4158            0                0
                                16509//2464     -4413//616      -104639//133056   191011//266112    72047//177408  118771//1064448    31501//2128896        3//78848          0
                                169651//66528  -104639//133056    -46861//22176    -16405//29568     58403//177408   45935//118272    121843//709632     31501//2128896       1//4158
                                2347//4752    191011//266112    -16405//29568     -3083//3168     -36415//88704    39203//177408     45935//118272    118771//1064448       7//1188
                                377//5544     72047//177408     58403//177408   -36415//88704     -5807//7392    -36415//88704      58403//177408     72047//177408      377//5544
                                7//1188    118771//1064448    45935//118272    39203//177408   -36415//88704    -3083//3168      -16405//29568     191011//266112     2347//4752
                                1//4158     31501//2128896   121843//709632    45935//118272    58403//177408  -16405//29568     -46861//22176    -104639//133056   169651//66528
                                0              3//78848      31501//2128896  118771//1064448   72047//177408  191011//266112   -104639//133056     -4413//616       16509//2464
                                0               0                1//4158          7//1188        377//5544      2347//4752      169651//66528      16509//2464       -108//11]

            @test CompactBases.diffop!(similar(B.S), B, B, 1, 2) ≈
                [    162//1     -49527//224     93901//2016    10723//1008     41//24        43//252        1//126          0             0
                     -23049//224      162//1      -78221//1344   -13759//2688   4773//1792    2023//1536    5039//21504       9//7168        0
                     -93901//2016   78221//1344        0//1       -9741//896    -905//256     4079//3584   10065//7168     5039//21504      1//126
                     -10723//1008   13759//2688     9741//896         0//1     -4621//896    -5055//1792    4079//3584     2023//1536      43//252
                     -41//24     -4773//1792      905//256      4621//896       0//1      -4621//896     -905//256      4773//1792      41//24
                     -43//252    -2023//1536    -4079//3584     5055//1792   4621//896        0//1      -9741//896    -13759//2688   10723//1008
                     -1//126    -5039//21504  -10065//7168    -4079//3584    905//256     9741//896        0//1      -78221//1344   93901//2016
                     0           -9//7168    -5039//21504   -2023//1536  -4773//1792   13759//2688   78221//1344     -162//1      23049//224
                     0             0            -1//126       -43//252     -41//24    -10723//1008  -93901//2016    49527//224     -162//1]

            @test CompactBases.diffop!(similar(B.S), B, B, 2, 2) ≈
                [    2700//1    -406215//112     232415//336      32465//168     1045//28       185//42          5//21            0               0
                     -406215//112   141615//28     -781435//672    -354385//1344  -16085//896    56255//5376    34625//10752      135//3584          0
                     232415//336  -781435//672      57055//112       7365//448   -45905//896   -25215//1792    25535//3584     34625//10752        5//21
                     32465//168  -354385//1344      7365//448       8255//112      595//64     -3815//128    -25215//1792     56255//5376       185//42
                     1045//28    -16085//896     -45905//896        595//64      5055//112      595//64     -45905//896     -16085//896       1045//28
                     185//42     56255//5376    -25215//1792     -3815//128      595//64      8255//112      7365//448    -354385//1344     32465//168
                     5//21     34625//10752    25535//3584    -25215//1792  -45905//896     7365//448     57055//112    -781435//672     232415//336
                     0          135//3584     34625//10752    56255//5376  -16085//896  -354385//1344  -781435//672     141615//28     -406215//112
                     0             0              5//21         185//42      1045//28     32465//168    232415//336    -406215//112       2700//1]

            @test CompactBases.diffop!(similar(B.S), B, B, 2, 3) ≈
                [  -36450//1    6885405//112   -3209175//112     158535//56      19305//28       1395//14          45//7             0               0
                   5361795//112  -164025//2      8973045//224   -2228985//448   -866295//896     -3375//256     127305//3584      3645//3584          0
                   -873225//112  3274155//224     -18225//2      1109295//448     26325//896   -421605//1792    -10125//512     127305//3584        45//7
                   -158535//56   2228985//448   -1109295//448          0//1      202095//448     26325//896    -421605//1792     -3375//256       1395//14
                   -19305//28    866295//896     -26325//896    -202095//448         0//1      202095//448      26325//896    -866295//896      19305//28
                   -1395//14      3375//256     421605//1792    -26325//896   -202095//448         0//1      1109295//448   -2228985//448     158535//56
                   -45//7    -127305//3584     10125//512     421605//1792   -26325//896  -1109295//448      18225//2     -3274155//224     873225//112
                   0        -3645//3584   -127305//3584      3375//256    866295//896   2228985//448   -8973045//224     164025//2     -5361795//112
                   0             0            -45//7        -1395//14     -19305//28    -158535//56     3209175//112   -6885405//112      36450//1]

            @test CompactBases.diffop!(similar(B.S), B, B, 3, 3) ≈
                [ -3499200//7   1673055//2     -759645//2     215595//7       72090//7      13500//7        1080//7            0              0
                  1673055//2  -9863370//7    18273735//28   -474525//8    -2080485//112   -79245//32     105435//448     10935//448          0
                  -759645//2  18273735//28   -4495905//14   2229525//56     157545//16   -176175//224   -419985//448    105435//448      1080//7
                  215595//7   -474525//8     2229525//56   -186705//14      26325//56    325215//112   -176175//224    -79245//32      13500//7
                  72090//7  -2080485//112    157545//16     26325//56     -57105//14     26325//56     157545//16   -2080485//112     72090//7
                  13500//7    -79245//32    -176175//224   325215//112     26325//56   -186705//14    2229525//56    -474525//8      215595//7
                  1080//7    105435//448   -419985//448  -176175//224    157545//16   2229525//56   -4495905//14   18273735//28    -759645//2
                  0       10935//448    105435//448   -79245//32   -2080485//112  -474525//8    18273735//28   -9863370//7     1673055//2
                  0            0          1080//7      13500//7       72090//7     215595//7     -759645//2     1673055//2    -3499200//7]
        end
    end

    @testset "Derivative accuracy" begin
        @info "B-splines derivative accuracy"
        f,g,h,a,b = derivative_test_functions(1)
        Ns = ceil.(Int, 10 .^ range(2, stop=3, length=10))
        @testset "k = $k" for k = 3:6
            @show k
            hs,ϵg,ϵh,ϵh′,pg,ph,ph′ = compute_derivative_errors(Ns, f, g, h, 1) do N
                t = LinearKnotSet(k, a, b, N)
                BSpline(t)
            end
            # These requirements are not derived, but empirical
            # results, mainly to ensure not regression occurs.
            @test pg > k + 3
            @test ph > k + (k > 3 && k < 6 ? 3 : 1)
            @test ph′ > k + (k > 3 && k < 6 ? 2 : 1)
        end
    end

    @testset "Hooke's law" begin
        function ode_hookes_law(xₘₐₓ, kspring, k, N)
            t = LinearKnotSet(k, 0, xₘₐₓ, N)
            # By omitting the first basis function, we enforce V(0) = 0
            B = BSpline(t,k′=1)[:,2:end]
            S = B'B

            x = axes(B,1)

            D = Derivative(x)
            ∇ = B'*D*B

            # Hooke's law
            F = x -> -kspring*x

            # Expand Hooke's law on B-splines
            cF = B \ F.(x)
            # Solve for expansion coefficients of potential
            cV = -∇ \ S*cF

            B,cF,cV
        end

        kspring = 0.1
        @testset "k=3" begin
            @test last(ode_hookes_law(3, kspring, 3, 1)) ≈ [0, 0.45] atol=1e-15
        end
        @testset "k=7" begin
            B,cF,cV = ode_hookes_law(3, kspring, 7, 30)
            @test cV ≈ (B \ (x -> kspring*x^2/2).(axes(B,1)))
        end
    end
end


text:
using FastGaussQuadrature, Test

@testset "Gauss–Radau" begin
    @testset "Check error" begin
        @test_throws DomainError gaussradau(0)
        @test_throws DomainError gaussradau(0,0,0)
    end

    @testset "Legendre" begin
        n = 1
        x, w = gaussradau(n)
        @test x[1] ≈ -1.
        @test w[1] ≈ 2.

        n = 2
        x, w = gaussradau(n)
        @test x[1] ≈ -1.
        @test x[2] ≈ 1/3
        @test w[1] ≈ .5
        @test w[2] ≈ 1.5

        tol = 1e-14
        n = 42
        ntests = 4
        x, w = gaussradau(n)
        @test length(x) == n && length(w) == n
        @test abs(x[37] - 0.908847278001044) < tol
        @test abs(w[37] - 0.031190846817016) < tol
        @test x[1] == -1

        @test dot( w,exp.(x)) ≈ exp(1)-exp(-1)
    end
    @testset "Jacobi" begin
        for (a,b) in ((0.1, 0.2), (0,-0.5))
            # compare with Gauss–Jacobi
            for n = 1:5
                x, w = gaussradau(n, a, b)
                x̃, w̃ = gaussjacobi(n, a, b)

                @test length(x) == n

                @test sum(w) ≈ sum(w̃)
                for j = 1:2n-2
                    @test dot(w,x.^j) ≈ dot(w̃,x̃.^j)
                end
            end
        end

        # compare with Gauss-Radau-Legendre
        @test all(gaussradau(2, 0, 0) .≈ gaussradau(2))
        @test all(gaussradau(3, 0, 0) .≈ gaussradau(3))
    end
end


text:
"""
cuboid ranges
convert given ranges and print overview on screen
used for gethydro, getparticles, getgravity, getclumps..., subregions...
"""
function prepranges(    dataobject::InfoType,
                        range_unit::Symbol,
                        verbose::Bool,
                        xrangem::Array{<:Any,1},
                        yrangem::Array{<:Any,1},
                        zrangem::Array{<:Any,1},
                        center::Array{<:Any,1};
                        dataranges::Array{<:Real,1}=[0.,1., 0.,1., 0.,1.] )

    xrange = zeros(Float64,2)
    yrange = zeros(Float64,2)
    zrange = zeros(Float64,2)

    selected_unit = 1. # :standard
    conv = 1. # :standard, variable used to convert to standard units
    if range_unit != :standard
        selected_unit = getunit(dataobject, range_unit)
        conv = dataobject.boxlen * selected_unit
    end

    center = prepboxcenter(dataobject, range_unit, center)

    # assign ranges to selected data range or missing ranges = full box
    if xrangem[1] === missing
        xmin = dataranges[1]
    else
        xrange[1]=xrangem[1]
        xmin = (xrange[1] + center[1]) / conv
    end

    if xrangem[2] === missing
        xmax = dataranges[2]
    else
        xrange[2]=xrangem[2]
        xmax = (xrange[2] + center[1]) / conv
    end

    if yrangem[1] === missing
        ymin = dataranges[3]
    else
        yrange[1]=yrangem[1]
        ymin = (yrange[1] + center[2]) / conv
    end

    if yrangem[2] === missing
        ymax = dataranges[4]
    else
        yrange[2] = yrangem[2]
        ymax = (yrange[2] + center[2]) / conv
    end

    if zrangem[1] === missing
        zmin = dataranges[5]
    else
        zrange[1] = zrangem[1]
        zmin = (zrange[1] + center[3]) / conv
    end

    if zrangem[2] === missing
        zmax = dataranges[6]
    else
        zrange[2]=zrangem[2]
        zmax = (zrange[2] + center[3]) / conv
    end

    center = center ./ conv

    # ensure that min-var is minimum and max-var is maximum of each dimension
    if xmin > xmax  error("[Mera]: xmin > xmax") end
    if ymin > ymax  error("[Mera]: ymin > ymax") end
    if zmin > zmax  error("[Mera]: zmin > zmax") end
    if xmax < xmin  error("[Mera]: xmax < xmin") end
    if ymax < ymin  error("[Mera]: ymax < ymin") end
    if zmax < zmin  error("[Mera]: zmax < zmin") end

    # ensure that ranges are inside the box
    xmin = maximum( [xmin, dataranges[1]] )
    ymin = maximum( [ymin, dataranges[3]] )
    zmin = maximum( [zmin, dataranges[5]] )
    xmax = minimum( [xmax, dataranges[2]] )
    ymax = minimum( [ymax, dataranges[4]] )
    zmax = minimum( [zmax, dataranges[6]] )

    if verbose
        if center != [0., 0., 0.]
            print("center: $(round.(center, digits=7)) ")
            center_val1, center_unit1= humanize(center[1] * dataobject.boxlen , dataobject.scale, 3, "length")
            center_val2, center_unit2= humanize(center[2] * dataobject.boxlen , dataobject.scale, 3, "length")
            center_val3, center_unit3= humanize(center[3] * dataobject.boxlen , dataobject.scale, 3, "length")
            if center_val1 == 0. center_unit1 = xmax_unit end
            if center_val2 == 0. center_unit2 = xmax_unit end
            if center_val3 == 0. center_unit3 = xmax_unit end
            print("==> [$center_val1 [$center_unit1] :: $center_val2 [$center_unit2] :: $center_val3 [$center_unit3]]\n")
            println()
        end

        #println("domain \t \t \t \thuman-readable units")
        println("domain:")

        print("xmin::xmax: $(round(xmin, digits=7)) :: $(round(xmax, digits=7))  \t")
        xmin_val, xmin_unit  = humanize(xmin * dataobject.boxlen, dataobject.scale, 3, "length")
        xmax_val, xmax_unit  = humanize(xmax * dataobject.boxlen, dataobject.scale, 3, "length")
        if xmin_val == 0. xmin_unit = xmax_unit end
        print("==> $xmin_val [$xmin_unit] :: $xmax_val [$xmax_unit]\n")

        print("ymin::ymax: $(round(ymin, digits=7)) :: $(round(ymax, digits=7))  \t")
        ymin_val, ymin_unit  = humanize(ymin * dataobject.boxlen, dataobject.scale, 3, "length")
        ymax_val, ymax_unit  = humanize(ymax * dataobject.boxlen, dataobject.scale, 3, "length")
        if ymin_val == 0. ymin_unit = xmax_unit end
        print("==> $ymin_val [$ymin_unit] :: $ymax_val [$ymax_unit]\n")

        print("zmin::zmax: $(round(zmin, digits=7)) :: $(round(zmax, digits=7))  \t")
        zmin_val, zmin_unit  = humanize(zmin * dataobject.boxlen, dataobject.scale, 3, "length")
        zmax_val, zmax_unit  = humanize(zmax * dataobject.boxlen, dataobject.scale, 3, "length")
        if zmin_val == 0. zmin_unit = xmax_unit end
        print("==> $zmin_val [$zmin_unit] :: $zmax_val [$zmax_unit]\n")
        println()
    end

    ranges = [xmin, xmax, ymin, ymax, zmin, zmax]
    return ranges
end

"""
cylinder ranges (height != 0.), sphere ranges (height==0.)
convert given ranges + radius and print overview on screen
used for subregions...
"""
function prepranges(    dataobject::InfoType,
                            center::Array{<:Any,1},
                            radius::Real,
                            height::Real,
                            range_unit::Symbol,
                            verbose::Bool)

    center = prepboxcenter(dataobject, range_unit, center)

    selected_unit = 1.
    if range_unit == :standard
        xmin = -radius + center[1]
        xmax =  radius + center[1]
        ymin = -radius + center[2]
        ymax =  radius + center[2]
        if height != 0.
            zmin = - height + center[3]
            zmax =   height + center[3]
        else
            zmin = - radius + center[3]
            zmax =   radius + center[3]
        end
        # given center relative to the data range in units: cell centers
        #todo

    else
        selected_unit = getunit(dataobject, range_unit)
        xmin = (-radius + center[1]) * selected_unit /dataobject.boxlen
        xmax = ( radius + center[1]) * selected_unit /dataobject.boxlen
        ymin = (-radius + center[2]) * selected_unit /dataobject.boxlen
        ymax = ( radius + center[2]) * selected_unit /dataobject.boxlen
        if height != 0.
            zmin = ( -height + center[3]) * selected_unit /dataobject.boxlen
            zmax = (  height + center[3]) * selected_unit /dataobject.boxlen
        else
            zmin = ( -radius + center[3]) * selected_unit /dataobject.boxlen
            zmax = (  radius + center[3]) * selected_unit /dataobject.boxlen
        end
        # given center relative to the data range in units: cell centers
        cx_shift = center[1] * selected_unit /dataobject.boxlen
        cy_shift = center[2] * selected_unit /dataobject.boxlen
        cz_shift = center[3] * selected_unit /dataobject.boxlen
        radius_shift = radius * selected_unit /dataobject.boxlen
        if height != 0. height_shift = height * selected_unit /dataobject.boxlen end
        center  =  center / (dataobject.boxlen * selected_unit )
    end

    if verbose
        if center != [0., 0., 0.]
            print("center: $(round.(center,digits=7)) ")
            center_val1, center_unit1= humanize(center[1] * dataobject.boxlen , dataobject.scale, 3, "length")
            center_val2, center_unit2= humanize(center[2] * dataobject.boxlen , dataobject.scale, 3, "length")
            center_val3, center_unit3= humanize(center[3] * dataobject.boxlen , dataobject.scale, 3, "length")
            if center_val1 == 0. center_unit1 = xmax_unit end
            if center_val2 == 0. center_unit2 = xmax_unit end
            if center_val3 == 0. center_unit3 = xmax_unit end
            print("==> [$center_val1 [$center_unit1] :: $center_val2 [$center_unit2] :: $center_val3 [$center_unit3]]\n")
            println()
        end

        #println("domain \t \t \t \thuman-readable units")
        println("domain:")

        print("xmin::xmax: $(round(xmin,digits=7)) :: $(round(xmax,digits=7))  \t")
        xmin_val, xmin_unit  = humanize(xmin * dataobject.boxlen , dataobject.scale, 3, "length")
        xmax_val, xmax_unit  = humanize(xmax * dataobject.boxlen, dataobject.scale, 3, "length")
        if xmin_val == 0. xmin_unit = xmax_unit end
        print("==> $xmin_val [$xmin_unit] :: $xmax_val [$xmax_unit]\n")

        print("ymin::ymax: $(round(ymin,digits=7)) :: $(round(ymax,digits=7))  \t")
        ymin_val, ymin_unit  = humanize(ymin * dataobject.boxlen, dataobject.scale, 3, "length")
        ymax_val, ymax_unit  = humanize(ymax * dataobject.boxlen, dataobject.scale, 3, "length")
        if ymin_val == 0. ymin_unit = xmax_unit end
        print("==> $ymin_val [$ymin_unit] :: $ymax_val [$ymax_unit]\n")

        print("zmin::zmax: $(round(zmin,digits=7)) :: $(round(zmax,digits=7))  \t")
        zmin_val, zmin_unit  = humanize(zmin * dataobject.boxlen, dataobject.scale, 3, "length")
        zmax_val, zmax_unit  = humanize(zmax * dataobject.boxlen, dataobject.scale, 3, "length")
        if zmin_val == 0. zmin_unit = xmax_unit end
        print("==> $zmin_val [$zmin_unit] :: $zmax_val [$zmax_unit]\n")
        println()

        R_val, R_unit  = humanize(radius_shift * dataobject.boxlen, dataobject.scale, 3, "length")
        println("Radius: $R_val [$R_unit]")
        if height != 0.
            h_val, h_unit  = humanize(height_shift * dataobject.boxlen, dataobject.scale, 3, "length")
            println("Height: $h_val [$h_unit]")
        end
    end

    # if xmin < 0. || ymin < 0. || zmin < 0. || xmax > 1. || ymax > 1. || zmax > 1.
    #     error("[Mera]: Given range(s) outside of box!")
    # end
    ranges = [xmin, xmax, ymin, ymax, zmin, zmax]

    if height != 0.
        return ranges, cx_shift, cy_shift, cz_shift, radius_shift, height_shift
    else
        return ranges, cx_shift, cy_shift, cz_shift, radius_shift
    end
end

"""
cylindrical shell ranges
convert given ranges + radius and print overview on screen
used for shellregions...
"""
function prep_cylindrical_shellranges(    dataobject::InfoType,
                            center::Array{<:Any,1},
                            radius_in::Real,
                            radius_out::Real,
                            height::Real,
                            range_unit::Symbol,
                            verbose::Bool)

    center = prepboxcenter(dataobject, range_unit, center)

    selected_unit = 1.
    if range_unit == :standard
        xmin = -radius_out + center[1]
        xmax =  radius_out + center[1]
        ymin = -radius_out + center[2]
        ymax =  radius_out + center[2]
        zmin = - height + center[3]
        zmax =   height + center[3]

        # given center relative to the data range in units: cell centers
        #todo

    else
        selected_unit = getunit(dataobject, range_unit)
        xmin = (-radius_out + center[1]) * selected_unit /dataobject.boxlen
        xmax = ( radius_out + center[1]) * selected_unit /dataobject.boxlen
        ymin = (-radius_out + center[2]) * selected_unit /dataobject.boxlen
        ymax = ( radius_out + center[2]) * selected_unit /dataobject.boxlen
        zmin = ( -height + center[3]) * selected_unit /dataobject.boxlen
        zmax = (  height + center[3]) * selected_unit /dataobject.boxlen

        # given center relative to the data range in units: cell centers
        cx_shift = center[1] * selected_unit /dataobject.boxlen
        cy_shift = center[2] * selected_unit /dataobject.boxlen
        cz_shift = center[3] * selected_unit /dataobject.boxlen
        radius_in_shift = radius_in * selected_unit /dataobject.boxlen
        radius_out_shift = radius_out * selected_unit /dataobject.boxlen
        height_shift = height * selected_unit /dataobject.boxlen
        center  =  center / (dataobject.boxlen * selected_unit )
    end

    if verbose
        if center != [0., 0., 0.]
            print("center: $(round.(center,digits=7)) ")
            center_val1, center_unit1= humanize(center[1] * dataobject.boxlen , dataobject.scale, 3, "length")
            center_val2, center_unit2= humanize(center[2] * dataobject.boxlen , dataobject.scale, 3, "length")
            center_val3, center_unit3= humanize(center[3] * dataobject.boxlen , dataobject.scale, 3, "length")
            if center_val1 == 0. center_unit1 = xmax_unit end
            if center_val2 == 0. center_unit2 = xmax_unit end
            if center_val3 == 0. center_unit3 = xmax_unit end
            print("==> [$center_val1 [$center_unit1] :: $center_val2 [$center_unit2] :: $center_val3 [$center_unit3]]\n")
            println()
        end

        #println("domain \t \t \t \thuman-readable units")
        println("domain:")

        print("xmin::xmax: $(round(xmin,digits=7)) :: $(round(xmax,digits=7))  \t")
        xmin_val, xmin_unit  = humanize(xmin * dataobject.boxlen , dataobject.scale, 3, "length")
        xmax_val, xmax_unit  = humanize(xmax * dataobject.boxlen, dataobject.scale, 3, "length")
        if xmin_val == 0. xmin_unit = xmax_unit end
        print("==> $xmin_val [$xmin_unit] :: $xmax_val [$xmax_unit]\n")

        print("ymin::ymax: $(round(ymin,digits=7)) :: $(round(ymax,digits=7))  \t")
        ymin_val, ymin_unit  = humanize(ymin * dataobject.boxlen, dataobject.scale, 3, "length")
        ymax_val, ymax_unit  = humanize(ymax * dataobject.boxlen, dataobject.scale, 3, "length")
        if ymin_val == 0. ymin_unit = xmax_unit end
        print("==> $ymin_val [$ymin_unit] :: $ymax_val [$ymax_unit]\n")

        print("zmin::zmax: $(round(zmin,digits=7)) :: $(round(zmax,digits=7))  \t")
        zmin_val, zmin_unit  = humanize(zmin * dataobject.boxlen, dataobject.scale, 3, "length")
        zmax_val, zmax_unit  = humanize(zmax * dataobject.boxlen, dataobject.scale, 3, "length")
        if zmin_val == 0. zmin_unit = xmax_unit end
        print("==> $zmin_val [$zmin_unit] :: $zmax_val [$zmax_unit]\n")
        println()

        Rin_val, Rin_unit  = humanize(radius_in_shift * dataobject.boxlen, dataobject.scale, 3, "length")
        println("Inner radius: $Rin_val [$Rin_unit]")
        Rout_val, Rout_unit  = humanize(radius_out_shift * dataobject.boxlen, dataobject.scale, 3, "length")
        println("Outer radius: $Rout_val [$Rout_unit]")
        radius_diff_shift = radius_out_shift - radius_in_shift
        Rdiff_val, Rdiff_unit  = humanize(radius_diff_shift * dataobject.boxlen, dataobject.scale, 3, "length")
        println("Radius diff: $Rdiff_val [$Rdiff_unit]")
        h_val, h_unit  = humanize(height_shift * dataobject.boxlen, dataobject.scale, 3, "length")
        println("Height: $h_val [$h_unit]")
    end

    # if xmin < 0. || ymin < 0. || zmin < 0. || xmax > 1. || ymax > 1. || zmax > 1.
    #     error("[Mera]: Given range(s) outside of box!")
    # end
    ranges = [xmin, xmax, ymin, ymax, zmin, zmax]
    return ranges, cx_shift, cy_shift, cz_shift, radius_in_shift, radius_out_shift, height_shift
end

"""
spherical shell ranges
convert given ranges + radius and print overview on screen
used for shellregions...
"""
function prep_spherical_shellranges(    dataobject::InfoType,
                            center::Array{<:Any,1},
                            radius_in::Real,
                            radius_out::Real,
                            range_unit::Symbol,
                            verbose::Bool)

    center = prepboxcenter(dataobject, range_unit, center)

    selected_unit = 1.
    if range_unit == :standard
        xmin = -radius_out + center[1]
        xmax =  radius_out + center[1]
        ymin = -radius_out + center[2]
        ymax =  radius_out + center[2]
        zmin = - radius_out + center[3]
        zmax =   radius_out + center[3]

        # given center relative to the data range in units: cell centers
        #todo

    else
        selected_unit = getunit(dataobject, range_unit)
        xmin = (-radius_out + center[1]) * selected_unit /dataobject.boxlen
        xmax = ( radius_out + center[1]) * selected_unit /dataobject.boxlen
        ymin = (-radius_out + center[2]) * selected_unit /dataobject.boxlen
        ymax = ( radius_out + center[2]) * selected_unit /dataobject.boxlen
        zmin = ( -radius_out + center[3]) * selected_unit /dataobject.boxlen
        zmax = (  radius_out + center[3]) * selected_unit /dataobject.boxlen

        # given center relative to the data range in units: cell centers
        cx_shift = center[1] * selected_unit /dataobject.boxlen
        cy_shift = center[2] * selected_unit /dataobject.boxlen
        cz_shift = center[3] * selected_unit /dataobject.boxlen
        radius_in_shift = radius_in * selected_unit /dataobject.boxlen
        radius_out_shift = radius_out * selected_unit /dataobject.boxlen

        center  =  center / (dataobject.boxlen * selected_unit )
    end

    if verbose
        if center != [0., 0., 0.]
            print("center: $(round.(center,digits=7)) ")
            center_val1, center_unit1= humanize(center[1] * dataobject.boxlen , dataobject.scale, 3, "length")
            center_val2, center_unit2= humanize(center[2] * dataobject.boxlen , dataobject.scale, 3, "length")
            center_val3, center_unit3= humanize(center[3] * dataobject.boxlen , dataobject.scale, 3, "length")
            if center_val1 == 0. center_unit1 = xmax_unit end
            if center_val2 == 0. center_unit2 = xmax_unit end
            if center_val3 == 0. center_unit3 = xmax_unit end
            print("==> [$center_val1 [$center_unit1] :: $center_val2 [$center_unit2] :: $center_val3 [$center_unit3]]\n")
            println()
        end

        #println("domain \t \t \t \thuman-readable units")
        println("domain:")

        print("xmin::xmax: $(round(xmin,digits=7)) :: $(round(xmax,digits=7))  \t")
        xmin_val, xmin_unit  = humanize(xmin * dataobject.boxlen , dataobject.scale, 3, "length")
        xmax_val, xmax_unit  = humanize(xmax * dataobject.boxlen, dataobject.scale, 3, "length")
        if xmin_val == 0. xmin_unit = xmax_unit end
        print("==> $xmin_val [$xmin_unit] :: $xmax_val [$xmax_unit]\n")

        print("ymin::ymax: $(round(ymin,digits=7)) :: $(round(ymax,digits=7))  \t")
        ymin_val, ymin_unit  = humanize(ymin * dataobject.boxlen, dataobject.scale, 3, "length")
        ymax_val, ymax_unit  = humanize(ymax * dataobject.boxlen, dataobject.scale, 3, "length")
        if ymin_val == 0. ymin_unit = xmax_unit end
        print("==> $ymin_val [$ymin_unit] :: $ymax_val [$ymax_unit]\n")

        print("zmin::zmax: $(round(zmin,digits=7)) :: $(round(zmax,digits=7))  \t")
        zmin_val, zmin_unit  = humanize(zmin * dataobject.boxlen, dataobject.scale, 3, "length")
        zmax_val, zmax_unit  = humanize(zmax * dataobject.boxlen, dataobject.scale, 3, "length")
        if zmin_val == 0. zmin_unit = xmax_unit end
        print("==> $zmin_val [$zmin_unit] :: $zmax_val [$zmax_unit]\n")
        println()

        Rin_val, Rin_unit  = humanize(radius_in_shift * dataobject.boxlen, dataobject.scale, 3, "length")
        println("Inner radius: $Rin_val [$Rin_unit]")
        Rout_val, Rout_unit  = humanize(radius_out_shift * dataobject.boxlen, dataobject.scale, 3, "length")
        println("Outer radius: $Rout_val [$Rout_unit]")
        radius_diff_shift = radius_out_shift - radius_in_shift
        Rdiff_val, Rdiff_unit  = humanize(radius_diff_shift * dataobject.boxlen, dataobject.scale, 3, "length")
        println("Radius diff: $Rdiff_val [$Rdiff_unit]")
    end

    # if xmin < 0. || ymin < 0. || zmin < 0. || xmax > 1. || ymax > 1. || zmax > 1.
    #     error("[Mera]: Given range(s) outside of box!")
    # end
    ranges = [xmin, xmax, ymin, ymax, zmin, zmax]

    return ranges, cx_shift, cy_shift, cz_shift, radius_in_shift, radius_out_shift
end

function prepboxcenter(dataobject::InfoType, range_unit::Symbol, center::Array{<:Any,1})

    selected_unit = 1.
    if range_unit != :standard
        selected_unit = getunit(dataobject, range_unit)
    end

    # check for :bc, :boxcenter
    Ncenter = length(center)
    if Ncenter  == 1
        if in(:bc, center) || in(:boxcenter, center)
            if range_unit == :standard
                bc = 0.5
            else
                bc = dataobject.boxlen * 0.5 * selected_unit # use range_unit
            end
            centerm = [bc, bc, bc]

        end
    else
        centerm = zeros(Float64, length(center))
        for i = 1:Ncenter
            if center[i] == :bc || center[i] == :boxcenter
                if range_unit == :standard
                    bc =0.5
                else
                    bc = dataobject.boxlen * 0.5 * selected_unit # use range_unit
                end
                centerm[i] = bc
            else
                centerm[i] = center[i]
            end
        end
    end

    return  centerm

end

function prepboxcenter(dataobject::InfoType, selected_unit::Real, centerm::Any)

    Nc = length(centerm)
    center = zeros(Float64, Nc)

    if centerm == :bc || centerm == :boxcenter
        center = dataobject.boxlen * 0.5 * selected_unit
    else
        center = centerm ./ dataobject.boxlen .* selected_unit
    end

    return center
end

function prepdatacenter(dataobject::InfoType, center::Array{<:Any,1}, range_unit::Symbol, data_centerm::Array{<:Any,1}, data_center_unit::Symbol)

    center = center_in_standardnotation(dataobject, center, range_unit)
    selected_unit = 1.
    if data_center_unit != :standard
        selected_unit = getunit(dataobject, data_center_unit)
    end

    Ndc = length(center)
    data_center = zeros(Float64, Ndc)
    if  Ndc == 3
        if data_centerm[1] === missing
            data_center[1] = center[1]
        else
            data_center[1] = prepboxcenter(dataobject, selected_unit, data_centerm[1])
        end

        if data_centerm[2] === missing
            data_center[2] = center[2]
        else
            data_center[2] = prepboxcenter(dataobject, selected_unit, data_centerm[2])
        end

        if data_centerm[3] === missing
            data_center[3] = center[3]
        else
            data_center[3] = prepboxcenter(dataobject, selected_unit, data_centerm[3])
        end
    elseif Ndc == 1
        data_center[1] = data_center[2] = data_center[3] = prepboxcenter(dataobject, selected_unit, data_centerm[1])
    end

    return data_center
end


text:
struct Qda
    Wi::AbstractVector
    ct::Array{Float64}
    wprior::AbstractVector
    lev::AbstractVector
    ni::AbstractVector
end

"""
    qda(X, y; prior = "unif")
Quadratic discriminant analysis  (QDA).
* `X` : X-data.
* `y` : y-data (class membership).
* `prior` : Type of prior probabilities for class membership
    (`unif`: uniform; `prop`: proportional).
"""
function qda(X, y; prior = "unif")
    X = ensure_mat(X)
    z = aggstat(X, y; fun = mean)
    ct = z.res
    lev = z.lev
    nlev = length(lev)
    ni = z.ni
    if isequal(prior, "unif")
        wprior = ones(nlev) / nlev
    elseif isequal(prior, "prop")
        wprior = mweights(ni)
    end
    res = matW(X, y)
    Qda(res.Wi, ct, wprior, lev, ni)
end

function predict(object::Qda, X)
    X = ensure_mat(X)
    m = size(X, 1)
    lev = object.lev
    nlev = length(lev)
    ds = similar(X, m, nlev)
    ni = object.ni
    for i = 1:nlev
        if ni[i] == 1
            S = object.Wi[i]
        else
            S = object.Wi[i] * ni[i] / (ni[i] - 1)
        end
        fm = dmnorm(; mu = object.ct[i, :], S = S)
        ds[:, i] .= vec(predict(fm, X).pred)
    end
    A = object.wprior' .* ds
    v = sum(A, dims = 2)
    posterior = scale(A', v)' # This could be replaced by code similar as in scale!
    z =  mapslices(argmax, posterior; dims = 2)  # if equal, argmax takes the first
    pred = reshape(replacebylev(z, object.lev), m, 1)
    (pred = pred, ds = ds, posterior = posterior)
end


text:
using Pkg
Pkg.activate(pwd())

using Merly

@page "/" HTTP.Response(200, "")

@page "/user/:id" HTTP.Response(200, string(request.params["id"]))

@route POST "/user" begin
    HTTP.Response(200, "")
end

start(host = "0.0.0.0", port = 3000, verbose = true)

max_stars_repo_name: mattiapenati/web-frameworks

max_issues_repo_name: mattiapenati/web-frameworks

max_forks_repo_name: mattiapenati/web-frameworks


text:
"""
    community_detection_nback(g, k)

Community detection using the spectral properties of
the non-backtracking matrix of graph `g` (see [Krzakala et al.](http://www.pnas.org/content/110/52/20935.short)).
`k` is the number of communities to detect.

See also [`community_detection_bethe`](@ref) for a related community ddetection
algorithm.

Returns a vector with the vertex assignments in the communities.
"""
function community_detection_nback(g::AGraph, k::Integer)
    #TODO insert check on connected_components
    ϕ = real(nonbacktrack_embedding(g, k))
    if k == 2
        c = community_detection_threshold(g, ϕ[1,:])
    else
        c = kmeans(ϕ, k).assignments
    end
    return c
end

function community_detection_threshold(g::AGraphOrDiGraph, coords::AbstractArray)
    # TODO use a more intelligent method to set the threshold
    # 0 based thresholds are highly sensitive to errors.
    c = ones(Int, nv(g))
    # idx = sortperm(λ, lt=(x,y)-> abs(x) > abs(y))[2:k] #the second eigenvector is the relevant one
    for i=1:nv(g)
        c[i] = coords[i] > 0 ?  1 : 2
    end
    return c
end

"""
    nonbacktrack_embedding(g::AGraph, k::Int)

Spectral embedding of the non-backtracking matrix of `g`
(see [Krzakala et al.](http://www.pnas.org/content/110/52/20935.short)).

    `g`: imput Graph
    `k`: number of dimensions in which to embed

Returns  a matrix `ϕ` where `ϕ[:,i]` are the coordinates for vertex `i`.

Note:  does not explicitly construct the [`nonbacktracking_matrix`](@ref).
See [`nonbacktracking_matrix`](@ref) for details.
"""
function nonbacktrack_embedding(g::AGraph, k::Int)
    B, edgeid = nonbacktracking_matrix(g)
    λ, eigv = Arpack.eigs(B, nev=k+1, v0=ones(Float64, size(B,1)))
    ϕ = zeros(ComplexF64, k-1, nv(g))
    # TODO decide what to do with the stationary distribution ϕ[:,1]
    # this code just throws it away in favor of eigv[:,2:k+1].
    # we might also use the degree distribution to scale these vectors as is
    # common with the laplacian/adjacency methods.
    E = edgetype(g)
    for n=1:k-1
        v = eigv[:,n+1]
        for i=1:nv(g)
            for e in in_edges(g, i)
                j = src(e)
                u = edgeid[e]
                ϕ[n,i] += v[u]
            end
        end
    end
    return ϕ
end

"""
    community_detection_bethe(g::AGraph, k=-1; kmax=15)

Community detection using the spectral properties of
the Bethe Hessian matrix associated to `g` (see [Saade et al.](http://papers.nips.cc/paper/5520-spectral-clustering-of-graphs-with-the-bethe-hessian)).
`k` is the number of communities to detect. If omitted or if `k < 1` the
optimal number of communities will be automatically selected.
In this case the maximum number of detectable communities is given by `kmax`.
Returns a vector containing the vertex assignments.
"""
function community_detection_bethe(g::AGraph, k::Int=-1; kmax::Int=15)
    A = adjacency_matrix(g)
    D = Diagonal(degree(g))
    r = (sum(degree(g)) / nv(g))^0.5

    Hr = (r^2-1)*Diagonal(ones(nv(g))) - r*A + D;
    # Hmr = (r^2-1)*eye(nv(g))+r*A+D;
    k >= 1 && (kmax = k)
    λ, eigv = Arpack.eigs(Hr, which=:SR, nev=min(kmax, nv(g)))
    q = findlast(x -> x<0, λ)
    k > q && @warn("Using eigenvectors with positive eigenvalues,
                    some communities could be meaningless. Try to reduce `k`.")
    k < 1 && (k = q)
    k < 1 && return fill(1, nv(g))
    m = Matrix(eigv[:,2:k]') # convert to Matrix since
                             # kmeans doesn't support Adjoint yet
    labels = kmeans(m, k).assignments
    return labels
end

"""
    label_propagation(g; maxiter=1000)

Community detection using the label propagation algorithm (see [Raghavan et al.](http://arxiv.org/abs/0709.2938)).
`g` is the input Graph, `maxiter` is the  maximum number of iterations.
Returns a vertex assignments and the convergence history
"""
function label_propagation(g::AGraphOrDiGraph; maxiter=1000)
    n = nv(g)
    label = collect(1:n)
    active_nodes = BitSet(vertices(g))
    c = NeighComm(collect(1:n), fill(-1,n), 1)
    convergence_hist = Vector{Int}()
    random_order = Vector{Int}(undef, n)
    i = 0
    while !isempty(active_nodes) && i < maxiter
        num_active = length(active_nodes)
        push!(convergence_hist, num_active)
        i += 1

        # processing nodes in random order
        for (j,node) in enumerate(active_nodes)
            random_order[j] = node
        end
        range_shuffle!(1:num_active, random_order)
        @inbounds for j=1:num_active
            u = random_order[j]
            old_comm = label[u]
            label[u] = vote!(g, label, c, u)
            if old_comm != label[u]
                for v in out_neighbors(g, u)
                    push!(active_nodes, v)
                end
            else
                delete!(active_nodes, u)
            end
        end
    end
    fill!(c.neigh_cnt, 0)
    renumber_labels!(label, c.neigh_cnt)
    label, convergence_hist
end

"""Type to record neighbor labels and their counts."""
mutable struct NeighComm
  neigh_pos::Vector{Int}
  neigh_cnt::Vector{Int}
  neigh_last::Int
end

"""Fast shuffle Array `a` in UnitRange `r` inplace."""
function range_shuffle!(r::UnitRange, a::AbstractVector)
    (r.start > 0 && r.stop <= length(a)) || error("out of bounds")
    @inbounds for i=length(r):-1:2
        j = rand(1:i)
        ii = i + r.start - 1
        jj = j + r.start - 1
        a[ii],a[jj] = a[jj],a[ii]
    end
end

"""Return the most frequency label."""
function vote!(g::AGraphOrDiGraph, m::Vector{Int}, c::NeighComm, u::Int)
    @inbounds for i=1:c.neigh_last-1
        c.neigh_cnt[c.neigh_pos[i]] = -1
    end
    c.neigh_last = 1
    c.neigh_pos[1] = m[u]
    c.neigh_cnt[c.neigh_pos[1]] = 0
    c.neigh_last = 2
    max_cnt = 0
    for neigh in out_neighbors(g,u)
        neigh_comm = m[neigh]
        if c.neigh_cnt[neigh_comm] < 0
            c.neigh_cnt[neigh_comm] = 0
            c.neigh_pos[c.neigh_last] = neigh_comm
            c.neigh_last += 1
        end
        c.neigh_cnt[neigh_comm] += 1
        if c.neigh_cnt[neigh_comm] > max_cnt
          max_cnt = c.neigh_cnt[neigh_comm]
        end
    end
    # ties breaking randomly
    range_shuffle!(1:c.neigh_last-1, c.neigh_pos)
    for lbl in c.neigh_pos
      if c.neigh_cnt[lbl] == max_cnt
        return lbl
      end
    end
end

function renumber_labels!(membership::Vector{Int}, label_counters::Vector{Int})
    N = length(membership)
    (maximum(membership) > N || minimum(membership) < 1) && error("Label must between 1 and |V|")
    j = 1
    @inbounds for i=1:length(membership)
        k = membership[i]
        if k >= 1
            if label_counters[k] == 0
                # We have seen this label for the first time
                label_counters[k] = j
                k = j
                j += 1
            else
                k = label_counters[k]
            end
        end
        membership[i] = k
    end
end


text:
"""
    post_godbolt(c.godbolt)

Post code snippet in [Compiler Explore](godbolt.org).

**WARNING**: When the code is posted to godbolt.org there is no way to delete
the code as of writing.

# Examples

```julia
c = @sc_llvm sum((1, 2, 3))
ce = post_godbolt(c.godbolt)         # post the code
string(ce)                           # get the URL
ce()                                 # post it and then open the URL
```
"""
post_godbolt

# ce - Compiler Explore
# cec - Compiler Explore Client

abstract type GodboltBase end

struct Godbolt <: GodboltBase
    code::Any
    cache::Dict{String,Any}
end

Godbolt(code) = Godbolt(code, Dict{String,Any}())

Base.propertynames(::Godbolt) = (:mca,)
Base.getproperty(ce::Godbolt, compiler::Symbol) = getproperty(ce, string(compiler))
Base.getproperty(ce::Godbolt, compiler::AbstractString) =
    get!(Fields(ce).cache, compiler) do
        _getproperty(ce, compiler)
    end::GodboltClient

function _getproperty(ce::Godbolt, compiler::AbstractString)
    compiler == "_default" && return default_client(ce)
    @unpack code = Fields(ce)
    compiler = get(Dict("mca" => "llvm-mcatrunk"), compiler, compiler)
    language = godbolt_language(code)
    return GodboltClient(; code, language, compiler)
end

godbolt_language(::CodeLLVM) = "llvm"
godbolt_language(::CodeNative) = "assembly"

function default_client(ce::Godbolt)
    @unpack code = Fields(ce)
    if code isa CodeNative
        return GodboltClient(; code, language = "analysis", compiler = "llvm-mcatrunk")
    elseif code isa CodeLLVM
        return ce.llctrunk
    else
        error("No default client for ", code)
    end
end

Base.@kwdef struct GodboltClient <: GodboltBase
    code::Any

    # `id` of https://godbolt.org/api/languages
    language::String

    # `id` of https://godbolt.org/api/compilers
    compiler::String

    options::Cmd = ``

    url::Base.RefValue{Union{Nothing,String}} = Ref{Union{Nothing,String}}(nothing)
end
# https://github.com/compiler-explorer/compiler-explorer/blob/master/docs/API.md

function clientstate(cec::GodboltClient)
    # cec - Compiler Explore Client
    @unpack code, language, compiler, options = cec
    options_str = join(options, " ")  # TODO: quote
    return Dict(
        "sessions" => [Dict(
            "language" => language,
            "source" => string(code),
            "compilers" => [Dict("id" => compiler, "options" => options_str)],
        )],
    )
end

function post_godbolt(ce::Godbolt)
    post_godbolt(GodboltClient(ce))
    return ce
end

function post_godbolt(cec::GodboltClient)
    @unpack url = Fields(cec)
    url[] === nothing || return cec
    url[] = godbolt_url(cec)
    return cec
end

get_url(cec::GodboltClient) = Fields(cec).url[]
Base.string(cec::GodboltClient) = something(get_url(cec), "?? not posted ??")
# TODO: maybe don't overload `string`?

function godbolt_url(cec::GodboltClient)
    url = godbolt_base64url(cec)
    # Use URL shortener for URL longer than 2000 characters
    # https://stackoverflow.com/a/417184
    if length(url) > 2000
        url = godbolt_shorturl(cec)
    end
    return url
end

function godbolt_base64url(cec::GodboltClient)
    s = base64encode(JSON.print, clientstate(cec))
    return "https://godbolt.org/clientstate/$s"
end

function godbolt_shorturl(cec::GodboltClient)
    response = HTTP.post(
        "https://godbolt.org/shortener",
        ["Content-Type" => "application/json", "Accept" => "application/json"],
        JSON.json(clientstate(cec)),
    )
    msg = JSON.parse(String(response.body))
    return msg["url"]
end

function (cec::GodboltClient)()
    url = get_url(post_godbolt(cec))::AbstractString
    DefaultApplication.open(url; wait = false)
    return
end

GodboltClient(cec::GodboltClient) = cec
GodboltClient(ce::Godbolt) = ce._default

(ce::GodboltBase)() = GodboltClient(ce)()
Base.string(ce::GodboltBase) = string(GodboltClient(ce))
Base.print(io::IO, ce::GodboltBase) = print(io, string(ce))

function Base.summary(io::IO, ce::GodboltBase)
    print(io, nameof(typeof(ce)), " (interface for opening URL at godbolt.org)")
end

function Base.show(io::IO, ::MIME"text/plain", ce::GodboltBase)
    summary(io, ce)
    println(io)
    print(io, GodboltClient(ce))
end


text:
const MeanOrKernelOrLikelihood = Union{Mean,Kernel,Likelihood}
const CompositeMeanOrKernel = Union{CompositeMean,CompositeKernel}

##########¤
# Display #
###########

function Base.show(io::IO, obj::MeanOrKernelOrLikelihood, depth::Int = 0)
    for _ in 1:depth
        print(io, "  ")
    end
    print(io, "Type: ", typeof(obj), ", Params: ")
    show(io, get_params(obj))
end

function Base.show(io::IO, obj::CompositeMeanOrKernel, depth::Int = 0)
    for _ in 1:depth
        print(io, "  ")
    end
    println(io, "Type: ", typeof(obj))
    for c in components(obj)
        show(io, c, depth + 1)
    end
end

##############################
# Parameter name definitions #
##############################

# This generates names like [:ll_1, :ll_2, ...] for parameter vectors
get_param_names(n::Int, prefix::Symbol) = [Symbol(prefix, :_, i) for i in 1:n]
get_param_names(v::Vector, prefix::Symbol) = get_param_names(length(v), prefix)

"""
    composite_param_names(objects, prefix)

Call `get_param_names` on each element of `objects` and prefix the returned name of the
element at index `i` with `prefix * i * '_'`.

# Examples
```jldoctest
julia> get_param_names(ProdKernel(Mat(1/2, 1/2, 1/2), SEArd([0.0, 1.0], 0.0)))
5-element Array{Symbol,1}:
 :pk1_ll
 :pk1_lσ
 :pk2_ll_1
 :pk2_ll_2
 :pk2_lσ
```
"""
function composite_param_names(objects, prefix)
    p = Symbol[]
    for (i, obj) in enumerate(objects)
        append!(p, [Symbol(prefix, i, :_, sym) for sym in get_param_names(obj)])
    end
    p
end

##############
# Parameters #
##############

function get_params(obj::CompositeMeanOrKernel)
    p = Array{Float64}(undef, 0)
    for c in components(obj)
        append!(p, get_params(c))
    end
    p
end

num_params(obj::CompositeMeanOrKernel) = sum(num_params, components(obj))

function set_params!(obj::CompositeMeanOrKernel, hyp::VecF64)
    length(hyp) == num_params(obj) ||
        throw(ArgumentError("$(typeof(obj)) object requires $(num_params(obj)) hyperparameters"))
    i = 1
    @inbounds for c in components(obj)
        j = i + num_params(c)
        set_params!(c, view(hyp, i:(j - 1)))
        i = j
    end
end

##########
# Priors #
##########

get_priors(obj::MeanOrKernelOrLikelihood) = obj.priors

function get_priors(obj::CompositeMeanOrKernel)
    p = []
    for c in components(obj)
        append!(p, get_priors(c))
    end
    p
end

function set_priors!(obj::MeanOrKernelOrLikelihood, priors::Array)
    length(priors) == num_params(obj) ||
        throw(ArgumentError("$(typeof(obj)) object requires $(num_params(obj)) priors"))

    obj.priors = priors
end

function set_priors!(obj::CompositeMeanOrKernel, priors::Array)
    length(priors) == num_params(obj) ||
        throw(ArgumentError("$(typeof(obj)) object requires $(num_params(obj)) priors"))
    i = 1
    @inbounds for c in components(obj)
        j = i + num_params(c)
        set_priors!(c, view(priors, i:(j - 1)))
        i = j
    end
end

function prior_logpdf(obj::MeanOrKernelOrLikelihood)
    num_params(obj) == 0 && return 0.0

    priors = get_priors(obj)
    isempty(priors) && return 0.0

    sum(logpdf(prior, param) for (prior, param) in zip(priors, get_params(obj)))
end

function prior_gradlogpdf(obj::MeanOrKernelOrLikelihood)
    num_params(obj) == 0 && return Float64[]

    priors = get_priors(obj)
    isempty(priors) && return zeros(num_params(obj))

    [gradlogpdf(prior, param) for (prior, param) in zip(priors, get_params(obj))]
end


text:
## complex-real case

##################################################

# after a quick fuse
# we must passthrough the diagonal rotator through the rest of Q
function absorb_Ut(QF::QFactorization{T, S, VV}, RF, storage, ctr) where {T, S <: Complex, VV}

    Ut = storage.VU[1]'
    fuse!(Ut, QF)

    return nothing
end

function passthrough_triu(QF::QFactorization{T, S, VV}, RF, storage, ctr, dir::Val{:right}) where {T, S <: Complex, VV}

    U = storage.VU[1]
    i = idx(U)

    flag = false

    if  i < ctr.tr

        flag = simple_passthrough!(RF, U)

    end

    if i > ctr.tr || !flag

        U = passthrough!(RF, U)
        storage.VU[1] = U

    end

     ctr.tr -= 1

    return nothing
end

## When Ct and B are identical, we can update just one and leave U,V alone
function simple_passthrough!(RF::RFactorizationRankOne{T, S}, U) where {T, S <: Complex}

    i = idx(U)
    N = length(RF)

    _ = passthrough!(RF.B, U)

    for k in 0:1

        a,b = vals(RF.B[i+k])
        ii = N + 1 - (i + k)
        RF.Ct[ii] = Rotator(conj(a), -b, i + k) # adjoint

    end

    true
end

function passthrough_Q(QF::QFactorization{T, S, VV}, RF, storage, ctr,  dir::Val{:right}) where {T, S <: Complex, VV}

    U = storage.VU[1]
    i = idx(U)

    if i < ctr.stop_index

        U = passthrough!(QF, U)

        storage.VU[1] = U

        return false

    else

        # handle details of knitting in
        D = QF.D
        U = passthrough!(D, U)
        fuse!(QF, U) # handles D bit

        return true

    end
end

##################################################
## Deflation
##
# deflate a term
# deflation for ComplexReal is different, as
# we replace Qi with I and move diagonal part into D
function deflate(QF::AbstractQFactorization{T, S}, k, cgr) where {T, S <: Complex}

    # when we deflate here we want to leave Q[k] = I and

    # move Dk matrix over to merge with D
    # Qi       Qi            Qi           Qi
    #   Qj   ->   I Dj  -->    I Dj   -->   I Dj
    #     Qk        QK           Qk Dj        Qk Dj    fuse Dj's with D
    #       QL        QL           QL           Ql Dj

    # then the Dk's are collected into D

    alpha, s = vals(QF.Q[k])
    i = idx(QF.Q[k])

    ## Make Q[k] an identity rotator
    QF.Q[k] = Rotator(one(Complex{T}), zero(T), i)

    # absorb Di into D
    Di =  DiagonalRotator(alpha, i)
    passthrough_phase!(Di, QF)

end


text:
using Test
using ForwardDiff: Dual, gradient, partials

using DiffEqBase: ODE_DEFAULT_NORM
const internalnorm = ODE_DEFAULT_NORM

val = rand(10)
par = rand(10)
u = Dual.(val, par)
reference = internalnorm(val, 1)
dual_real = internalnorm(u, 1)
dual_dual = internalnorm(u, u[1])
@test reference === dual_real
@test reference == dual_dual
@test partials(dual_dual, 1) ≈ gradient(x->internalnorm(x, x[1]), val)'par