text:
#=
  https://dtai.cs.kuleuven.be/problog/tutorial/basic/10_inhibitioneffects.html
  """
  ...
  In some cases, however, if both causes are present this reduces the likelihood of the effect happening.
  For example when the presence of a second drug blocks the impact of the original drug. In such a case we
  can use negated heads. Depending on the chosen probabilities the effect can be lower than the joint occurance
  of c1 and c2 or be lower then either p1 or p2:
  """

  cf ~/blog/inhibition1.blog
     ~/webppl/inhibition1.wppl

=#

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

@model function inhibition1()
    c1 ~ flip(0.5)
    c2 ~ flip(0.5)

    e1 ~ c1 ? flip(0.3) : c2 ? flip(0.2) : flip(0.0)
    e2 ~ c2 ? flip(0.4) : c1 ? flip(0.2) : flip(0.0)
    # ee ~ flip(0.5)
    ee ~ Dirac(e1 || e2)

    true ~ Dirac(c1 == true)
    true ~ Dirac(c2 == true)
end

model = inhibition1()

num_chns = 4

# chns = sample(model, Prior(), 10_000)

# chns = sample(model, MH(), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, MH(), 10_000)

# chns = sample(model, PG(15), MCMCThreads(), 1_000, num_chns)
# chns = sample(model, PG(15), 1_000)

# chns = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chns)
chns = sample(model, SMC(1000), 10_000)

# chns = sample(model, IS(), 1_000)

#
display(chns)
show_var_dist_pct(chns,:c1,20)
show_var_dist_pct(chns,:c2,20)
show_var_dist_pct(chns,:e1,20)
show_var_dist_pct(chns,:e2,20)
show_var_dist_pct(chns,:ee,20)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
#
###
# AbstractDomain interface #TODO
###

"""
args:
    AbstractDomain{T,D}
    direction < D
ret:
    Bool
"""
function isperiodic end

"""
args:
    AbstractDomain
    direction
ret:
    Tuple of end points
"""
function end_points end

"""
check if domain end points match

args:
    AbstractDomain
    direction
ret:
    Bool
"""
function domains_match end

###
# IntervalDomain
###

""" 1D interval """
struct IntervalDomain{T<:Number,Ttag} <: AbstractDomain{T,1}
    x0::T
    x1::T
    isperiodic::Bool
    tags::Ttag

    function IntervalDomain(x0::Number, x1::Number, isperiodic, tags)
        T = promote_type(eltype.((x0, x1))...)
        new{T,typeof(tags)}(T(x0), T(x1), isperiodic, tags)
    end
end

function IntervalDomain(x0 = -1e0, x1 = 1e0,
                        isperiodic=false, tags = (nothing, nothing)
                       )
    IntervalDomain(x0, x1, isperiodic, tags)
end

isperiodic(dom::IntervalDomain) = dom.isperiodic
end_points(dom::IntervalDomain) = (dom.x0, dom.x1)

function domains_match(int1::IntervalDomain, int2::IntervalDomain)
    bools = isapprox.(end_points.((int1, int2))...)
    *(bools...)
end

###
# Box Domain
###

""" D-dimensional logically reectangular domain """
struct BoxDomain{T,D,Ti} <: AbstractDomain{T,D}
    intervals::Ti

    function BoxDomain(intervals::IntervalDomain...)
        T = promote_type(eltype.(intervals)...)
        D = length(intervals)
        new{T,D,typeof(intervals)}(intervals)
    end
end

isperiodic(box::BoxDomain, dir::Integer) = box.intervals[dir] |> isperiodic
end_points(box::BoxDomain, dir::Integer) = box.intervals[dir] |> end_points
end_points(box::BoxDomain) = end_points.(box.intervals)

Base.:*(int1::IntervalDomain, int2::IntervalDomain) = BoxDomain(int1, int2)
Base.:*(box1::BoxDomain, int2::IntervalDomain) = BoxDomain(box1.intervals..., int2)
Base.:*(int1::IntervalDomain, box2::BoxDomain) = BoxDomain(int1, box2.intervals...)
Base.:*(box1::BoxDomain, box2::BoxDomain) = BoxDomain(box1.intervals..., box2.intervals...)

###
# DeformedDomain
###

"""
Deform D-dimensional domain via mapping

args:
   -reference domain
   -mapping function
        (x1,...,xD) = map(r1, ..., rD)
   -isrescaling - do optimizations if
    mapping is a simple rescaling
        x1 = a1 + λ1 * x1(r1), ...,
        xD = aD + λD * xD(rD)
   -isseparable - do optimizations if
    mapping is separable
        x1 = x1(r1), ...,
        xD = xD(rD)
"""
struct DeformedDomain{T,D,Tdom<:AbstractDomain{T,D}, Tm} <: AbstractDomain{T,D}
    domain::Tdom
    mapping::Tm
#   isrescaling::Bool
    isseparable::Bool
end

function deform(domain, mapping = nothing, isseparable = false)
    DeformedDomain(domain, mapping, isseparable)
end

function domains_match(box1::BoxDomain{<:Number, D1}, box2::BoxDomain{<:Number, D2}) where{D1,D2}
    if D1 != D2
        @warn "dimension mismatch between $box1, $box2"
        return false
    end

    ret = true
    for i=1:D1
        ret * domains_match(box1.intervals[i], box2.intervals[i])
    end

    ret
end

# TODO - make a struct for mappings
# struct Deformations{D,Tmap} <: AbstractDeformation{D}
#   mapping::Tmap
#   isrescaling::Bool
#   isseparable::Bool
# end

function map_from_ref(domain, ref_domain;D=D) # TODO

    if domains_match(domain, ref_domain)
        return domain
    end

    xe = domain.end_points

    mapping = domain.mapping ==! nothing ? domain.mapping : (r -> r)
    mapping = mapping ∘ map
end

###
# Conveniences
###

function unit_box(D, args...) # TODO not exactly unit haha
    interval = IntervalDomain(-true, true, args...)
    BoxDomain((interval for i=1:D)...)
end

function polar2D(r, θ)
    x = @. r * cos(θ)
    y = @. r * sin(θ)
    x, y
end

function annulus_2D(r0, r1)
    intR = IntervalDomain(r0, r1, false, (:Inner, :Outer))
    intθ = IntervalDomain(-π,  π, true , (:Periodic, :Periodic))

    dom = intR * intθ

    deform(dom, polar2D)
end

###
# Transfinite Interpolation
###

function gordon_hall end # TODO

#

ext: jl


text:
module RocketFaultedObservableTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "FaultedObservable" begin

    println("Testing: faulted")

    struct DummyScheduler <: Rocket.AbstractScheduler end

    @testset begin
        @test faulted(1)      == FaultedObservable{Any, AsapScheduler}(1, AsapScheduler())
        @test faulted(Int, 1) == FaultedObservable{Int, AsapScheduler}(1, AsapScheduler())

        @test faulted("Error")      == faulted("Error")
        @test faulted(Int, "Error") != faulted(Float64, "Error")

        @test faulted(1; scheduler = DummyScheduler()) == FaultedObservable{Any, DummyScheduler}(1, DummyScheduler())

        @test_throws ErrorException faulted(Int)
        @test_throws ErrorException faulted(Float64)
    end

    @testset begin
        source = faulted(Int, "Error")
        io = IOBuffer()

        show(io, source)

        printed = String(take!(io))

        @test occursin("FaultedObservable", printed)
        @test occursin(string(eltype(source)), printed)
    end

    run_testset([
        (
            source = faulted(0),
            values = @ts(e(0))
        ),
        (
            source = faulted(0, scheduler = AsyncScheduler(0)),
            values = @ts(0 ~ e(0))
        )
    ])

end

end


text:
# Run after including ../src/CGP.jl
#  or uncomment the following line
include("../src/CGP.jl")

# Note that the parameters are set in the file ../src/SetParams.jl
# The default functions are defined in ../src/Func.jl
# The type  MyInt is defined in CGP.jl

#using Base.Test
using Printf
using Main.CGP

funcs = default_funcs()
const context = construct_contexts(Main.CGP.numinputs)[Main.CGP.numinputs]

# Generate random circuits and count the logic functions that are computed.
# Similar to results of Raman and Wagner 2011 to generate their Figure 2b.
#num_iterations = 50
num_iterations = 50000
println("num_iterations: ",num_iterations)
counts = create_count_function_hash( p.numinputs, p.numoutputs )
for _ = 1:num_iterations
    c = random_chromosome(p, funcs)
    #print_chromosome( c )
    # Executing on these inputs tests all possible bit combinations for inputs
    outputs = execute_chromosome(c, context)
  #=
    print("outputs: ",outputs)
    if length(outputs) > 1
      println("   hamming: ",hamming(outputs[1],outputs[2]))
    else
      println()
    end
  =#
    concat_outputs = concatenate_outputs(p.numinputs, outputs)
    #print("concat:  ")
    #Printf.@printf("0x%2x\n",concat_outputs)
    increment_count_function_hash(concat_outputs, counts )
    outputs
end

#if num_iterations <= 50
  Main.CGP.print_count_function_hash(counts)
#end
#Main.CGP.print_count_function_summary(counts)


text:
function picard(L,z;maxiter=10,tol=0.,relax=1.,n_eig_val=3,v0=[],output=true)
    if output
        println("Launching Picard...")
        println("Iter    Res:     dz:     z:")
        println("----------------------------------")
    end
    z0=complex(Inf)
    n=0
    flag=itsol_converged
    if v0==[]
        v0=ones(ComplexF64,size(L(0))[1])
    end
    M=L(1,oplist=["M"],in_or_ex=true) #get mass matrix
    try
        while abs(z-z0)>tol && n< maxiter
            if output; println(n,"\t\t",abs(z-z0),"\t", z );flush(stdout); end
            z0=z
            X=L(z0,oplist=["M","__aux__"]) #TODO put aux operator as default to linopfamily call and definition
            z,v0=Arpack.eigs(-X,M,sigma=z0,v0=v0,nev=n_eig_val)
            z=sqrt.(z)
            indexing=sortperm(z.-z0,by=abs)
            z,v0=z[indexing[1]],v0[:,indexing[1]]
            z=z0+relax*(z-z0)
            #TODO consider relaxation in v0
            n+=1
        end
    catch excp
        if output
            println("Error occured:")
            println(excp)
            println(typeof(excp))
            println("...aborted Picard!")
        end

        flag=itsol_unknown
        if typeof(excp) <: Arpack.ARPACKException
           flag=itsol_arpack_exception
           if excp==Arpack.ARPACKException(-9999)
               flag=itsol_arpack_9999
           end
        elseif excp== LinearAlgebra.SingularException(0) #This means that the solution is so good that L(z) cannot be LU factorized...TODO: implement other strategy into perturb
           flag=itsol_singular_exception
           L.params[L.eigval]=z
        end
    end
    if flag==itsol_converged
        L.params[L.eigval]=z
        if output;  println(n,"\t\t",abs(z-z0),"\t", z ); end

        if n>=maxiter
            flag=itsol_maxiter
            if output; println("Warning: Maximum number of iterations has been reached!");end

        elseif abs(z-z0)<=tol
            flag=itsol_converged
            if output; println("Solution has converged!"); end
        elseif isnan(z)
            flag=itsol_isnan
            if output; println("Warning: computer arithmetics problem. Eigenvalue is NaN"); end
        else
            if output; println("Warning: This should not be possible....\n If you can read this contact GAM!");end
            flag=itsol_impossible
        end
        if output
            println("...finished Picard!")
            println("#####################")
            println(" Picard results ")
            println("#####################")
            println("Number of steps: ",n)
            println("Last step parameter variation:",abs(z0-z))
            println("Eigenvalue:",z)
            println("Eigenvalue/(2*pi):",z/2/pi)
        end
    end
    return Solution(L.params,v0,[],L.eigval), n, flag
end


text:
module CLocales

using ResourceBundles
using ResourceBundles.LocaleIdTranslations
using ResourceBundles.LC

import ResourceBundles: LocaleId, S0, locale, CLocaleType

export current_locale, newlocale, duplocale, freelocale, nl_langinfo, clocale_id
export NlItem, category, offset, typ
export strcoll, localize_isless
export LC_GLOBAL_LOCALE

"""
    newlocale(cat::LC.CategorySet, loc::LocaleId[, base::CLocaleType])::CLocaleType

Return a pointer to an opaque locale-object of glibc.
Modifies and invalidate the previously returned `base`.
`categories` may be on or a tuple of symbols `:CTYPE, :NUMERIC, :TIME, :COLLATE, :MONETARY,
:MESSAGES, :ALL, :PAPER, :NAME, :ADDRESS, :TELEPHONE, :MEASUREMENT, :IDENTIFICATION`.
The locale id is converted to a string of a locale name in POSIX-form.
This name must exist in the output of `locale -a` on a system using glibc.
If no such object exists, return `Ptr{Nothing}(0)`.
"""
newlocale(cat::CategorySet, loc::LocaleId) = newlocale(cat, loc, CL0)

function newlocale(cat::CategorySet, loc::LocaleId, base::CLocaleType)
    cloc = loc_to_cloc(loc)
    newlocale_c(LC.mask(cat), cloc, base)
end

"""
    current_clocale()

Return the task-specific current glibc-locale.
Any call to `set_locale!` or `newlocale` invalidates the returned pointer.
"""
function current_clocale()
    locale().cloc
end

"""
    strcoll(a::AbstractString, b::AbstractString[, loc::LocaleId]) -> -1, 0, 1

Compare two strings a and b using the locale specified by `loc` or the default-locale
for category `:MESSAGES`.
Return `a < b ? -1 : a > b ? 1 : 0`.
"""
strcoll(a::AbstractString, b::AbstractString) = strcoll_c(a, b, current_clocale())
function strcoll(a::AbstractString, b::AbstractString, loc::LocaleId)
    ploc = newlocale(LC.ALL, loc)
    if ploc != CL0
        res = strcoll_c(a, b, ploc)
        freelocale(ploc)
        res
    else
        error("no posix locale found for $loc")
    end
end

"""
    localized_isless([loc::LocaleId])

Return isless function comparing two strings accoding to locale specific collation rules.
If loc is not specified, use category `COLLATE` of current locale.
"""
localized_isless() = (a::AbstractString, b::AbstractString) -> strcoll(a, b) < 0
localized_isless(loc::LocaleId) = (a::AbstractString, b::AbstractString) -> strcoll(a, b, loc) < 0

"""
Mask composed of output-type, category, and offset
"""
struct NlItem
    mask::Cint
end

"""
    nl_item(cat, offset, typ=0) -> NlItem(::Cint)

Create an NlItem object, which wraps Cint mask.
Result is:
typ == 0: char*
typ == 1: int of same size as pointer
typ == 2: *uint8
"""
function nl_item(cat::Category, offset::Integer, typ::Integer=0)
    NlItem(Cint(cat.id & 0x00fff)<<16 + Cint(offset) & 0xffff + Cint(typ)<<28)
end

category(nli::NlItem) = LC.ALL_CATS[((nli.mask>>16) & 0x0fff)+1]
offset(nli::NlItem) = Int(nli.mask & 0xffff)
typ(nli::NlItem) = Int((nli.mask>>>28) & 0xffff)

nl_cmask(cat::NlItem) = cat.mask & Cint(0x0ffffff)
nl_ctype(cat::NlItem) = Int(cat.mask>>>28)

"""
    nl_langinfo(cat::NlItem[, loc::CLocaleType] )
        -> string value | integer value

Provide information about the locale as stored in the glibc implementation.
"""
function nl_langinfo(cat::NlItem, loc::CLocaleType=current_clocale())
    res = nl_langinfo_c(nl_cmask(cat), loc)
    nl_convert(Val(nl_ctype(cat)), res)
end

const PTYPE = sizeof(Ptr) == 8 ? Int64 : Int32
const PU0 = Ptr{UInt8}(0)

nl_convert(::Val{0}, res::Ptr{UInt8}) = res === PU0 ? "" : unsafe_string(res)
nl_convert(::Val{1}, res::Ptr{UInt8}) = Int(Base.bitcast(PTYPE, res))
nl_convert(::Val{2}, res::Ptr{UInt8}) = res === PU0 ? 0 : Int(unsafe_wrap(Array, res, 1)[1])

err_day(i::Integer) = throw(ArgumentError("day of week ($i) not between 1 and 7"))
err_month(i::Integer) = throw(ArgumentError("month number ($i) not between 1 and 12"))

"""
    clocale_id(category[, cloc::CLocaleType])

Return the name string of a clocale or current clocale.
"""
clocale_id(cat::Category, loc::CLocaleType) = nl_langinfo(nl_item(cat, -1), loc)
clocale_id(cat::Category) = nl_langinfo(nl_item(cat, -1), current_clocale())

include("libc.jl")

## Interface constants derived from /usr/include/langinfo.h

const CTYPE_CODESET = nl_item(LC.CTYPE, 14)
const NUM_LC_CTYPE = 86

const RADIXCHAR = nl_item(LC.NUMERIC, 0)
const THOUSEP = nl_item(LC.NUMERIC, 1)
const THOUSANDS_SEP = nl_item(LC.NUMERIC, 1)
const GROUPING = nl_item(LC.NUMERIC, 2, 2)
const NUMERIC_CODESET = nl_item(LC.NUMERIC, 5)
const NUM_LC_NUMERIC = 6

ABDAY(i::Integer) = 1 <= i <= 7 ? nl_item(LC.TIME, i-1) : err_day(i)
DAY(i::Integer) = 1 <= i <= 7 ? nl_item(LC.TIME, i+6) : err_day(i)
ABMON(i::Integer) = 1 <= i <= 12 ? nl_item(LC.TIME, 13+i) : err_month(i)
MON(i::Integer) = 1 <= i <= 12 ? nl_item(LC.TIME, 25+i) : err_month(i)
const AM_STR = nl_item(LC.TIME, 38)
const PM_STR = nl_item(LC.TIME, 39)
const D_T_FMT = nl_item(LC.TIME, 40)
const D_FMT = nl_item(LC.TIME, 41)
const T_FMT = nl_item(LC.TIME, 42)
const T_FMT_AMPM = nl_item(LC.TIME, 43)
const ERA = nl_item(LC.TIME, 44)
const ERA_YEAR = nl_item(LC.TIME, 45)
const ERA_D_FMT = nl_item(LC.TIME, 46)
const ALT_DIGITS = nl_item(LC.TIME, 47)
const ERA_D_T_FMT = nl_item(LC.TIME, 48)
const ERA_T_FMT = nl_item(LC.TIME, 49)
const NUM_ERA_ENTRIES = nl_item(LC.TIME, 50, 1)
const WEEK_NDAYS = nl_item(LC.TIME, 101, 2)
const WEEK_1STDAY = nl_item(LC.TIME, 102, 1)
const WEEK_1STWEEK = nl_item(LC.TIME, 103, 2)
const FIRST_WEEKDAY = nl_item(LC.TIME, 104, 2)
const FIRST_WORKDAY = nl_item(LC.TIME, 105, 2)
const CAL_DIRECTION = nl_item(LC.TIME, 106, 2)
const TIMEZONE = nl_item(LC.TIME, 107)
const DATE_FMT = nl_item(LC.TIME, 108)
const TIME_CODESET = nl_item(LC.TIME, 110)
const NUM_LC_TIME = 111

const COLLATE_CODESET = nl_item(LC.COLLATE, 18)
const NUM_LC_COLLATE = 19

const YESEXPR = nl_item(LC.MESSAGES, 0)
const NOEXPR = nl_item(LC.MESSAGES, 1)
const YESSTR = nl_item(LC.MESSAGES, 2)
const NOSTR = nl_item(LC.MESSAGES, 3)
const MESSAGES_CODESET = nl_item(LC.MESSAGES, 4)
const NUM_LC_MESSAGES = 5

const INT_CURR_SYMBOL = nl_item(LC.MONETARY, 0)
const CURRENCY_SYMBOL = nl_item(LC.MONETARY, 1)
const MON_DECIMAL_POINT = nl_item(LC.MONETARY, 2)
const MON_THOUSANDS_SEP = nl_item(LC.MONETARY, 3)
const MON_GROUPING = nl_item(LC.MONETARY, 4, 2)
const POSITIVE_SIGN = nl_item(LC.MONETARY, 5)
const NEGATIVE_SIGN = nl_item(LC.MONETARY, 6)
const INT_FRAC_DIGITS = nl_item(LC.MONETARY, 7, 2)
const FRAC_DIGITS = nl_item(LC.MONETARY, 8, 2)
const P_CS_PRECEDES = nl_item(LC.MONETARY, 9, 2)
const P_SEP_BY_SPACE = nl_item(LC.MONETARY, 10, 2)
const N_CS_PRECEDES = nl_item(LC.MONETARY, 11, 2)
const N_SEP_BY_SPACE = nl_item(LC.MONETARY, 12, 2)
const P_SIGN_POSN = nl_item(LC.MONETARY, 13, 2)
const N_SIGN_POSN = nl_item(LC.MONETARY, 14, 2)
const CRNCYSTR = nl_item(LC.MONETARY, 15)
const INT_P_CS_PRECEDES = nl_item(LC.MONETARY, 16, 2)
const INT_P_SEP_BY_SPACE = nl_item(LC.MONETARY, 17, 2)
const INT_N_CS_PRECEDES = nl_item(LC.MONETARY, 18, 2)
const INT_N_SEP_BY_SPACE = nl_item(LC.MONETARY, 19, 2)
const INT_P_SIGN_POSN = nl_item(LC.MONETARY, 20, 2)
const INT_N_SIGN_POSN = nl_item(LC.MONETARY, 21, 2)
const MONETARY_CODESET = nl_item(LC.MONETARY, 45)
const NUM_LC_MONETARY = 46

const PAPER_HEIGHT = nl_item(LC.PAPER, 0, 1)
const PAPER_WIDTH = nl_item(LC.PAPER, 1, 1)
const PAPER_CODESET = nl_item(LC.PAPER, 2)
const NUM_LC_PAPER = 3

const NAME_FMT = nl_item(LC.NAME, 0)
const NAME_GEN = nl_item(LC.NAME, 1)
const NAME_MR = nl_item(LC.NAME, 2)
const NAME_MRS = nl_item(LC.NAME, 3)
const NAME_MISS = nl_item(LC.NAME, 4)
const NAME_MS = nl_item(LC.NAME, 5)
const NAME_CODESET = nl_item(LC.NAME, 6)
const NUM_LC_NAME = 7

const ADDRESS_POSTAL_FMT = nl_item(LC.ADDRESS, 0)
const ADDRESS_COUNTRY_NAME = nl_item(LC.ADDRESS, 1)
const ADDRESS_COUNTRY_POST = nl_item(LC.ADDRESS, 2)
const ADDRESS_COUNTRY_AB2 = nl_item(LC.ADDRESS, 3)
const ADDRESS_COUNTRY_AB3 = nl_item(LC.ADDRESS, 4)
const ADDRESS_COUNTRY_CAR = nl_item(LC.ADDRESS, 5)
const ADDRESS_COUNTRY_NUM = nl_item(LC.ADDRESS, 6, 1)
const ADDRESS_COUNTRY_ISBN = nl_item(LC.ADDRESS, 7)
const ADDRESS_LANG_NAME = nl_item(LC.ADDRESS, 8)
const ADDRESS_LANG_AB = nl_item(LC.ADDRESS, 9)
const ADDRESS_LANG_LIB = nl_item(LC.ADDRESS, 10)
const ADDRESS_LANG_TERM = nl_item(LC.ADDRESS, 11)
const ADDRESS_CODESET = nl_item(LC.ADDRESS, 12)
const NUM_LC_ADDRESS = 13

const TELEPHONE_TEL_INT_FMT = nl_item(LC.TELEPHONE, 0)
const TELEPHONE_TEL_DOM_FMT = nl_item(LC.TELEPHONE, 1)
const TELEPHONE_INT_SELECT = nl_item(LC.TELEPHONE, 2)
const TELEPHONE_INT_PREFIX = nl_item(LC.TELEPHONE, 3)
const TELEPHONE_CODESET = nl_item(LC.TELEPHONE, 4)
const NUM_LC_TELEPHONE = 5

const MEASUREMENT = nl_item(LC.MEASUREMENT, 0, 2)
const MEASUREMENT_CODESET = nl_item(LC.MEASUREMENT, 1)
const NUM_LC_MEASUREMENT = 2

const IDENTIFICATION_TITLE = nl_item(LC.IDENTIFICATION, 0)
const IDENTIFICATION_SOURCE = nl_item(LC.IDENTIFICATION, 1)
const IDENTIFICATION_ADDRESS = nl_item(LC.IDENTIFICATION, 2)
const IDENTIFICATION_CONTACT = nl_item(LC.IDENTIFICATION, 3)
const IDENTIFICATION_EMAIL = nl_item(LC.IDENTIFICATION, 4)
const IDENTIFICATION_TEL = nl_item(LC.IDENTIFICATION, 5)
const IDENTIFICATION_FAX = nl_item(LC.IDENTIFICATION, 6)
const IDENTIFICATION_LANGUAGE = nl_item(LC.IDENTIFICATION, 7)
const IDENTIFICATION_TERRITORY = nl_item(LC.IDENTIFICATION, 8)
const IDENTIFICATION_AUDIENCE = nl_item(LC.IDENTIFICATION, 9)
const IDENTIFICATION_APPLICATION = nl_item(LC.IDENTIFICATION, 10)
const IDENTIFICATION_ABBREVIATION = nl_item(LC.IDENTIFICATION, 11)
const IDENTIFICATION_REVISION = nl_item(LC.IDENTIFICATION, 12)
const IDENTIFICATION_DATE = nl_item(LC.IDENTIFICATION, 13)
const IDENTIFICATION_CATEGORY = nl_item(LC.IDENTIFICATION, 14)
const IDENTIFICATION_CODESET = nl_item(LC.IDENTIFICATION, 15)
const NUM_LC_IDENTIFICATION = 16

LC_GLOBAL_LOCALE = CLocaleType(-1)
CL0 = CLocaleType(0)

end # module


text:
using Base.Broadcast: BroadcastFunction, broadcasted, materialize

ChainRulesCore.@non_differentiable Base.getindex(m::AbstractAttenMask, I::Integer...)
ChainRulesCore.@non_differentiable Base.getindex(m::MaskIndexer, I::Integer...)
ChainRulesCore.@non_differentiable Base.getindex(m::AbstractAttenMask, I::Tuple)
ChainRulesCore.@non_differentiable Base.getindex(m::MaskIndexer, I::Tuple)
ChainRulesCore.@non_differentiable (::Type{<:AbstractAttenMask})(args...)
ChainRulesCore.@non_differentiable (::Type{<:AbstractAttenMaskOp})(args...)
ChainRulesCore.@non_differentiable getmask(arg...)

function ChainRulesCore.rrule(::typeof(apply_mask), op::NaiveAttenMaskOp, mask, score)
    naive_apply_mask_pullback(Ȳ) = (NoTangent(), NoTangent(), NoTangent(), @thunk unthunk(Ȳ) .* mask)
    return score .* mask, naive_apply_mask_pullback
end

function ChainRulesCore.rrule(config::RuleConfig, ::typeof(apply_broadcast_mask), f, mask, score, scale)
    tmp = getmask(mask, score, scale)
    Y, back = rrule_via_ad(config, broadcast, f, tmp, score)
    function fallback_apply_broadcast_mask_pullback(Ȳ)
        Ȳs = (back∘unthunk)(Ȳ)
        return (NoTangent(), Ȳs[2], NoTangent(), Ȳs[4], NoTangent())
    end
    Y, fallback_apply_broadcast_mask_pullback
end

function ChainRulesCore.rrule(config::RuleConfig, ::typeof(apply_broadcast_mask), ::typeof(+), mask, score, scale)
    apply_broadcast_mask_pullback(Ȳ) = (NoTangent(), NoTangent(), NoTangent(), Ȳ, NoTangent())
    apply_broadcast_mask(+, mask, score, scale), apply_broadcast_mask_pullback
end

function ChainRulesCore.rrule(config::RuleConfig, ::typeof(apply_broadcast_mask), ::typeof(*), mask, score, scale)
    function apply_broadcast_mask_pullback(Ȳ)
        thk = @thunk begin
            tmp = unthunk(Ȳ)
            @. tmp * (mask * scale)
        end
        return (NoTangent(), NoTangent(), NoTangent(), thk, NoTangent())
    end
    apply_broadcast_mask(*, mask, score, scale), apply_broadcast_mask_pullback
end

@init @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" begin
    Zygote.unbroadcast(x::AbstractAttenMask, _) = nothing
end


text:
@testset "Schur Complement" begin
    T = Float64
    n = 11
    m = 16
    M = rand(n+m,n+m)
    S = ContactImplicitMPC.Schur(M, n=n, m=m)

    u = rand(T,n)
    v = rand(T,m)
    D = rand(m,m)
    ContactImplicitMPC.schur_factorize!(S, D)
    # @benchmark ContactImplicitMPC.schur_factorize!(S, D)
    ContactImplicitMPC.schur_solve!(S, u, v)
    # @benchmark ContactImplicitMPC.schur_solve!(S, u, v)
    M1 = [S.A S.B; S.C D]
    @test norm(M1*[S.x; S.y] - [u; v], Inf) < 1e-9
end

@testset "Schur Complement: ill-conditioning" begin
    T = Float64
    s = get_simulation("flamingo", "flat_2D_lc", "flat")
    ref_traj = deepcopy(get_trajectory(s.model, s.env,
        joinpath(module_dir(), "src/dynamics/flamingo/gaits/gait_forward_36_4.jld2"),
        load_type = :split_traj_alt))

    ix, iy1, iy2 = linearization_var_index(s.model, s.env)
    idyn, irst, ibil = linearization_term_index(s.model, s.env)
    nz = num_var(s.model, s.env)
    nθ = num_data(s.model)
    nx = length(ix)
    ny = length(iy1)

    r0 = zeros(T, nz)
    rz0 = zeros(T, nz, nz)
    rθ0 = zeros(T, nz, nθ)
    z0 = zeros(T, nz)
    θ0 = zeros(T, nθ)
    z0 .= max.(1e-6, ref_traj.z[10])
    θ0 .= ref_traj.θ[10]
    κ0 = 0.0
    s.res.r!(r0, z0, θ0, κ0)
    s.res.rz!(rz0, z0, θ0)
    s.res.rθ!(rθ0, z0, θ0)

    A = rz0[idyn, ix]
    B = rz0[idyn, iy1]
    C = rz0[irst, ix]
    D = rz0[irst, iy1] - rz0[irst, iy2] * Diagonal(z0[iy2] ./ z0[iy1])
    M = [A B ; C D]
    S = ContactImplicitMPC.Schur(M, n=nx, m=ny)

    u = r0[idyn]
    v = r0[irst] - r0[ibil] ./ z0[iy1]

    # @benchmark schur_factorize!(S, D)
    ContactImplicitMPC.schur_factorize!(S, D)
    # @benchmark ContactImplicitMPC.schur_solve!(S, u, v)
    ContactImplicitMPC.schur_solve!(S, u, v)
    M1 = [S.A S.B; S.C D]
    norm(M1*[S.x; S.y] - [u; v], Inf)
    @test norm(M1*[S.x; S.y] - [u; v], Inf) < 1e-7
end


text:
function is_complex_datatype(datatype::BLPDataType) :: Bool
    return (
            datatype == BLPAPI_DATATYPE_SEQUENCE ||
            datatype == BLPAPI_DATATYPE_CHOICE
        )
end

is_simple_datatype(datatype::BLPDataType) = !is_complex_datatype(datatype)

# an enumeration is also a simple datatype
function is_enumeration_datatype(datatype::BLPDataType) :: Bool
    return datatype == BLPAPI_DATATYPE_ENUMERATION
end

is_choice_datatype(datatype::BLPDataType) = datatype == BLPAPI_DATATYPE_CHOICE

function AbstractSchemaTypeDefinition(handle::Ptr{Cvoid})

    ptr_check(handle, "Failed to create AbstractSchemaTypeDefinition")

    if blpapi_SchemaTypeDefinition_isComplex(handle) != 0

        if blpapi_SchemaTypeDefinition_numElementDefinitions(handle) > 0
            elements = [ SchemaElementDefinition(blpapi_SchemaTypeDefinition_getElementDefinitionAt(handle, i)) for i in 0:(blpapi_SchemaTypeDefinition_numElementDefinitions(handle) - 1) ]
        else
            elements = Vector{SchemaElementDefinition}()
        end

        return ComplexSchemaTypeDefinition(
                BLPName(blpapi_SchemaTypeDefinition_name(handle)),
                unsafe_string(blpapi_SchemaTypeDefinition_description(handle)),
                SchemaStatus(blpapi_SchemaTypeDefinition_status(handle)),
                BLPDataType(blpapi_SchemaTypeDefinition_datatype(handle)),
                elements)

    elseif blpapi_SchemaTypeDefinition_isEnumeration(handle) != 0

        return EnumerationSchemaTypeDefinition(
                BLPName(blpapi_SchemaTypeDefinition_name(handle)),
                unsafe_string(blpapi_SchemaTypeDefinition_description(handle)),
                SchemaStatus(blpapi_SchemaTypeDefinition_status(handle)),
                BLPDataType(blpapi_SchemaTypeDefinition_datatype(handle)),
                BLPConstantList(blpapi_SchemaTypeDefinition_enumeration(handle))
            )

    elseif blpapi_SchemaTypeDefinition_isSimple(handle) != 0
        return SimpleSchemaTypeDefinition(
                BLPName(blpapi_SchemaTypeDefinition_name(handle)),
                unsafe_string(blpapi_SchemaTypeDefinition_description(handle)),
                SchemaStatus(blpapi_SchemaTypeDefinition_status(handle)),
                BLPDataType(blpapi_SchemaTypeDefinition_datatype(handle))
            )
    end

    error("Failed to create AbstractSchemaTypeDefinition: unsupported datatype $(blpapi_SchemaTypeDefinition_datatype(handle)).")
end

function SchemaElementDefinition(handle::Ptr{Cvoid})

    ptr_check(handle, "Failed to create SchemaElementDefinition")

    # alternate_names
    if blpapi_SchemaElementDefinition_numAlternateNames(handle) > 0
        alternate_names = [ BLPName(blpapi_SchemaElementDefinition_getAlternateName(handle, i)) for i in 0:(blpapi_SchemaElementDefinition_numAlternateNames(handle)-1) ]
    else
        alternate_names = Vector{BLPName}()
    end

    return SchemaElementDefinition(
            BLPName(blpapi_SchemaElementDefinition_name(handle)),
            SchemaStatus(blpapi_SchemaElementDefinition_status(handle)),
            AbstractSchemaTypeDefinition(blpapi_SchemaElementDefinition_type(handle)),
            alternate_names,
            UInt64(blpapi_SchemaElementDefinition_minValues(handle)),
            UInt64(blpapi_SchemaElementDefinition_maxValues(handle))
        )
end

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
# CZAS NORMALNEGO ALGORYTMU 4-1024 (1 PROBA)
indeksyNormal = [4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 535, 536, 537, 538, 539, 540, 541, 542, 543, 544, 545, 546, 547, 548, 549, 550, 551, 552, 553, 554, 555, 556, 557, 558, 559, 560, 561, 562, 563, 564, 565, 566, 567, 568, 569, 570, 571, 572, 573, 574, 575, 576, 577, 578, 579, 580, 581, 582, 583, 584, 585, 586, 587, 588, 589, 590, 591, 592, 593, 594, 595, 596, 597, 598, 599, 600, 601, 602, 603, 604, 605, 606, 607, 608, 609, 610, 611, 612, 613, 614, 615, 616, 617, 618, 619, 620, 621, 622, 623, 624, 625, 626, 627, 628, 629, 630, 631, 632, 633, 634, 635, 636, 637, 638, 639, 640, 641, 642, 643, 644, 645, 646, 647, 648, 649, 650, 651, 652, 653, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 666, 667, 668, 669, 670, 671, 672, 673, 674, 675, 676, 677, 678, 679, 680, 681, 682, 683, 684, 685, 686, 687, 688, 689, 690, 691, 692, 693, 694, 695, 696, 697, 698, 699, 700, 701, 702, 703, 704, 705, 706, 707, 708, 709, 710, 711, 712, 713, 714, 715, 716, 717, 718, 719, 720, 721, 722, 723, 724, 725, 726, 727, 728, 729, 730, 731, 732, 733, 734, 735, 736, 737, 738, 739, 740, 741, 742, 743, 744, 745, 746, 747, 748, 749, 750, 751, 752, 753, 754, 755, 756, 757, 758, 759, 760, 761, 762, 763, 764, 765, 766, 767, 768, 769, 770, 771, 772, 773, 774, 775, 776, 777, 778, 779, 780, 781, 782, 783, 784, 785, 786, 787, 788, 789, 790, 791, 792, 793, 794, 795, 796, 797, 798, 799, 800, 801, 802, 803, 804, 805, 806, 807, 808, 809, 810, 811, 812, 813, 814, 815, 816, 817, 818, 819, 820, 821, 822, 823, 824, 825, 826, 827, 828, 829, 830, 831, 832, 833, 834, 835, 836, 837, 838, 839, 840, 841, 842, 843, 844, 845, 846, 847, 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 858, 859, 860, 861, 862, 863, 864, 865, 866, 867, 868, 869, 870, 871, 872, 873, 874, 875, 876, 877, 878, 879, 880, 881, 882, 883, 884, 885, 886, 887, 888, 889, 890, 891, 892, 893, 894, 895, 896, 897, 898, 899, 900, 901, 902, 903, 904, 905, 906, 907, 908, 909, 910, 911, 912, 913, 914, 915, 916, 917, 918, 919, 920, 921, 922, 923, 924, 925, 926, 927, 928, 929, 930, 931, 932, 933, 934, 935, 936, 937, 938, 939, 940, 941, 942, 943, 944, 945, 946, 947, 948, 949, 950, 951, 952, 953, 954, 955, 956, 957, 958, 959, 960, 961, 962, 963, 964, 965, 966, 967, 968, 969, 970, 971, 972, 973, 974, 975, 976, 977, 978, 979, 980, 981, 982, 983, 984, 985, 986, 987, 988, 989, 990, 991, 992, 993, 994, 995, 996, 997, 998, 999, 1000, 1001, 1002, 1003, 1004, 1005, 1006, 1007, 1008, 1009, 1010, 1011, 1012, 1013, 1014, 1015, 1016, 1017, 1018, 1019, 1020, 1021, 1022, 1023, 1024]
wynikiNormal = [5.14e-7, 1.028e-6, 5.14e-7, 1.028e-6, 5.14e-7, 1.028e-6, 1.542e-6, 2.571e-6, 3.085e-6, 3.599e-6, 4.626e-6, 5.655e-6, 7.196e-6, 8.739e-6, 1.0281e-5, 1.2338e-5, 4.4723e-5, 1.7478e-5, 2.2618e-5, 2.5702e-5, 3.3413e-5, 6.4257e-5, 3.9582e-5, 3.8554e-5, 4.1125e-5, 6.6827e-5, 6.9398e-5, 7.608e-5, 8.1221e-5, 8.6876e-5, 9.1502e-5, 0.000100241, 9.6643e-5, 9.9213e-5, 0.000122345, 0.000167582, 0.00013057, 0.000158329, 0.000152161, 0.000176835, 0.000171181, 0.000199968, 0.000250345, 0.000212819, 0.000222072, 0.000244176, 0.000254972, 0.000279647, 0.000286843, 0.000312032, 0.000404562, 0.000342361, 0.000397879, 0.000387084, 0.000410731, 0.000431293, 0.000489381, 0.0004791, 0.000505831, 0.00056392, 0.00056392, 0.000581911, 0.000612755, 0.00087081, 0.000804498, 0.000880578, 0.000963855, 0.000926842, 0.00083277, 0.000836369, 0.000869268, 0.000888288, 0.00103994, 0.000975678, 0.00100549, 0.00109905, 0.00108723, 0.00126715, 0.00156222, 0.00149539, 0.00146711, 0.00136225, 0.00141571, 0.00146609, 0.00146403, 0.00151852, 0.00156787, 0.00173854, 0.0019714, 0.00191846, 0.00195444, 0.00201613, 0.00204646, 0.00213744, 0.00226647, 0.00250499, 0.00221404, 0.0023004, 0.00234255, 0.00263351, 0.00301699, 0.00287254, 0.0027759, 0.00288848, 0.0029697, 0.00315116, 0.00330281, 0.00300877, 0.00308177, 0.00322519, 0.00391299, 0.00335267, 0.00380196, 0.003712, 0.00888905, 0.00396388, 0.00429082, 0.00429802, 0.00400655, 0.00422297, 0.00493031, 0.00437719, 0.00437975, 0.00538268, 0.00507733, 0.00505882, 0.00578364, 0.00509892, 0.00511743, 0.00607563, 0.00537959, 0.00625657, 0.00612549, 0.00604684, 0.0065208, 0.00604016, 0.00697779, 0.00642621, 0.00729291, 0.00664828, 0.00701326, 0.00723534, 0.00735048, 0.00757616, 0.00747437, 0.00835906, 0.00829326, 0.0079324, 0.00861558, 0.00810409, 0.00951826, 0.00932189, 0.00922833, 0.00962004, 0.00902168, 0.0100585, 0.0101048, 0.00957994, 0.0106769, 0.0109345, 0.0102708, 0.0138857, 0.0110897, 0.011137, 0.0113159, 0.0123024, 0.0130462, 0.0125589, 0.0126175, 0.0125908, 0.0128242, 0.0131161, 0.0135865, 0.0130462, 0.0137813, 0.0137191, 0.013916, 0.0142923, 0.0150125, 0.0139355, 0.0144969, 0.0144995, 0.0151343, 0.0153728, 0.0152011, 0.0154027, 0.0160863, 0.016255, 0.0161506, 0.0166061, 0.0168369, 0.0169541, 0.017903, 0.0185286, 0.0188242, 0.0197171, 0.0193105, 0.0199952, 0.0197552, 0.0205607, 0.0205571, 0.0207365, 0.0215682, 0.023462, 0.0235386, 0.0220885, 0.022364, 0.0222288, 0.0237905, 0.0228827, 0.024266, 0.0244094, 0.0250463, 0.0268311, 0.0253954, 0.0255537, 0.0259172, 0.026277, 0.0267417, 0.0274413, 0.0270697, 0.0276434, 0.0277575, 0.0286051, 0.0295382, 0.0291177, 0.0299232, 0.0294852, 0.0303828, 0.031123, 0.0305169, 0.0306084, 0.0314176, 0.0326811, 0.0315908, 0.0316556, 0.0329669, 0.0342315, 0.0354637, 0.0357217, 0.036281, 0.0367581, 0.0368701, 0.0383892, 0.0381758, 0.0386503, 0.0389865, 0.038785, 0.0402917, 0.0410432, 0.0404012, 0.041528, 0.0437133, 0.0417018, 0.042734, 0.0436007, 0.044436, 0.0422569, 0.0434285, 0.0439944, 0.043722, 0.0440906, 0.0468469, 0.0437148, 0.0463293, 0.0473846, 0.0487469, 0.048256, 0.0501251, 0.0483223, 0.0508103, 0.0506232, 0.0519777, 0.0519942, 0.0536685, 0.0551767, 0.0545871, 0.0543028, 0.0542961, 0.0557411, 0.0570684, 0.0562865, 0.0566865, 0.0552965, 0.0587108, 0.0573481, 0.0597096, 0.0637866, 0.0613161, 0.0624917, 0.0623246, 0.0639789, 0.063729, 0.0651766, 0.0641537, 0.0655354, 0.0662587, 0.0666669, 0.0703974, 0.0688151, 0.0668812, 0.065409, 0.066926, 0.0691169, 0.0725538, 0.0727816, 0.0733265, 0.0743386, 0.0750177, 0.0746959, 0.0771017, 0.0776476, 0.0782028, 0.0788557, 0.0792952, 0.076511, 0.0782198, 0.0778769, 0.0814007, 0.0832513, 0.0850953, 0.0860725, 0.0863999, 0.0854237, 0.0881359, 0.099311, 0.0972568, 0.102221, 0.0954082, 0.0933978, 0.0993768, 0.0929269, 0.0938964, 0.0960092, 0.0967453, 0.0979138, 0.098727, 0.0996374, 0.100889, 0.101569, 0.101402, 0.105395, 0.101822, 0.105141, 0.107467, 0.106581, 0.109128, 0.110695, 0.107076, 0.122094, 0.113382, 0.118433, 0.121925, 0.13772, 0.122154, 0.123836, 0.124641, 0.128001, 0.124951, 0.122037, 0.136801, 0.149063, 0.131214, 0.13052, 0.142691, 0.153619, 0.136026, 0.15122, 0.14342, 0.154906, 0.143233, 0.140714, 0.150502, 0.144756, 0.162844, 0.157602, 0.150509, 0.164092, 0.145797, 0.148571, 0.15148, 0.147333, 0.149221, 0.15173, 0.146198, 0.173646, 0.157258, 0.152108, 0.17138, 0.151935, 0.151878, 0.14948, 0.154634, 0.168015, 0.16078, 0.162169, 0.167383, 0.167535, 0.163324, 0.16555, 0.170367, 0.170004, 0.169016, 0.174337, 0.175648, 0.174929, 0.179601, 0.179029, 0.176425, 0.178938, 0.180928, 0.176575, 0.176232, 0.183423, 0.181434, 0.184428, 0.184556, 0.183646, 0.181007, 0.188763, 0.192984, 0.191903, 0.195225, 0.194671, 0.188872, 0.200806, 0.204544, 0.215357, 0.218004, 0.222678, 0.218146, 0.220267, 0.221524, 0.225612, 0.229129, 0.22377, 0.229627, 0.233997, 0.245438, 0.250883, 0.246389, 0.267191, 0.250421, 0.248201, 0.26506, 0.266477, 0.255528, 0.254468, 0.251468, 0.243596, 0.244948, 0.257266, 0.257016, 0.267592, 0.267461, 0.268949, 0.261414, 0.265142, 0.260219, 0.262508, 0.267136, 0.264619, 0.283398, 0.29317, 0.291515, 0.286036, 0.278449, 0.302661, 0.426915, 0.379697, 0.35911, 0.332429, 0.317078, 0.311217, 0.318142, 0.310883, 0.301983, 0.294825, 0.308267, 0.304089, 0.316754, 0.308043, 0.31422, 0.323442, 0.312821, 0.313458, 0.326075, 0.32027, 0.343862, 0.341855, 0.332053, 0.324593, 0.320684, 0.330552, 0.328432, 0.325064, 0.333403, 0.340747, 0.338037, 0.338648, 0.345018, 0.39581, 0.381063, 0.39292, 0.363885, 0.454182, 0.4314, 0.534177, 0.467974, 0.39375, 0.371925, 0.383573, 0.391666, 0.367958, 0.366461, 0.385287, 0.369164, 0.377458, 0.374833, 0.380357, 0.438789, 0.474939, 0.498723, 0.463478, 0.414744, 0.42826, 0.459133, 0.439587, 0.441646, 0.443758, 0.521191, 0.490236, 0.457921, 0.447679, 0.441414, 0.448829, 0.463271, 0.45182, 0.440822, 0.466739, 0.452837, 0.4554, 0.46665, 0.469842, 0.450062, 0.441317, 0.462834, 0.448949, 0.442331, 0.461824, 0.448548, 0.529146, 0.524394, 0.48534, 0.464375, 0.455645, 0.476671, 0.48406, 0.470089, 0.477547, 0.493352, 0.55745, 0.540892, 0.617659, 0.541162, 0.512501, 0.507709, 0.511819, 0.516729, 0.571563, 0.715315, 0.752682, 0.703535, 0.687922, 0.624454, 0.604883, 0.720002, 0.607014, 0.651741, 0.630744, 0.57852, 0.572346, 0.54857, 0.556822, 0.556545, 0.546309, 0.585392, 0.550879, 0.544821, 0.576733, 0.565506, 0.559287, 0.585784, 0.593912, 0.595862, 0.5665, 0.577297, 0.611184, 0.617884, 0.589468, 0.592561, 0.625631, 0.592532, 0.586158, 0.593459, 0.589012, 0.613548, 0.629008, 0.617611, 0.635149, 0.641955, 0.609838, 0.615657, 0.669637, 0.661501, 0.627022, 0.642316, 0.662542, 0.665231, 0.661377, 0.665768, 0.683028, 0.64944, 0.652104, 0.666069, 0.671978, 0.701906, 0.664813, 0.706375, 0.678683, 0.705303, 0.71854, 0.703028, 0.715178, 0.705109, 0.714977, 0.768602, 0.732617, 0.748559, 0.757922, 0.729265, 0.757838, 0.75264, 0.750407, 0.754484, 0.754359, 0.751514, 0.786787, 0.756142, 0.763772, 0.781275, 0.755479, 0.782195, 0.80318, 0.805854, 0.793695, 0.808595, 0.811548, 0.811839, 0.824553, 0.818959, 0.827369, 0.852773, 0.830636, 0.836004, 0.806208, 0.851084, 0.893533, 0.809894, 0.859965, 0.838532, 0.84628, 0.842111, 0.886989, 0.862851, 0.858341, 0.855688, 0.899716, 0.889327, 0.868092, 0.905066, 0.872358, 0.870096, 0.859366, 0.871247, 0.886268, 0.881831, 0.877771, 0.887588, 0.889106, 0.903536, 0.904155, 0.912371, 0.909833, 0.938653, 0.910469, 0.918731, 0.919646, 0.928429, 0.928636, 0.975046, 0.943006, 0.941074, 0.94698, 0.954384, 0.948486, 0.94945, 0.964386, 0.972441, 0.975272, 0.971686, 0.970482, 0.975734, 0.991236, 0.998943, 1.00289, 0.998847, 1.003, 0.996035, 1.00595, 1.01327, 1.01645, 1.02379, 1.02575, 1.03968, 1.04201, 1.06053, 1.04194, 1.0495, 1.06204, 1.06883, 1.06686, 1.06557, 1.07623, 1.11139, 1.08334, 1.09024, 1.09181, 1.08498, 1.09915, 1.09357, 1.10975, 1.10661, 1.11873, 1.11975, 1.12468, 1.13458, 1.14283, 1.13566, 1.14559, 1.14207, 1.17234, 1.16304, 1.16586, 1.17995, 1.17339, 1.19186, 1.22418, 1.26995, 1.20195, 1.20714, 1.21935, 1.20635, 1.20671, 1.35243, 1.23342, 1.24468, 1.32899, 1.25477, 1.23931, 1.25691, 1.26389, 1.24944, 1.25299, 1.26664, 1.27004, 1.28787, 1.29314, 1.29858, 1.31107, 1.31271, 1.30543, 1.31402, 1.30775, 1.30722, 1.3269, 1.34099, 1.34337, 1.35419, 1.36398, 1.36321, 1.38821, 1.3702, 1.3786, 1.37707, 1.38045, 1.37877, 1.41017, 1.39667, 1.39985, 1.41066, 1.47349, 1.41461, 1.41776, 1.43156, 1.44022, 1.45012, 1.51846, 1.45592, 1.444, 1.45931, 1.49137, 1.46669, 1.46519, 1.48058, 1.58388, 1.51114, 1.4954, 1.50232, 1.51139, 1.50919, 1.52006, 1.5318, 1.51591, 1.5271, 1.54203, 1.55806, 1.55311, 1.70687, 1.56267, 1.55396, 1.59739, 1.58317, 1.59509, 1.6076, 1.5982, 1.63075, 1.61585, 1.6115, 1.64811, 1.71577, 1.65892, 1.62869, 1.64102, 1.66309, 1.664, 1.65984, 1.67378, 1.68167, 1.70118, 1.67333, 1.69758, 1.68997, 1.71424, 1.71161, 1.71653, 1.72393, 1.73689, 1.73327, 1.73588, 1.74776, 1.82853, 1.75715, 1.77829, 2.74038, 2.75825, 2.43404, 2.39458, 2.36478, 2.35195, 2.12313, 2.06044, 2.03302, 2.05849, 2.09263, 2.07677, 2.08951, 2.10338, 2.14197, 2.11046, 2.12169, 2.1464, 2.12369, 2.14594, 2.15349, 2.16745, 2.11186, 2.06566, 2.34357, 2.367, 2.3537, 2.33694, 2.55914, 2.38232, 2.36434, 2.11141, 1.97219, 2.01552, 2.00781, 2.00139, 2.02354, 2.05682, 2.0446, 2.04537, 2.04182, 2.07909, 2.20516, 2.07363, 2.08177, 2.08623, 2.06895, 2.07549, 2.09201, 2.19106, 2.11866, 2.11132, 2.12169, 2.24747, 2.13486, 2.1414, 2.14225, 2.16881, 2.15069, 2.18262, 2.18526, 2.41744, 2.20639, 2.22079, 2.22011, 2.23561, 2.21661, 2.22428, 2.23001, 2.24267, 2.32486, 2.27308, 2.25277, 2.28587, 2.27437, 2.31592, 2.2978, 2.32042, 2.30438, 2.31054, 2.31654, 2.33804, 2.32388, 2.36031, 2.33414, 2.39056, 2.3562, 2.3781, 2.3917, 2.45733, 2.44022, 2.40223, 2.66146, 2.44358, 2.44035, 2.44318, 2.45293, 2.47395, 2.4582, 2.45835, 2.47562, 2.4868, 2.51302, 2.53959, 2.51411, 2.54618, 2.51901, 2.52062, 2.53212, 2.56471, 2.58369, 2.57237, 2.5705, 2.577, 2.5833, 2.84643, 2.60912, 2.68037, 2.67995, 2.86548, 2.62414, 2.64238, 2.65793, 2.66739, 2.70008, 2.95411, 2.74937, 2.68106, 2.69243, 2.70925, 2.71207, 3.00298, 2.81568, 2.84134, 2.73066, 2.76289, 2.75086, 2.77529, 2.95398, 2.79265, 2.79618, 2.78259, 3.00443, 2.80783, 2.81067, 2.85485, 2.83892, 2.95514, 3.76249, 2.97783, 2.88877, 2.90692, 2.89399, 2.90611, 2.92184, 2.94588, 6.71817]

# CZASY STRASSENA(STOP = 64) 4-1024 (20 PROB)
indeksyStrass = [4, 8, 16, 32, 64, 128, 256, 512, 1024]
wynikiStrass = [6.1685e-7, 2.0305e-6, 1.56529e-5, 0.000134323, 0.000807479, 0.00482459, 0.0336755, 0.248749, 1.78561]
wynikiPrzeindeksowaneStrass = [6.1685e-7, 2.0305e-6, 2.0305e-6, 2.0305e-6, 2.0305e-6, 1.56529e-5, 1.56529e-5, 1.56529e-5, 1.56529e-5, 1.56529e-5, 1.56529e-5, 1.56529e-5, 1.56529e-5, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000134323, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.000807479, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.00482459, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.0336755, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 0.248749, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561, 1.78561]

# Potęgi dwjójki dla 4096- 1, dla 2048 - 2
reszta = [2048, 4096]
ts = [12.6364, 95.074557115]
tn = [82.5707, 834.604167245]

# 20 powtórzeń
ind = [4, 8, 16, 32, 64, 128, 256, 512, 1024]
vs = [2.827e-7, 9.7675e-7, 1.13864e-5, 6.7521e-5, 0.0005998, 0.00421873, 0.0304357, 0.258856, 1.82804]
vn = [2.313e-7, 7.9675e-7, 1.12321e-5, 6.88323e-5, 0.000578544, 0.00521803, 0.0416395, 0.39834, 6.66768]

max_stars_repo_path: II rok/I semestr/ANM/Pracownia 1/prog/wynikiczas.jl

max_stars_repo_name: antonitomaszewski/Studia

max_issues_repo_name: antonitomaszewski/Studia

max_forks_repo_name: antonitomaszewski/Studia


text:
## --------------  Guess

export Guess
export make_guess, set_params_from_guess!, random_guess, n_values
export validate_guess

Guess{T}() where T = Guess(OrderedDict{ObjectId, ValueVector{T}}());

Base.isapprox(g1 :: Guess{T1}, g2 :: Guess{T2}; atol = 1e-8) where {T1, T2} = false;

function Base.isapprox(g1 :: Guess{T}, g2 :: Guess{T}; atol = 1e-8) where T
    if length(g1.d) != length(g2.d)
        return false
    end
    for (oId, vv1) in g1.d
        if haskey(g2.d, oId)
            vv2 = g2.d[oId]
            if !isapprox(vv1, vv2; atol)
                return false
            end
        else
            return false
        end
    end
    return true
end

lb(g :: Guess{F1}) where F1 = fill(F1(TransformationLb), (n_values(g), ));
ub(g :: Guess{F1}) where F1 = fill(F1(TransformationUb), (n_values(g), ));

function validate_guess(x :: ModelObject, g :: Guess{F1}) where F1
    isValid = validate_guess(g);
    isValid = isValid  &&  validate_param_match(x, g);
    return isValid
end

function validate_guess(g :: Guess{F1}) where F1
    isValid = true;
    n = n_values(g);
    if n > 0
        takenV = zeros(Int, n);
        for (_, vv) in g.d
            if !isempty(vv)
                isValid = isValid && validate_vv(vv);
                takenV[start_index(vv) : last_index(vv)] .+= 1;
            end
        end
        isValid = isValid  &&  all(takenV .== 1);
    end
    return isValid
end

# Check that parameters in Guess match calibrated parameters in ModelObject
function validate_param_match(x :: ModelObject, g :: Guess{F1}) where F1
    isValid = true;
    pvecV = collect_pvectors(x);
    if n_value_vectors(g) != length(pvecV)
        @warn "No of ParamVectors differs: $x";
        isValid = false;
    end
    if isValid
        for (_, pvec) in pvecV
            objId = get_object_id(pvec);
            nCal, _ = n_calibrated_params(pvec);
            if haskey(g.d, objId)
                vVec = g.d[objId];
                isValid = isValid  &&  validate_param_match(pvec, vVec);
            elseif (nCal > 0)
                @warn """
                    Guess has no entry for ParamVector with calibrated params
                    $objId
                    """;
                isValid = false;
            end
        end
    end
    return isValid
end

get_value_vector(guess :: Guess{T}, objId :: ObjectId) where T =
    guess.d[objId];

n_value_vectors(g :: Guess{T1}) where T1 = length(g.d);

"""
	$(SIGNATURES)

Total number of scalar values (elements in all arrays) in all parameters.
"""
function n_values(g :: Guess{T}) where T
    n = 0;
    for (_, vv) in g.d
        n += n_values(vv);
    end
    return n
end

"""
	$(SIGNATURES)

Values as a vector. Transformed.
"""
function get_values(o :: ModelObject, guess :: Guess{F1}) where F1
    pvecV = collect_pvectors(o);
    return get_values(pvecV, guess);
end

function get_values(pvecV :: PVectorCollection, guess :: Guess{F1}) where F1
    valueV = Vector{F1}();
    for (objId, vv) in guess.d
        pvec = find_pvector(pvecV, objId);
        @assert !isnothing(pvec);
        pValueV = get_values(pvec, vv);
        if !isempty(pValueV)
            append!(valueV, pValueV);
            @assert length(valueV) == last_index(vv);
        end
    end
    return valueV
end

# """
# 	$(SIGNATURES)

# Set values from a Vector{Real}. These are usually the transformed values.
# """
# function set_values!(guess :: Guess{F1}, valueV :: AbstractVector{F1}) where F1
#     for (_, vv) in guess.d
#         set_values!(vv, valueV)
#     end
#     return nothing
# end

"""
    $(SIGNATURES)

Make vector of parameters and bounds for an object.
Including nested objects.
This should be done only once to guarantee that order never changes.
"""
function make_guess(m :: ModelObject)
    pvecV = collect_pvectors(m);
    @assert !isempty(pvecV)  "$m contains no ParamVectors"
    vv = make_guess(pvecV);
    return vv
end

"""
    $(SIGNATURES)

Make vector from a list of param vectors.
Output contains values, lower bounds, upper bounds.
"""
function make_guess(pvv :: PVectorCollection)
    vv = Guess{ValueType}();
    startIdx = 1;
    for (objId, pv) in pvv
        vVec = make_value_vector(pv, startIdx);
        vv.d[objId] = vVec;
        startIdx += n_values(vVec);
    end
    @assert validate_guess(vv);
    return vv
end

"""
    $(SIGNATURES)

Make vector of values, lb, ub for optimization algorithm.
Vectors are transformed using the `ParameterTransformation` specified in the `ParamVector`.
Returns a `ValueVector`.
"""
function make_value_vector(pvec :: ParamVector,  startIdx :: Integer)
    T1 = ValueType;
    pList = calibrated_params(pvec);
    n, nElem = n_calibrated_params(pvec);

    vv = ValueVector{T1}();
    idxLast = startIdx - 1;
    for p in pList
        pValue = transform_param(pvec.pTransform,  p);
        pLen = length(pValue);
        # pIdxV = idxLast .+ (1 : pLen);
        vv.d[p.name] = ParamInfo(p.name, idxLast + 1,
            make_bounds(pvec.pTransform, size(pValue))...);
        idxLast += pLen;
    end
    @assert (idxLast - startIdx + 1) == nElem;
    return vv
end

function make_bounds(pTransform, sz)
    lbnd = lb(pTransform);
    ubnd = ub(pTransform);
    if !isempty(sz)
        lbnd = fill(lbnd, sz);
        ubnd = fill(ubnd, sz);
    end
    return lbnd, ubnd
end

"""
    $(SIGNATURES)

Make vector of guesses into model parameters. For object and children.
This changes the values in `m` and in its `pvector`.
"""
function set_params_from_guess!(m :: ModelObject, guess :: Guess{F1},
    guessV :: AbstractVector{F1}) where F1
    objV = collect_model_objects(m);
    # Copy param vectors into model
    for obj in objV
        set_own_params_from_guess!(obj, guess, guessV);
    end
end

"""
    $(SIGNATURES)

Copy values from a *vector* of `ParamVector` into a *vector* of `ModelObject`s.
The order of the objects must match the order of the `ParamVector`s.
The order of the values in `vAllInV` must match the order of `ParamVector`s.
Returns `true` if all values used; `false` otherwise.
Also ensures that all fixed parameters match `ParamVector`.
"""
function set_own_params_from_guess!(o, g :: Guess{F1}, guessV :: AbstractVector{F1}) where F1
    if has_pvector(o)
        objId = get_object_id(o);
        vv = get_value_vector(g, objId);
        set_own_params_from_vector!(o, vv, guessV);
        @assert check_own_calibrated_params(o, get_pvector(o))
        @assert check_own_fixed_params(o, get_pvector(o))
        # startIdx = idxLast + 1;
    end

    # if startIdx == length(vAll) + 1
    #     success = true;
    # else
    #     success = false;
    #     @warn "Not all values used: $(startIdx - 1)  vs  $(length(vAll))"
    # end
    # return success
end

"""
    $(SIGNATURES)

Copy values from vector into param vector and object.
Calibrated parameters.
Also ensures that fixed parameters are set according to `ParamVector`.
Uses the values in `vAll` starting from (optional) `startIdx`.
Returns index of last value used.
Order in `vAll` must match order in `pvec`. E.g., because `vAll` is generated by `make_guess`.
"""
function set_own_params_from_vector!(x :: ModelObject, vAll :: ValueVector{F1},
    guessV :: AbstractVector{F1}) where F1
    pvec = get_pvector(x);
    set_params_from_vector!(pvec, vAll, guessV);
    # d11, idxEnd = vector_to_dict(pvec, vAll, true; startIdx = startIdx);
    # set_values_from_dict!(pvec, d11);
    # This is key! It copies the values (calibrated and fixed) from the ParamVector into the ModelObject.
    sync_own_values!(x);
    return nothing
end

function set_params_from_vector!(pvec :: ParamVector, vv :: ValueVector{F1},
    guessV :: AbstractVector{F1}) where F1
    for (pName, pInfo) in vv.d
        set_value_from_pinfo!(pvec, pInfo, guessV; isCalibrated = true);
    end
end

function set_value_from_pinfo!(pvec :: ParamVector,
    pInfo :: ParamInfo{T},
    guessV :: AbstractVector{F1};
    isCalibrated = true) where {F1, T}

    p = retrieve(pvec, pInfo.pName);
    @assert size(calibrated_value(p)) == size(lb(pInfo));
    @assert !isnothing(p)  "Param $(pInfo.pName) not found in $pvec";
    if is_calibrated(p) == isCalibrated
        valV = reshape_vector(pInfo, guessV[indices(pInfo)]);
        v = untransform_param(pvec.pTransform, p, valV);
        set_calibrated_value!(p, v);
    end
end

function random_guess(g :: Guess{F1}, rng :: AbstractRNG) where F1
    lb(g) .+ (ub(g) .- lb(g)) .* rand(rng, n_values(g));
end

# ----------


text:
# Packages required for testing
using Test
using Random
using Logging
using LightGraphs
using TOML
# Package Under Test
using GraphUtils

# Set logging level
global_logger(SimpleLogger(stderr, Logging.Debug))

# Fix randomness during tests
Random.seed!(0)

# Check equality of two arrays
@inline function array_isapprox(x::AbstractArray{F},
                  y::AbstractArray{F};
                  rtol::F=sqrt(eps(F)),
                  atol::F=zero(F)) where {F<:AbstractFloat}

    # Easy check on matching size
    if length(x) != length(y)
        return false
    end

    for (a,b) in zip(x,y)
        @test isapprox(a,b, rtol=rtol, atol=atol)
    end
end

# Check if array equals a single value
@inline function array_isapprox(x::AbstractArray{F},
                  y::F;
                  rtol::F=sqrt(eps(F)),
                  atol::F=zero(F)) where {F<:AbstractFloat}

    for a in x
        @test isapprox(a, y, rtol=rtol, atol=atol)
    end
end

# Define package tests
@time @testset "GraphUtils Package Tests" begin
    testdir = joinpath(dirname(@__DIR__), "test")
    @time @testset "GraphUtils.GraphSorting" begin
        include(joinpath(testdir, "test_sorting.jl"))
    end
    @time @testset "GraphUtils.Angles" begin
        include(joinpath(testdir, "test_angles.jl"))
    end
    @time @testset "GraphUtils.Construction" begin
        include(joinpath(testdir, "test_construction.jl"))
    end
    @time @testset "GraphUtils.Connectivity" begin
        include(joinpath(testdir, "test_connectivity.jl"))
    end
    @time @testset "GraphUtils.CubicSplines" begin
        include(joinpath(testdir, "test_cubic_splines.jl"))
    end
    @time @testset "GraphUtils.Graphs" begin
        include(joinpath(testdir, "test_graphs.jl"))
    end
    @time @testset "GraphUtils.Trees" begin
        include(joinpath(testdir, "test_abstract_trees.jl"))
    end
    @time @testset "GraphUtils.CachedElements" begin
        include(joinpath(testdir, "test_cached_element.jl"))
    end
    @time @testset "GraphUtils.AbstractIDs" begin
        include(joinpath(testdir, "test_abstract_ids.jl"))
    end
    @time @testset "GraphUtils.FactoryWorlds" begin
        include(joinpath(testdir, "test_factory_worlds.jl"))
    end
    @time @testset "GraphUtils.FileSystem" begin
        include(joinpath(testdir, "test_filesystem.jl"))
    end
    @time @testset "GraphUtils.Printing" begin
        include(joinpath(testdir, "test_printing.jl"))
    end
    @time @testset "GraphUtils.Logging" begin
        include(joinpath(testdir, "test_logging.jl"))
    end
end


text:
import Robotlib: T2R, Rt2T
"""
    T_TF_S, meanRMS,RMS, norms = calibNAXP(points_S, lines_S, POSES, T_TF_S, planes::AbstractVector{Int},  iters = 50; doplot=false)

This functions implements the algorithm from the paper
"Six DOF eye-to-hand calibration from 2D measurements using planar constraints"
which solves the problem `Prb = n'A(X*Ps)`
If result is bad, check if you send data in correct form\n
`POSES ∈ R(4,4,N)` is always the position of the tool frame from the robot FK\n
`points_S ∈ R(3,N)` are measurements from a line laser scanner. Given
in the sensor frame. This script assumes that the laser plane is the
sensor XY plane with y-axis pointing outwards from the sensor.\n
`lines_S ∈ R(3,N)` are vectors corresponding to the direction of the measured laser line
`planes ∈ Z(N)` is a vector of indices corresponding to which plane a mesurement
comes from. It must be same length as `points_S` and enumerate the planes starting at 1. Example (N=15, number of planes = 3): [1,1,1,1,1,2,2,2,2,2,3,3,3,3,3]\n
`iters` determines the number of iterations\n
@inproceedings{carlson2015six,
  title={Six DOF eye-to-hand calibration from 2D measurements using planar constraints},
  author={Carlson, Fredrik Bagge and Johansson, Rolf and Robertsson, Anders},
  booktitle={Intelligent Robots and Systems (IROS), 2015 IEEE/RSJ International Conference on},
  pages={3628--3632},
  year={2015},
  organization={IEEE}
}
"""
function calibNAXP(points_S, lines_S, POSES, T_TF_S, planes::AbstractVector{Int},  iters::Integer = 50; doplot=false, trueT_TF_S=nothing, variance=nothing)
    N_poses   = size(POSES,3)
    N_planes  = maximum(planes)
    RMScalibs = zeros(iters)
    ALLcalibs = zeros(iters, N_planes)
    norms     = trueT_TF_S == nothing ? nothing : zeros(iters)
    local points, lines, N_RB
    for c = 1:iters
        # Convert points to RB cordinate system
        points, lines = transform_data(POSES, points_S, lines_S, T_TF_S)
        N_RB = findplanes(planes, points)

        # Estimation
        if variance == nothing
            A,y = get_matrices(N_RB, POSES, points_S, lines_S)
            w   = A\y
        else
            A,y,Σy = get_matrices(N_RB, POSES, points_S, lines_S, variance)
            w   = wls(A,y,Σy)
        end
        # w   = tls(A,y)
        R,t = w2Rt(w)

         # Re-estimation of translation
        w₂     = R[1:3,1:2][:]
        y₂     = y - A[:,1:6]*w₂
        t₂     = A[:,7:9]\y₂
        T_TF_S = Rt2T(R,t₂)

        # Calculate iteration errors
        RMSi = map(1:N_planes) do j
            ind     = planes .== j
            ind     = findall(ind)
            sqrt(pointDiff(T_TF_S,POSES[:,:,ind],points_S[1:3,ind])[1])
        end
        RMScalibs[c]   = mean(RMSi)
        ALLcalibs[c,:] = RMSi
        trueT_TF_S == nothing || (norms[c] = norm(T_TF_S-trueT_TF_S))
        #  any(abs(RMSi) > 1e-1) && warn("Points does not seem to lie on a plane")
    end
    if doplot
        try
            plotPlanes(unique(N_RB,dims=2))
            plotLines(points,lines)
        catch ex
            if ex isa UndefVarError
                @error "Plotting failed. Plotting is only available if you have manually run `using Plots`"
            else
                rethrow(ex)
            end
        end
    end
    T_TF_S, RMScalibs,ALLcalibs, norms
end

function w2Rt(w)
    H   = reshape(w,3,3)
    Rx  = H[1:3,1]
    Ry  = H[1:3,2]
    Rz  = cross(Rx,Ry) # These lines must be changed if laser plane is changed
    Rnn = [Rx Ry Rz]
    # @show Rnn
    R   = orthonormal(Rnn)
    t   = H[1:3,3]
    R,t
end

function get_matrices(N_RB, POSES, points_S, lines_S,variance=nothing)
    N_poses = size(POSES,3)
    A = zeros(2*N_poses,9)
    y = zeros(2*N_poses)
    𝜸 = 1
    if variance != nothing
        α = mapslices(lines_S, dims=1) do l
            acos(l[2]) # = acos(eʸ'l)
        end[:]
        Σ(α) = [sin(α)^2 -sin(α)*cos(α); -sin(α)*cos(α) cos(α)^2]
        Σp,σα = variance
        Σp2 = [Σp .+ 𝜸^2*σα^2*abs.(Σ(α)) for α in α]
        Σy = spzeros(2*N_poses,2*N_poses)
    end
    for i = 1:N_poses
        Ni     = N_RB[:,i]
        Ra     = POSES[1:3,1:3,i]
        Ta     = POSES[1:3,4,i]
        Pi     = points_S[:,i]
        y[i]   = norm(Ni)^2-Ni'Ta
        A[i,:] = (reshape(repeat(Ni'Ra,3,1)',9,1).*[reshape(repeat(Pi[1:2],1,3)',6,1);1;1;1])' #TODO: rewrite
        # Sample one additional point along line
        Pi = points_S[:,i] + 𝜸*lines_S[:,i]
        y[i+N_poses] = norm(N_RB[:,i])^2-Ni'Ta
        A[i+N_poses,:] = (reshape(repeat(Ni'Ra,3,1)',9,1) .*[reshape(repeat(Pi[1:2],1,3)',6,1);1;1;1])'
        if variance != nothing
            Σy[i,i] = Σp
            Σy[i+N_poses, i+N_poses] = tr(Σp2[i])
            Σy[i+N_poses, i] = Σy[i, i+N_poses] = Σp
        end
    end
    variance == nothing ? (A,y) : (A,y,Symmetric(Σy))
end

function transform_data(POSES, points_S, lines_S, T_TF_S)
    N_poses = size(POSES,3)
    points  = zeros(4,N_poses)
    lines   = similar(lines_S)
    for i = 1:N_poses
        T_RB_S      = POSES[:,:,i] * T_TF_S
        points[:,i] = T_RB_S*[points_S[:,i]; 1]
        lines[:,i]  = T2R(T_RB_S)*lines_S[:,i]
    end
    points, lines
end

function findplanes(planes, points, rep=true)
    N_poses  = length(planes)
    N_planes = maximum(planes)
    N_out    = rep ? N_poses : N_planes
    μ_RB     = zeros(3,N_out)
    N_RB     = μ_RB
    for j = 1:N_planes
        pointind           = planes .== j
        planeind           = rep ? pointind : j

        μᵢ                 = mean(points[1:3,pointind],dims = 2)[:]
        μ_RB[:,planeind]  .= μᵢ
        D,V                = eigen(cov(points[1:3,pointind]')) # PCA
        N_RBi              = V[:,1] # Eigenvector of smalles eigval is the normal of the plane
        if μᵢ'N_RBi < 0 # Make sure all normals point in same direction
            N_RBi = -1*N_RBi
        end
        N_RBi             = N_RBi*(N_RBi'μᵢ)
        N_RB[:,planeind] .= N_RBi # Repeat plane normal for all planes with same ind
    end
    N_RB
end

function pointDiff(T_TF_S, T_RB_TF, points_S)
    N = size(points_S,2)
    @assert size(T_RB_TF,3) == N
    points_RB = zeros(4,N)
    for i = 1:N
        T_RB_S = T_RB_TF[:,:,i]*T_TF_S
        points_RB[:,i] = T_RB_S*[points_S[:,i];1]
    end
    points_RB = points_RB[1:3,:]
    S = eigvals(cov(points_RB'))
    d = minimum(S)
end

## The code below is used to verify the statements made in the appendix to chapter 13 in the PhD thesis cited in the README.

# function bootstrap(fun::Function, estimate, N, k=100)
#     k < 1 && return (fun(1:N), nothing)
#     results = map(1:k) do _
#         sample  = rand(1:N, N)
#         result  = fun(sample)
#     end
#     fun(1:N), estimate(results)
# end
#
#
# function columnwise_cov(data)
#     numcols = size(data[1], 2)
#     map(1:numcols) do c
#         columndata = hcat([d[:,c] for d in data]...)'
#         cov(columndata)
#     end
# end
#
# """
# This function is the same as `calibNAXP` except it estimates the covariance of the estimated normals using bootstrapping, and solves the optimization problem using the weighted total least-squares algorithm instead of standard least-square. It does not really provide any additional benefints over the standard algorithm, but is considerable slower.
# """
# function calibNAXP_bootstrap(points_S, lines_S, POSES, T_TF_S, planes::AbstractVector{Int},  iters::Integer = 50; trueT_TF_S=nothing, nbootstrap=500)
#     N_poses   = size(POSES,3)
#     N_planes  = maximum(planes)
#     RMScalibs = zeros(iters)
#     ALLcalibs = zeros(iters, N_planes)
#     norms     = trueT_TF_S == nothing ? nothing : zeros(iters)
#     for c = 1:iters
#         # Convert points to RB cordinate system
#         points, lines = transform_data(POSES, points_S, lines_S, T_TF_S)
#         N_RB, Σn = Robotlib.Calibration.bootstrap(columnwise_cov, N_poses, nbootstrap) do sample
#             findplanes(planes[sample], points[:,sample], false)
#         end
#         N_RB = N_RB[:,planes]
#         Σn = Σn[planes]
#         rowQ = map(1:N_poses) do i
#             Q,T,p,n = Σn[i], POSES[:,:,i], points_S[:,i], N_RB[:,i]
#             R = T[1:3,1:3]
#             t = T[1:3,4]
#             Qaa = R'Q*R
#             σ = [p[1], p[2], 1]*[p[1], p[2], 1]'
#             QAA = kron(ones(3,3), Qaa) .* kron(σ, ones(3,3))
#             Qyy = t'Q*t # TODO: + cov of norm(n)^2 Q
#             Qay = kron(ones(3), R'Q*(n-t)) .* kron([p[1], p[2], 1], ones(3))
#
#             [QAA Qay; [Qay' Qyy]]
#         end
#         Qaa,Qay,Qyy = rowcovariance(rowQ)
#         # Estimation
#         A,y = get_matrices(N_RB, POSES, points_S, lines_S)
#         w   = wtls(A,y,kron([1 0.5;0.5 1],Qaa),kron([1 0.5;0.5 1],Qay),kron([1 0.5;0.5 1],Qyy), iters=15)
#         # w   = tls(A,y)
#         R,t = w2Rt(w)
#
#          # Re-estimation of translation
#         w₂     = R[1:3,1:2][:]
#         y₂     = y - A[:,1:6]*w₂
#         t₂     = c <= 5 ? A[:,7:9]\y₂ : tls(A[:,7:9],y₂)
#         T_TF_S = Rt2T(R,t₂)
#
#         # Calculate iteration errors
#         RMSi = map(1:N_planes) do j
#             ind     = planes .== j
#             ind     = findall(ind)
#             sqrt(pointDiff(T_TF_S,POSES[:,:,ind],points_S[1:3,ind])[1])
#         end
#         RMScalibs[c]   = mean(RMSi)
#         ALLcalibs[c,:] = RMSi
#         trueT_TF_S == nothing || (norms[c] = norm(T_TF_S-trueT_TF_S))
#         #  any(abs(RMSi) > 1e-1) && warn("Points does not seem to lie on a plane")
#     end
#     T_TF_S, RMScalibs,ALLcalibs, norms
# end


text:
push!(LOAD_PATH, "../")

using Documenter, NumerAPI

makedocs(
    sitename = "NumerAPI.jl",
    format = Documenter.HTML(),
    modules = [NumerAPI],
    pages = [
      "Home" => "index.md",
			"tournament.md",
			"signals.md"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
#=deploydocs(
    repo = "<repository url>"
)=#


text:
module TestBangBang
using Test

@testset "$file" for file in sort([file for file in readdir(@__DIR__) if match(
    r"^test_.*\.jl$",
    file,
) !== nothing])

    if VERSION < v"1.2-" && file == "test_zygote.jl"
        @info "Skip $file for Julia $VERSION to avoid segfault"
        # For example:
        # https://travis-ci.com/tkf/BangBang.jl/builds/134791701
        # https://travis-ci.com/tkf/BangBang.jl/builds/134795161
        continue
    end

    include(file)
end

end  # module


text:
using DataDeps

register(DataDep(
  "wiki4he",
  "https://archive.ics.uci.edu/ml/datasets/wiki4HE",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/00334/wiki4HE.csv",
  "2270c9b4cc6a9492c740f2bd42a2fee23c7afab383ec31364b08365cc1d43ae2",
  post_fetch_method=(path -> begin
    UCIData.process_dataset(path,
        target_index=1,
        feature_indices=2:53,
        categoric_indices=[2:4; 6:10],
        delim=';',
        skipto=2,
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data


text:
import FastGaussQuadrature
import FastTransforms

using DecFP: Dec128

using ..CommonFunctions

struct Quadrature{T,N}
    order::Int
    nodes::Vector{T}
    weights::Vector{T}
end

function Quadrature(order, nodes::Vector{T}, weights::Vector{T}) where {T}
    @assert length(nodes) == length(weights)
    Quadrature{T, length(nodes)}(order, nodes, weights)
end

CommonFunctions.nnodes(Q::Quadrature{T,N}) where {T,N} = N
CommonFunctions.nodes(Q::Quadrature) = Q.nodes
CommonFunctions.order(Q::Quadrature) = Q.order
weights(Q::Quadrature) = Q.weights
Base.eachindex(Q::Quadrature) = eachindex(Q.nodes, Q.weights)

"""
Integrate a function `f(x)` over the interval [0,1] using the quadrature `quad`.
"""
function quadrature(quad::Quadrature, f::Function)
    sum(quad.weights[i] * f(quad.nodes[i]) for i in eachindex(quad))
end

"Scale nodes and weights from the interval [-1,+1] to the interval [0,1]."
function shift!(b,c)
    b .= b ./ 2
    c .= (c .+ 1) ./ 2
end

"Scale nodes and weights from the interval [0,1] to the interval [-1,+1]."
function unshift!(b,c)
    b .= b .* 2
    c .= 2 .* c .- 1
end

function RiemannQuadratureLeft(T=Float64)
    Quadrature(1, [0.0,], [1.0,])
end

function RiemannQuadratureRight(T=Float64)
    Quadrature(1, [1.0,], [1.0,])
end

function GaussLegendreQuadrature(s, T=Float64; shift=true)
    c, b = FastGaussQuadrature.gausslegendre(s)
    if shift shift!(b,c) end
    Quadrature(2s, c, b)
end

function LobattoLegendreQuadrature(s, T=Float64; shift=true)
    c, b = FastGaussQuadrature.gausslobatto(s)
    if shift shift!(b,c) end
    Quadrature(2s-2, c, b)
end

function get_nodes_gauss_chebyshev(s, k=1, T=Float64)
    x = zeros(T, s)
    if k == 1
        for i in eachindex(x)
            x[i] = @dec128 ( sin( π*(s-2i+1) / (2s) ) + 1 ) / 2
        end
        return x[end:-1:1]
    elseif k == 2
        for i in eachindex(x)
            x[i] = @dec128 ( cos( π*i / (s-1) ) + 1 ) / 2
        end
        return x
    end
end

function GaussChebyshevQuadrature(s, k=1, T=Float64; shift=true)
    c = get_nodes_gauss_chebyshev(s,k,T)
    b = zeros(c)

    tj::Dec128 = 0
    th::Dec128 = 0

    if k == 1
        for i in eachindex(b)
            tj = 0
            th = @dec128 π * (2i-1) / (2s)
            for j in 1:div(s,2)
                tj += @dec128 cos(2j*th) / (4j^2 - 1)
            end
            b[i] = @dec128 (1 - 2tj) / s
        end
    elseif k == 2
        for i in eachindex(b)
            tj = 0
            th = @dec128 π * i / (s+1)
            for j in 1:div(s+1,2)
                tj += @dec128 sin((2j-1) * th) / (2j - 1)
            end
            b[i] = @dec128 (2tj * sin(th)) / (s+1)
        end
    else
        error("GaussChebyshevQuadrature: k=", k, " not supported.")
    end

    if !shift unshift!(b,c) end

    Quadrature(2s-2,c,b)
end

function get_nodes_lobatto_chebyshev(s, T=Float64)
    x = zeros(T, s)
    for i in eachindex(x)
        x[i] = @dec128 ( cos( (i-1)*π / (s-1) ) + 1 ) / 2
    end
    return x[end:-1:1]
end

function LobattoChebyshevQuadrature(s, T=Float64)
    c = get_nodes_lobatto_chebyshev(s)
    b = zeros(c)

    # TODO compute weights

    error("LobattoChebyshevQuadrature not implemented, yet!")
end

function ClenshawCurtisQuadrature(s, T=Float64; shift=true)
    c, b = FastTransforms.clenshawcurtis(s, zero(T), zero(T))
    c = flipdim(c,1)
    if shift shift!(b,c) end
    Quadrature(s^2, c, b)
end


text:
function adagram_isblank(c::Char)
  return c == ' ' || c == '\t'
end

function adagram_isblank(s::AbstractString)
  return all((c->begin
            c == ' ' || c == '\t'
        end),s)
end

@resumable function word_iterator(f::IO, end_pos::Int64=-1) :: AbstractString
  while (end_pos < 0 || position(f) < end_pos) && !eof(f)
    w = readuntil(f, ' ')
    if length(w) < 1 break end
    if !adagram_isblank(w)
        @yield w
    end
  end
end

@resumable function looped_word_iterator(f::IO, start_pos::Int64,
  end_pos::Int64) :: AbstractString
  while !eof(f)
    w = readuntil(f, ' ')
    if length(w) < 1 break end
    if !adagram_isblank(w)
      @yield w
    end
    if position(f) >= end_pos seek(f, start_pos) end
  end
end

function count_words(f::IOStream, min_freq::Int=5)
  counts = Dict{AbstractString, Int64}()

  for word in word_iterator(f)
    if get(counts, word, 0) == 0
      counts[word] = 1
    else
      counts[word] += 1
    end
  end

  for word in [keys(counts)...]
    if counts[word] < min_freq
      delete!(counts, word)
    end
  end

  V = length(counts)
  id2word = Array(AbstractString, V)
  freqs = zeros(Int64, V)
  i = 1
  for (word, count) in counts
    id2word[i] = word
    freqs[i] = count
    i += 1
  end

  return freqs, id2word
end

function align(f::IO)
  while !adagram_isblank(read(f, Char))
    continue
  end

  while adagram_isblank(read(f, Char))
    continue
  end

  seek(f, position(f)-1)
end

function read_words(f::IO,
    dict::Dictionary, doc::DenseArray{Int32},
    batch::Int, last_pos::Int)
  words = Stateful(word_iterator(f, last_pos))
  i = 1
  for word in words
    if i > batch break end

    id = get(dict.word2id, word, -1)
    if id == -1 continue end

    doc[i] = id
    i += 1
  end

  return view(doc, 1:i-1)
end

function read_words(str::AbstractString,
    dict::Dictionary, doc::DenseArray{Int32},
    batch::Int, last_pos::Int)
  i = 1
  for word in split(str, ' ')
    id = get(dict.word2id, word, -1)
    if id == -1
      continue
    end

    doc[i] = id
    i += 1
  end

  return view(doc, 1:i-1)
end

function read_words(f::IOStream, start_pos::Int64, end_pos::Int64,
    dict::Dictionary, doc::DenseArray{Int32},
    freqs::DenseArray{Int64}, threshold::Float64,
    words_read::DenseArray{Int64}, total_words::Float64)
  words = Stateful(looped_word_iterator(f, start_pos, end_pos))
  i = 1
  while i <= length(doc) && words_read[1] < total_words && !isempty(words)
    word = ""
    try
      word = popfirst!(words)
    catch e
      println("EOF exception")
      break
    end
    id = get(dict.word2id, word, -1)
    if id == -1
      continue
    elseif rand() < 1. - sqrt(threshold / (freqs[id] / total_words))
      words_read[1] += 1
      continue
    end

    doc[i] = id
    i += 1
  end

  return view(doc, 1:i-1)
end


text:
#####################
# Generic transforms
#####################

# A function set can have several associated transforms. The default transform is
# associated with the grid of the set, e.g. the FFT and the DCTII for Chebyshev expansions
# which convert between coefficient space and value space. In this case, the
# transform maps coefficients to or from a GridBasis.
#
# transform_operator takes two arguments, a source and destination set, in order
# to allow for different transforms.
#
# An important case is where the source or destination is a GridBasis. In that case,
# the routines transform_from_grid and transform_to_grid are invoked,
# with the grid as a third argument.
# Dictionary's can intercept these functions and define their coefficient-to-value
# and value-to-coefficient transforms.

# The default transform space is the space associated with the grid of the set
transform_dict(Φ::Dictionary; options...) = GridBasis(Φ)

transform(src::Dictionary; options...) = transform(src, transform_dict(src); options...)

transform(src::Dictionary, dest::Dictionary; options...) =
    transform(operatoreltype(src, dest), src, dest; options...)

# If source or destination is a GridBasis, specialize
transform(T, src::Dictionary, dest::GridBasis; options...) =
    transform_to_grid(T, src, dest, grid(dest); options...)
transform(T, src::GridBasis, dest::Dictionary; options...) =
    transform_from_grid(T, src, dest, grid(src); options...)

transform_to_grid(T, src::Dictionary, dest::GridBasis; options...) =
    transform_to_grid(T, src, dest, grid(dest); options...)
transform_from_grid(T, src::GridBasis, dest::Dictionary; options...) =
    transform_from_grid(T, src, dest, grid(src); options...)

# Convert grid to GridBasis
transform(T, src::Dictionary, grid::AbstractGrid; options...) =
    transform_to_grid(T, src, GridBasis{T}(grid); options...)
transform(T, grid::AbstractGrid, dest::Dictionary; options...) =
    transform_from_grid(T, GridBasis{T}(grid), dest; options...)


text:
"""
symmetrize matrix
"""
sym(M::Matrix{Float64}) = (M' + M) / 2

# link functions:
invLogit(xb::Float64) = 1 / (1 + exp(-xb))
logit(p::Float64) = log(p / (1-p))
# end of Link functions

const quants = [0., .025, .25, .5, .75, .975, 1.]

rig{A <: Real, B <: Real}(a::A,b::B) = 1 / rand( Gamma(a,1/b) )

function predict(X0::Matrix{Float64},β::Matrix{Float64}, inv_link)
  inv_link.(X0 * β)
end

function dic{T}(post_samples::Vector{T},loglikelihood)
  const D = -2 * loglikelihood.(post_samples)
  return mean(D) + var(D) / 2
end

"""
cpo = conditional posterior ordinate.
      which is the leave-one-out density evaluated at the
      actual observed y's. used for *comparing* models. Not
      for prediction.
"""
function cpo{T}(post_params::Vector{T}, den,
                y::Vector{Float64},
                X::Matrix{Float64})

  const N = size(X,1)
  const J = length(post_params)
  const CPO = Matrix{Float64}(N,J)

  for i in 1:N
    for j in 1:J
      CPO[i,j] = 1 / den(post_params[j], y[i], X[i,:])
    end
  end

  return 1 ./ mean(CPO,2)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 125fdd06-d330-11eb-0cf4-17cfac8cc4f7
begin
	using Plots,DataFrames,CSV,Dates,GLM
	# gr()
	# plotly()
	pyplot()
	# pgfplotsx()
	default(fontfamily="Gill Sans",
		titlefont = (12,:grey),
		legendfont = (11),
		guidefont = (10),
		tickfont = (9),
		annotationfontsize=(8),
		annotationcolor=:blue

	  )
end

# ╔═╡ a1341efb-8d58-4b9a-8071-6f3e665273ec
md"# PIP and DLA receipts
source: [Stat-Xplore](https://stat-xplore.dwp.gov.uk/);
retrieved 02/July/2021
 "

# ╔═╡ 6ee3c8bd-5dca-4c7d-9e41-be59c2ba2b7e
md"## DLA "

# ╔═╡ 923d961b-0912-4d1e-976e-20b89150b563
begin
	cd("/home/graham_s/julia/vw/ScottishTaxBenefitModel/")
	dla=CSV.File("data/receipts/dla_2002-2020_from_stat_explore.csv")|>DataFrame
	dla.Date = Date.( dla.Date, dateformat"u-yy" ) .+Year(2000)
	dla
end

# ╔═╡ fd0658ac-4fef-471a-af49-db01751ef2a9
begin
	p = plot(dla.Date,[dla.Scotland dla.Wales dla.England], ylims = (0, 3_000_000), labels=["DLA-Sco" "DLA-Wal" "DLA-Eng"], formatter=:plain )

end

# ╔═╡ 1686c58b-ebb0-4502-8021-7842e369eef8
md"## Personal Independence Payment"

# ╔═╡ a01723f5-e22a-4d36-9afa-6675d591bfeb
begin
	pip=CSV.File( "data/receipts/pip_2002-2020_from_stat_explore.csv",
		missingstrings=[".."],
		types=Dict([:Date=>String]))|>DataFrame
	pip.Date = Date.( pip.Date, dateformat"yyyymm" )
	pip
end

# ╔═╡ 8dabffad-137c-4c91-8406-96127b960898
begin
	pipp = plot!(p,pip.Date,[pip.Scotland,pip.Wales,pip.England],title="DLA/PIP Cases UK 2002-2020",labels=["PIP-Sco" "PIP-Wal" "PIP-Eng"])
	# svg( pipp, "docs/pip-dla.svg" )
	savefig(pipp, "docs/pip-dla-uk.svg" )
	pipp
end

# ╔═╡ c9c0cb68-9bd6-4044-a92b-e83532dce766
begin
	psco = plot(dla.Date,dla.Scotland, labels="DLA", formatter=:plain )
	pscop = plot!( psco, pip.Date, pip.Scotland, title="PIP/DLA Scotland", labels="PIP" )
	savefig(pscop, "docs/pip-dla-scotland.svg" )
	pscop
end

# ╔═╡ Cell order:
# ╠═125fdd06-d330-11eb-0cf4-17cfac8cc4f7
# ╟─a1341efb-8d58-4b9a-8071-6f3e665273ec
# ╠═6ee3c8bd-5dca-4c7d-9e41-be59c2ba2b7e
# ╠═923d961b-0912-4d1e-976e-20b89150b563
# ╠═fd0658ac-4fef-471a-af49-db01751ef2a9
# ╠═1686c58b-ebb0-4502-8021-7842e369eef8
# ╠═a01723f5-e22a-4d36-9afa-6675d591bfeb
# ╠═8dabffad-137c-4c91-8406-96127b960898
# ╠═c9c0cb68-9bd6-4044-a92b-e83532dce766


text:
module GeneralizedNullspaces
# copyright (c) 2017-2019 Ralph Smith
# SPDX-License-Identifier: MIT-Expat
# License-Filename: LICENSE.md

export gnsd, gnsd!, GNSD

using LinearAlgebra

include("QRDowndate.jl")
using .QRDowndate

"""
    GNSD <: Factorization

A generalized nullspace decomposition of a matrix, typically obtained from
[`gnsd`](@ref). If `A` is a square matrix, then the GNSD is `A == V*B*V'`
where `V` is a unitary matrix and `B` is block upper triangular.
The diagonal blocks of `B` are square and null, with orders listed in
a vector `μ`, except for a possible last block which
will represent the range of `A^ν` if it is not empty, where `ν == length(μ)`.
The order `μ[j]` is the number of zero-diagonal Jordan blocks of order `≥ j`
in the Jordan decomposition of a matrix close to `A`.
The blocks on the first superdiagonal of `B` have full column rank.
Columns of `V` span the corresponding subspaces of `A`.
"""
struct GNSD{T,M<:AbstractMatrix} <: Factorization{T}
    B::M
    V::M
    μ::Vector{Int}
    GNSD{T,M}(B::AbstractMatrix{T},V::AbstractMatrix{T},μ::Vector{Int}) where {T,M} =
        new(B,V,μ)
end
GNSD(B::AbstractMatrix{T},V::AbstractMatrix{T},μ::Vector{Int}) where {T} = GNSD{T,typeof(B)}(B,V,μ)

const pverbosity = Ref{Int}(0)
function set_verbosity(k)
    pverbosity[] = k
end

const useviews = Ref{Bool}(true)

"""
    gnsd(A) -> F::GNSD

Compute the generalized nullspace decomposition of a matrix `A`.

Returns a [`GNSD`](@ref) object.
`B` and `V` can be obtained from the factorization `F` with `F[:B]` and
`F[:V]`, such that `A = V * B * V'`.

Iterating the decomposition produces the tuple `(B, V, ν, μ)`
where `ν` is the index of `A` and `μ` holds the diagonal block orders.
See [`GNSD`](@ref) for the structure of `B`.

Uses the QR/update algorithm from [^Guglielmi2015].

[^Guglielmi2015]: N.Guglielmi, M.Overton, & G.Stewart,
"An efficient algorithm for computing the generalized null space decomposition",
SIAM J. Matrix Anal. Appl. 36 (1), 38-54 (2015).
"""
gnsd(A::Matrix{T},tol=sqrt(eps(T))) where T = gnsd!(copy(A),tol)

"""
    gnsd!(A) -> F::GNSD

`gnsd!` is the same as [`gnsd`](@ref), but saves space by
overwriting the input `A` instead of creating a copy.
"""
function gnsd!(B::Matrix{T},tol=sqrt(eps(T))) where T
    verbosity = pverbosity[]
    m,n = size(B)
    F = qr(B)
    Q = Matrix(F.Q)
    R = F.R
    μ = Vector{Int}(undef,0) # sigh.
    # we will multiply rotations into V
    V = Matrix{T}(I,size(B))

    k = 1 # block number
    idxk = 1 # starting index of block k
    local rxnorm

    while true
        (verbosity > 0) && println("starting block $k")
        j = idxk-1 # pointer to working row
        while j<n
            j += 1
            # R22 = R[j:n,j:n]
            if useviews[]
                x,rxnorm = nullvector(view(R,j:n,j:n))
            else
                x,rxnorm = nullvector(R[j:n,j:n])
            end
            (verbosity > 0) && println("row j=$j rxnorm: $rxnorm")
            (rxnorm > tol) && break
            if j>1
                x = vcat(zeros(T,j-1),x)
            end

            # Transposition is reversed w.r.t. paper because
            # of Julia Givens convention.

            # Find unitary Vi s.t. Vi * x = e1 (up to factor)
            # (Vi stands for $V_{i,i+1}$ in the paper.)
            # Find unitary U s.t. R = U * R22 * Vi' where R22 is Up.Tri.,
            # Save R22 in R
            # Update Q <- Vi * Q * U'
            # Accumulate Vi into V
            # Update B <- Vi * B * Vi'

            for i=n:-1:j+1
                Vi,r = givens(x[i-1],x[i],i-1,i)
                x[i-1] = r
                x[i] = 0
                if useviews[]
                    rmul!(view(R,1:i,:),adjoint(Vi)) # slow
                else
                    R[1:i,:] .= rmul!(R[1:i,:],adjoint(Vi)) # allocates
                end
                rmul!(V,adjoint(Vi))
                rmul!(B,adjoint(Vi))
                lmul!(Vi,B)
                lmul!(Vi,Q)
                Ui,r = givens(R[i-1,i-1],R[i,i-1],i-1,i)
                R[i-1,i-1] = r
                R[i,i-1] = 0
                if useviews[]
                    lmul!(Ui,view(R,:,i:n)) # slow
                else
                    R[:,i:n] = lmul!(Ui,R[:,i:n]) # allocates
                end
                rmul!(view(Q,idxk:n,:),adjoint(Ui))

            end
            # Find unitary W s.t. Rtilde = W' * Rhat is Up.Tri.,
            #    w/ zero first row
            # Update Q <- Q * W
            for i=j+1:n
                W,r = givens(R[i,i],R[j,i],i,j)
                R[i,i] = r
                R[j,i] = 0
                if useviews[]
                    lmul!(W,view(R,:,i+1:n)) # slow
                else
                    R[:,i+1:n] = lmul!(W,R[:,i+1:n]) # slow
                end
                rmul!(view(Q,idxk:n,:),adjoint(W))
            end
        end

        if j==n && rxnorm < tol
            # done (nilpotent)
            push!(μ,j-idxk+1)
            break
        elseif j==idxk
            # done (nonsingular last block, not counted in μ)
            break
        else
            # downdate to prepare for next block
            push!(μ,j-idxk)
            for i=idxk:j-1
                if useviews[]
                    qrdowndate!(view(Q,i:n,j:n),view(R,j:n,j:n))
                else
                    Qxx = Q[i:n,j:n]
                    Rxx = R[j:n,j:n]
                    qrdowndate!(Qxx,Rxx)
                    Q[i:n,j:n] = Qxx
                    R[j:n,j:n] = Rxx
                end
            end
            # ready for next block
            idxk=j
            k += 1
        end
    end

    # Actually null out the (sub)diagonal blocks
    i=1 # first index of current block
    for j in eachindex(μ)
        B[i:n, i:i+μ[j]-1] .= 0
        i += μ[j]
    end
    GNSD(B,V,μ)
end

"""
    nullvector(R, tol=0) -> x, nRx

Compute an approximate null vector `x` for upper triangular matrix `R`.
Also returns the norm of `R * x`, which is a rough approximation to
the smallest singular value of `R`.
"""
function nullvector(R::AbstractMatrix{T}, tol=zero(real(T))) where T
    # For now, we allow R to be a view.
    all(x->isa(x, Base.OneTo), axes(R)) ||
        throw(ArgumentError("only implemented for 1-based indexing"))
    n = size(R,1)
    x = zeros(T,n)

    # Case 1: tiny/zero diagonal element
    # Do a partial upper triangular back-solve, but skip the initial division.

    # Caveat: If `tol` is too small, and it is met at some index, but
    # should have been met earlier, the backsolve is unstable.
    # Consider using a better default or rescaling.
    idxsmall = findfirst(x -> (abs(x) <= tol), diag(R))
    if (idxsmall != nothing) && (idxsmall > 0)
        b = zeros(T,idxsmall); b[end] = 1
        x[idxsmall] = 1
        @inbounds for i in idxsmall-1:-1:1
            b[i] -= R[i,idxsmall]
        end
        @inbounds for j in idxsmall-1:-1:1
            xj = x[j] = R[j,j] \ b[j]
            for i in j-1:-1:1
                b[i] -= R[i,j] * xj
            end
        end
        xnrm = norm(x)
        rmul!(x,1/xnrm)
        rxnorm = norm(R*x)
    else
        # Case 2:
        # The paper suggests a LINPACK norm-estimation scheme.
        # This is from GVL4, p.141, but loop is reversed to give a 1-norm.
        p = zeros(T,n)
        y = zeros(T,n)
        pp = zeros(T,n)
        pm = zeros(T,n)
        for k=1:n
            yp = (1-p[k])/R[k,k]
            ym = (-1-p[k])/R[k,k]
            pp[k:end] .= p[k:end] .+ yp .* R[k:end,k]
            pm[k:end] .= p[k:end] .+ ym .* R[k:end,k]
            if abs(yp)+norm(pp,1) >= abs(ym)+norm(pm,1)
                y[k] = yp
                p[k:end] = pp[k:end]
            else
                y[k] = ym
                p[k:end] = pm[k:end]
            end
        end
        normalize!(y)

        x = UpperTriangular(R) \ y
        # We often need a test for nullity later, viz. ‖R x‖;
        # return it from here to avoid an extra MxV.
        rxnorm = 1/norm(x)
        rmul!(x,rxnorm)
    end
    x,rxnorm
end

# iteration for destructuring into components
Base.iterate(S::GNSD) = (S.B, Val(:V))
Base.iterate(S::GNSD, ::Val{:V}) = (S.V, Val(:ν))
Base.iterate(S::GNSD, ::Val{:ν}) = (S.ν, Val(:μ))
Base.iterate(S::GNSD, ::Val{:μ}) = (S.μ, Val(:done))
Base.iterate(S::GNSD, ::Val{:done}) = nothing

function Base.getproperty(F::GNSD, d::Symbol)
    if d == :ν
        return length(F.μ)
    else
        return getfield(F, d)
    end
end

end # module


text:
# ---
# title: 203. Remove Linked List Elements
# id: problem203
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Linked List
# link: <https://leetcode.com/problems/remove-linked-list-elements/description/>
# hidden: true
# ---
#
# Remove all elements from a linked list of integers that have value **_val_**.
#
# **Example:**
#
#
#
#     Input:  1->2->6->3->4->5->6, _**val**_ = 6
#     Output: 1->2->3->4->5
#
#
#
## @lc code=start
using LeetCode

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