text:
cd("File Path")
wdir     = pwd()
# # Write API Key to environment

ENV["ALPHA_VANTAGE_API_KEY"] = "YOUR KEY HERE"

# Core Packages needed
using AlphaVantage
using DataFrames
using DataFramesMeta
using Dates
using CSV

# For making nice plots
using Plots
using LaTeXStrings
using Statistics
using ColorSchemes
using TimeSeries

include("Project_Functions.jl")

################################################################################
                                    # Main Part #
################################################################################

## Prerequisites necessary to use this code

# (1) Create an AlphaVantage account for free at  https://www.alphavantage.co/support/#api-key and request an API
#
# (2) Construct a spreadsheet (.csv) with  all sales/purchase information in it:
#
# The structure of the spreadsheet should be like this:
#
# Date       | Quantity_Stock1 | Price_Stock1    | Payments_Stock1| Cost_Stock1 |
# 20xx-xx-xx |      Float      |      Float      |      Float     |      Float  |
#
# "Stock1" should be the string-handle used to pull the security-series from AlphaVantage.
# E.g. if you want to pull Apple's stock price, use "Quantity_AAPL", "Price_AAPL", etc
#
# For symbol search define SEARCHWORD and find ticker:
# https://www.alphavantage.co/query?function=SYMBOL_SEARCH&keywords=SEARCHWORD&apikey=6M3D8EV1HA6WKQFI
#
# Adding new infomation just requires adding  additional rows. Adding a new security
# means adding four new columns

## Procedures

# Set start/end date of of data sample:

start_date          = "2019-01-04";  # Trading information available starts on "yyyy-mm-dd".
end_date            = [];  # if "end_date" is empty, the last available value will be set.

# (1) Load Trading Information avialable: - Which titles and in which quantity have been traded on this account over the entire horizon?

raw_trades, trading_start, header_types, full_name_list = get_trading_info(wdir);

# (2) Use available trading information and match it to desired sample horizon.
#     Only load security data of securities traded in respective interval.
#
#     This requires the following steps:

# (2a) - Which titles and which quantity did I trade within selected start/end_date (today)?
#      - Find securities in "full_name_list" that are not in trading sample within selected time span.
#      - Create a subset that only pulls, analyzes and prints securities available in in
#        the selected time span.
# (2b) - Determine initial/already exisiting positions of a security for any date
#        starting after 'trading_start'. Neded for accounting/valuation purposes later.

if Date(start_date) < Date(trading_start)

    # Date to start portfolio analysis: should not set to values before 'trading_start'
    start_date           = trading_start;
    println("Warning: 'start_date' set to "*trading_start*" , as quantities held are zero before.")

    # Initial quantities are zero by definition
    initial_quantities, traded_securities  =  get_initial_quantities(raw_trades,start_date, header_types);
    init_name_quantities = Dict(full_name_list[i] => initial_quantities[i] for i = 1:length(initial_quantities));

    # Find which funds where traded in regarded time period
    name_list = full_name_list[traded_securities.> 0];

    # Extract securities actually traded in interval from raw_trades: Save in sub-DataFrame 'trades'.
    # If full horizon is selected, all securities and all trades are selected. 'trade' will be identical to 'raw_trades'.

    traded_securities = findall(x -> x in name_list,full_name_list);
    traded_pos_start  = convert.(Int64,hcat([(j-1)*length(header_types)+1 for (i,j) in enumerate(traded_securities)].+1));
    traded_pos_end    = traded_pos_start.+(length(header_types)-1);

    index_to_rearrange_raw_trades = rearrange_trade_data(traded_pos_start,traded_pos_end, header_types, name_list);
    trades = raw_trades[:,index_to_rearrange_raw_trades]

elseif Date(start_date) >= Date(trading_start)

    # Find initial quantities, update raw_trades to have only those securites that were traded
    initial_quantities, traded_securities    =  get_initial_quantities(raw_trades,start_date, header_types);
    init_name_quantities = Dict(full_name_list[i] => initial_quantities[i] for i = 1:length(initial_quantities))

    # Find which funds where traded in regarded time period
    name_list = full_name_list[traded_securities.> 0];

    # Extract securities actually traded in interval from raw_trades: Save in sub-DataFrame 'trades'.

    traded_securities = findall(x -> x in name_list,full_name_list);
    traded_pos_start  = convert.(Int64,hcat([(j-1)*length(header_types)+1 for (i,j) in enumerate(traded_securities)].+1));
    traded_pos_end    = traded_pos_start.+(length(header_types)-1);

    index_to_rearrange_raw_trades = rearrange_trade_data(traded_pos_start,traded_pos_end, header_types, name_list);
    trades = raw_trades[:,index_to_rearrange_raw_trades]

else
    println("Warning: Something went wrong when defining the date value. Check the spreadsheet.")

end

# (3) Obtain price information and clean raw data

# (3a) Preparations

# Labels of reduced security list in symbol notation (used for meta-programming)
labels      = dynamic_tuple(name_list)

# Labels of full security list in symbol notation (used for meta-programming)
full_labels = dynamic_tuple(full_name_list);

# String version of labels: Used as identifier in loops and to access list of variables
df_list     = string.(labels);
var_list    = [df_list[i]*"_Stock" for i = 1:length(labels)]

# Produce header row in symbol notation
csv_headers = reduce(vcat,[header_types.*df_list[i] for i = 1:length(labels)]);
# Add first column header entry
csv_headers = ["Date"; csv_headers];

csv_headers_full = reduce(vcat,[header_types.*string.(full_labels)[i] for i = 1:length(full_labels)]);
# Add first column header entry
csv_headers_full = ["Date"; csv_headers_full];

# Exchange headers in dataframes such that they are accessible with meta-programming using symbol notation.
trades        = names!(trades, Symbol.(csv_headers));
raw_trades    = names!(raw_trades, Symbol.(csv_headers_full))

# (3b) Download Market Data

# Mapping that applies the "draw_series"-function to every variable in name_list

# Raw data: It can happen that dataframes are of different length (holidays).
# We will take care of that in following steps.
vals          = map(draw_series,name_list);

# Construct NamedTuples of raw data, which makes data moer accessible.
# i.e. untrimmed.'label[1]' returns a DataFrame of the first security.
untrimmed     = NamedTuple{Tuple(labels)}(vals);

# Get all stocks to the same size (Differences due to different trading days)

# Trim data to same size
data          = trim_dataframe!(untrimmed);

# Export securities as individual DataFrames

for (i,j) in enumerate(labels)
    rhs = "data["*string(i)*"]";
    lhs = String(j)*"_Stock";
    eval(Meta.parse(lhs*"="*rhs))
end

# (4) Adjust portfolio information for trades and check if any day at which trading took place
# has been deleted in trimming process.

# Define Timeline: Choose whether stock/ETF you like as all df should be of the same length
timeline      = data[1].Date

start_index   = minimum(findall(x -> x >= Date(start_date), raw_trades.Date));
if isempty(end_date)
 end_index = maximum(findall(x -> x >= Date(start_date), raw_trades.Date));
else
    end_index     = maximum(findall(x -> x <= Date(end_date), raw_trades.Date));
end

# # Further adjust dataframe with trading info accordingly
trades        = trades[start_index:end_index,:];

# Find all days in which sales / purchases took place
trading_days  = trades[!, :Date]
trade_index   = sort(findall(x -> x in trading_days, timeline), rev = true)
#
trading_days_in_timeline = timeline[trade_index]

##
if trading_days == trading_days_in_timeline

else

    println("Warning: There are days at which trades took place that have been thrown out of the sample. Think about altering the trimming routine")

    missing_days = findall(x -> x == 1, trading_days .∉ [trading_days_in_timeline])
    missing_days = trading_days[missing_days]

    data = recover_missing_data!(untrimmed, data, missing_days,trading_days,labels)

    # Export (again) as individual DataFrames
    for (i,j) in enumerate(labels)
        rhs = "data["*string(i)*"]";
        lhs = String(j)*"_Stock";
        eval(Meta.parse(lhs*"="*rhs))
    end

    # Reset Timeline
    timeline      = data[1].Date

    trade_index   = sort(findall(x -> x in trading_days, timeline), rev = true)

    trading_days_in_timeline = timeline[trade_index]

    ##
    if trading_days == trading_days_in_timeline
    println("Discarded information on trading day data has been succesfully recovered. Missing data has been set to the last available trading day.")
    else
    println("Something went wrong. Check the algorithm for potential undetected bugs.")
    end
end

##

for (i,v) in enumerate(df_list)

    # add quantities
    df_new_entry = Meta.parse("Cum_Quantity_"*v)
    df_entry     = Meta.parse("Quantity_"*v)
    #raw_trades[!,df_new_entry].= initial_quantities[i].+cumsum(raw_trades[!,df_entry])
    raw_trades[!,df_new_entry].= cumsum(raw_trades[!,df_entry])

    # add information about invested amount, payments and cost, price_tag
    df_invested_amount = Meta.parse("Invested_Amount_"*v);
    df_takeouts        = Meta.parse("Payouts_"*v);
    df_payments_entry  = Meta.parse("Payments_"*v);
    df_cost_entry      = Meta.parse("Cost_in_Euro_"*v);
    df_price_tag       = Meta.parse("Price_"*v);

    raw_trades[!,df_payments_entry] .= raw_trades[!,df_payments_entry];
    raw_trades[!,df_cost_entry]     .= raw_trades[!,df_cost_entry];
    raw_trades[!,df_price_tag]      .= raw_trades[!,df_price_tag];
    raw_trades[!,df_invested_amount].= 0.0;

    # Calculate the amount invested and adjust accordingly it after a sale occurs (Fifo method)
    for i_trade = 1:size(raw_trades,1)

    # special condition of first entry of Invested Amount (nothing to cumulate)

        if i_trade == 1
            raw_trades[i_trade, df_invested_amount] = raw_trades[i_trade,df_entry]*raw_trades[i_trade,df_price_tag];
        else

            if raw_trades[i_trade,df_entry] >= 0
                raw_trades[i_trade, df_invested_amount] = raw_trades[i_trade-1, df_invested_amount]+raw_trades[i_trade,df_entry]*raw_trades[i_trade,df_price_tag];
            elseif raw_trades[i_trade,df_entry] < 0

                amount_to_sell         = raw_trades[i_trade,df_entry];
                #cumulated_quantity_aux = raw_trades[!,df_new_entry];
                cumulated_quantity_aux = cumsum(raw_trades[!,df_entry]);
                remaining_quantity     =  cumulated_quantity_aux[i_trade-1] - abs(amount_to_sell); # cumulated quantity up to current sale

                if remaining_quantity > 0 # not everything has been sold -> calculate remaining invested amount according to FIFO accounting

                    # Find ALL purchases and constrain index to all purchases BEFORE current sale (and even time horizon considered)

                    all_purchases =  reverse(findall( x -> x > 0, raw_trades[!,df_entry]));
                    buy_index     =  all_purchases[all_purchases.<i_trade];

                    if isempty(buy_index)
                    all_purchases =  reverse(findall( x -> x > 0, raw_trades[!,df_entry]));
                    buy_index     =  all_purchases[all_purchases.<(i_trade)];

                    else
                    end

                    check_aux = remaining_quantity;
                    i_sell    = 0;

                    while check_aux > 0

                        i_sell = i_sell +1;

                        if raw_trades[buy_index[i_sell],df_entry] < check_aux #trades.Quantity_Tech(buy_index(i_sell)) < check_aux

                            # Case: Last purchase was smaller than overall
                            # quantity remaining. Remaining amount has to be
                            # allotted the the previous purchases until every
                            # part of the remaining amount has been allocated.

                            raw_trades[i_trade, df_invested_amount] = raw_trades[i_trade, df_invested_amount] + raw_trades[buy_index[i_sell], df_entry]*raw_trades[buy_index[i_sell], df_price_tag];
                            check_aux = cumulated_quantity_aux[buy_index[i_sell+1]] - abs(raw_trades[i_trade, df_entry]);

                        elseif raw_trades[buy_index[i_sell],df_entry] > check_aux

                            # Case: Last purchase was higher than overall
                            # quantity remaining. Invested amount has to be
                            # priced using the price of the last purchase and
                            # the remaining quantity.
                            raw_trades[i_trade, df_invested_amount] = raw_trades[i_trade, df_invested_amount] + check_aux*raw_trades[buy_index[i_sell], df_price_tag];
                            check_aux = 0; # used to exit while-condition - has to hold in this case as all remaining quantity is priced with price of the last purchase

                        end

                    end

                else # everything has been sold -> remaining invested amount has to be zero
                end
            else
            end
        end
    end

end

##

trades        = raw_trades[start_index:end_index,:];

for (i,v) in enumerate(df_list)

    # add quantities
    df_new_entry = Meta.parse("Cum_Quantity_"*v)

    # add information about invested amount, payments and cost, price_tag
    df_invested_amount = Meta.parse("Invested_Amount_"*v);

    # Broadcast Variables to full dataframe over all dates
    for (j,k) in enumerate(trade_index)
        # Quantity
        lhs = var_list[i]*"["*string(1:k) *",:Quantity_held]"
        rhs = " .= trades["*string(j)*",:"*string(df_new_entry)*"]"
        iter_string = Meta.parse(lhs*rhs)
        eval(iter_string)

        # Invested_Amount
        lhs_invs_ = var_list[i]*"["*string(1:k) *",:Invested_Amount]"
        rhs_invs_ = " .= trades["*string(j)*",:"*string(df_invested_amount)*"]"
        iter_string_invs = Meta.parse(lhs_invs_*rhs_invs_)
        eval(iter_string_invs)
    end

    lhs = var_list[i]*"[!,:Current_value]"
    rhs = " .="*var_list[i]*"[!,:adjusted_close].*"*var_list[i]*"[!,:Quantity_held]"
    iter_string    = Meta.parse(lhs*rhs)
    eval(iter_string)

end

# (5) Calculate some portfolio stats: Total Value, Returns, Portfolio Weights etc
total_value, total_amount_invested, total_return, portfolio_weights = calc_portfolio_stats(var_list, df_list)


text:
using Oceananigans.Fields: TracerFields

using Oceananigans.Buoyancy:
    required_tracers, ρ′, ∂x_b, ∂y_b,
    thermal_expansionᶜᶜᶜ, thermal_expansionᶠᶜᶜ, thermal_expansionᶜᶠᶜ, thermal_expansionᶜᶜᶠ,
    haline_contractionᶜᶜᶜ, haline_contractionᶠᶜᶜ, haline_contractionᶜᶠᶜ, haline_contractionᶜᶜᶠ

function instantiate_linear_equation_of_state(FT, α, β)
    eos = LinearEquationOfState(FT, α=α, β=β)
    return eos.α == FT(α) && eos.β == FT(β)
end

function instantiate_seawater_buoyancy(FT, EquationOfState; kwargs...)
    buoyancy = SeawaterBuoyancy(FT, equation_of_state=EquationOfState(FT); kwargs...)
    return typeof(buoyancy.gravitational_acceleration) == FT
end

function density_perturbation_works(arch, FT, eos)
    grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(1, 1, 1))
    C = datatuple(TracerFields((:T, :S), arch, grid))
    density_anomaly = ρ′(2, 2, 2, grid, eos, C.T, C.S)
    return true
end

function ∂x_b_works(arch, FT, buoyancy)
    grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(1, 1, 1))
    C = datatuple(TracerFields(required_tracers(buoyancy), arch, grid))
    dbdx = ∂x_b(2, 2, 2, grid, buoyancy, C)
    return true
end

function ∂y_b_works(arch, FT, buoyancy)
    grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(1, 1, 1))
    C = datatuple(TracerFields(required_tracers(buoyancy), arch, grid))
    dbdy = ∂y_b(2, 2, 2, grid, buoyancy, C)
    return true
end

function ∂z_b_works(arch, FT, buoyancy)
    grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(1, 1, 1))
    C = datatuple(TracerFields(required_tracers(buoyancy), arch, grid))
    dbdz = ∂z_b(2, 2, 2, grid, buoyancy, C)
    return true
end

function thermal_expansion_works(arch, FT, eos)
    grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(1, 1, 1))
    C = datatuple(TracerFields((:T, :S), arch, grid))
    α = thermal_expansionᶜᶜᶜ(2, 2, 2, grid, eos, C.T, C.S)
    α = thermal_expansionᶠᶜᶜ(2, 2, 2, grid, eos, C.T, C.S)
    α = thermal_expansionᶜᶠᶜ(2, 2, 2, grid, eos, C.T, C.S)
    α = thermal_expansionᶜᶜᶠ(2, 2, 2, grid, eos, C.T, C.S)
    return true
end

function haline_contraction_works(arch, FT, eos)
    grid = RegularCartesianGrid(FT, size=(3, 3, 3), extent=(1, 1, 1))
    C = datatuple(TracerFields((:T, :S), arch, grid))
    β = haline_contractionᶜᶜᶜ(2, 2, 2, grid, eos, C.T, C.S)
    β = haline_contractionᶠᶜᶜ(2, 2, 2, grid, eos, C.T, C.S)
    β = haline_contractionᶜᶠᶜ(2, 2, 2, grid, eos, C.T, C.S)
    β = haline_contractionᶜᶜᶠ(2, 2, 2, grid, eos, C.T, C.S)
    return true
end

EquationsOfState = (LinearEquationOfState, SeawaterPolynomials.RoquetEquationOfState, SeawaterPolynomials.TEOS10EquationOfState)
buoyancy_kwargs = (Dict(), Dict(:constant_salinity=>35.0), Dict(:constant_temperature=>20.0))

@testset "Buoyancy" begin
    @info "Testing buoyancy..."

    @testset "Equations of State" begin
        @info "  Testing equations of state..."
        for FT in float_types
            @test instantiate_linear_equation_of_state(FT, 0.1, 0.3)

            for EOS in EquationsOfState
                for kwargs in buoyancy_kwargs
                    @test instantiate_seawater_buoyancy(FT, EOS; kwargs...)
                end
            end

            for arch in archs
                @test density_perturbation_works(arch, FT, SeawaterPolynomials.RoquetEquationOfState())
            end

            buoyancies =
                (nothing, BuoyancyTracer(), SeawaterBuoyancy(FT),
                 (SeawaterBuoyancy(FT, equation_of_state=eos(FT)) for eos in EquationsOfState)...)

            for arch in archs
                for buoyancy in buoyancies
                    @test ∂x_b_works(arch, FT, buoyancy)
                    @test ∂y_b_works(arch, FT, buoyancy)
                    @test ∂z_b_works(arch, FT, buoyancy)
                end
            end

            for arch in archs
                for EOS in EquationsOfState
                    @test thermal_expansion_works(arch, FT, EOS())
                    @test haline_contraction_works(arch, FT, EOS())
                end
            end
        end
    end
end

ext: jl


text:
import Base: IEEEFloat

uint(x::Float64) = reinterpret(UInt64, x)
uint(x::Float32) = reinterpret(UInt32, x)
uint(x::Float16) = reinterpret(UInt16, x)

const nanbit64 = xor(uint(NaN), uint(Inf))
const nanbit32 = xor(uint(NaN32), uint(Inf32))
const nanbit16 = xor(uint(NaN16), uint(Inf16))

@inline function isolate_nan(a::T, b::T) where {T<:IEEEFloat}
    return uint(a) | uint(b)
end

@inline function isolate_nan(a::U, b::T) where {U<:Unsigned, T<:IEEEFloat}
    return a | uint(b)
end

@inline function isolate_nan(a::U, b::U) where {U<:Unsigned}
    return a | b
end

@inline function isolate_nan(vec::AbstractVector{T}) where {T<:IEEEFloat}
    reduce( isolate_nan, vec )
end
#=
@inline function isolate_nan(vec::AbstractVector{T}) where {T<:IEEEFloat}
    reduce( isolate_nan, vec )
end
=#

function hasnan(vec::AbstractVector{T}) where {T<:Float64}
    nanbit64 === nanbit64 & isolate_nan(vec)
end
function hasnan(vec::AbstractVector{T}) where {T<:Float32}
    nanbit32 === nanbit32 & isolate_nan(vec)
end
function hasnan(vec::AbstractVector{T}) where {T<:Float16}
    nanbit16 === nanbit16 & isolate_nan(vec)
end


text:
export
    ContinuedFraction, NonPeriodicContinuedFraction, PeriodicContinuedFraction,
    confrac,
    ContinuedFractionConvergentsIterator, NonPeriodicContinuedFractionConvergentsIterator, PeriodicContinuedFractionConvergentsIterator,
    convergents

abstract type ContinuedFraction{T<:Integer} end

@auto_hash_equals struct NonPeriodicContinuedFraction{T<:Integer} <: ContinuedFraction{T}
    denominators::Array{T,1}
end

@auto_hash_equals struct PeriodicContinuedFraction{T<:Integer} <: ContinuedFraction{T}
    denominators::Array{T,1}
    periodstart::Int
end

Base.start(::ContinuedFraction) = 1
Base.next(it::NonPeriodicContinuedFraction, state) = it.denominators[state], state + 1
Base.next(it::PeriodicContinuedFraction, state) = it.denominators[state], (state == length(it.denominators) ? it.periodstart : state + 1)
Base.done(it::NonPeriodicContinuedFraction, state) = state > length(it.denominators)
Base.done(::PeriodicContinuedFraction, _) = false

Base.eltype(it::ContinuedFraction) = Base.eltype(typeof(it))
Base.eltype(::Type{NonPeriodicContinuedFraction{T}}) where {T} = T
Base.eltype(::Type{PeriodicContinuedFraction{T}}) where {T} = T

Base.length(it::NonPeriodicContinuedFraction) = length(it.denominators)
Base.iteratorsize(::PeriodicContinuedFraction) = Base.IsInfinite()

Base.show(io::IO, cf::ContinuedFraction) = begin
    print(io, "[")
    for (i, d) in enumerate(cf.denominators)
        i == 2 && print(io, ";")
        i >= 3 && print(io, ",")
        isa(cf, PeriodicContinuedFraction) && i == cf.periodstart && print(io, "(")
        print(io, d)
    end
    isa(cf, PeriodicContinuedFraction) && print(io, ")")
    print(io, "]")
end

Base.rationalize(cf::ContinuedFraction{T}, len::Int = -1) where {T<:Integer} = begin
    len >= 0 || (len = length(cf.denominators))
    len > 0 || return one(T) // zero(T)   # [] = +Inf;  [a; b] = a + 1 / [b] = a + 1 / (b + 1 / []) = a + 1 / (b + 1 / +Inf) = a + 1 / b

    @pipe cf |> convergents |> drop(_, len-1) |> first |> Rational(_...)
end

confrac(d::Array{T,1}) where {T<:Integer} = NonPeriodicContinuedFraction(d)
confrac(d::Array{T,1}, p::Int) where {T<:Integer} = PeriodicContinuedFraction(d, p)

confrac(rat::Rational{T}) where {T<:Integer} = begin
    n, d = rat.num, rat.den

    denominators = T[]
    while d != 0
        a = fld(n, d)
        push!(denominators, a)
        n, d = d, n - a * d
    end

    confrac(denominators)
end

# == (n + sqrt(δ)) / d
confrac(n::T, δ::T, d::T) where {T<:Integer} = begin
    d == 0 && throw(Error("The denominator is zero."))
    δ < 0 && throw(Error("The number is not real, so it does not have a continued fraction expansion."))
    δ == 0 && return confrac(n // d)

    if (δ - n^2) % d != 0
        δ *= d^2
        n *= abs(d)
        d *= abs(d)
    end

    sqrtδ = isqrt(δ)
    sqrtδ^2 == δ && return confrac((n + sqrtδ) // d)

    biP = d
    biK = sqrtδ + (biP > 0 ? 0: 1) + n
    if biK > 0
        biM = d > 0 ? biK ÷ d : (d + 1 - biK) ÷ (-d)
    else
        biM = d > 0 ? (biK + 1 - d) ÷ d : (-biK) ÷ (-d)
    end

    denominators = T[]
    periodstart = 0
    push!(denominators, biM)
    biM = biM * d  - n

    cont = -1
    K = P = L = M = -1

    while (cont < 0) || (K != P) || (L != M)
        if (cont < 0) && (biP > 0) && (biP <= sqrtδ + biM) && (biM > 0) && (biM <= sqrtδ)
            # period starts
            periodstart = length(denominators) + 1
            K = P = biP
            L = M = biM
            cont = 0
        end

        if cont >= 0
            # both numerator and denominator are positive
            P = (δ - M^2) ÷ P
            Z = (sqrtδ + M) ÷ P
            M = Z * P - M
            cont += 1
        else
            biP = (δ - biM^2) ÷ biP
            Z = (sqrtδ + biM + (biP > 0 ? 0 : 1)) ÷ biP
            biM = Z * biP - biM
        end

        push!(denominators, Z)
        #if (cont > 100000)
        #    # too many convergents
        #    break
        #end
    end

    confrac(denominators, periodstart)
end

abstract type ContinuedFractionConvergentsIterator{T<:Integer} end

struct NonPeriodicContinuedFractionConvergentsIterator{T<:Integer} <: ContinuedFractionConvergentsIterator{T}
    cf::NonPeriodicContinuedFraction{T}
end

struct PeriodicContinuedFractionConvergentsIterator{T<:Integer} <: ContinuedFractionConvergentsIterator{T}
    cf::PeriodicContinuedFraction{T}
end

convergents(cf::NonPeriodicContinuedFraction) = NonPeriodicContinuedFractionConvergentsIterator(cf)
convergents(cf::PeriodicContinuedFraction) = PeriodicContinuedFractionConvergentsIterator(cf)

Base.start(it::ContinuedFractionConvergentsIterator{T}) where {T<:Integer} = start(it.cf), one(T), zero(T), zero(T), one(T)
Base.next(it::ContinuedFractionConvergentsIterator, state) = begin
    cfs, p, q, r, s = state

    c, cfs = next(it.cf, cfs)
    p, q, r, s = c*p+r, c*q+s, p, q

    (p, q), (cfs, p, q, r, s)
end
Base.done(it::ContinuedFractionConvergentsIterator, state) = done(it.cf, state[1])

Base.eltype(it::ContinuedFractionConvergentsIterator) = Base.eltype(typeof(it))
Base.eltype(::Type{NonPeriodicContinuedFractionConvergentsIterator{T}}) where {T} = Tuple{T,T}
Base.eltype(::Type{PeriodicContinuedFractionConvergentsIterator{T}}) where {T} = Tuple{T,T}

Base.length(it::NonPeriodicContinuedFractionConvergentsIterator) = length(it.cf.denominators)
Base.iteratorsize(::PeriodicContinuedFractionConvergentsIterator) = Base.IsInfinite()


text:
"""
Atomistic solution for an ODE.

### Fields

* `t`: time of current time step
* `t̅`: time of previous time step
* `q`: current solution of q
* `q̅`: previous solution of q
* `q̃`: compensated summation error of q
* `v`: vector field of q
* `v̅`: vector field of q̅
"""
mutable struct AtomicSolutionODE{DT,TT} <: AtomicSolution{DT,TT}
    t::TT
    t̅::TT

    q::Vector{DT}
    q̅::Vector{DT}
    q̃::Vector{DT}

    v::Vector{DT}
    v̅::Vector{DT}

    function AtomicSolutionODE{DT, TT}(nd) where {DT <: Number, TT <: Real}
        new(zero(TT), zero(TT), zeros(DT, nd), zeros(DT, nd), zeros(DT, nd),
                                zeros(DT, nd), zeros(DT, nd))
    end
end

AtomicSolutionODE(DT, TT, nd) = AtomicSolutionODE{DT, TT}(nd)

function CommonFunctions.set_solution!(asol::AtomicSolutionODE, sol)
    t, q = sol
    asol.t  = t
    asol.q .= q
    asol.v .= 0
end

function CommonFunctions.get_solution(asol::AtomicSolutionODE)
    (asol.t, asol.q)
end

function CommonFunctions.reset!(asol::AtomicSolutionODE, Δt)
    asol.t̅  = asol.t
    asol.q̅ .= asol.q
    asol.v̅ .= asol.v
    asol.t += Δt
end

function update!(asol::AtomicSolutionODE{DT}, v::Vector{DT}) where {DT}
    for k in eachindex(v)
        update!(asol, v[k])
    end
end

function update!(asol::AtomicSolutionODE{DT}, v::DT, k::Int) where {DT}
    asol.q[k], asol.q̃[k] = compensated_summation(v, asol.q[k], asol.q̃[k])
end


text:
@inline getbody(mechanism::Mechanism, id::Int64) = mechanism.bodies[id]
@inline getbody(mechanism::Mechanism, id::Nothing) = mechanism.origin
function getbody(mechanism::Mechanism, name::String)
    if mechanism.origin.name == name
        return mechanism.origin
    else
        for body in mechanism.bodies
            if body.name == name
                return body
            end
        end
    end
    return nothing
end
@inline geteqconstraint(mechanism::Mechanism, id::Int64) = mechanism.eqconstraints[id]
function geteqconstraint(mechanism::Mechanism, name::String)
    for eqc in mechanism.eqconstraints
        if eqc.name == name
            return eqc
        end
    end
    return nothing
end
@inline getineqconstraint(mechanism::Mechanism, id::Int64) = mechanism.ineqconstraints[id]
function getineqconstraint(mechanism::Mechanism, name::String)
    for ineqc in mechanism.eqconstraints
        if ineqc.name == name
            return ineqc
        end
    end
    return nothing
end

getcomponent(mechanism::Mechanism, id::Nothing) = mechanism.origin
function getcomponent(mechanism::Mechanism, id::Int64)
    if haskey(mechanism.bodies, id)
        return getbody(mechanism, id)
    elseif haskey(mechanism.eqconstraints, id)
        return geteqconstraint(mechanism, id)
    elseif haskey(mechanism.ineqconstraints, id)
        return getineqconstraint(mechanism, id)
    else
        return nothing
    end
end

function getcomponent(mechanism::Mechanism, name::String)
    component = getbody(mechanism,name)
    if component == nothing
        component = geteqconstraint(mechanism,name)
    end
    if component == nothing
        component = getineqconstraint(mechanism,name)
    end
    return component
end

function getshape(mechanism::Mechanism, id)
    for shape in mechanism.shapes
        for bodyid in shape.bodyids
            if bodyid == id
                return shape
            end
        end
    end

    return nothing
end

function setentries!(mechanism::Mechanism)
    graph = mechanism.graph
    ldu = mechanism.ldu

    for (id, body) in pairs(mechanism.bodies)
        for cid in directchildren(graph, id)
            setLU!(mechanism, getentry(ldu, (id, cid)), id, geteqconstraint(mechanism, cid))
        end

        diagonal = getentry(ldu, id)
        setDandΔs!(mechanism, diagonal, body)
        for cid in ineqchildren(graph, id)
            extendDandΔs!(mechanism, diagonal, body, getineqconstraint(mechanism, cid))
        end
    end

    for node in mechanism.eqconstraints
        id = node.id

        for cid in directchildren(graph, id)
            setLU!(mechanism, getentry(ldu, (id, cid)), node, cid)
        end

        for cid in loopchildren(graph, id)
            setLU!(getentry(ldu, (id, cid)))
        end

        diagonal = getentry(ldu, id)
        setDandΔs!(mechanism, diagonal, node)
    end
end

@inline function normf(mechanism::Mechanism, body::Body{T}) where T
    f = dynamics(mechanism, body)
    return dot(f, f)
end

@inline function normf(mechanism::Mechanism, c::EqualityConstraint)
    f = g(mechanism, c)
    return dot(f, f)
end

@inline function normf(mechanism::Mechanism, ineqc::InequalityConstraint)
    f = gs(mechanism, ineqc)
    d = h(ineqc)
    return dot(f, f) + dot(d, d)
end

@inline function normfμ(mechanism::Mechanism, ineqc::InequalityConstraint)
    f = gs(mechanism, ineqc)
    d = hμ(ineqc, mechanism.μ)
    return dot(f, f) + dot(d, d)
end

@inline function GtλTof!(mechanism::Mechanism, body::Body, eqc::EqualityConstraint)
    body.f -= ∂g∂pos(mechanism, eqc, body.id)' * eqc.s1
    return
end

@inline function NtγTof!(mechanism::Mechanism, body::Body, ineqc::InequalityConstraint)
    body.f -= ∂g∂pos(mechanism, ineqc, body)' * ineqc.γ1
    return
end

@inline function normf(mechanism::Mechanism)
    mechanism.normf = 0

    # Allocates otherwise
    for body in mechanism.bodies
        mechanism.normf += normf(mechanism, body)
    end
    foreach(addNormf!, mechanism.eqconstraints, mechanism)
    foreach(addNormf!, mechanism.ineqconstraints, mechanism)

    return sqrt(mechanism.normf)
end

@inline function meritf(mechanism::Mechanism)
    mechanism.normf = 0

    # Allocates otherwise
    for body in mechanism.bodies
        mechanism.normf += normf(mechanism, body)
    end
    foreach(addNormf!, mechanism.eqconstraints, mechanism)
    foreach(addNormfμ!, mechanism.ineqconstraints, mechanism)

    return sqrt(mechanism.normf)
end

@inline function normΔs(mechanism::Mechanism)
    mechanism.normΔs = 0

    # Allocates otherwise
    mechanism.normΔs += mapreduce(normΔs, +, mechanism.bodies)
    foreach(addNormΔs!, mechanism.eqconstraints, mechanism)
    foreach(addNormΔs!, mechanism.ineqconstraints, mechanism)

    return sqrt(mechanism.normΔs)
end

@inline function addNormf!(ineqc::InequalityConstraint, mechanism::Mechanism)
    mechanism.normf += normf(mechanism, ineqc)
    return
end

@inline function addNormfμ!(ineqc::InequalityConstraint, mechanism::Mechanism)
    mechanism.normf += normfμ(mechanism, ineqc)
    return
end

@inline function addNormf!(eqc::EqualityConstraint, mechanism::Mechanism)
    mechanism.normf += normf(mechanism, eqc)
    return
end

@inline function addNormΔs!(component::Component, mechanism::Mechanism)
    mechanism.normΔs += normΔs(component)
    return
end

function feasibilityStepLength!(mechanism::Mechanism)
    ldu = mechanism.ldu

    τ = 0.995
    mechanism.α = 1.

    for ineqc in mechanism.ineqconstraints
        feasibilityStepLength!(mechanism, ineqc, getineqentry(ldu, ineqc.id), τ)
    end

    return
end

function feasibilityStepLength!(mechanism, ineqc::InequalityConstraint{T,N}, ineqentry::InequalityEntry, τ) where {T,N}
    s1 = ineqc.s1
    γ1 = ineqc.γ1
    Δs = ineqentry.Δs
    Δγ = ineqentry.Δγ

    for i = 1:N
        αmax = τ * s1[i] / Δs[i]
        (αmax > 0) && (αmax < mechanism.α) && (mechanism.α = αmax)
        αmax = τ * γ1[i] / Δγ[i]
        (αmax > 0) && (αmax < mechanism.α) && (mechanism.α = αmax)
    end
    return
end

function setFrictionForce!(mechanism)
    foreach(setFrictionForce!, mechanism.ineqconstraints, mechanism)
end

function saveToStorage!(mechanism::Mechanism, t)
    No = mechanism.No
    for (ind, body) in enumerate(mechanism.bodies)
        mechanism.storage.x[ind][t] = body.x[No]
        mechanism.storage.q[ind][t] = body.q[No]
    end
    for (ind, eqc) in enumerate(mechanism.eqconstraints)
        mechanism.storage.eqmultipliers[ind][t] = eqc.s1
    end
end

@inline function updatePos!(body::Body, Δt)
    body.x[1] = body.x[2]
    body.x[2] = getx3(body, Δt)
    body.q[1] = body.q[2]
    body.q[2] = getq3(body, Δt)
    return
end

function verifyConstraints!(mechanism::Mechanism)
    for eqc in mechanism.eqconstraints
        if norm(g(mechanism, eqc)) > 1e-3
            @info string("Probably disconnected bodies at constraint: ", eqc.id)
        end
    end
end

function simulate!(mechanism::Mechanism;save::Bool = false,debug::Bool = false)
    debug && verifyConstraints!(mechanism)
    bodies = mechanism.bodies
    eqcs = mechanism.eqconstraints
    ineqcs = mechanism.ineqconstraints
    Δt = mechanism.Δt
    foreach(s0tos1!, bodies)
    foreach(s0tos1!, eqcs)
    foreach(s0tos1!, ineqcs)

    for i = mechanism.steps
        newton!(mechanism, warning = debug)
        save && saveToStorage!(mechanism, i)
        foreach(updatePos!, bodies, Δt)
    end
    return
end

function simulate!(mechanism::Mechanism, control!::Function;save::Bool = false,debug::Bool = false)
    debug && verifyConstraints!(mechanism)
    bodies = mechanism.bodies
    eqcs = mechanism.eqconstraints
    ineqcs = mechanism.ineqconstraints
    Δt = mechanism.Δt
    foreach(s0tos1!, bodies)
    foreach(s0tos1!, eqcs)
    foreach(s0tos1!, ineqcs)

    for i = mechanism.steps
        control!(mechanism, (i - 1) * Δt)
        newton!(mechanism, warning = debug)
        save && saveToStorage!(mechanism, i)
        foreach(updatePos!, bodies, Δt)
    end
    return
end

function simulate!(mechanism::Mechanism, controller::Controller;save::Bool = false,debug::Bool = false)
    debug && verifyConstraints!(mechanism)
    bodies = mechanism.bodies
    eqcs = mechanism.eqconstraints
    ineqcs = mechanism.ineqconstraints
    Δt = mechanism.Δt
    foreach(s0tos1!, bodies)
    foreach(s0tos1!, eqcs)
    foreach(s0tos1!, ineqcs)

    for i = mechanism.steps
        control!(mechanism, controller)
        newton!(mechanism, warning = debug)
        save && saveToStorage!(mechanism, i)
        foreach(updatePos!, bodies, Δt)
    end
    return
end

function plotθ(mechanism::Mechanism{T}, id) where T
    n = length(mechanism.bodies)
    θ = zeros(T, n, length(mechanism.steps))
    for i = 1:n
        qs = mechanism.storage.q[i]
        for (t, q) in enumerate(qs)
            θ[i,t] = angleaxis(q)[1] * sign(angleaxis(q)[2][1])
        end
    end

    p = plot(collect(0:mechanism.Δt:mechanism.tend - mechanism.Δt), θ[id[1],:])
    for ind in Iterators.rest(id, 2)
        plot!(collect(0:mechanism.Δt:mechanism.tend - mechanism.Δt), θ[ind,:])
    end
    return p
end

function plotλ(mechanism::Mechanism{T}, id) where T
    n = sum(length.(mechanism.eqconstraints))
    λ = zeros(T, n, length(mechanism.steps))
    startpos = 1
    endpos = 0
    for i = 1:length(mechanism.eqconstraints)
        endpos = startpos + length(mechanism.eqconstraints[i]) - 1

        λs = mechanism.storage.λ[i]
        for (t, val) in enumerate(λs)
            λ[startpos:endpos,t] = val
        end

        startpos = endpos + 1
    end

    p = plot(collect(0:mechanism.Δt:mechanism.tend - mechanism.Δt), λ[id[1],:])
    for ind in Iterators.rest(id, 2)
        plot!(collect(0:mechanism.Δt:mechanism.tend - mechanism.Δt), λ[ind,:])
    end
    return p
end


text:
export by, aggregate

by(args...; kwargs...) = throw(ArgumentError("by function was removed from DataFrames.jl. " *
                                             "Use the `combine(groupby(...), ...)` or `combine(f, groupby(...))` instead."))

aggregate(args...; kwargs...) = throw(ArgumentError("aggregate function was removed from DataFrames.jl. " *
                                                    "Use the `combine` function instead."))


text:
#############################################################################
# Joulia
# A Large-Scale Spatial Power System Model for Julia
# See https://github.com/JuliaEnergy/Joulia.jl
#############################################################################
# Example: ELMOD-DE

using Joulia

using DataFrames, CSV
using Gurobi

data_path = "data"

pp_df = CSV.read(joinpath(data_path, "power_plants.csv"))
avail_con_df = CSV.read(joinpath(data_path, "avail_con.csv"))
prices_df = CSV.read(joinpath(data_path, "prices.csv"))

storages_df = CSV.read(joinpath(data_path, "storages.csv"))

lines_df = CSV.read(joinpath(data_path, "lines.csv"))

load_df = CSV.read(joinpath(data_path, "load.csv"))
nodes_df = CSV.read(joinpath(data_path, "nodes.csv"))
exchange_df = CSV.read(joinpath(data_path, "exchange.csv"))

res_df = CSV.read(joinpath(data_path, "res.csv"))
avail_pv = CSV.read(joinpath(data_path, "avail_pv.csv"))
avail_windon = CSV.read(joinpath(data_path, "avail_windon.csv"))
avail_windoff = CSV.read(joinpath(data_path, "avail_windoff.csv"))

avail = Dict(:PV => avail_pv,
	:WindOnshore => avail_windon,
	:WindOffshore => avail_windoff,
	:global => avail_con_df)

# generation of Joulia data types
pp = PowerPlants(pp_df, avail=avail_con_df, prices=prices_df)
storages = Storages(storages_df)
lines = Lines(lines_df)
nodes = Nodes(nodes_df, load_df, exchange_df)
res = RenewableEnergySource(res_df, avail)

# generation of the Joulia model
elmod = JouliaModel(pp, res, storages, nodes, lines)

# slicing the data in weeks with the first full week starting at hour 49
slices = week_slices(1)[2:end]

# running the Joulia model for week 30 using the Gurobi solver
results = run_model(elmod, 1:672, solver=Gurobi.Optimizer)

# elmod.results[:LOST_LOAD].value |> sum

# sum(nodes.load[n][t] for n in keys(nodes.load), t in 1:672)

# asdf = filter(x -> x.value > 0, elmod.results[:LOST_LOAD])

# test = combine(groupby(asdf, :Node), :value => sum => :value)

# sort!(test, :value, rev=true)

for (k, df) in elmod.results
	filename = string(k) * ".csv"
	CSV.write("results/$filename", df)
end


text:
# ---
# title: 446. Arithmetic Slices II - Subsequence
# id: problem446
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/arithmetic-slices-ii-subsequence/description/>
# hidden: true
# ---
#
# A sequence of numbers is called arithmetic if it consists of at least three
# elements and if the difference between any two consecutive elements is the
# same.
#
# For example, these are arithmetic sequences:
#
#
#
#     1, 3, 5, 7, 9
#     7, 7, 7, 7
#     3, -1, -5, -9
#
# The following sequence is not arithmetic.
#
#
#
#     1, 1, 2, 5, 7
#
#
#
# A zero-indexed array A consisting of N numbers is given. A **subsequence**
# slice of that array is any sequence of integers (P0, P1, ..., Pk) such that 0
# ≤ P0 < P1 < ... < Pk < N.
#
# A **subsequence** slice (P0, P1, ..., Pk) of array A is called arithmetic if
# the sequence A[P0], A[P1], ..., A[Pk-1], A[Pk] is arithmetic. In particular,
# this means that k ≥ 2.
#
# The function should return the number of arithmetic subsequence slices in the
# array A.
#
# The input contains N integers. Every integer is in the range of -231 and 231-1
# and 0 ≤ N ≤ 1000\. The output is guaranteed to be less than 231-1.
#
#
#
# **Example:**
#
#
#
#     Input: [2, 4, 6, 8, 10]
#
#     Output: 7
#
#     Explanation:
#     All arithmetic subsequence slices are:
#     [2,4,6]
#     [4,6,8]
#     [6,8,10]
#     [2,4,6,8]
#     [4,6,8,10]
#     [2,4,6,8,10]
#     [2,6,10]
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Test
using LuxurySparse
using LinearAlgebra
using SparseArrays

@testset "broadcast *" begin

    @testset "Diagonal .* $(nameof(typeof(M)))" for M in [pmrand(3)]
        M1 = Diagonal(rand(3))
        out = M1 .* M
        @test typeof(out) <: Diagonal
        @test out ≈ Matrix(M1) .* M

        out = M .* M1
        @test typeof(out) <: Diagonal
        @test out ≈ Matrix(M) .* M1
    end

    @testset "PermMatrix .* $(nameof(typeof(M)))" for M in [
        rand(3, 3),
        pmrand(3),
        sprand(3, 3, 0.5),
    ]
        M1 = pmrand(3)
        out = M1 .* M
        @test typeof(out) <: PermMatrix
        @test out ≈ Matrix(M1) .* M

        out = M .* M1
        @test typeof(out) <: PermMatrix
        @test out ≈ M .* Matrix(M1)
    end

    @testset "IMatrix .* $(nameof(typeof(M)))" for M in [
        rand(3, 3),
        pmrand(3),
        sprand(3, 3, 0.5),
    ]
        eye = IMatrix{3}()
        out = eye .* M
        @test typeof(out) <: Diagonal
        @test out ≈ Matrix(eye) .* M

        out = M .* eye
        @test typeof(out) <: Diagonal
        @test out ≈ M .* Matrix(eye)
    end

    @test IMatrix{3}() .* IMatrix{3}() === IMatrix{3}()
    d = Diagonal(rand(3))
    sp = sprand(3, 3, 0.5)
    @test d .* sp ≈ Matrix(d) .* Matrix(sp)
    @test typeof(d .* sp) <: Diagonal

end

@testset "broadcast -" begin

    @testset "Diagonal .- $(nameof(typeof(M)))" for M in [pmrand(3)]
        M1 = Diagonal(rand(3))
        @test M1 .- M ≈ Matrix(M1) .- M
        @test M .- M1 ≈ Matrix(M) .- M1
    end

    @testset "PermMatrix .* $(nameof(typeof(M)))" for M in [
        2.0,
        rand(3, 3),
        pmrand(3),
        sprand(3, 0.5),
        sprand(3, 3, 0.5),
    ]
        M1 = pmrand(3)
        @test M1 .- M ≈ Matrix(M1) .- M
        @test M .- M1 ≈ M .- Matrix(M1)
    end

    @testset "IMatrix .* $(nameof(typeof(M)))" for M in [
        2.0,
        rand(3, 3),
        pmrand(3),
        sprand(3, 0.5),
        sprand(3, 3, 0.5),
    ]
        eye = IMatrix{3}()
        @test eye .- M ≈ Matrix(eye) .- M

        @test M .- eye ≈ M .- Matrix(eye)
    end

    @test IMatrix{3}() .- IMatrix{3}() ≈ zeros(3, 3)
end


text:
#=

  Number lock problem in Turing.jl

  From Presh Talwalkar (MindYourDecisions)
  """
  Puzzles like this have been shared with the dubious claim that "only a
  genius can solve" them. But they are still fun problems so let's work one
  out.

  A number lock requires a 3 digit code. Based on these hints, can you crack
  the code?

   682 - one number is correct and in the correct position
   645 - one number is correct but in the wrong position
   206 - two numbers are correct but in the wrong positions
   738 - nothing is correct
   780 - one number is correct but in the wrong position

  Video:  https://youtu.be/-etLb-8sHBc
  """

  This model use the same approach as
  - MiniZinc model: http://hakank.org/minizinc/number_lock.mzn
  - Picat model: http://hakank.org/picat/number_lock.pi

  Summary Statistics
  parameters      mean       std   naive_se      mcse       ess      rhat   ess_per_sec
      Symbol   Float64   Float64    Float64   Float64   Float64   Float64       Float64

        x[1]    0.0000    0.0000     0.0000    0.0000       NaN       NaN           NaN
        x[2]    5.0000    0.0000     0.0000    0.0000       NaN       NaN           NaN
        x[3]    2.0000    0.0000     0.0000    0.0000       NaN       NaN           NaN

   Distributions of variable (num:0)
   [0, 5, 2]	=>	10000 (1.0)

  Cf ~/webppl/number_lock2.wppl

=#

using Turing, StatsPlots
include("jl_utils.jl")

@model function number_lock2(y, hints)
    # m = 5 # number of rows (i.e. number of hints)
    # n = 3 # number of columns (i.e. length of the numbers)
    m,n = size(y)

    # The unknown number
    x ~ filldist(DiscreteUniform(0,9),n)

    #=
    # Number of digits that has correct position and value.
    num_correct_positions = tzeros(m)
    for r in 1:m
        num_correct_positions[r] ~ Dirac(sum([y[r,i] == x[i] ? 1 : 0 for i in 1:n]))
    end

    # Number of digits that has correct value (wether in correct position or not)
    num_correct_values = tzeros(m)
    for r in 1:m
        num_correct_values[r] ~ Dirac(sum([y[r,j] == x[i] ? 1 : 0  for i in 1:n for j in 1:n]))
    end

    for r in 1:m
        hints[r,1] ~ Dirac(num_correct_positions[r])
        hints[r,2] ~ Dirac(num_correct_values[r])
    end
    =#

    # It's faster if we skip the two intermediate arrays num_correct_position and num_correct_values
    for r in 1:m
        hints[r,1] ~ Dirac(sum([y[r,i] == x[i] ? 1 : 0 for i in 1:n]))
        hints[r,2] ~ Dirac(sum([y[r,j] == x[i] ? 1 : 0  for i in 1:n for j in 1:n]))
    end

    return x
end

function run_number_lock2(y, hints)

    model = number_lock2(y,hints)

    # chns = sample(model, Prior(), 10_000)
    # chns = sample(model, MH(), 10_000)
    # chns = sample(model, PG(5), 10_000)
    chns = sample(model, SMC(), 10_000)
    # chns = sample(model, IS(), 10_000)

    display(chns)
    # display(plot(chns))

    chains_params = Turing.MCMCChains.get_sections(chns, :parameters)
    genq = generated_quantities(model, chains_params)
    show_var_dist_pct(genq)
    println()

end

# This instance: [0,5,2]
# The number combinations
y = [6 8 2
     6 4 5
     2 0 6
     7 3 8
     7 8 0]

# Hints
hints = [1 1
         0 1
         0 2
         0 0
         0 1]

run_number_lock2(y, hints)

# Another instance: https://twitter.com/natebrix/status/1418256269052321793
# The only changed is that the second row is
#   6 1 4
# instead of
#   6 4 5

# Solution:
# [0, 4, 2]	=>	10000 (1.0)#
#
y = [6 8 2
     6 1 4 # <- changed from 6 4 5
     2 0 6
     7 3 8
     7 8 0]

# Hints
hints = [1 1
         0 1
         0 2
         0 0
         0 1]
run_number_lock2(y, hints)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
module Enplot

using Colors: RGB, weighted_color_mean
using RecipesBase
using ColorSchemes
import Pixell: Enmap, pix2sky

# convert 0:255 integer to 0:1 Float64
rgb(c...) = RGB{Float64}((c./255)...)

# build a ColorScheme based on a list of colors and their locations in
# 0:1 colormap
function build_colorscheme(colors, locs)
    xs = collect(0:0.01:1)
    map(xs) do x
        i = min(findlast((<=)(x), locs), length(locs)-1)
        r = (x-locs[i])/(locs[i+1]-locs[i])
        weighted_color_mean(1-r, colors[i], colors[i+1])
    end |> ColorScheme
end

# here we define some common color schemes
cschemes = Dict{Symbol,ColorScheme}()
cschemes[:planck] = build_colorscheme([rgb(0, 0, 255), rgb(0, 215, 255), rgb(255, 237, 217), rgb(255, 180, 0), rgb(255, 75, 0), rgb(100, 0, 0)], [0, 0.332, 0.5, 0.664, 0.828, 1])

function register_colors!()
    for (k, v) in cschemes
        ColorSchemes.colorschemes[k] = v
    end
end

export register_colors!

# here we define common plot recipes
@recipe function f(imap::Enmap)
    seriestype   := :heatmap
    aspect_ratio := :equal
    xformatter   := x -> pix2sky_formatter(x, imap)
    yformatter   := x -> pix2sky_formatter(x, imap; ind=2)
    color          --> :planck
    xlim           --> (1, size(imap.data,1))  # not nice
    ylim           --> (1, size(imap.data,2))
    colorbar       --> true
    minorticks     --> true
    tick_direction --> :out
    grid           --> false
    size           --> (500,
                        500 * size(imap.data,2) / size(imap.data,1))
    framestyle     --> :box

    imap.data'
end

# format pixel index with its sky coordinate. ra: ind=1, dec: ind=2
function pix2sky_formatter(x, imap::Enmap; ind=1, digits=2)
    res = round(rad2deg(pix2sky(imap, float.([x, x]))[ind]), digits=digits)
    string(res)*"°"
end

end  # module


text:
using StaticArrays, LinearAlgebra
using Base.Iterators: flatten

import RecipesBase

export Dataset, AbstractDataset, minima, maxima
export minmaxima, columns

abstract type AbstractDataset{D, T} end

# Size:
@inline Base.length(d::AbstractDataset) = length(d.data)
@inline Base.size(d::AbstractDataset{D,T}) where {D,T} = (length(d.data), D)
@inline Base.size(d::AbstractDataset, i) = size(d)[i]
@inline Base.IteratorSize(d::AbstractDataset) = Base.HasLength()

# Itereting interface:
@inline Base.eachindex(D::AbstractDataset) = Base.OneTo(length(D.data))
@inline Base.iterate(d::AbstractDataset, state = 1) = iterate(d.data, state)
@inline Base.eltype(::Type{<:AbstractDataset{D, T}}) where {D, T} = SVector{D, T}

# 1D indexing  over the container elements:
@inline Base.getindex(d::AbstractDataset, i) = d.data[i]
@inline Base.lastindex(d::AbstractDataset) = length(d)
@inline Base.lastindex(d::AbstractDataset, k) = size(d)[k]
@inline Base.firstindex(d::AbstractDataset) = 1

# 2D indexing exactly like if the dataset was a matrix
# with each column a dynamic variable
@inline Base.getindex(d::AbstractDataset, i::Int, j::Int) = d.data[i][j]
@inline Base.getindex(d::AbstractDataset, i::Colon, j::Int) =
[d.data[k][j] for k in 1:length(d)]
@inline Base.getindex(d::AbstractDataset, i::Int, j::Colon) = d.data[i]
@inline Base.getindex(d::AbstractDataset, r::AbstractRange) = d.data[r]
# Indexing with ranges
@inline Base.getindex(d::AbstractDataset, i::AbstractRange, j::Int) =
[d.data[k][j] for k in i]
@inline Base.getindex(d::AbstractDataset, i::AbstractRange, ::Colon) =
Dataset([d[k] for k in i])
@inline Base.getindex(d::AbstractDataset, i::Int, j::AbstractRange) =
[d.data[i][k] for k in j]
# Indexing with boolean vectors
Base.getindex(d::AbstractDataset, i::AbstractVector{Bool}, j::Colon) =
Dataset(d[i])

# this function could be done generated
function Base.getindex(d::AbstractDataset{D,T}, i::AbstractVector{Int},
    j::AbstractVector{Int}) where {D,T}
    I = length(i)
    J = length(j)
    ret = zeros(T, J, I)
    for k=1:I
        for l=1:J
            ret[l,k] = d[i[k],j[l]]
        end
    end
    return Dataset(transpose(ret))
end

function mygetindex(d::AbstractDataset{D,T}, I::AbstractVector{Int},
    J::AbstractVector{Int}) where {D,T}

    L = length(J)
    sind::SVector{L, Int} = SVector{L, Int}(J)

    return Base.getindex(d, I, sind)
end

function Base.getindex(d::AbstractDataset{D,T},
    ::Colon, j::AbstractVector{Int}) where {D, T}
    return Base.getindex(d, 1:length(d), j)
end

"""
    columns(dataset) -> x, y, z, ...
Return the individual columns of the dataset.
"""
function columns end
@generated function columns(data::AbstractDataset{D, T}) where {D, T}
    gens = [:(data[:, $k]) for k=1:D]
    quote tuple($(gens...)) end
end

# Other commonly used functions:
Base.append!(d1::AbstractDataset, d2::AbstractDataset) = append!(d1.data, d2.data)
Base.push!(d::AbstractDataset, new_item) = push!(d.data, new_item)
@inline dimension(::AbstractDataset{D,T}) where {D,T} = D
@inline Base.eltype(d::AbstractDataset{D,T}) where {D,T} = T
import Base: ==
==(d1::AbstractDataset, d2::AbstractDataset) = d1.data == d2.data

"""
    Dataset{D, T} <: AbstractDataset{D,T}
A dedicated interface for datasets.
It contains *equally-sized datapoints* of length `D`, represented by `SVector{D, T}`.

When indexed with 1 index, a `dataset` is like a vector of datapoints.

When indexed with 2 indices it behaves like a matrix that has each of the columns be the
timeseries of each of the dynamic variables.

## Description of indexing
In the following let `i, j` be integers,  `typeof(data) <: AbstractDataset`
and `v1, v2` be `<: AbstractVector{Int}` (`v1, v2` could also be ranges).

* `data[i]` gives the `i`th datapoint (returns an `SVector`)
* `data[v1]` will return a vector of datapoints
* `data[v1, :]` using a `Colon` as a second index will return a `Dataset` of
  these points
* `data[:, j]` gives the `j`th variable timeseries, as `Vector`
* `data[v1, v2]` returns a `Dataset` with the appropriate entries (first indices
  being "time"/point index, while second being dynamic variables)
* `data[i, j]` value of the `j`th variable, at the `i`th timepoint

Use `Matrix(dataset)` or `Dataset(matrix)` to convert. It is assumed
that each *column* of the `matrix` is one dynamic variable.
If you have various timeseries vectors `x, y, z, ...` pass them like
`Dataset(x, y, z, ...)`. You can use `columns(dataset)` to obtain the reverse,
i.e. all columns of the dataset in a tuple.
"""
struct Dataset{D, T} <: AbstractDataset{D,T}
    data::Vector{SVector{D,T}}
end
# Empty dataset:
Dataset{D, T}() where {D,T} = Dataset(SVector{D,T}[])

# Identity constructor:
Dataset{D, T}(s::Dataset{D, T}) where {D,T} = s
Dataset(s::Dataset) = s

###########################################################################
# Dataset(Vectors of stuff)
###########################################################################
Dataset{1, T}(s::AbstractVector{T}) where {T<:Number} =
Dataset{1, T}(reinterpret(SVector{1, T}, s))

function Dataset(v::Vector{<:AbstractArray{T}}) where {T<:Number}
    D = length(v[1])
    L = length(v)
    data = Vector{SVector{D, T}}(undef, L)
    for i in 1:length(v)
        D != length(v[i]) && throw(ArgumentError(
        "All data-points in a Dataset must have same size"
        ))
        @inbounds data[i] = SVector{D,T}(v[i])
    end
    return Dataset{D, T}(data)
end

@generated function _dataset(vecs::Vararg{<:AbstractVector{T},D}) where {D, T}
    gens = [:(vecs[$k][i]) for k=1:D]

    quote
        L = length(vecs[1])
        data = Vector{SVector{$D, T}}(undef, L)
        for i in 1:L
            @inbounds data[i] = SVector{$D, T}($(gens...))
        end
        data
    end
end

function Dataset(vecs::Vararg{<:AbstractVector{T}}) where {T}
    return Dataset(_dataset(vecs...))
end

#####################################################################################
#                                Dataset <-> Matrix                                 #
#####################################################################################
function Base.Matrix{S}(d::AbstractDataset{D,T}) where {S, D, T}
    mat = Matrix{S}(undef, length(d), D)
    for j in 1:D
        for i in 1:length(d)
            @inbounds mat[i,j] = d.data[i][j]
        end
    end
    mat
end
Base.Matrix(d::AbstractDataset{D,T}) where {D, T} = Matrix{T}(d)

function Dataset(mat::AbstractMatrix{T}) where {T}
    N, D = size(mat)
    Dataset{D,T}(reshape(reinterpret(SVector{D,T}, vec(transpose(mat))), (N,)))
end

#####################################################################################
#                                   Pretty Printing                                 #
#####################################################################################
function Base.summary(d::Dataset{D, T}) where {D, T}
    N = length(d)
    return "$D-dimensional Dataset{$(T)} with $N points"
end

function matstring(d::AbstractDataset{D, T}) where {D, T}
    N = length(d)
    if N > 50
        mat = zeros(eltype(d), 50, D)
        for (i, a) in enumerate(flatten((1:25, N-24:N)))
            mat[i, :] .= d[a]
        end
    else
        mat = Matrix(d)
    end
    s = sprint(io -> show(IOContext(io, :limit=>true), MIME"text/plain"(), mat))
    s = join(split(s, '\n')[2:end], '\n')
    tos = summary(d)*"\n"*s
    return tos
end

Base.show(io::IO, ::MIME"text/plain", d::AbstractDataset) = print(io, matstring(d))
Base.show(io::IO, d::AbstractDataset) = print(io, summary(d))

#####################################################################################
#                                     Plotting                                      #
#####################################################################################

# Support the plotting of Datasets-as-Matrices
RecipesBase.@recipe f(::Type{T}, vals::T) where {T<:AbstractDataset} = Matrix(vals)

#####################################################################################
#                                 Minima and Maxima                                 #
#####################################################################################
"""
    minima(dataset)
Return an `SVector` that contains the minimum elements of each timeseries of the
dataset.
"""
function minima(data::AbstractDataset{D, T}) where {D, T<:Real}
    m = Vector(data[1])
    for point in data
        for i in 1:D
            if point[i] < m[i]
                m[i] = point[i]
            end
        end
    end
    return SVector{D,T}(m)
end

"""
    maxima(dataset)
Return an `SVector` that contains the maximum elements of each timeseries of the
dataset.
"""
function maxima(data::AbstractDataset{D, T}) where {D, T<:Real}
    m = Vector(data[1])
    for point in data
        for i in 1:D
            if point[i] > m[i]
                m[i] = point[i]
            end
        end
    end
    return SVector{D, T}(m)
end

"""
    minmaxima(dataset)
Return `minima(dataset), maxima(dataset)` without doing the computation twice.
"""
function minmaxima(data::AbstractDataset{D, T}) where {D, T<:Real}
    mi = Vector(data[1])
    ma = Vector(data[1])
    for point in data
        for i in 1:D
            if point[i] > ma[i]
                ma[i] = point[i]
            elseif point[i] < mi[i]
                mi[i] = point[i]
            end
        end
    end
    return SVector{D, T}(mi), SVector{D, T}(ma)
end

#####################################################################################
#                                     SVD                                           #
#####################################################################################
using LinearAlgebra
# SVD of Base seems to be much faster when the "long" dimension of the matrix
# is the first one, probably due to Julia's column major structure.
# This does not depend on using `svd` or `svdfact`, both give same timings.
# In fact it is so much faster, that it is *much* more worth it to
# use `Matrix(data)` instead of `reinterpret` in order to preserve the
# long dimension being the first.
"""
    svd(d::AbstractDataset) -> U, S, Vtr
Perform singular value decomposition on the dataset.
"""
function svd(d::AbstractDataset)
    F = svd(Matrix(d))
    return F[:U], F[:S], F[:Vt]
end


text:
import Nabla: ∇, compute_Dv, approximate_Dv, compute_Dv_update
@testset "Finite-difference estimates of sensitivities" begin
    # Define a dummy test function.
    foo(x::∇Scalar) = 5x
    foo(x::∇Scalar, y::∇Scalar) = 5x + 6y
    foo(x::Vector{<:∇Scalar}, y::Vector{<:∇Scalar}) = 10x + 11y
    foo(x::Matrix{<:∇Scalar}, y::Matrix{<:∇Scalar}) = 15x + 16y

    # Create sensitivity intercepts.
    @explicit_intercepts foo Tuple{∇Scalar}
    @explicit_intercepts foo Tuple{∇Scalar, ∇Scalar}
    @explicit_intercepts foo Tuple{Vector{<:∇Scalar}, Vector{<:∇Scalar}}
    @explicit_intercepts foo Tuple{Matrix{<:∇Scalar}, Matrix{<:∇Scalar}}

    # Define sensitivity implementations.
    _Vec = Vector{<:∇Scalar}
    _Mat = Matrix{<:∇Scalar}
    Nabla.∇(::typeof(foo), ::Type{Arg{1}}, p, z, z̄, x::∇Scalar) = 5z̄
    Nabla.∇(::typeof(foo), ::Type{Arg{1}}, p, z, z̄, x::∇Scalar, y::∇Scalar) = 5z̄
    Nabla.∇(::typeof(foo), ::Type{Arg{1}}, p, z, z̄, x::_Vec, y::_Vec) = 10z̄
    Nabla.∇(::typeof(foo), ::Type{Arg{1}}, p, z, z̄, x::_Mat, y::_Mat) = 15z̄
    Nabla.∇(::typeof(foo), ::Type{Arg{2}}, p, z, z̄, x::∇Scalar, y::∇Scalar) = 6z̄
    Nabla.∇(::typeof(foo), ::Type{Arg{2}}, p, z, z̄, x::_Vec, y::_Vec) = 11z̄
    Nabla.∇(::typeof(foo), ::Type{Arg{2}}, p, z, z̄, x::_Mat, y::_Mat) = 16z̄

    # Check that Nabla and finite differencing yield the correct results for scalars.
    @test compute_Dv(foo, 1.0, 1.0, 1.0) ≈ 5
    @test compute_Dv_update(foo, 1.0, 1.0, 1.0) ≈ 5
    @test approximate_Dv(foo, 1.0, 1.0, 1.0) ≈ 5
    @test compute_Dv(foo, 1.0, (1.0, 1.0), (1.0, 1.0)) ≈ 5 + 6
    @test compute_Dv_update(foo, 1.0, (1.0, 1.0), (1.0, 1.0)) ≈ 5 + 6
    @test approximate_Dv(foo, 1.0, (1.0, 1.0), (1.0, 1.0)) ≈ 5 + 6
    @test compute_Dv(foo, 5.0, (1.0, 1.0), (1.0, 1.0)) ≈ (5 + 6) * 5.0
    @test compute_Dv_update(foo, 5.0, (1.0, 1.0), (1.0, 1.0)) ≈ (5 + 6) * 5.0
    @test approximate_Dv(foo, 5.0, (1.0, 1.0), (1.0, 1.0)) ≈ (5 + 6) * 5.0
    @test compute_Dv(foo, 1.0, (10.0, 5.0), (1.0, 1.0)) ≈ 5 + 6
    @test compute_Dv_update(foo, 1.0, (10.0, 5.0), (1.0, 1.0)) ≈ 5 + 6
    @test approximate_Dv(foo, 1.0, (10.0, 5.0), (1.0, 1.0)) ≈ 5 + 6
    @test compute_Dv(foo, 1.0, (1.0, 1.0), (7.5, 6.3)) ≈ 5 * 7.5 + 6 * 6.3
    @test compute_Dv_update(foo, 1.0, (1.0, 1.0), (7.5, 6.3)) ≈ 5 * 7.5 + 6 * 6.3
    @test approximate_Dv(foo, 1.0, (1.0, 1.0), (7.5, 6.3)) ≈ 5 * 7.5 + 6 * 6.3
    @test compute_Dv(foo, 5.0, (10.0, 5.0), (7.5, 6.3)) ≈ (5 * 7.5 + 6 * 6.3) * 5
    @test compute_Dv_update(foo, 5.0, (10.0, 5.0), (7.5, 6.3)) ≈ (5 * 7.5 + 6 * 6.3) * 5
    @test approximate_Dv(foo, 5.0, (10.0, 5.0), (7.5, 6.3)) ≈ (5 * 7.5 + 6 * 6.3) * 5

    # Check that finite differencing yields approximately correct results for vectors.
    N = 3
    z̄, x, y, vx, vy = ones.((N, N, N, N, N))
    @test compute_Dv(foo, z̄, (x, y), (vx, vy)) ≈
        sum(z̄ .* (UniformScaling(10) * vx)) + sum(z̄ .* (UniformScaling(11) * vy))
    @test compute_Dv_update(foo, z̄, (x, y), (vx, vy)) ≈
        sum(z̄ .* (UniformScaling(10) * vx)) + sum(z̄ .* (UniformScaling(11) * vy))
    @test approximate_Dv(foo, z̄, (x, y), (vx, vy)) ≈
        sum(z̄ .* (UniformScaling(10) * vx)) + sum(z̄ .* (UniformScaling(11) * vy))
    z̄ = randn(N)
    @test compute_Dv(foo, z̄, (x, y), (vx, vy)) ≈
        sum(z̄ .* (UniformScaling(10) * vx)) + sum(z̄ .* (UniformScaling(11) * vy))
    @test compute_Dv_update(foo, z̄, (x, y), (vx, vy)) ≈
        sum(z̄ .* (UniformScaling(10) * vx)) + sum(z̄ .* (UniformScaling(11) * vy))
    @test approximate_Dv(foo, z̄, (x, y), (vx, vy)) ≈
        sum(z̄ .* (UniformScaling(10) * vx)) + sum(z̄ .* (UniformScaling(11) * vy))
    v = randn(N)
    @test compute_Dv(foo, z̄, (x, y), (vx, vy)) ≈
        sum(z̄ .* (UniformScaling(10) * vx)) + sum(z̄ .* (UniformScaling(11) * vy))
    @test compute_Dv_update(foo, z̄, (x, y), (vx, vy)) ≈
        sum(z̄ .* (UniformScaling(10) * vx)) + sum(z̄ .* (UniformScaling(11) * vy))
    @test approximate_Dv(foo, z̄, (x, y), (vx, vy)) ≈
        sum(z̄ .* (UniformScaling(10) * vx)) + sum(z̄ .* (UniformScaling(11) * vy))

    # Check that finite differencing yields approximately correct results for matrices.
    M, N = 3, 2
    Z̄, X, Y, VX, VY = ones.((M, M, M, M, M), (N, N, N, N, N))
    @test approximate_Dv(foo, Z̄, (X, Y), (VX, VY)) ≈
        compute_Dv(foo, Z̄, (X, Y), (VX, VY))
    @test compute_Dv(foo, Z̄, (X, Y), (VX, VY)) ≈
        compute_Dv_update(foo, Z̄, (X, Y), (VX, VY))
    Ȳ = randn(M, N)
    @test approximate_Dv(foo, Z̄, (X, Y), (VX, VY)) ≈
        compute_Dv(foo, Z̄, (X, Y), (VX, VY))
    @test compute_Dv(foo, Z̄, (X, Y), (VX, VY)) ≈
        compute_Dv_update(foo, Z̄, (X, Y), (VX, VY))
end


text:
println("tell me sth")
test_me=100
typeof(test_me)
tsstr=""""abba" test chars in string"""
println(tsstr)
typeof('a')

myphonebook = Dict("Jenny" => "867-5309", "Ghostbusters" => "555-2368")
myphonebook["Jenny"]
rand(4,3)
n=10
m=10
for i in 1:n
    println("testing $i")
end

C=[ i+j for i in 1:n,j in 1:m]

max_stars_repo_name: stepinski/julia-in

max_issues_repo_name: stepinski/julia-init

max_forks_repo_name: stepinski/julia-init


text:
using Gadfly
using Distributions

"""
Bimodal Spike Demo
"""
function drawBimodalDemo()
    μ = [.0 2.]
    Σ = [1. .05]
    w = [.5 .5]
    xs = collect(-2:.01:μ[2]*2)

    n1 = Normal(μ[1], Σ[1])
    n2 = Normal(μ[2], Σ[2])
    p = w[1]*pdf(n1, xs) + w[2]*pdf(n2, xs)

    μ = mean(xs .* p)

    pp = plot(x=xs, y=p, Geom.line, Theme(default_color=colorant"black", line_width=3px),
              xintercept=[μ], Geom.vline(color=colorant"red", size=3px))
end


text:
data_test_cases = [
    [0x42, 0x43, 0x02, 0x00, 0xa1, 0x4c],
    [0x02, 0x03, 0x00, 0x00],
    [0x00, 0x01, 0x04, 0x00, 0x01, 0x02, 0x03, 0x04,
    0xff, 0xf1, 0x01, 0x00, 0xff],
]

@testset "Is valid data" begin
    test_valid(v) = is_valid_extra_data(pointer(v), UInt16(length(v)))
    for data in data_test_cases
        @test test_valid(data)
        data[2] = 0x00
        @test !test_valid(data)
        data[2] = 0xff
        push!(data, 0x00)
        @test !test_valid(data)
        pop!(data)
        @test !test_valid(data[1:end-1])
        data = empty!(copy(data))
        @test test_valid(data)
    end
end

@testset "Parse fields" begin
    test_parse(v) = LibDeflate.parse_fields(pointer(v), UInt32(1), UInt16(length(v)))
    for data in data_test_cases
        # We merely test it doesn't fail
        @test test_parse(data) !== nothing
        data[2] = 0x00
        @test test_parse(data) == LibDeflateErrors.gzip_bad_extra
        data[2] = 0xa0
        push!(data, 0x00)
        @test test_parse(data) == LibDeflateErrors.gzip_extra_too_long
        pop!(data)
        @test test_parse(data[1:end-1]) == LibDeflateErrors.gzip_extra_too_long
        data = empty!(copy(data))
        @test test_parse(data) !== nothing
    end
end

header_data = UInt8[
    # header
    0x1f, 0x8b, 0x08, 0x1e, 0xb3, 0x2c, 0x51, 0x60, 0xff, 0x00,

    # Extra data
    0x0a, 0x00, 0x42, 0x43, 0x02, 0x00, 0xa1, 0x4c,
    0x02, 0x03, 0x00, 0x00,

    # Filename: "filename.fna"
    0x66, 0x69, 0x6c, 0x65, 0x6e, 0x61, 0x6d, 0x65, 0x2e, 0x66, 0x6e, 0x61, 0x00,

    # Complicated unicode comment "αβ学中文"
    0xce, 0xb1, 0xce, 0xb2, 0xe5, 0xad, 0xa6, 0xe4, 0xb8, 0xad, 0xe6, 0x96, 0x87, 0x00,

    # CRC16
    0x78, 0x18
]

function test_header_example(data::Vector{UInt8}, header::LibDeflate.GzipHeader)
    @test header.mtime == 0x60512cb3
    @test length(header.extra) == 2
    @test first(header.extra).tag == (0x42, 0x43)
    @test first(header.extra).data == 0x00000011:0x00000012
    @test last(header.extra).tag == (0x02, 0x03)
    @test last(header.extra).data === nothing # empty field
    @test String(data[header.filename]) == "filename.fna"
    @test String(data[header.comment]) == "αβ学中文"
    true
end

@testset "Parse header" begin
    header = unsafe_parse_gzip_header(pointer(header_data), UInt(51))[2]
    test_header_example(header_data, header)

    header = parse_gzip_header(header_data)[2]
    test_header_example(header_data, header)

    header = unsafe_parse_gzip_header(pointer(header_data), UInt(51), LibDeflate.GzipExtraField[])[2]
    test_header_example(header_data, header)

    header_data[end-2] = 0x01
    @test unsafe_parse_gzip_header(pointer(header_data), UInt(51)) == LibDeflateErrors.gzip_string_not_null_terminated
    header_data[end-2] = 0x00
end

test_data = [
    "",
    "Abracadabra!",
    "En af dem der red med fane",
    rand(UInt8, 1000)
]

test_comment = "This is a comment"
test_filename = "testfile.foo"

@testset "Compression" begin
    outdata = zeros(UInt8, 1250)
    compressor = Compressor()
    for data in test_data
        n_bytes = unsafe_gzip_compress!(
            compressor, pointer(outdata), UInt(length(outdata)),
            pointer(data), UInt(sizeof(data)),
            LibDeflate.ReadableMemory(test_comment), LibDeflate.ReadableMemory(test_filename),
            nothing, true
        )
        decompressed = transcode(GzipDecompressor, outdata[1:n_bytes])
        @test decompressed == Vector{UInt8}(data)

        gzip_compress!(
            compressor, outdata, data;
            comment=test_comment, filename=test_filename, extra=data_test_cases[1], header_crc=false
        )
        decompressed = transcode(GzipDecompressor, outdata)
        @test decompressed == Vector{UInt8}(data)

        # Resize for next iteration
        resize!(outdata, 1250)
    end
end

complex_test_case = vcat(header_data, UInt8[
    # Data: compressed "Abracadabra"
    0x01, 0x0b, 0x00, 0xf4, 0xff, 0x41, 0x62, 0x72, 0x61, 0x63, 0x61, 0x64, 0x61, 0x62, 0x72, 0x61,

    # CRC32
    0x60, 0x76, 0x76, 0x91,

    # isize
    0x0b, 0x00, 0x00, 0x00
])

@testset "Decompression" begin
    decompressor = Decompressor()
    outdata = zeros(UInt8, 5) # begin with small buffer, let it resize
    for data in test_data
        compressed = transcode(GzipCompressor, data)

        result = unsafe_gzip_decompress!(
            decompressor, outdata, UInt(1001),
            pointer(compressed), UInt(length(compressed)),
        )
        @test outdata[1:result.len] == Vector{UInt8}(data)

        empty!(outdata)
        result = gzip_decompress!(decompressor, outdata, compressed)
        @test outdata[1:result.len] == Vector{UInt8}(data)

        resize!(outdata, 5)
    end

    # Hard test case
    res = gzip_decompress!(decompressor, outdata, complex_test_case)
    test_header_example(complex_test_case, res.header)
    @test res.len == 11
end


text:
struct KrylovJL{AT,bT,A,K} <: SciMLLinearSolveAlgorithm
    solver::Function
    A::AT
    b::bT
    args::A
    kwargs::K
end

function KrylovJL(A, b, args...; solver = gmres, kwargs...)
    return KrylovJL(solver, A, b, args, kwargs)
end

function SciMLBase.solve(prob::LinearProblem, alg::KrylovJL, args...; kwargs...)
    @unpack A, b, p = prob
    x, stats = alg.solver(A, b, args...; kwargs...)
    resid = mul!(similar(b), A, x) - b
    retcode = stats.solved ? :Success : :Failure
    return SciMLBase.build_solution(prob, alg, x, resid; retcode = retcode)
end


text:
using Test, ANSIColoredPrinters

function repr_html(printer::HTMLPrinter)
    result = repr("text/html", printer)
    take!(printer.buf)
    return result
end

@testset "plain" begin
    buf = IOBuffer()
    printer = HTMLPrinter(buf)

    print(buf, "This is a plain text.")
    @test repr_html(printer) == "
This is a plain text.
"
end

@testset "escape" begin
    buf = IOBuffer()
    printer = HTMLPrinter(buf)

    print(buf, "\"HTMLWriter\" uses '
' & '' elements.")
    result = repr_html(printer)
    @test result == "
"HTMLWriter" uses '<pre>' & " *
                    "'<span>' elements.
"
end

@testset "single modification" begin
    buf = IOBuffer()
    printer = HTMLPrinter(buf)

    @testset "bold/faint" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[1m", " Bold ")
        print(buf, "\e[2m", " Faint ")
        print(buf, "\e[22m", " Normal ")
        result = repr_html(printer)
        @test result == """
 Normal  Bold """ *
                        """ Faint  Normal
"""
    end

    @testset "fg/bg" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[36m", " CyanFG ")
        print(buf, "\e[39;44m", " BlueBG ")
        print(buf, "\e[49m", " Normal ")
        result = repr_html(printer)
        @test result == """
 Normal  CyanFG """ *
                        """ BlueBG  Normal
"""
    end

    @testset "blink" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[5m", " Blink ")
        print(buf, "\e[6m", " RapidBlink ")
        print(buf, "\e[25m", " Normal ")
        result = repr_html(printer)
        @test result == """
 Normal  Blink """ *
                        """ RapidBlink  Normal
"""
    end
end

@testset "nested modification" begin
    buf = IOBuffer()
    printer = HTMLPrinter(buf)

    @testset "bold/italic" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[1m", " Bold ")
        print(buf, "\e[3m", " Bold-Italic ")
        print(buf, "\e[23m", " Bold ")
        print(buf, "\e[0m", " Normal ")
        result = repr_html(printer)
        @test result == """
 Normal  Bold """ *
                        """ Bold-Italic """ *
                        """ Bold  Normal
"""

        print(buf, "\e[3m", " Italic ")
        print(buf, "\e[1m", " Bold-Italic ")
        print(buf, "\e[22m", " Italic ")
        print(buf, "\e[0m", " Normal ")
        result = repr_html(printer)
        @test result == """
 Italic """ *
                        """ Bold-Italic """ *
                        """ Italic  Normal
"""
    end

    @testset "bold/fg" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[1m", " Bold ")
        print(buf, "\e[36m", " Bold-CyanFG ")
        print(buf, "\e[39m", " Bold ")
        print(buf, "\e[0m", " Normal ")
        result = repr_html(printer)
        @test result == """
 Normal  Bold """ *
                        """ Bold-CyanFG """ *
                        """ Bold  Normal
"""
    end

    @testset "strike/blink" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[9m", " Strike ")
        print(buf, "\e[5m", " Strike-Blink ")
        print(buf, "\e[25m", " Strike ")
        print(buf, "\e[29m", " Normal ")
        result = repr_html(printer)
        @test result == """
 Normal  Strike """ *
                        """ Strike-Blink """ *
                        """ Strike  Normal
"""
    end

    @testset "fg/bg" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[96m", " LightCyanFG ")
        print(buf, "\e[45m", " LightCyanFG-MagentaBG ")
        print(buf, "\e[49m", " LightCyanFG ")
        print(buf, "\e[39m", " Normal ")
        result = repr_html(printer)
        @test result == """
 Normal  LightCyanFG """ *
                        """ LightCyanFG-MagentaBG """ *
                        """ LightCyanFG  Normal
"""
    end
end

@testset "overlapped modification" begin
    buf = IOBuffer()
    printer = HTMLPrinter(buf)

    @testset "bold/fg" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[1m", " Bold ")
        print(buf, "\e[36m", " Bold-CyanFG ")
        print(buf, "\e[22;36m", " CyanFG ")
        print(buf, "\e[39m", " Normal ")
        result = repr_html(printer)
        @test result == """
 Normal  Bold """ *
                        """ Bold-CyanFG """ *
                        """ CyanFG  Normal
"""
    end

    @testset "fg/bg" begin
        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[96m", " LightCyanFG ")
        print(buf, "\e[45m", " LightCyanFG-MagentaBG ")
        print(buf, "\e[39m", " MagentaBG ")
        print(buf, "\e[49m", " Normal ")
        result = repr_html(printer)
        @test result == """
 Normal  LightCyanFG """ *
                        """ LightCyanFG-MagentaBG """ *
                        """ MagentaBG  Normal
"""
    end
end

@testset "force reset" begin
    buf = IOBuffer()
    printer = HTMLPrinter(buf)

    print(buf, "\e[7m", " Invert ")
    print(buf, "\e[8m", " Conceal ")
    result = repr_html(printer)
    @test result == """
 Invert  Conceal """ *
                    """
"""
end

@testset "256 colors" begin
    buf = IOBuffer()
    printer = HTMLPrinter(buf)

    print(buf, "\e[38;5;4m", " Blue(FG) ")
    print(buf, "\e[48;5;16m", " #000(BG) ")
    print(buf, "\e[38;5;110m", " #87afd7(FG) ")
    print(buf, "\e[48;5;255m", " #eee(BG) ")

    result = repr_html(printer)
    @test result == """
 Blue(FG) """ *
                    """ #000(BG) """ *
                    """""" *
                    """ #87afd7(FG) """ *
                    """""" *
                    """ #eee(BG) """ *
                    """
"""
end

@testset "24-bit colors" begin
    buf = IOBuffer()
    printer = HTMLPrinter(buf)

    print(buf, "\e[38;2;0;128;255m", " #0080ff(FG) ")
    print(buf, "\e[48;2;238;238;238m", " #eee(BG) ")
    print(buf, "\e[38;2;170;187;204m", " #abc(FG) ")

    result = repr_html(printer)
    @test result == """
 #0080ff(FG) """ *
                    """ #eee(BG) """ *
                    """""" *
                    """ #abc(FG) """ *
                    """
"""
end

@testset "callback" begin
    @testset "use default" begin
        counter = 0
        function cb(io, printer, tag, attrs)
            startswith(tag, "/") && return nothing
            push!(attrs, :id => tag * string(counter))
            counter += 1
            return nothing
        end

        buf = IOBuffer()
        printer = HTMLPrinter(buf, root_class = "root test", root_tag = "code", callback = cb)

        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[38;5;255m", " #eee(FG) ")
        print(buf, "\e[0m", " Normal ")

        result = repr_html(printer)
        @test result == """ Normal """ *
                        """""" *
                        """ #eee(FG)  Normal """
    end

    @testset "prevent default" begin
        dom = Tuple[(:rootnode, Tuple[])]
        function cb(io, printer, tag, attrs)
            text = String(take!(io))
            parent = dom[end]
            children = parent[end]
            isempty(text) || push!(children, (:textnode, text))

            if startswith(tag, "/")
                pop!(dom)
            else
                parent = (Symbol(tag), attrs, Tuple[])
                push!(children, parent)
                push!(dom, parent)
            end
            return true
        end

        buf = IOBuffer()
        printer = HTMLPrinter(buf, callback = cb)

        print(buf, "\e[0m", " Normal ")
        print(buf, "\e[38;5;255m", " #eee(FG) ")
        print(buf, "\e[0m", " Normal ")

        result = repr_html(printer)
        @test result == ""
        @test dom[1] ==
                (:rootnode, Tuple[
                    (:pre, Dict{Symbol,String}(), Tuple[
                        (:textnode, " Normal "),
                        (:span, Dict(:class => "sgr38_5", :style => "color:#eee"), Tuple[
                            (:textnode, " #eee(FG) ")
                        ]),
                        (:textnode, " Normal ")
                    ])
                ])
    end
end


text:
include("Include.jl")

function updateCholesky(A,zV,CCOV,AVG_P_SUCC,PTHRESH)

  if (AVG_P_SUCC<PTHRESH)
		CA = sqrt(1-CCOV);
		F = 1 + ((1-CA^2)*(norm(zV))^2)/(CA^2);
		A = CA*A+CA/(norm(zV)^2)*(sqrt(F)-1)*A*zV*transpose(zV);
	else
		A = A;
	end

  return A;
end

function updateStepSize(LAMBDA_SUCC,AVG_P_SUCC,P_SUCC_TARGET,CP,D,SIGMA)

	# Compute AVG_P_SUCC -
	AVG_P_SUCC = (1-CP)*AVG_P_SUCC+CP*LAMBDA_SUCC;

	# Compute the step-size -
	F = (1/D)*(AVG_P_SUCC-(P_SUCC_TARGET/(1-P_SUCC_TARGET))*(1-AVG_P_SUCC));
	SIGMA = SIGMA*exp(F);

	if (SIGMA<= 1e-5)
		SIGMA = 1e-5;
	end

  return (SIGMA,AVG_P_SUCC);
end

function sample_function(parameter_guess,kVI,lower_bound_array,upper_bound_array,NPARAMETERS,A,SIGMA)

  # Generate a new vector -
  zV = randn(NPARAMETERS);

  # Compute the new perturbed parameter vector -
  perturbed_parameter_array = parameter_guess.*(1 + SIGMA*A*zV);
  corrected_parameter_array = parameter_bounds_function(perturbed_parameter_array,lower_bound_array,upper_bound_array);

  # return -
  return (corrected_parameter_array,zV)
end

function bounds_function(number_of_parameters)

  # Setup my upper bound, and lower bounds on parameters -
  lower_bound_array = zeros(number_of_parameters)
  upper_bound_array = [
    4.0   ;   # "n_gene_AP1_gene_AhR"	;	# 1
    500.0 ;   # "K_gene_AP1_gene_AhR"	;	# 2
    4.0   ;   # "n_gene_AP1_gene_PU1"	;	# 3
    500.0 ;   # "K_gene_AP1_gene_PU1"	;	# 4
    4.0   ;   # "n_gene_AP1_gene_PPARg"	;	# 5
    500.0 ;   # "K_gene_AP1_gene_PPARg"	;	# 6
    4.0   ;   # "n_gene_AhR_gene_Trigger"	;	# 7
    500.0 ;   # "K_gene_AhR_gene_Trigger"	;	# 8
    4.0   ;  # "n_gene_CD11b_gene_PU1_gene_cRAF"	;	# 9
    500.0 ;  # "K_gene_CD11b_gene_PU1_gene_cRAF"	;	# 10
    4.0   ;  # "n_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"	;	# 11
    500.0 ;  # "K_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"	;	# 12
    4.0   ;  # "n_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"	;	# 13
    500.0 ;  # "K_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"	;	# 14
    4.0   ;  # "n_gene_CEBPa_gene_Trigger"	;	# 15
    500.0 ;  # "K_gene_CEBPa_gene_Trigger"	;	# 16
    4.0   ;  # "n_gene_CEBPa_gene_PPARg"	;	# 17
    500.0 ;  # "K_gene_CEBPa_gene_PPARg"	;	# 18
    4.0   ;  # "n_gene_CEBPa_gene_CEBPa"	;	# 19
    500.0 ;  # "K_gene_CEBPa_gene_CEBPa"	;	# 20
    4.0   ;  # "n_gene_CEBPa_gene_GFI1"	;	# 21
    1000.0 ;  # "K_gene_CEBPa_gene_GFI1"	;	# 22
    4.0   ;  # "n_gene_E2F_gene_E2F"	;	# 23
    500.0 ;  # "K_gene_E2F_gene_E2F"	;	# 24
    4.0   ;  # "n_gene_E2F_gene_PPARg"	;	# 25
    500.0 ;  # "K_gene_E2F_gene_PPARg"	;	# 26
    4.0   ;  # "n_gene_E2F_gene_CEBPa"	;	# 27
    500.0 ;  # "K_gene_E2F_gene_CEBPa"	;	# 28
    4.0   ;  # "n_gene_E2F_gene_GFI1"	;	# 29
    500.0 ;  # "K_gene_E2F_gene_GFI1"	;	# 30
    4.0   ;  # "n_gene_E2F_gene_cRAF"	;	# 31
    500.0 ;  # "K_gene_E2F_gene_cRAF"	;	# 32
    4.0   ;  # "n_gene_EGR1_gene_Trigger"	;	# 33
    500.0 ;  # "K_gene_EGR1_gene_Trigger"	;	# 34
    4.0   ;  # "n_gene_EGR1_gene_PU1"	;	# 35
    500.0 ;  # "K_gene_EGR1_gene_PU1"	;	# 36
    4.0   ;  # "n_gene_EGR1_gene_PPARg"	;	# 37
    500.0 ;  # "K_gene_EGR1_gene_PPARg"	;	# 38
    4.0   ;  # "n_gene_EGR1_gene_GFI1"	;	# 39
    500.0 ;  # "K_gene_EGR1_gene_GFI1"	;	# 40
    4.0   ;  # "n_gene_GFI1_gene_CEBPa"	;	# 41
    500.0 ;  # "K_gene_GFI1_gene_CEBPa"	;	# 42
    4.0   ;  # "n_gene_GFI1_gene_EGR1"	;	# 43
    500.0 ;  # "K_gene_GFI1_gene_EGR1"	;	# 44
    4.0   ;  # "n_gene_IRF1_gene_Trigger"	;	# 45
    500.0 ;  # "K_gene_IRF1_gene_Trigger"	;	# 46
    4.0   ;  # "n_gene_IRF1_gene_AhR"	;	# 47
    500.0 ;  # "K_gene_IRF1_gene_AhR"	;	# 48
    4.0   ;  # "n_gene_IRF1_gene_PPARg"	;	# 49
    500.0 ;  # "K_gene_IRF1_gene_PPARg"	;	# 50
    4.0   ;  # "n_gene_OCT1_gene_PPARg"	;	# 51
    500.0 ;  # "K_gene_OCT1_gene_PPARg"	;	# 52
    4.0   ;  # "n_gene_OCT4_gene_Trigger"	;	# 53
    500.0 ;  # "K_gene_OCT4_gene_Trigger"	;	# 54
    4.0   ;  # "n_gene_OCT4_gene_AhR"	;	# 55
    500.0 ;  # "K_gene_OCT4_gene_AhR"	;	# 56
    4.0   ;  # "n_gene_OCT4_gene_cRAF"	;	# 57
    500.0 ;  # "K_gene_OCT4_gene_cRAF"	;	# 58
    4.0   ;  # "n_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"	;	# 59
    500.0 ;  # "K_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"	;	# 60
    4.0   ;  # "n_gene_P21_gene_GFI1"	;	# 61
    500.0 ;  # "K_gene_P21_gene_GFI1"	;	# 62
    4.0   ;  # "n_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"	;	# 63
    500.0 ;  # "K_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"	;	# 64
    4.0   ;  # "n_gene_P47Phox_gene_PPARg"	;	# 65
    500.0 ;  # "K_gene_P47Phox_gene_PPARg"	;	# 66
    4.0   ;  # "n_gene_PPARg_gene_Trigger"	;	# 67
    500.0 ;  # "K_gene_PPARg_gene_Trigger"	;	# 68
    4.0   ;  # "n_gene_PPARg_gene_CEBPa"	;	# 69
    500.0 ;  # "K_gene_PPARg_gene_CEBPa"	;	# 70
    4.0   ;  # "n_gene_PPARg_gene_EGR1"	;	# 71
    500.0 ;  # "K_gene_PPARg_gene_EGR1"	;	# 72
    4.0   ;  # "n_gene_PPARg_gene_PU1"	;	# 73
    500.0 ;  # "K_gene_PPARg_gene_PU1"	;	# 74
    4.0   ;  # "n_gene_PPARg_gene_AP1"	;	# 75
    500.0 ;  # "K_gene_PPARg_gene_AP1"	;	# 76
    4.0   ;  # "n_gene_PU1_gene_Trigger"	;	# 77
    500.0 ;  # "K_gene_PU1_gene_Trigger"	;	# 78
    4.0   ;  # "n_gene_PU1_gene_CEBPa"	;	# 79
    500.0 ;  # "K_gene_PU1_gene_CEBPa"	;	# 80
    4.0   ;  # "n_gene_PU1_gene_PU1"	;	# 81
    1000.0 ;  # "K_gene_PU1_gene_PU1"	;	# 82
    4.0   ;  # "n_gene_PU1_gene_AP1"	;	# 83
    500.0 ;  # "K_gene_PU1_gene_AP1"	;	# 84
    4.0   ;  # "n_gene_PU1_gene_OCT1"	;	# 85
    500.0 ;  # "K_gene_PU1_gene_OCT1"	;	# 86
    4.0   ;  # "n_gene_PU1_gene_AhR"	;	# 87
    500.0 ;  # "K_gene_PU1_gene_AhR"	;	# 88
    4.0   ;  # "n_gene_PU1_gene_GFI1"	;	# 89
    500.0 ;  # "K_gene_PU1_gene_GFI1"	;	# 90
    100.0 ; # "W_gene_AP1_RNAP"	;	# 91
    100.0 ;  # "W_gene_AP1_gene_AhR"	;	# 92
    100.0 ;  # "W_gene_AP1_gene_PU1"	;	# 93
    100.0 ;  # "W_gene_AP1_gene_PPARg"	;	# 94
    100.0 ;  # "W_gene_AhR_RNAP"	;	# 95
    100.0 ;  # "W_gene_AhR_gene_Trigger"	;	# 96
    100.0 ;  # "W_gene_CD11b_RNAP"	;	# 97
    100.0 ;  # "W_gene_CD11b_gene_PU1_gene_cRAF"	;	# 98
    100.0 ;  # "W_gene_CD14_RNAP"	;	# 99
    0.0   ;  # "W_gene_CD14_gene_PPARg_gene_CEBPa_gene_EGR1_gene_cRAF"	;	# 100
    100.0 ;  # "W_gene_CD38_RNAP"	;	# 101
    100.0 ;  # "W_gene_CD38_gene_IRF1_gene_PPARg_gene_Trigger_gene_cRAF"	;	# 102
    100.0 ;  # "W_gene_CEBPa_RNAP"	;	# 103
    100.0 ;  # "W_gene_CEBPa_gene_Trigger"	;	# 104
    100.0 ;  # "W_gene_CEBPa_gene_PPARg"	;	# 105
    100.0 ;  # "W_gene_CEBPa_gene_CEBPa"	;	# 106
    100.0 ;  # "W_gene_CEBPa_gene_GFI1"	;	# 107
    100.0 ;  # "W_gene_E2F_RNAP"	;	# 108
    100.0 ;  # "W_gene_E2F_gene_E2F"	;	# 109
    100.0 ;  # "W_gene_E2F_gene_PPARg"	;	# 110
    100.0 ;  # "W_gene_E2F_gene_CEBPa"	;	# 111
    100.0 ;  # "W_gene_E2F_gene_GFI1"	;	# 112
    100.0 ;  # "W_gene_E2F_gene_cRAF"	;	# 113
    100.0 ;  # "W_gene_EGR1_RNAP"	;	# 114
    100.0 ;  # "W_gene_EGR1_gene_Trigger"	;	# 115
    100.0 ;  # "W_gene_EGR1_gene_PU1"	;	# 116
    100.0 ;  # "W_gene_EGR1_gene_PPARg"	;	# 117
    100.0 ;  # "W_gene_EGR1_gene_GFI1"	;	# 118
    100.0 ;  # "W_gene_GFI1_RNAP"	;	# 119
    100.0 ;  # "W_gene_GFI1_gene_CEBPa"	;	# 120
    100.0 ;  # "W_gene_GFI1_gene_EGR1"	;	# 121
    100.0 ;  # "W_gene_IRF1_RNAP"	;	# 122
    100.0 ;  # "W_gene_IRF1_gene_Trigger"	;	# 123
    100.0 ;  # "W_gene_IRF1_gene_AhR"	;	# 124
    100.0 ;  # "W_gene_IRF1_gene_PPARg"	;	# 125
    100.0 ;  # "W_gene_OCT1_RNAP"	;	# 126
    100.0 ;  # "W_gene_OCT1_gene_PPARg"	;	# 127
    100.0 ;  # "W_gene_OCT4_RNAP"	;	# 128
    100.0 ;  # "W_gene_OCT4_gene_Trigger"	;	# 129
    100.0 ;  # "W_gene_OCT4_gene_AhR"	;	# 130
    100.0 ;  # "W_gene_OCT4_gene_cRAF"	;	# 131
    100.0 ;  # "W_gene_P21_RNAP"	;	# 132
    100.0 ;  # "W_gene_P21_gene_Trigger_gene_AP1_gene_PPARg_gene_PU1_gene_IRF1_gene_CEBPa_gene_cRAF"	;	# 133
    100.0 ;  # "W_gene_P21_gene_GFI1"	;	# 134
    100.0 ;  # "W_gene_P47Phox_RNAP"	;	# 135
    100.0 ;  # "W_gene_P47Phox_gene_PU1_gene_CEBPa_gene_cRAF"	;	# 136
    100.0 ;  # "W_gene_P47Phox_gene_PPARg"	;	# 137
    100.0 ;  # "W_gene_PPARg_RNAP"	;	# 138
    100.0 ;  # "W_gene_PPARg_gene_Trigger"	;	# 139
    100.0 ;  # "W_gene_PPARg_gene_CEBPa"	;	# 140
    100.0 ;  # "W_gene_PPARg_gene_EGR1"	;	# 141
    100.0 ;  # "W_gene_PPARg_gene_PU1"	;	# 142
    100.0 ;  # "W_gene_PPARg_gene_AP1"	;	# 143
    100.0 ;  # "W_gene_PU1_RNAP"	;	# 144
    100.0 ;  # "W_gene_PU1_gene_Trigger"	;	# 145
    100.0 ;  # "W_gene_PU1_gene_CEBPa"	;	# 146
    100.0 ;  # "W_gene_PU1_gene_PU1"	;	# 147
    100.0 ;  # "W_gene_PU1_gene_AP1"	;	# 148
    100.0 ;  # "W_gene_PU1_gene_OCT1"	;	# 149
    100.0 ;  # "W_gene_PU1_gene_AhR"	;	# 150
    100.0 ;  # "W_gene_PU1_gene_GFI1"	;	# 151
    100.0 ;  # "W_gene_Trigger_RNAP"	;	# 152
    100.0 ;  # "W_gene_cRAF_RNAP"	;	# 153
  ];

  return (lower_bound_array,upper_bound_array)
end

function estimate_parameters(pObjectiveFunction,initial_parameter_array,data_dictionary,number_of_iterations,error_target)

  # Get the initial parameter set -
  kVI = initial_parameter_array;
  N = length(kVI);

  # Initialize the bounds -
  (lower_bound_array,upper_bound_array) = bounds_function(N);

  # Initialize the algorithm parameters -
  D = 1+N/2;
  P_SUCC_TARGET = 2/11;
  CP = 1/12;
  CC = 2/(N+2);
  CCOV = 2/(N^2+6);
  PTHRESH = 0.99;

  # Set the initial search directions - apply -
  A = eye(N,N);
  AVG_P_SUCC = P_SUCC_TARGET;

  # Setup the run -
  should_loop_continue = true
  parameter_array = kVI;
  SIGMA = 0.10;

  # Initialize -
  error_array = zeros(2)
  error_array[1] = Inf
  counter = 1;
  while (should_loop_continue == true)

    # Generate new parameter set -
  	(perturbed_parameter_array,zV) = sample_function(parameter_array,initial_parameter_array,lower_bound_array,upper_bound_array,N,A,SIGMA);

    # Evaluate the parameter array -
    objective_function_array = pObjectiveFunction(perturbed_parameter_array);
    error_array[2] = sum(objective_function_array);

    # Was the step good or bad?
    delta_error = error_array[2] - error_array[1];
  	if (delta_error<0)
  		LAMBDA_SUCC = 1;
  	else
  		LAMBDA_SUCC = 0;
  	end

  	# Update the step size -
  	(SIGMA,AVG_P_SUCC) = updateStepSize(LAMBDA_SUCC,AVG_P_SUCC,P_SUCC_TARGET,CP,D,SIGMA);

    # Ok, recompute the direction array -
    if (LAMBDA_SUCC == 1)

      # Message -
      msg = "Step improved the objective function. New error = "*string(error_array[2])*" Old error = "*string(error_array[1])*" counter = "*string(counter)
      println(msg)

  		# Keep the peturbed parameters -
  		parameter_array = perturbed_parameter_array;

  		# Update the covariance matrix -
  		A = updateCholesky(A,zV,CCOV,AVG_P_SUCC,PTHRESH);

  		# Update the error -
  		error_array[1] = error_array[2];

  		# If we have a successful step and we are stuck at the min, reset this
  		if (SIGMA<=1e-5)
  			SIGMA = 0.05;
  		end;

  		# Update the counter =
  		counter = counter + 1;

      # Write parameters to disk ...
      writedlm("./parameter_best_v1.dat",parameter_array);

  	else
  		# If we have a successful step and we are stuck at the min, reset this
  		if (SIGMA<=1e-5)
  			SIGMA = 0.05;
  		end;

      msg = "No improvement. Current error = "*string(error_array[1])*" Step error = "*string(error_array[2])*" counter = "*string(counter)
      println(msg)
  	end

  	# Check to see if we need to set FLAG -
  	if (counter>=number_of_iterations || error_array[1] <= error_target)
  		should_loop_continue = false;
  	end
  end

  return parameter_array
end

function main()

  # Initialize -
  number_of_parameters = 153

  # Setup initial parameter guess -
  initial_parameter_array = zeros(number_of_parameters);

  # Load the data dictionary -
  data_dictionary = DataDictionary(0.0,0.0,0.0);

  # Load the parameter array -
  previous_parameter_array = readdlm("./parameter_best_v1.dat.5")

  # Update data dictionary to match new parameters before calculating obj
  parameter_mapping_array = data_dictionary["parameter_name_mapping_array"]
  for index = 1:length(parameter_mapping_array)
    if index <= data_dictionary["number_of_binding"]
      #value = data_dictionary["binding_parameter_dictionary"][parameter_mapping_array[index]]
      value = previous_parameter_array[index]
      initial_parameter_array[index] = value;
    else
      #value = data_dictionary["control_parameter_dictionary"][parameter_mapping_array[index]]
      value = previous_parameter_array[index]
      initial_parameter_array[index] = value;
    end
  end

  # Setup stoping criteria -
  NUMBER_OF_ITERATIONS = 100;
  ERR_TARGET = 10;

  # Setup function pointers -
  pObjectiveFunction = objective_function;

  # estimate parameters -
  parameter_array = estimate_parameters(pObjectiveFunction,initial_parameter_array,data_dictionary,NUMBER_OF_ITERATIONS,ERR_TARGET);

  # Write parameters to disk ...
  writedlm("./parameter_best_v1.dat",parameter_array);
end

# Run me ...
main()


text:
function setEu(aimCode)

    euVal = zeros(2, 1)
    if aimCode == 1
        euVal[1, 1] = 1
        euVal[2, 1] = 1
    elseif aimCode == 2
        euVal[1, 1] = -2
        euVal[2, 1] = -2
    elseif aimCode == 3
        euVal[1, 1] = 0
        euVal[2, 1] = 0
    elseif aimCode == 4
        euVal[1, 1] = 1
        euVal[2, 1] = 0
    elseif aimCode == 5
        euVal[1, 1] = -5
        euVal[2, 1] = -5
    end

    return euVal
end


text:
# Interpreter-based implementations of the PDDL.jl interface #

abstract type Interpreter end

include("utils.jl")
include("diff.jl")
include("evaluate.jl")
include("transition.jl")
include("available.jl")
include("execute.jl")
include("relevant.jl")
include("regress.jl")
include("update.jl")

include("concrete/concrete.jl")
include("abstract/abstract.jl")


text:
using FiniteDifferences: j′vp

# Default tolerances for testing Zygote against FiniteDifferences.
const _rtol = 1e-10
const _atol = 1e-10

#
# Print stuff for debugging
#

function print_adjoints(adjoint_ad, adjoint_fd, rtol, atol)
    @show typeof(adjoint_ad), typeof(adjoint_fd)
    adjoint_ad, adjoint_fd = to_vec(adjoint_ad)[1], to_vec(adjoint_fd)[1]
    println("atol is $atol, rtol is $rtol")
    println("ad, fd, abs, rel")
    abs_err = abs.(adjoint_ad .- adjoint_fd)
    rel_err = abs.((adjoint_ad .- adjoint_fd) ./ adjoint_ad)
    display([adjoint_ad adjoint_fd abs_err rel_err])
    println()
end

#
# Version of isapprox that works for lots of types.
#

function fd_isapprox(x_ad::Nothing, x_fd, rtol, atol)
    return fd_isapprox(x_fd, zero(x_fd), rtol, atol)
end
function fd_isapprox(x_ad::AbstractArray, x_fd::AbstractArray, rtol, atol)
    return all(fd_isapprox.(x_ad, x_fd, rtol, atol))
end
function fd_isapprox(x_ad::Real, x_fd::Real, rtol, atol)
    return isapprox(x_ad, x_fd; rtol=rtol, atol=atol)
end
function fd_isapprox(x_ad::NamedTuple, x_fd, rtol, atol)
    f = (x_ad, x_fd)->fd_isapprox(x_ad, x_fd, rtol, atol)
    return all([f(getfield(x_ad, key), getfield(x_fd, key)) for key in keys(x_ad)])
end
function fd_isapprox(x_ad::Tuple, x_fd::Tuple, rtol, atol)
    return all(map((x, x′)->fd_isapprox(x, x′, rtol, atol), x_ad, x_fd))
end
function fd_isapprox(x_ad::Dict, x_fd::Dict, rtol, atol)
    return all([fd_isapprox(get(()->nothing, x_ad, key), x_fd[key], rtol, atol) for
        key in keys(x_fd)])
end

#
# Check Zygote against FiniteDifferences.
#

function adjoint_test(
    f, ȳ, x...;
    rtol=_rtol,
    atol=_atol,
    fdm=FiniteDifferences.Central(5, 1),
    print_results=false,
)

    # Compute forwards-pass and j′vp.
    y, back = Zygote.forward(f, x...)
    adj_ad = back(ȳ)
    adj_fd = j′vp(fdm, f, ȳ, x...)

    # If unary, pull out first thing from ad.
    adj_ad = length(x) == 1 ? first(adj_ad) : adj_ad

    # Check that forwards-pass agrees with plain forwards-pass.
    @test y ≈ f(x...)

    # Check that ad and fd adjoints (approximately) agree.
    print_results && print_adjoints(adj_ad, adj_fd, rtol, atol)
    @test fd_isapprox(adj_ad, adj_fd, rtol, atol)
end


text:
using Downloads

const AWS_LAMBDA_RUNTIME_API = ENV["AWS_LAMBDA_RUNTIME_API"]
const HANDLER = ENV["_HANDLER"]  # from Dockerfile's CMD command
const LAMBDA_TASK_ROOT = ENV["LAMBDA_TASK_ROOT"]
const RUNTIME_URL = "http://$AWS_LAMBDA_RUNTIME_API/2018-06-01"

abstract type State end
struct Start <: State end
struct Received <: State end
struct Handled <: State end

function error(::Start, exception, request_id)
    @error "Unable to receive request" exception
end

function error(::Received, exception, request_id)
    @error "Unable to handle event" exception request_id
    post_error("runtime/invocation/$request_id/error", exception)
end

function error(::Handled, exception, request_id)
    @error "Unable to notify lambda runtime (invocation response)" exception request_id
end

function http_get(url)
    io = IOBuffer()
    response = request(url; output=io)

    seekstart(io)
    body = read(io, String)

    return response.headers, body
end

function http_post(url, headers, body)
    io = IOBuffer()
    write(io, body)
    seekstart(io)  # We need to seekstart for the runtime to properly read the body

    return request(url, method="POST", headers=headers, input=io)
end

function post_error(path::AbstractString, ex::Exception)
    headers = [
        "Content-type" => "application/json",
        "Lambda-Runtime-Function-Error-Type" => typeof(ex),
    ]

    url = "$RUNTIME_URL/$path"

    try
        http_post(url, headers, """{"errorType": "$(typeof(ex))", "errorMessage": "$(ex)"}""")
        @info "Notified lambda runtime about the error" url ex
    catch failure
        @error "Unable to notify lambda runtime about the error" url ex failure
    end
end

# Initialize Lambda function by loading user module
try
    @info "Loading user module"
    global mod, func = Symbol.(split(HANDLER, "."))
    @eval using $(mod)
catch ex
    @error "Initializtion error" ex
    post_error("runtime/init/error", ex)
    @info "Shutting down container"
    exit(1)
end

const MY_MODULE = Base.eval(Main, mod)
const MY_HANDLER = Base.eval(MY_MODULE, func)

@info "Start processing requests"
while true
    local state, request_id
    try
        state = Start()
        request_headers, request_body = http_get("$RUNTIME_URL/runtime/invocation/next")
        request_id_idx = findfirst(x -> lowercase(x[1]) == "lambda-runtime-aws-request-id", request_headers)
        request_id = request_headers[request_id_idx][2]
        @info "Received event" request_id
        state = Received()

        response = MY_HANDLER(request_body, request_headers)

        @info "Got response from handler" response
        state = Handled()

        http_post("$RUNTIME_URL/runtime/invocation/$request_id/response", [], response)
        @info "Notified Lambda Runtime"
    catch ex
        error(state, ex, request_id)
    end
end

@info "All done"


text:
function REPARATION_basic!(offspring::Sol)

    debug && DEBUG_correct_solution(offspring)

    if !SOLUTION_is_feasible(offspring)
        # repair solution
        for var in 1:offspring.prob.n
            if offspring.x[var] == 0
                if SOLUTION_is_feasible_after_adding_item(offspring, var)
                    SOLUTION_add_item!(offspring, var)
                end
            else
                if !SOLUTION_is_feasible(offspring)
                    SOLUTION_remove_item!(offspring, var)
                end
            end
        end
    end

    debug && DEBUG_feasible_solution(offspring)

end

function REPARATION_advanced!(offspring::Sol)

    debug && DEBUG_correct_solution(offspring)

    if !SOLUTION_is_feasible(offspring)

        λ = rand() # we pick a random direction for the reparation
        momkp = offspring.prob

        # count the number of ones (thus also the number of zeros) in the solution
        nb_ones = 0
        for component in offspring.x
            component == 1 ? nb_ones += 1 : nothing
        end

        idx_of_ones = Vector{Int}(undef, nb_ones) # vars to one
        idx_ones = 0

        idx_of_zeros = Vector{Int}(undef, momkp.n-nb_ones) # vars to zero
        idx_zeros = 0

        # fill array idx_of_ones and idx_of_zeros
        for var in 1:momkp.n
            if offspring.x[var] == 1
                idx_ones += 1
                idx_of_ones[idx_ones] = var
            else
                idx_zeros += 1
                idx_of_zeros[idx_zeros] = var
            end
        end

        # compute the ratio of each variable : ratios[idx_var]
        ratios = Vector{Float64}(undef, momkp.n)
        for var in 1:momkp.n
            profit = λ*momkp.P[1, var] + (1-λ)*momkp.P[2, var]
            weight = 0
            for dim in 1:momkp.m
                weight += momkp.W[dim, var]
            end
            ratios[var] = profit/weight
        end

        # sort the ratios of ones (worst to best)
        sort!(idx_of_ones, by = idx_var -> ratios[idx_var])

        # sort the ratios of zeros (best to worst)
        sort!(idx_of_zeros, by = idx_var -> ratios[idx_var], rev = true)

        # remove items to repair the solution starting from the worst variables
        idx = 0
        while idx < length(idx_of_ones) && !SOLUTION_is_feasible(offspring)
            idx += 1
            SOLUTION_remove_item!(offspring, idx_of_ones[idx])
        end

        # add items in the solution (while remaining feasible) starting from the best variables
        idx = 1
        while idx <= length(idx_of_zeros) && SOLUTION_is_feasible_after_adding_item(offspring, idx_of_zeros[idx])
            SOLUTION_add_item!(offspring, idx_of_zeros[idx])
            idx += 1
        end
    end

    debug && DEBUG_feasible_solution(offspring)

end


text: compute(n::Int, k::Int)::Float64 = round(7(1 - prod(1 .- k ./ (6n ÷ 7 + 1:n))), digits=9)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler