text:
using MsgPack
using MsgPack: pack, unpack, Ext
using Test
using Random

ck_pack(a, b) = pack(a) == b && unpack(b) == a
round_trip(x) = unpack(pack(x))

# nil
@test ck_pack(nothing, [0xc0])

# bool
@test ck_pack(false, [0xc2])
@test ck_pack(true, [0xc3])

# positive fixint
@test ck_pack(5, [0x05])
# UInt8
@test ck_pack(128, [0xcc,0x80])
# UInt16
@test ck_pack(2^8, [0xcd,0x01,0x00])
# UInt32
@test ck_pack(2^16, [0xce,0x00,0x01,0x00,0x00])
# UInt64
@test ck_pack(2^32, [0xcf,0x00,0x00,0x00,0x01,0x00,0x00,0x00,0x00])
@test ck_pack((UInt64(2)^64-1), [0xcf,0xff,0xff,0xff,0xff,0xff,0xff,0xff,0xff])

# negative fixint
@test ck_pack(-5, [0xfb])
@test ck_pack(-32, [0xe0])
# Int8
@test ck_pack(-33, [0xd0,0xdf])
@test ck_pack(-2^7, [0xd0,0x80])
# Int16
@test ck_pack(-2^8, [0xd1,0xff,0x00])
@test ck_pack(-2^15, [0xd1,0x80,0x00])
# Int32
@test ck_pack(-2^16, [0xd2,0xff,0xff,0x00,0x00])
@test ck_pack(-2^31, [0xd2,0x80,0x00,0x00,0x00])
# Int64
@test ck_pack(-2^32, [0xd3,0xff,0xff,0xff,0xff,0x00,0x00,0x00,0x00])
@test ck_pack(-2^63, [0xd3,0x80,0x00,0x00,0x00,0x00,0x00,0x00,0x00])

# Float32
@test ck_pack((Float32(5.876)), [0xca, 0x40, 0xbc, 0x08, 0x31])
# Float64
@test ck_pack(0.654789321, [0xcb, 0x3f, 0xe4, 0xf4, 0x08, 0xbb, 0xee, 0xe1, 0xa8])

# fixstr
@test ck_pack("Меня зовут Шон",
              [0xba,0xd0,0x9c,0xd0,0xb5,0xd0,0xbd,0xd1,0x8f,
               0x20,0xd0,0xb7,0xd0,0xbe,0xd0,0xb2,0xd1,0x83,
               0xd1,0x82,0x20,0xd0,0xa8,0xd0,0xbe,0xd0,0xbd])
# str8 - currently unimplemented
# str16
s = randstring(2^8)
@test ck_pack(s, vcat(0xda, 0x01, 0x00, codeunits(s)))
s = randstring(2^16 - 1)
@test ck_pack(s, vcat(0xda, 0xff, 0xff, codeunits(s)))
# str32
s = randstring(2^16)
@test ck_pack(s, vcat(0xdb, 0x00, 0x01, 0x00, 0x00, codeunits(s)))
s = ""

# bin8
b = rand(UInt8, 121)
@test ck_pack(b, vcat(0xc4, 0x79, b))
# bin16
b = rand(UInt8, 2^11)
@test ck_pack(b, vcat(0xc5, 0x08, 0x00, b))
# bin32
b = rand(UInt8, 2^20)
@test ck_pack(b, vcat(0xc6, 0x00, 0x10, 0x00, 0x00, b))
b = []

# fixarray
@test ck_pack(Any[1, (Float32(2)), UInt8[0x22, 0x54]],
              [0x93, 0x01, 0xca, 0x40, 0x00, 0x00, 0x00, 0xc4, 0x02, 0x22, 0x54])
@test ck_pack(Any[nothing, "heya", true, false],
              [0x94, 0xc0, 0xa4, 0x68, 0x65, 0x79, 0x61, 0xc3, 0xc2])
@test ck_pack(Any[true, UInt8[0xff, 0xde, 0x11], [1, 2, 3]],
              [0x93, 0xc3, 0xc4, 0x03, 0xff, 0xde, 0x11, 0x93, 0x01, 0x02, 0x03])

# tuple
@test pack(Any[(1, [true, false], UInt8[0xff, 0xde])]) ==
              [0x91, 0x93, 0x01, 0x92, 0xc3, 0xc2, 0xc4, 0x02, 0xff, 0xde]
@test pack(("hi", ((Float32(2)), 0xff))) ==
              [0x92, 0xa2, 0x68, 0x69, 0x92, 0xca, 0x40, 0x00, 0x00, 0x00, 0xcc, 0xff]

# fixmap
@test ck_pack((Dict(1=>2, "hi"=>"mom")),
              [0x82,0xa2,0x68,0x69,0xa3,0x6d,0x6f,0x6d,0x01,0x02])

# Ext

# fixext 1
@test ck_pack(Ext(22, [0x61]), [0xd4, 0x16, 0x61])
# fixext 2
@test ck_pack(Ext(-32, [0x61, 0x61], impltype=true), [0xd5, 0xe0, 0x61, 0x61])
# fixext 4
@test ck_pack(
    Ext(22, [0x61, 0x61, 0x61, 0x61], impltype=true),
    [0xd6, 0x16, 0x61, 0x61, 0x61, 0x61]
)
# fixext 8
@test ck_pack(
    Ext(22, [0x61, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61]),
    [0xd7, 0x16, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61]
)
# fixext 16
@test ck_pack(
    Ext(22, [
        0x61, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61,
        0x61, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61
    ]),
    [
        0xd8, 0x16,
        0x61, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61,
        0x61, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61
    ]
)

# no elements
@test ck_pack(Ext(-32, UInt8[], impltype=true), [0xc7, 0x00, 0xe0])
# ext 8
@test ck_pack(
    Ext(22, [0x61, 0x61, 0x61, 0x61, 0x61, 0x61]),
    [0xc7, 0x06, 0x16, 0x61, 0x61, 0x61, 0x61, 0x61, 0x61]
)
# ext 16
b = rand(UInt8, 2^14)
@test ck_pack(Ext(78, b), vcat(0xc8, 0x40, 0x00, 0x4e, b))
# ext 32
b = rand(UInt8, 2^19)
@test ck_pack(
    Ext(-123, b, impltype=true),
    vcat(0xc9, 0x00, 0x08, 0x00, 0x00, 0x85, b)
)

# Custom type

struct A
    a::Int
    b::String
end

Base.:(==)(x::A, y::A) = x.a == y.a && x.b == y.b

MsgPack.register(A, 4)

struct B{T}
    a::T
end

Base.:(==)(x::B{T}, y::B{T}) where T = x.a == y.a

MsgPack.register(B, 5)

# Test generic encoding.

let
    x = [A(2, "hi"), A(3, "you")]
    @test round_trip(x) == x
end

let
    x = Dict(1 => A(3, "you"), 2 => A(2, "hi"))
    @test round_trip(x) == x
end

let
    x = [B(2), B(3)]
    @test round_trip(x) == x
end

let
    x = Dict(1 => B("you"), 2 => B("hi"))
    @test round_trip(x) == x
end

# Test custom encoding.

function MsgPack.encode(x::A)::Vector{UInt8}
    println("Use custom encoder for A")
    tmp = x.a, x.b
    MsgPack.pack(tmp)
end

function MsgPack.decode(::Type{A}, x::Vector{UInt8})::A
    println("Use custom decoder for A")
    a, b = MsgPack.unpack(x)
    A(a, b)
end

let
    x = [A(2, "hi"), A(3, "you")]
    @test round_trip(x) == x
end

let
    x = Dict(1 => A(3, "you"), 2 => A(2, "hi"))
    @test round_trip(x) == x
end

function MsgPack.encode(x::B)::Vector{UInt8}
    println("Use custom encoder for B")
    MsgPack.pack(x.a)
end

function MsgPack.decode(::Type{B}, x::Vector{UInt8})::B
    println("Use custom decoder for B")
    a = MsgPack.unpack(x)
    B(a)
end

let
    x = [B(2), B(3)]
    @test round_trip(x) == x
end

let
    x = Dict(1 => B("you"), 2 => B("hi"))
    @test round_trip(x) == x
end


text:
# Knet array to Float32 array conversion
karr2arr(ws) = map(w -> convert(Array{Float32}, w), ws)

function extract_file_from_tar(tar, from, to)
    sc = count(x->x=='/', from[1:end-1])
    run(`tar -xf $tar --strip-components=$sc -C $to --directory $from`)
end

function extract_file_from_zip(zip, from, to)
    filecontent = readstring(`unzip -p $zip $from`)
    write(to, filecontent)
end

# FIXME: pretty elapsed time conversion
function pretty_time(t)
    t = round(t)
    h = Int(div(t, 3600))
    m = Int(div(mod(t, 3600), 60))
    s = Int(mod(mod(t, 3600), 60))

    retval = ""
    h != 0 && (retval = string(retval, h, " hours "))
    m != 0 && (retval = string(retval, m, " minutes "))
    s != 0 && (retval = string(retval, s, " seconds"))

    return retval
end

function convert_weight(atype, w::Dict)
    w = convert(Dict{Any,Any}, w)
    for k in keys(w)
        w[k] = convert_weight(atype, w[k])
    end
    return w
end

function convert_weight(atype, w::Array)
    map(i->convert_weight(atype,w[i]), [1:length(w)...])
end

function convert_weight{T<:Number}(atype, w::Union{KnetArray{T},Array{T}})
    return convert(atype, w)
end

function get_weights(o)
    if o[:loadfile] == nothing
        w, srnn = initweights(o)
    else
        w = load(o[:loadfile], "w")
        w = convert_weight(o[:atype], w)
        srnn = nothing
    end
    return w, srnn
end

function get_wcnn(o)
    if o[:finetune] && o[:cnnfile]
        return get_vgg_weights(o[:cnnfile]; last_layer=o[:lastlayer])
    elseif o[:finetune] && !o[:cnnfile]
        return map(i -> convert(o[:atype], i), load(o[:loadfile], "wcnn"))
    else
        return nothing
    end
end

function get_opts(o,w)
    if o[:loadfile] == nothing || o[:newoptimizer]
        return init_opts(o,w)
    elseif o[:loadfile] != nothing
        opts = load(o[:loadfile], "opts")
        return copy_opts(opts,w,false)
    end
end

function copy_opts(opt::Knet.Sgd, w, save)
    Sgd(;lr=opt.lr,gclip=opt.gclip)
end

function copy_opts(opt::Knet.Adam, w, save)
    optcopy = Adam()

    # scalar elements
    optcopy.lr = opt.lr
    optcopy.beta1 = opt.beta1
    optcopy.beta2 = opt.beta2
    optcopy.t = opt.t
    optcopy.eps = opt.eps
    optcopy.gclip = opt.gclip

    # array elements
    if save
        optcopy.fstm  = Array(opt.fstm)
        optcopy.scndm = Array(opt.scndm)
    else
        optcopy.fstm = convert(typeof(w), opt.fstm)
        optcopy.scndm = convert(typeof(w), opt.scndm)
    end

    return optcopy
end

function copy_opts(opt::Knet.Adagrad, w, save)
    optcopy = Adagrad()
    optcopy.lr = opt.lr
    optcopy.eps = opt.eps
    optcopy.gclip = opt.gclip
    optcopy.G = opt.G
    return optcopy
end

function copy_opts(opts::Dict, w::Dict, save)
    optcopy = Dict()
    for k in keys(w)
        optcopy[k] = copy_opts(opts[k], w[k], save)
    end
    return optcopy
end

function copy_opts(opts::Array, w::Array, save)
    map(i -> copy_opts(opts[i], w[i], save), [1:length(w)...])
end

function init_opts(o, w::Dict)
    opts = Dict()
    for (k,v) in w
        opts[k] = init_opts(o,v)
    end
    return opts
end

function init_opts(o, w::Array)
    map(i->init_opts(o,w[i]), [1:length(w)...])
end

function init_opts{T<:Number}(o, w::Union{KnetArray{T},Array{T}})
    (eval(parse(o[:optim])))(;lr=o[:lr],gclip=o[:gclip])
end

function copy_weights(w::Dict)
    Dict(k => copy_weights(v) for (k,v) in w)
end

function copy_weights(w::Array)
    map(i->copy_weights(w[i]), [1:length(w)...])
end

function copy_weights{T<:Number}(w::Union{KnetArray{T}, Array{T}})
    Array(w)
end

ext: jl

max_stars_repo_name: ilkerkesen/kusrp

max_issues_repo_name: ilkerkesen/kusrp

max_forks_repo_name: ilkerkesen/kusrp


text:
using HDF5
using Test

@testset "properties" begin

fn = tempname()
h5open(fn, "w";
       userblock = 1024,
       alignment = (0, sizeof(Int)),
       libver_bounds = (:earliest, :latest),
      ) do hfile
    # generic
    g = create_group(hfile, "group")
    d = create_dataset(g, "dataset", datatype(Int), dataspace((500,500)),
                 alloc_time = HDF5.API.H5D_ALLOC_TIME_EARLY,
                 chunk = (5, 10),
                 obj_track_times = false)
    attributes(d)["metadata"] = "test"

    flush(hfile)

    fcpl = HDF5.get_create_properties(hfile)
    fapl = HDF5.get_access_properties(hfile)
    gcpl = HDF5.get_create_properties(hfile["group"])
    dcpl = HDF5.get_create_properties(hfile["group/dataset"])
    acpl = HDF5.get_create_properties(attributes(hfile["group/dataset"])["metadata"])

    # Retrievability of properties
    @test isvalid(fcpl)
    @test isvalid(fapl)
    @test isvalid(gcpl)
    @test isvalid(dcpl)
    @test isvalid(acpl)

    # Retrieving property values:
    @test fcpl.userblock == 1024
    @test fcpl.obj_track_times

    @test fapl.alignment == (0, sizeof(Int))
    @test fapl.driver == Drivers.POSIX()
    @test_throws HDF5.API.H5Error fapl.driver_info
    @test fapl.fclose_degree == :strong
    @test fapl.libver_bounds == (:earliest, Base.thisminor(HDF5.libversion))

    @test gcpl.local_heap_size_hint == 0
    @test gcpl.obj_track_times

    @test HDF5.UTF8_LINK_PROPERTIES.char_encoding == :utf8
    @test HDF5.UTF8_LINK_PROPERTIES.create_intermediate_group

    @test dcpl.alloc_time == :early
    @test dcpl.chunk == (5, 10)
    @test dcpl.layout == :chunked
    @test !dcpl.obj_track_times

    @test acpl.char_encoding == :utf8

    nothing
end

rm(fn, force=true)

end


text:
using OSQP
using SparseArrays
using LinearAlgebra

#Convert model from continuous to discrete time
function c2d(A, B, Δt)
    (nx, nu) = size(B)
    M = exp([Array(A) * Δt  Array(B) * Δt;
    zeros(nu, nx + nu)])
    Ad = M[1:nx, 1:nx]
    Bd = M[1:nx, nx + 1:nx + nu]
   	return  sparse(Ad), sparse(Bd)
end

const M = 10;
const θ = π / 4;
const d = 0.5;
A = sparse([
    [0 0 0 1. 0 0];
    [0 0 0 0 1. 0];
    [0 0 0 0 0 1.];
    [0 0 0 0 0 0];
    [0 0 0 0 0 0];
    [0 0 0 0 0 0]
])
B = sparse([
    [0 0 0 0];
    [0 0 0 0];
    [0 0 0 0];
    [(sqrt(2) / M) (-sqrt(2) / M) (-sqrt(2) / M) (sqrt(2) / M)];
    [(sqrt(2) / M)  (sqrt(2) / M)  (sqrt(2) / M) (sqrt(2) / M)];
    [(cos(π / 4 - θ) * d)  (cos(π / 4 - θ) * d) (cos(π / 4 - θ) * d) (cos(π / 4 - θ) * d)]
])
Ad, Bd = c2d(A, B, 0.05);
(nx, nu) = size(Bd)

# Constraints
u0 = 10.5916
umin = [9.6, 9.6, 9.6, 9.6] .- u0
umax = [13., 13., 13., 13.] .- u0
xmin = [-π / 6, -pi / 6, -Inf,  - Inf,  - Inf,  -1]
xmax = [π / 6,  π / 6,	Inf,  Inf,  Inf,  Inf]

# Objective function
Q = spdiagm(0 => [0., 0., 10., 10., 10., 10.])
QN = Q
R =  sparse(0.1I, 4, 4)

# Initial and reference states
x0 = spzeros(6, 1)
xr = [0.,0.,1.,0.,0.,0.]

# Prediction horizon
N = 10
# Cast MPC problem to a QP: x = (x(0),x(1),...,x(N),u(0),...,u(N-1))
# - quadratic objective
P = blockdiag(kron(sparse(1.0I, N, N), Q), QN,
                       kron(sparse(1.0I, N, N), R))
# - linear objective
q = vcat(kron(ones(N, 1), (-Q * xr)), (-QN * xr), spzeros(N * nu, 1))
# - linear dynamics
Ax = kron(sparse(1.0I, N + 1, N + 1), -sparse(1.0I, nx, nx)) + kron(spdiagm(-1 => ones(N + 1 - 1)), Ad)
Bu = kron(vcat(spzeros(1, N), sparse(1.0I, N, N)), Bd)
Aeq = hcat(Ax, Bu)
leq = hcat(-x0', spzeros(1, N * nx))
ueq = leq
# - input and state constraints
Aineq = spdiagm(0 => ones((N + 1) * nx + N * nu))
lineq = hcat(kron(ones(1, N + 1), xmin'), kron(ones(1, N), umin'))
uineq = hcat(kron(ones(1, N + 1), xmax'), kron(ones(1, N), umax'))
# - OSQP constraints
A = vcat(Aeq, Aineq)
l = hcat(leq, lineq)
u = hcat(ueq, uineq)

# Create an OSQP object

prob = OSQP.Model()
# Setup workspace
q = reshape(Array(q), length(q))
l = reshape(Array(l), length(l))
u = reshape(Array(u), length(u))
OSQP.setup!(prob;P = P, q = q, A = A, l = l, u = u)

results = OSQP.solve!(prob)
# Simulate in closed loop
nsim = 5000
for i = 1:nsim
    results = OSQP.solve!(prob)    # Check solver status

    # Apply first control input to the plant
    ctrl = results.x[(N + 1) * nx + 1:(N + 1) * nx + nu]
    global x0 = Ad * x0 + Bd * ctrl

    # Update initial state
    global l[1:nx] = -x0
    global u[1:nx] = -x0
    OSQP.update!(prob, l = l, u = u)
end


text:
using Documenter, Rasters, Plots, Logging, Statistics, Dates

using Rasters.LookupArrays, Rasters.Dimensions

ENV["GKSwstype"] = "100"

# Plots warnings are brWarn doctests. They dont warn the second time.
# Downloads also show op in doctests. So download everything first.
function flush_info_and_warnings()
    # RasterStack(AWAP, (:tmin, :tmax); date=DateTime(2001, 1, 1))
    # RasterStack(EarthEnv{HabitatHeterogeneity}, (:evenness, :range, :contrast, :correlation))
    # RasterStack(WorldClim{BioClim}, (1, 3, 5, 7, 12))
    st = RasterStack(WorldClim{Climate}; month=1);
    plot(st)
end
flush_info_and_warnings()

Logging.disable_logging(Logging.Warn)

# Make the docs, without running the tests again
makedocs(
    modules = [Rasters],
    sitename = "Rasters.jl",
    strict = true,
    clean = false,
)

# Enable logging to console again
Logging.disable_logging(Logging.BelowMinLevel)

deploydocs(
    repo = "github.com/rafaqz/Rasters.jl.git",
)


text:
#
# Brent.jl --
#
# Find a local root or a local minimum of an univariate function by Brent's
# methods described in:
#
# [1] Richard Brent, "Algorithms for minimization without derivatives,"
#     Prentice-Hall, inc. (1973).
#
# Global minimizer `fmin_global` is described in:
#
# [2] Ferréol Soulez, Éric Thiébaut, Michel Tallon, Isabelle Tallon-Bosc
#     and Paulo Garcia, "Optimal a posteriori fringe tracking in optical
#     interferometry", Proc. SPIE 9146, Optical and Infrared Interferometry
#     IV, 91462Y (July 24, 2014); doi:10.1117/12.2056590
#
#-----------------------------------------------------------------------------
#
# This file is part of OptimPack.jl which is licensed under the MIT
# "Expat" License:
#
# Copyright (C) 1973, Richard Brent.
# Copyright (C) 2015-2019, Éric Thiébaut.
#
#-----------------------------------------------------------------------------

# FIXME: make high level drivers fmin1, fmin2, ... which
#        check their arguments and set default tolerances
#        while low level functions do not do that?

# The following functions compute the default and minimal tolerances
# for the tolearances of the Brent's methods.
fmin_min_rtol(::Type{T}) where {T<:AbstractFloat} = convert(T,2*eps(T))
fmin_def_rtol(::Type{T}) where {T<:AbstractFloat} = convert(T,sqrt(eps(T)))
fmin_min_atol(::Type{T}) where {T<:AbstractFloat} = zero(T)
fmin_def_atol(::Type{T}) where {T<:AbstractFloat} = convert(T,3*floatmin(T))
fzero_min_rtol(::Type{T}) where {T<:AbstractFloat} = zero(T)
fzero_def_rtol(::Type{T}) where {T<:AbstractFloat} = convert(T,4*eps(T))
fzero_min_atol(::Type{T}) where {T<:AbstractFloat} = zero(T)
fzero_def_atol(::Type{T}) where {T<:AbstractFloat} = convert(T,2*floatmin(T))

# Declare specialized versions of these functions for the common
# floating-point types and manage to make them return a precomputed
# constant value.
# `goldstp` is the square of the inverse of the golden ratio.
let goldstp = (3.0 - sqrt(5.0))/2.0
    for T in (Float16, Float32, Float64)
        for func in (:fmin_min_rtol, :fmin_min_atol,
                     :fmin_def_rtol, :fmin_def_atol,
                     :fzero_min_rtol, :fzero_min_atol,
                     :fzero_def_rtol, :fzero_def_atol)
            @eval ($func)(::Type{$T}) = $(eval(func)(T))
        end
        @eval fmin_goldstp(::Type{$T}) = $(convert(T, goldstp))
    end
end

fmin_get_atol(T, x) = (x == nothing ? fmin_def_atol(T) : convert(T, x))
fmin_get_rtol(T, x) = (x == nothing ? fmin_def_rtol(T) : convert(T, x))
fzero_get_atol(T, x, a, b) = convert(T, (x == nothing ? eps(T)*abs(a - b) : x))
fzero_get_rtol(T, x) = (x == nothing ? fzero_def_rtol(T) : convert(T, x))

# _flt(T,x) method is to fast convert x to a given floating point type T
_flt(::Type{T}, x::T) where {T<:AbstractFloat} = x
_flt(::Type{T}, x::Real) where {T<:AbstractFloat} = T(x)

_flt(::Type{T}, x::AbstractRange{T}) where {T<:AbstractFloat} = x
_flt(::Type{T}, x::AbstractRange{<:Real}) where {T<:AbstractFloat} =
    convert(StepRangeLen{T}, x)

_flt(::Type{T}, x::AbstractVector{T}) where {T<:AbstractFloat} = x
_flt(::Type{T}, x::AbstractVector{<:Real}) where {T<:AbstractFloat} =
    convert(AbstractVector{T}, x)

#
# `fzero()` seeks a local root of a function F(X) in an interval [A,B].
#
# It is assumed that F(A) and F(B) have opposite signs.  This is checked,
# and an error is raised if this is not satisfied.  FZERO returns a zero X
# in the given interval [A,B] to within a tolerance: RTOL*abs(X) + ATOL.
#
# This function has been derived from Richard Brent's F77 code ZEROIN which
# itself is a slightly modified translation of the Algol 60 procedure ZERO
# given in:
#
# [1] Richard Brent, "Algorithms for minimization without derivatives,"
#     Prentice-Hall, inc. (1973).
#
# Arguments:
#
#   f    - The user-supplied function whose local root is being sought.
#          Called as F(X) to evaluates the function at any X in the
#          interval [A,B].
#
#   a, b - The endpoints of the initial search interval.
#
# Keywords:
#
#   atol - The absolute tolerance for the solution.
#
#   rtol - The relative tolerance for the solution.  The recommended (and
#          default) value for RTOL is 4*EPSILON where EPSILON is the
#          relative machine precision defined as the smallest representable
#          number such that 1 + EPSILON > 1.
#
#   T    - The floating-point data type to use for computations.  By
#          default, T=Float64.  The function f(x) must return a result
#          of type T or which can be automatically promoted to type T.
#
# Result:
#
#    The returned value is a tuple of 2 values:
#
#    (X, FX) - Where X is the estimated value of an abscissa for which
#           F is approximately zero in [A,B]; FX is the function value
#           at X.
#
function fzero(f::Function, a::Real, b::Real;
               T::Type=Float64, atol=nothing, rtol=nothing)
    # Enforce floating point type and get tolerances.
    @assert T <: AbstractFloat
    _a = T(a)
    _b = T(b)
    _atol = fzero_get_atol(T, atol, _a, _b)
    _rtol = fzero_get_rtol(T, rtol)
    return fzero(f, _a, _b, _atol, _rtol)
end

function fzero(f::Function, a::Real, b::Real, atol::Real, rtol::Real)
    T = float(promote_type(typeof(a), typeof(b), typeof(atol), typeof(rtol)))
    return fzero(T, f, a, b, atol, rtol)
end

function fzero(::Type{T}, f::Function, a::Real, b::Real,
              atol::Real, rtol::Real) where {T<:AbstractFloat}
    return fzero(f, _flt(T,a), _flt(T,b), _flt(T,atol), _flt(T,rtol))
end

function fzero(f::Function, a::T, b::T, atol::T, rtol::T) where {T<:AbstractFloat}

    # Explicitly declare types of variables to prevent accidental changes
    # of type during the execution of the code.  We also re-assign the
    # value of the passed arguments to force conversion (FIXME: there should
    # be a macro for that).
    local fa::T, fb::T, c::T, fc::T, tol::T, m::T, e::T, d::T

    # Some constants.
    ZERO::T = zero(T)
    ONE::T = one(T)
    TWO::T = ONE + ONE
    HALF::T = ONE/TWO
    THREE::T = ONE + TWO

    # Compute the function value at the endpoints and check the
    # assumptions.
    fa = f(a)
    fa == ZERO && return (a, fa)
    fb = f(b)
    fb == ZERO && return (b, fb)
    if (fa > ZERO) == (fb > ZERO)
        error("f(a) and f(b) must have different signs")
    end
    fc = fb # to trigger bound update below
    while true
        if (fb > ZERO) == (fc > ZERO)
            # Drop point C (make it coincident with point A)
            # and adjust bounds of interval.
            c, fc = a, fa
            e = d = b - a
        end

        # Make sure B is the best point so far.
        if abs(fc) < abs(fb)
            a, fa = b, fb
            b, fb = c, fc
            c, fc = a, fa
        end

        # Compute tolerance.  In original Brent's method, the precision and
        # the computed tolerance are given by:
        #    PREC = 4*EPS*abs(X) + 2*T
        #    TOL = 2*EPS*abs(b) + T = PREC/2
        # and we want:
        #    PREC = RTOL*abs(X) + ATOL
        # thus the expression of the tolerance parameter becomes:
        tol = HALF*(rtol*abs(b) + atol)

        # Check for convergence.
        m = (c - b)*HALF
        if abs(m) <= tol || fb == ZERO
	    return (b, fb)
        end

        # See if a bisection is forced.
        if abs(e) < tol || abs(fa) <= abs(fb)
            # Bisection.
            d = e = m
        else
            local p::T
            local q::T
            local r::T
            local s::T = fb/fa
            local two_p::T
            if a == c
                # Linear interpolation.
                p = TWO*m*s
                q = ONE - s
            else
                # Inverse quadratic interpolation.
                q = fa/fc
                r = fb/fc
                p = (TWO*m*q*(q - r) - (b - a)*(r - ONE))*s
                q = (q - ONE)*(r - ONE)*(s - ONE)
            end
            if p > ZERO
                q = -q
            else
                p = -p
            end
            two_p = p + p
            if two_p < THREE*m*q - tol*abs(q) && two_p < abs(e*q)
                # Take the interpolation point.
                e = d
                d = p/q
            else
                # Force a bisection.
                d = e = m
            end
        end
        a, fa = b, fb
        if abs(d) > tol
            b += d
        elseif m > ZERO
            b += tol
        else
            b -= tol
        end
        fb = f(b)
        fb == ZERO && return (b, fb)
    end
end

# FIXME: AUTOMATIC BRACKETING
#
# If the interval to consider is not bounded or only left/right bounded,
# the idea is to find a suitable interval (A,B) where at least one
# minimum must exists (if the function is continue) and start Brent's
# algorithm with correct values for X, FX, ... (in order to save some
# function evaluations).

#
# `fmin()` seeks a local minimum of a function F(X) in an interval [A,B].
#
# The method used is a combination of golden section search and successive
# parabolic interpolation.  Convergence is never much slower than that for
# a Fibonacci search.  If F has a continuous second derivative which is
# positive at the minimum (which is not at A or B), then convergence is
# superlinear, and usually of the order of about 1.324....
#
# The values RTOL and ATOL define a tolerance TOL = RTOL*abs(X) + ATOL.  F
# is never evaluated at two points closer than TOL.
#
# If F is a unimodal function and the computed values of F are always
# unimodal when separated by at least SQEPS*abs(X) + (ATOL/3), then FMIN
# approximates the abscissa of the global minimum of F on the interval
# [A,B] with an error less than 3*SQEPS*abs(FMIN) + ATOL.
#
# If F is not unimodal, then FMIN may approximate a local, but perhaps
# non-global, minimum to the same accuracy.
#
# This function has been derived from Richard Brent's FORTRAN77 code FMIN
# which itself is a slightly modified translation of the Algol 60 procedure
# LOCALMIN given in:
#
# [1] Richard Brent, "Algorithms for minimization without derivatives,"
#     Prentice-Hall, inc. (1973).
#
# Arguments:
#
#   a, b - The endpoints of the interval.
#
#   f    - The user-supplied function whose local minimum is being
#          sought.  Called as F(X) to evaluates the function at X.
#
# Keywords:
#
#   atol - A positive absolute error tolerance.
#
#   rtol - A positive relative error tolerance.  RTOL should be no smaller
#          than twice the relative machine precision, and preferably not
#          much less than the square root of the relative machine
#          precision.
#
#   T    - The floating-point data type to use for computations.  By
#          default, T=Float64.  The function f(x) must return a result
#          of type T or which can be automatically promoted to type T.
#
# Result:
#
#   The returned value is a tuple of 4 values:
#
#   (X, FX, XLO, XHI) - X is the estimated value of an abscissa for which F
#          attains a local minimum value in [A,B]; FX is the function value
#          at X; XLO and XHI are the bounds for the position of the local
#          minimum.
#
function fmin(f::Function, a::Real, b::Real;
              T::Type=Float64, atol=nothing, rtol=nothing)
    # Enforce floating point type and get tolerances.
    @assert T <: AbstractFloat
    return fmin(f, _flt(T,a), _flt(T,b), fmin_get_atol(T, atol), fmin_get_rtol(T, rtol))
end

function fmin(f::Function, a::Real, b::Real, atol::Real, rtol::Real)
    T = float(promote_type(typeof(a), typeof(b), typeof(atol), typeof(rtol)))
    return fmin(T, f, a, b, atol, rtol)
end

function fmin(::Type{T}, f::Function, a::Real, b::Real,
              atol::Real, rtol::Real) where {T<:AbstractFloat}
    return fmin(f, _flt(T,a), _flt(T,b), _flt(T,atol), _flt(T,rtol))
end

function fmin(f::Function, a::T, b::T, atol::T, rtol::T) where {T<:AbstractFloat}

    # Enforce floating point type.
    a::T = _a
    b::T = _b
    atol::T = fmin_get_atol(T, atol)
    rtol::T = fmin_get_rtol(T, rtol)

    # Make sure A and B are properly ordered.
    if a > b
        a, b = b, a
    end

    # Initialize the search.
    c::T = fmin_goldstp(T) # square of the inverse of the golden ratio
    x::T = a + c*(b - a)
    fx::T = f(x)
    _fmin(f, a, b, x, fx, x, fx, x, fx, atol, rtol)
end

# This function is called to start Brent's algorithm with a single given
# point, x, inside the search interval [a,b] and with known function value
# fx = f(x).
function _fmin1(f::Function, a::T, b::T, x::T, fx::T,
                atol::T, rtol::T) where {T<:AbstractFloat}
    # Make sure A and B are properly ordered.
    if a > b
        a, b = b, a
    end
    if x < a || x > b
        error("given point outside search interval")
    end
    _fmin(f, a, b, x, fx, x, fx, x, fx, atol, rtol)
end

# This function is called to start Brent's algorithm with two given points,
# (x and w not necessarily distinct) inside the search interval [a,b] and
# with known function values (fx and fw).
function _fmin2(f::Function, a::T, b::T, x::T, fx::T, w::T, fw::T,
                atol::T, rtol::T) where {T<:AbstractFloat}
    # Make sure A and B are properly ordered and check that given
    # points are in the interval.
    if a > b
        a, b = b, a
    end
    if w < a || w > b || x < a || x > b
        error("given point(s) outside search interval")
    end

    # Reorder the points as assumed by Brent's algorithm.
    if fw < fx
        x, fx, w, fw = w, fw, x, fx
    end
    _fmin(f, a, b, x, fx, w, fw, w, fw, atol, rtol)
end

# This function is called to start Brent's algorithm with 3 given
# points (x, w, and v not necessarily distinct) inside the search
# interval [a,b] and with known function values (fx, fw and fv).
function _fmin3(f::Function, a::T, b::T, x::T, fx::T,
                w::T, fw::T, v::T, fv::T,
                atol::T, rtol::T) where {T<:AbstractFloat}
    # Make sure A and B are properly ordered.
    if a > b
        a, b = b, a
    end
    if v < a || v > b || w < a || w > b || x < a || x > b
        error("given point(s) outside search interval")
    end

    # Reorder the points as assumed by Brent's algorithm.
    if fw < fx
        x, fx, w, fw = w, fw, x, fx
    end
    if fv < fx
        x, fx, v, fv = v, fv, x, fx
    end
    if fv < fw
        w, fw, v, fv = v, fv, w, fw
    end
    _fmin(f, a, b, x, fx, w, fw, v, fv, atol, rtol)
end

#
# Original Brent's algorithm assumes that the minimum is in (a,b) with
# a <= b and keeps track of the following variables:
#
#   x, fx = f(x) - least function value found so far;
#   w, fw = f(w) - previous value of x, fx;
#   v, fv = f(v) - previous value of w, fw;
#   d - computed step (new try is: u = x + d, unless d too small);
#   e - the previous value of d, if a parabolic step is taken; the
#       difference between the most distant current endpoint and x,
#       if a golden step is taken.
#
# Other variables need not be saved.
#   u, fu = f(u) - next point to try and its function value.
#
# Thus the main loop of Brent's algorithm can be entered with any x, w, v
# (not necessarily distinct) which are in [a,b] and such that:
#
#    fx = f(x)  <=  fw = f(w)  <=  fv = f(v)
#
# with:
#
#    d = x - w
#    e = w - v
#
# The following function is the main loop of Brent's algorithm.  It assumes
# that all parameters are properly set (as explained above).
#
function _fmin(f::Function, a::T, b::T, x::T, fx::T,
               w::T, fw::T, v::T, fv::T,
               atol::T, rtol::T) where {T<:AbstractFloat}
    # Explicitly declare types of variables to prevent accidental changes
    # of type during the execution of the code.  We also re-assign the
    # value of the passed arguments to force conversion (FIXME: there should
    # be a macro for that).
    a::T = a; b::T = b
    x::T = x; fx::T = fx
    w::T = w; fw::T = fw
    v::T = v; fv::T = fv
    atol::T = atol; rtol::T = rtol
    ZERO::T = zero(T)
    ONE::T = one(T)
    TWO::T = ONE + ONE
    HALF::T = ONE/TWO
    c::T = fmin_goldstp(T)
    local d::T = x - w
    local e::T = w - v
    local m::T, tol::T, t2::T, u::T, fu::T
    local take_golden_step::Bool

    while true

        # Compute mid-point and check the stopping criterion.
        m = HALF*(a + b)
        tol = rtol*abs(x) + atol
        t2 = tol + tol
        if abs(x - m) <= t2 - HALF*(b - a)
            return (x, fx, a, b)
        end

        # Determine next step to take.
        take_golden_step = true
        if abs(e) > tol
            # Fit a parabola (make sure final Q >= 0).
            local p::T, q::T, r::T
            r = (x - w)*(fx - fv)
            q = (x - v)*(fx - fw)
            if q > r
                p = (x - w)*r - (x - v)*q
                q = (q - r)*TWO
            else
                p = (x - v)*q - (x - w)*r
                q = (r - q)*TWO
            end
            if abs(p) < HALF*q*abs(e) && q*(a - x) < p < q*(b - x)
                # Take the parabolic interpolation step.
                take_golden_step = false
                e = d
                d = p/q
                u = x + d
                # F must not be evaluated too close to A or B.
                if u - a < t2 || b - u < t2
                    d = (x < m ? tol : -tol)
                end
            end
        end
        if take_golden_step
            # Take a golden-section step.
            e = (x < m ? b : a) - x
            d = c*e
        end

        # F must not be evaluated too close to X.
        if abs(d) >= tol
            u = x + d
        elseif d > ZERO
            u = x + tol
        else
            u = x - tol
        end
        fu = f(u)

        # Update A, B, V, W, and X.
        if fu <= fx
            if u < x
                b = x
            else
                a = x
            end
            v, fv = w, fw
            w, fw = x, fx
            x, fx = u, fu
        else
            if u < x
                a = u
            else
                b = u
            end
            if fu <= fw || w == x
                v, fv = w, fw
	        w, fw = u, fu
            elseif fu <= fv || v == x || v == w
                v, fv = u, fu
            end
        end
    end
end

# This function is called to start Brent's algorithm with a bracket of the
# minimum defined by 3 points (x, w and v) with known function values (fx,
# fw and fv) and such that the least function value is at x which is inside
# the interval [v,w].
function _fminbrkt(f::Function, x::T, fx::T, w::T, fw::T, v::T, fv::T,
                   atol::T, rtol::T) where {T<:AbstractFloat}
    # Reorder the points as assumed by Brent's algorithm.
    if fv < fw
        v, fv, w, fw = w, fw, v, fv
    end
    if v < w
        a, b = v, w
    else
        a, b = w, v
    end
    if x < a || x > b || fx > fw
        error("illegal bracket")
    end
    _fmin(f, a, b, x, fx, w, fw, v, fv, atol, rtol)
end

#
# `fmin_global()` finds the global minimum of an univariate function F.
# The argument X is a vector of coordinates in monotonic order; X[1] and
# X[end] are the endpoints of the global search interval and the other
# values of X are such that no more than a single local minimum lies in any
# subinterval [X(i),X(i+2)].
#
# X = range(A,stop=B,length=N) if these arguments are supplied instead;
# i.e. the global search is performed in the (included) interval [A,B] which is
# cut in N pieces of equal lengths.
#
# If specified, keywords `atol` and `rtol` set the absolute and relative
# tolerances for the precision.
#
# This function implements the BraDi ("Bracket" then "Dig") algorithm
# described in [1].
#
# [1] Ferréol Soulez, Éric Thiébaut, Michel Tallon, Isabelle Tallon-Bosc
#     and Paulo Garcia, "Optimal a posteriori fringe tracking in optical
#     interferometry", Proc. SPIE 9146, Optical and Infrared Interferometry
#     IV, 91462Y (July 24, 2014); doi:10.1117/12.2056590
#
# SEE ALSO: fmin.
function fmin_global(f::Function, a::Real, b::Real, n::Integer;
                     T::Type=Float64, kwds...)
    return fmin_global(f, range(convert(T,a), stop=convert(T,b), length=n);
                       T=T, kwds...)
end

function fmin_global(::Type{T}, f::Function, x::AbstractVector{<:Real},
                     atol=Real, rtol=Real) where {T<:AbstractFloat}
    return fmin_global(f, _flt(T, x), _flt(T, atol), _flt(T, rtol))
end

function fmin_global(f::Function, x::AbstractVector{T},
                     atol::T, rtol::T) where {T<:AbstractFloat}
    atol::T = fmin_get_atol(T, atol)
    rtol::T = fmin_get_rtol(T, rtol)

    xbest = xa = xb = xc = x[1]
    fbest = fa = fb = fc = f(xc)
    n = length(x)
    for j in 2 : n + 1
        xa = xb
        fa = fb
        xb = xc
        fb = fc
        if j <= n
            xc = x[j]
            fc = f(xc)
            if fc < fbest
                xbest = xc
                fbest = fc
            end
        end
        if fa >= fb <= fc
            # A minimum has been bracketed in [XA,XC].
            xm, fm = _fminbrkt(f, xb, fb, xa, fa, xc, fc, atol, rtol)
            if fm < fbest
                xbest = xm
                fbest = fm
            end
        end
    end
    return (xbest, fbest)
end


text:
#An example of nGMCA
using CompressedSensing
using PyPlot
srand(1)

#Generate some sparse positive data, S and A
n=1000
m=500
r=20
k_a=2
k_s=1

A=zeros(n,r)
S=zeros(r,m)

for i =1:n
    A[i,int(rand(k_a)*(r-1))+1]=rand(k_a)
end
for i =1:m
    S[int(rand(k_s)*(r-1))+1,i]=rand(k_s)
end

#calculate the Gini Index of the data
print(string("GI(S): ",GI(S),"    GI(A): ",GI(A)))

#Find the response matrix
Y=A*S

#attempt to rebuild the inputs from only the output
(A0,S0)=nGMCA(Y,r;verbose=true,maxIter=100,threshold=1e-7,phaseRatio=.15)

#See how well it correlates with the inputs
figure(figsize=(12,5))
subplot(121)
title("Correlation for A")
pcolormesh(cor(A0,A).^2)
colorbar()
subplot(122)
title("Correlation for S")
pcolormesh(cor(S0',S').^2)
colorbar();
savefig("Example_4_Fig_1.png")

max_stars_repo_path: examples/Example 4.jl

max_stars_repo_name: dahlend/CompressedSensing

max_issues_repo_path: examples/Example 4.jl

max_forks_repo_path: examples/Example 4.jl


text:
function iterate!(job::BasicMCJob, ::Type{UnvAMWG}, ::Type{Univariate})
  job.sstate.tune.proposed += 1

  if mod(job.sstate.tune.proposed, job.tuner.period) == 0
    set_batch!(job.sstate.tune, job.tuner)
    set_delta!(job.sstate.tune, job.tuner)
  end

  if isinf(job.sampler.lower) && isinf(job.sampler.upper)
    job.sstate.pstate.value = job.pstate.value+sqrt(exp(job.sstate.tune.logσ))*randn()
    job.parameter.logtarget!(job.sstate.pstate)

    job.sstate.ratio = job.sstate.pstate.logtarget-job.pstate.logtarget
  else
    job.sstate.proposal = Truncated(Normal(
      job.pstate.value, sqrt(exp(job.sstate.tune.logσ))), job.sampler.lower, job.sampler.upper
    )
    job.sstate.pstate.value = rand(job.sstate.proposal)
    job.parameter.logtarget!(job.sstate.pstate)

    job.sstate.ratio = job.sstate.pstate.logtarget-job.pstate.logtarget
    job.sstate.ratio -= logpdf(job.sstate.proposal, job.sstate.pstate.value)

    job.sstate.proposal = Truncated(Normal(
      job.sstate.pstate.value, sqrt(exp(job.sstate.tune.logσ))), job.sampler.lower, job.sampler.upper
    )

    job.sstate.ratio += logpdf(job.sstate.proposal, job.pstate.value)
  end

  if (job.sstate.ratio > 0 || (job.sstate.ratio > log(rand())))
    job.pstate.value = job.sstate.pstate.value

    job.pstate.logtarget = job.sstate.pstate.logtarget

    if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
      job.pstate.loglikelihood = job.sstate.pstate.loglikelihood
    end
    if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
      job.pstate.logprior = job.sstate.pstate.logprior
    end

    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = true
      end
    end

    job.sstate.tune.accepted += 1
  else
    if !isempty(job.sstate.diagnosticindices)
      if haskey(job.sstate.diagnosticindices, :accept)
        job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]] = false
      end
    end
  end

  if haskey(job.sstate.diagnosticindices, :logσ)
    job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:logσ]] = job.sstate.tune.logσ
  end

  if mod(job.sstate.tune.proposed, job.tuner.period) == 0
    rate!(job.sstate.tune)

    if job.tuner.verbose && job.sstate.tune.totproposed <= job.range.burnin
      println(
        "Burnin iteration ",
        job.fmt_iter(job.sstate.tune.totproposed),
        " of ",
        job.range.burnin,
        ": ",
        job.fmt_perc(100*job.sstate.tune.rate),
        " % acceptance rate"
      )
    end

    tune!(job.sstate.tune, job.tuner)
    reset!(job.sstate.tune)
  end
end

function iterate!(job::BasicMCJob, ::Type{MuvAMWG}, ::Type{Multivariate})
  job.sstate.tune.proposed += 1

  if mod(job.sstate.tune.proposed, job.tuner.period) == 0
    set_batch!(job.sstate.tune, job.tuner)
    set_delta!(job.sstate.tune, job.tuner)

    if job.tuner.verbose && job.sstate.tune.totproposed <= job.range.burnin
      println("Burnin iteration ", job.fmt_iter(job.sstate.tune.totproposed), " of ", job.range.burnin, ":")
    end
  end

  for i in 1:job.pstate.size
    if isinf(job.sampler.lower[i]) && isinf(job.sampler.upper[i])
      job.sstate.pstate.value[i] = job.pstate.value[i]+sqrt(exp(job.sstate.tune.logσ[i]))*randn()
      job.parameter.logtarget!(job.sstate.pstate)

      job.sstate.ratio = job.sstate.pstate.logtarget-job.pstate.logtarget
    else
      job.sstate.proposal = Truncated(Normal(
        job.pstate.value[i], sqrt(exp(job.sstate.tune.logσ[i]))), job.sampler.lower[i], job.sampler.upper[i]
      )
      job.sstate.pstate.value[i] = rand(job.sstate.proposal)
      job.parameter.logtarget!(job.sstate.pstate)

      job.sstate.ratio = job.sstate.pstate.logtarget-job.pstate.logtarget
      job.sstate.ratio -= logpdf(job.sstate.proposal, job.sstate.pstate.value)

      job.sstate.proposal = Truncated(Normal(
        job.sstate.pstate.value[i], sqrt(exp(job.sstate.tune.logσ[i]))), job.sampler.lower[i], job.sampler.upper[i]
      )

      job.sstate.ratio += logpdf(job.sstate.proposal, job.pstate.value)
    end

    if (job.sstate.ratio > 0 || (job.sstate.ratio > log(rand())))
      job.pstate.value[i] = job.sstate.pstate.value[i]

      job.pstate.logtarget = job.sstate.pstate.logtarget

      if in(:loglikelihood, job.outopts[:monitor]) && job.parameter.loglikelihood! != nothing
        job.pstate.loglikelihood = job.sstate.pstate.loglikelihood
      end
      if in(:logprior, job.outopts[:monitor]) && job.parameter.logprior! != nothing
        job.pstate.logprior = job.sstate.pstate.logprior
      end

      if !isempty(job.sstate.diagnosticindices)
        if haskey(job.sstate.diagnosticindices, :accept)
          job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]][i] = true
        end
      end

      job.sstate.tune.accepted[i] += 1
    else
      if !isempty(job.sstate.diagnosticindices)
        if haskey(job.sstate.diagnosticindices, :accept)
          job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:accept]][i] = false
        end
      end
    end

    if haskey(job.sstate.diagnosticindices, :logσ)
      job.pstate.diagnosticvalues[job.sstate.diagnosticindices[:logσ]][i] = job.sstate.tune.logσ[i]
    end

    if mod(job.sstate.tune.proposed, job.tuner.period) == 0
      rate!(job.sstate.tune, i)

      if job.tuner.verbose && job.sstate.tune.totproposed <= job.range.burnin
        println(
          "  parameter ",
          i,
          " of ",
          job.pstate.size,
          ": ",
          job.fmt_perc(100*job.sstate.tune.rate[i]),
          " % acceptance rate"
        )
      end

      tune!(job.sstate.tune, job.tuner, i)
      reset!(job.sstate.tune, i)
    end
  end

  if mod(job.sstate.tune.proposed, job.tuner.period) == 0
    reset!(job.sstate.tune)
  end
end


text:
using Documenter, Integrals

makedocs(;
    modules=[Integrals],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
    ],
    repo="https://github.com/jw3126/Integrals.jl/blob/{commit}{path}#L{line}",
    sitename="Integrals.jl",
    authors="Jan Weidner",
    assets=[],
)

deploydocs(;
    repo="github.com/jw3126/Integrals.jl",
)


text:
function coerceH(H)
  n = length(H[1, :]) # n – число столбцов
  r = length(H[:, 1]) # r – число строк

  # проxод вниз
  for i = 1:r-1
    col = n - r + i # текущий столбец

    row = -1 # первая строка с 1 в текущем столбце
    for j = i:r
      if H[j, col] == 1
        row = j
        break
      end
    end

    findCollumn = false
    if row == -1 # вырожденная подматрица, нужно менять столбцы
      for ii = i:r
        for jj = n-r:1
          if H[ii, jj] == 1
            println("меняем столбцы $col и $jj")
            H[:, col], H[:, jj] = H[:, jj], H[:, col]
            row = ii

            findCollumn = true
            break
          end
        end

        if findCollumn break end
      end
    end

    if row != i # нужно менять строки
      H[i, :], H[row, :] = H[row, :], H[i, :] # меняем строки местами
    end

    for j = i+1:r
      if H[j, col] == 1
        H[j, :] $= H[i, :] # $ — xor
      end
    end
  end

  # проxод вверx
  for i = r:-1:2
    col = n - r + i # текущий столбец

    for j = i-1:-1:1
      if H[j, col] == 1
        H[j, :] $= H[i, :] # $ — xor
      end
    end
  end

  P = H[:, 1:n-r]'  # первые n-k столбцов, ' — транспонирование
  E = int(eye(n-r)) # создаём единичную матрицу размера n-k
  G = [E P]

  return (G, H)
end

function minDistance(G)
  k = length(G[:, 1]) # число строк
  minWeight = 1000000 # какое-то большое число

  for i = 1:2^k-1 # перебираем все вектора кроме 0
    b = collect(reverse(bits(i))[1:k]) # k символов бит
    m = map(parseint, b) # bits возвращает 64 символа, конвертируем в int

    d = m' * G
    d = map(x -> x % 2, d) # берём по модулю 2

    w = sum(d)
    minWeight = min(w, minWeight)
  end

  return minWeight
end

function syndrome(H)
  n = length(H[1, :]) # n – число столбцов
  r = length(H[:, 1]) # r – число строк

  syndromeWeights = fill!(Array(Int64, 2^r), n + 1) # массив весов, заполняем значением n + 1
  S = Array(Int64, 2^r - 1, n)

  for i = 1:2^n-1 # перебираем все вектора кроме 0
    b = collect(reverse(bits(i))[1:n]) # n символов бит
    e = map(parseint, b) # текущий вектор ошибки

    s = e' * H'
    s = map(x -> x % 2, s) # берём по модулю 2
    w = sum(e)
    sind = parseint(join(s, ""), 2) # конвертируем битовый вектор в число

    if w != 0 && sind != 0 && syndromeWeights[sind] > w
      S[sind, :] = e
      syndromeWeights[sind] = w
    end
  end

  return S
end

function toTexTable(M)
  rows    = [M[i, :] for i in 1:length(M[:, 1])]
  texRows = [join(row, " & ") for row in rows]
  endRows = join(texRows, " \\\\ \n")
end

function toSyndromeTable(S)
  tableRows = ["\$0000\$ & \$0000000000\$", ["\$$(bits(i)[61:64])\$ & \$$(join(S[i, :], ""))\$" for i in 1:length(S[:, 1])]]

  join(tableRows, " \\\\ \n")
end

HS = [[1 0 0 1],
     [0 1 0 1],
     [0 0 1 1]]
HK = [[1 0 1 0 1 0 0 1 0 1],
      [0 1 0 0 0 0 1 1 0 1],
      [1 1 1 1 0 0 0 0 1 0],
      [0 1 1 1 0 1 0 1 0 0]]

HTEST = [[1 1 0 0 1 1],
         [0 0 1 0 1 1],
         [0 1 0 1 0 1]]

println(toTexTable(HK))

#S = syndrome(HK)
#println(S)
#println("S = ")
#println(toTexTable(S))
#println(toSyndromeTable(S))

G, H = coerceH(HK)
#w = minDistance(G)

#println("H = ")
#println(H)

println("G = ")
println(G)

#println("H = ")
#println(H)
println(toTexTable(G))

#println("w = $w")

max_stars_repo_name: ChShersh/university-courses

max_issues_repo_name: ChShersh/university-courses

max_forks_repo_name: ChShersh/university-courses


text:
## consumer

TopicPartition = Tuple{String, Int}

mutable struct KafkaConsumer
    client::KafkaClient
    rkparlist::Ptr{Cvoid}
end

function Base.setproperty!(c::KafkaConsumer, f::Symbol, v)
    f == :rkparlist || error("KafkaConsumer.$f unassignable")
    c.rkparlist == C_NULL || kafka_topic_partition_list_destroy(c.rkparlist)

    v isa Ptr{Cvoid} && return setfield!(c, f, v)

    rkparlist = kafka_topic_partition_list_new(length(v))
    for (t, p) in v
        kafka_topic_partition_list_add(rkparlist, t, p)
    end
    setfield!(c, f, rkparlist)
end

function KafkaConsumer(conf::Dict)
    @assert haskey(conf, "bootstrap.servers") "`bootstrap.servers` should be specified in conf"
    @assert haskey(conf, "group.id") "`group.id` should be specified in conf"
    client = KafkaClient(KAFKA_TYPE_CONSUMER, conf)
    consumer = KafkaConsumer(client, C_NULL)
    finalizer(c -> c.rkparlist = C_NULL, consumer)
    return consumer
end

function KafkaConsumer(bootstrap_servers::String, group_id::String, conf::Dict=Dict())
    conf["bootstrap.servers"] = bootstrap_servers
    conf["group.id"] = group_id
    return KafkaConsumer(conf)
end

function Base.show(io::IO, c::KafkaConsumer)
    group_id = c.client.conf["group.id"]
    bootstrap_servers = c.client.conf["bootstrap.servers"]
    print(io, "KafkaConsumer($group_id @ $bootstrap_servers)")
end

function subscribe(c::KafkaConsumer, tpars::Vector{TopicPartition})
    c.rkparlist = tpars
    kafka_subscribe(c.client.rk, c.rkparlist)
end

function assign(c::KafkaConsumer, tpars::Vector{TopicPartition})
    c.rkparlist = tpars
    kafka_assign(c.client.rk, c.rkparlist)
end

function poll(::Type{K}, ::Type{P}, c::KafkaConsumer, timeout::Int=1000) where {K,P}
    c_msg_ptr = kafka_consumer_poll(c.client.rk, timeout)
    if c_msg_ptr != nothing
        c_msg = unsafe_load(c_msg_ptr)
        msg = Message{K,P}(c_msg)
        kafka_message_destroy(c_msg_ptr)
        return msg
    else
        return nothing
    end
end

poll(c::KafkaConsumer, timeout::Int=1000) = poll(Vector{UInt8}, Vector{UInt8}, c, timeout)

function Base.seek(c::KafkaConsumer,
                   topic_partition::TopicPartition,
                   offset::Integer,
                   timeout::Integer=1000)
    topic, par = topic_partition
    kafka_topic_partition_list_find(c.rkparlist, topic, par) ||
        error("Seek on an assigned/subscribed topic partition $topic_partition")

    kt = KafkaTopic(c.client, topic)
    kafka_seek(kt.rkt, par, offset, timeout)
end


text:
#---------------------------------------------------------------------#
#This routine advances the solution in time using Split-Explicit using
#a General Order RK method in Butcher tableau form following Schlegel et al. GMD 2012.
#Written by F.X. Giraldo / P.R. Mugg on 9/22/19
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey, CA 93943-5216
#---------------------------------------------------------------------#
function multirate_erk_butcher_substepping(q0,c,αI,βI,I,αJ,βJ,J,M,Δt,time)

    #Compute Coefficients
    #Outer Loop
    α̃I = zeros(I+1,I)
    c̃I = zeros(I+1,1)
    cI = zeros(I,1)
    for i = 2:I
        cI[i] = sum(αI[i,:])
        α̃I[i,:] = αI[i,:] - αI[i-1,:]
        c̃I[i] = cI[i] - cI[i-1]
    end
    for j = 1:I
        α̃I[I+1,j] = βI[j] - αI[I,j]
    end
    c̃I[I+1] = 1.0 - cI[I]

    #Inner Loop
    α̃J = zeros(J+1,J)
    c̃J = zeros(J+1,1)
    cJ = zeros(J,1)
    for i = 2:J
        cJ[i] = sum(αJ[i,:])
        α̃J[i,:] = αJ[i,:] - αJ[i-1,:]
        c̃J[i] = cJ[i] - cJ[i-1]
    end
    for j = 1:J
        α̃J[J+1,j] = βJ[j] - αJ[J,j]
    end
    c̃J[J+1] = 1.0 - cJ[J]

    #Initialize Arrays
    Q = zeros(I+1,1)
    V = zeros(I+1,J+1)

    #Outer Loop
    Q[1] = q0
    for i = 2:I+1
        r_i = 0
        for j = 1:i-1
            t = time + Δt*cI[j]
            (rhs,f𝑗,g𝑗) = rhs_function(Q[j],c,t)
            r_i += α̃I[i,j]*g𝑗
        end

        #Substepping Loop
        qm = Q[i-1]
        for m = 1:M
            V[i,1] = qm

            #Inner Loop
            for j = 2:J+1
                R_sum = 0
                for k = 1:j-1
                    t = time + Δt*( cI[i-1] + c̃I[i]*cJ[k]/M*(m-1) ) #not correct so am losing accuracy
                    (rhs,f𝑘,g𝑘) = rhs_function(V[i,k],c,t)
                    R_sum += α̃J[j,k]*( r_i + c̃I[i]*f𝑘 )
                end #k
                V[i,j] = V[i,j-1] + Δt/M*R_sum
            end #j
            qm = V[i,J+1]
        end #m
        Q[i] = qm
    end #i

    return qp = Q[I+1]
end #end function


text:
## Be careful, some of this functionality depends on the storage type provided.
## Use higher level function calls for Whittle and de-biased Whitle if possible!

# Low level function to compute the Whittle summand given a spectral quantity S (the sdf or EI) and periodogram ordinate I.
generalwhittle_summand(S::SHermitianCompact, I) = real(log(det(S)) + tr(I*inv(S)))
generalwhittle_summand(S::Matrix{T}, I) where {T} = real(log(det(S)) + tr(I/S))
generalwhittle_summand(S::T, I) where {T<:Number} = real(log(S) + I/S)

# Compute the Whittle likelihood from the Periodogram and some spectral quantity (sdf or EI).
generalwhittle(S, I::Vector{T}) where {T} = sum(generalwhittle_summand(Sᵢ,Iᵢ) for (Sᵢ,Iᵢ) ∈ zip(S,I))
generalwhittle(funcmemory, data::GenWhittleData) = @views generalwhittle(unpack(funcmemory)[data.Ω_used_index], data.I)
generalwhittle(store::TimeSeriesModelStorageFunction, data::GenWhittleData) = generalwhittle(store.funcmemory, data)
generalwhittle(store::AdditiveStorage, data::GenWhittleData) = generalwhittle(store.store1, data) # function value will always be stored in leftmost storage.

## gradient
# Compute the gradient of the Whittle likelihood from the Periodogram and some spectral quantity (sdf or EI).
function grad_generalwhittle!(∇ℓ, S, ∇S, I::Vector{T}) where {T<:Number} # univariate case
    length(S) == length(I) == size(∇S,2) || throw(ArgumentError("S and I should be same length as second dim of ∇S."))
    ∇ℓ .= zero(eltype(∇ℓ))
    @inbounds for iω = 1:length(S)
        Ipart = (1.0/real(S[iω]) - I[iω]/(real(S[iω])^2))
        for jpar = 1:size(∇S, 1)
            ∇ℓ[jpar] += Ipart * real(∇S[jpar, iω])
        end
    end
    return nothing
end
function grad_generalwhittle!(∇ℓ, S, ∇S, I::Vector{SMatrix{D,D,T,n}}) where {D,T,n} # multivariate case
    length(S) == length(I) == size(∇S,2) || throw(ArgumentError("S and I should be same length as second dim of ∇S."))
    ∇ℓ .= zero(eltype(∇ℓ))
    Id = diagm(@SVector ones(D))
    @inbounds for iω = 1:length(S)
        invS = inv(S[iω])
        Ipart = (Id - I[iω] * invS)
        for jpar = 1:size(∇S, 1)
            ∇ℓ[jpar] += real(tr(Ipart * ∇S[jpar, iω] * invS))
        end
    end
    return nothing
end
grad_generalwhittle!(∇ℓ, funcmemory, gradmemory, data::GenWhittleData) = @views grad_generalwhittle!(∇ℓ, unpack(funcmemory)[data.Ω_used_index], unpack(gradmemory)[:,data.Ω_used_index], data.I)
grad_generalwhittle!(∇ℓ, store::TimeSeriesModelStorageGradient, data::GenWhittleData) = grad_generalwhittle!(∇ℓ, store.funcmemory, store.gradmemory, data)
grad_generalwhittle!(∇ℓ, funcmemory, store::TimeSeriesModelStorageGradient, data::GenWhittleData) = grad_generalwhittle!(∇ℓ, funcmemory, store.gradmemory, data)
function grad_generalwhittle!(∇ℓ, objS, store::AdditiveStorage, data::GenWhittleData)
    @views grad_generalwhittle!(∇ℓ[1:store.npar1], objS, store.store1, data)
    @views grad_generalwhittle!(∇ℓ[store.npar1+1:end], objS, store.store2, data)
end
grad_generalwhittle!(∇ℓ, store::AdditiveStorage, data::GenWhittleData) = grad_generalwhittle!(∇ℓ, extract_S(store), store::AdditiveStorage, data::GenWhittleData)

## hessian
# Compute the hessian of the Whittle likelihood from the Periodogram and some spectral quantity (sdf or EI).
function hess_generalwhittle!(∇²ℓ, S, ∇S, ∇²S, I::Vector{T}) where {T <: Number} # univariate
    length(S) == length(I) == size(∇S,2) == size(∇²S,2) || throw(ArgumentError("S and I should be same length as second dim of ∇S and ∇²S."))
    ∇²ℓ .= zero(eltype(∇²ℓ))
    @inbounds for iω = 1:length(S) # fill lower triangle
        invS = inv(S[iω])
        otherIpart = (I[iω] * invS)
        Ipart = (1 - otherIpart)
        for jpar = 1:size(∇S, 1)
            ∂SⱼinvS = ∇S[jpar, iω] * invS
            for kpar = 1:jpar
                ∂SₖinvS = ∇S[kpar, iω] * invS
                ∇²ℓ[jpar, kpar] += real((otherIpart * ∂SⱼinvS) * ∂SₖinvS + Ipart * (∇²S[indexLT(jpar,kpar,size(∇S, 1)), iω] * invS - ∂SⱼinvS * ∂SₖinvS) )
            end
        end
    end
    for jpar = 1:size(∇S, 1)-1, kpar = jpar+1:size(∇S, 1) # fill upper triangle from symmetry
        @views ∇²ℓ[jpar, kpar] = ∇²ℓ[kpar, jpar]
    end
    return nothing
end
function hess_generalwhittle!(∇²ℓ, S, ∇S, ∇²S, I::Vector{SMatrix{D,D,T,n}}) where {D,T,n}
    length(S) == length(I) == size(∇S,2) == size(∇²S,2) || throw(ArgumentError("S and I should be same length as second dim of ∇S and ∇²S."))
    ∇²ℓ .= zero(eltype(∇²ℓ))
    Id = diagm(@SVector ones(D))
    @inbounds for iω = 1:length(S) # fill lower triangle
        invS = inv(S[iω])
        otherIpart = (I[iω] * invS)
        Ipart = (Id - otherIpart)
        for jpar = 1:size(∇S, 1)
            ∂SⱼinvS = ∇S[jpar, iω] * invS
            for kpar = 1:jpar
                ∂SₖinvS = ∇S[kpar, iω] * invS
                ∇²ℓ[jpar, kpar] += real(tr( (otherIpart * ∂SⱼinvS) * ∂SₖinvS + Ipart * (∇²S[indexLT(jpar,kpar,size(∇S, 1)), iω] * invS - ∂SⱼinvS * ∂SₖinvS) ))
            end
        end
    end
    @inbounds for jpar = 1:size(∇²ℓ, 1)-1, kpar = jpar+1:size(∇²ℓ, 1) # fill upper triangle from symmetry
        @views ∇²ℓ[jpar, kpar] = ∇²ℓ[kpar, jpar]
    end
    return nothing
end
hess_generalwhittle!(∇²ℓ, funcmemory, gradmemory, hessmemory, data::GenWhittleData) = @views hess_generalwhittle!(∇²ℓ, unpack(funcmemory)[data.Ω_used_index], unpack(gradmemory)[:,data.Ω_used_index], unpack(hessmemory)[:,data.Ω_used_index], data.I)
hess_generalwhittle!(∇²ℓ, store::TimeSeriesModelStorageHessian, data::GenWhittleData) = hess_generalwhittle!(∇²ℓ, store.funcmemory, store.gradmemory, store.hessmemory, data)
hess_generalwhittle!(∇²ℓ, funcmemory, store::TimeSeriesModelStorageHessian, data::GenWhittleData) = hess_generalwhittle!(∇²ℓ, funcmemory, store.gradmemory, store.hessmemory, data)
function hess_generalwhittle!(∇²ℓ, objS, store::AdditiveStorage, data)
    @views hess_generalwhittle!(∇²ℓ[1:store.npar1,1:store.npar1], objS, store.store1, data)
    @views hess_generalwhittle!(∇²ℓ[store.npar1+1:end,store.npar1+1:end], objS, store.store2, data)
    @views offdiag_add_hess_generalwhittle!(∇²ℓ[1:store.npar1,store.npar1+1:end], objS, store.store1, store.store2, data)
    for i = store.npar1+1:size(∇²ℓ,2), j = 1:store.npar1
        @views ∇²ℓ[i,j] = ∇²ℓ[j,i]
    end
    return nothing
end
hess_generalwhittle!(∇ℓ, store::AdditiveStorage, data::GenWhittleData) = hess_generalwhittle!(∇ℓ, extract_S(store), store, data)

# Compute the offdiagonal parts of the hessian when the model is additive.
function offdiag_add_hess_generalwhittle!(∇²ℓ, S, ∇S1, ∇S2, I::Vector{T}) where {T<:Number} # univariate
    length(S) == length(I) == size(∇S1,2) == size(∇S2,2) || throw(ArgumentError("S and I should be same length as second dim of ∇S1 and ∇S2."))
    ∇²ℓ .= zero(eltype(∇²ℓ))
    @inbounds for iω = 1:length(S) # fill lower triangle
        invS = inv(S[iω])
        otherIpart = (I[iω] * invS)
        Ipart = (1 - otherIpart)
        for jpar = 1:size(∇S1, 1)
            ∂SⱼinvS = ∇S1[jpar, iω] * invS
            for kpar = 1:size(∇S2, 1)
                ∂SₖinvS = ∇S2[kpar, iω] * invS
                ∇²ℓ[jpar, kpar] += real( (otherIpart * ∂SⱼinvS)*∂SₖinvS - Ipart*(∂SⱼinvS * ∂SₖinvS) )
            end
        end
    end
    return nothing
end
function offdiag_add_hess_generalwhittle!(∇²ℓ, S, ∇S1, ∇S2, I::Vector{SMatrix{D,D,T,n}}) where {D,T,n}
    length(S) == length(I) == size(∇S1,2) == size(∇S2,2) || throw(ArgumentError("S and I should be same length as second dim of ∇S1 and ∇S2."))
    ∇²ℓ .= zero(eltype(∇²ℓ))
    Id = diagm(@SVector ones(D))
    for iω = 1:length(S) # fill lower triangle
        invS = inv(S[iω])
        otherIpart = (I[iω] * invS)
        Ipart = (Id - otherIpart)
        for jpar = 1:size(∇S1, 1)
            ∂SⱼinvS = ∇S1[jpar, iω] * invS
            for kpar = 1:size(∇S2, 1)
                ∂SₖinvS = ∇S2[kpar, iω] * invS
                ∇²ℓ[jpar, kpar] += real(tr( (otherIpart * ∂SⱼinvS)*∂SₖinvS - Ipart*(∂SⱼinvS * ∂SₖinvS) ))
            end
        end
    end
    return nothing
end
offdiag_add_hess_generalwhittle!(∇²ℓ, funcmemory, gradmemory1, gradmemory2, data::GenWhittleData) = @views offdiag_add_hess_generalwhittle!(∇²ℓ, unpack(funcmemory)[data.Ω_used_index], unpack(gradmemory1)[:,data.Ω_used_index], unpack(gradmemory2)[:,data.Ω_used_index], data.I)
offdiag_add_hess_generalwhittle!(∇²ℓ, funcmemory, store1::TimeSeriesModelStorageHessian, store2::TimeSeriesModelStorageHessian, data::GenWhittleData) = offdiag_add_hess_generalwhittle!(∇²ℓ, funcmemory, store1.gradmemory, store2.gradmemory, data)

function offdiag_add_hess_generalwhittle!(∇²ℓ, objS, store1::AdditiveStorage, store2::AdditiveStorage, data::GenWhittleData)
    @views offdiag_add_hess_generalwhittle!(∇²ℓ[1:store1.npar1,:], objS, store1.store1, store2, data)
    @views offdiag_add_hess_generalwhittle!(∇²ℓ[store1.npar1+1:end,:], objS, store1.store2, store2, data)
end

function offdiag_add_hess_generalwhittle!(∇²ℓ, objS, store1::AdditiveStorage, store2, data::GenWhittleData)
    @views offdiag_add_hess_generalwhittle!(∇²ℓ[1:store1.npar1,:], objS, store1.store1, store2, data)
    @views offdiag_add_hess_generalwhittle!(∇²ℓ[store1.npar1+1:end,:], objS, store1.store2, store2, data)
end

function offdiag_add_hess_generalwhittle!(∇²ℓ, objS, store1, store2::AdditiveStorage, data::GenWhittleData)
    @views offdiag_add_hess_generalwhittle!(∇²ℓ[:,1:store1.npar1], objS, store1, store2.store1, data)
    @views offdiag_add_hess_generalwhittle!(∇²ℓ[:,store1.npar1+1:end], objS, store1, store2.store2, data)
end


text:
using BenchmarkTools
const SUITE = BenchmarkGroup()
for file in sort(readdir(@__DIR__))
    if startswith(file, "bench_") && endswith(file, ".jl")
        SUITE[chop(file, head = length("bench_"), tail = length(".jl"))] =
            include(file)
    end
end


text:
function pushing!_animate(tumor, root, cellbox, dimv::Val)
    queue = [root]
    N = length(tumor)

    points = [Point2f0(cell.position...) for cell in tumor]
    pointsNode = Node(points)
    colors = fill(:lightgrey, length(tumor))
    colorsNode = Node(colors)

    scene = meshscatter(pointsNode, markersize = 1.0, color = colorsNode, scale_plot = false, shading=false)
    xlims!(-30,30)
    ylims!(-30,30)
    scene[Axis][:ticks][:ranges] = (-20:4:20, -20:4:20)
    display(scene)

    #record(scene, "test.mp4") do io
        while !isempty(queue)
            root = popfirst!(queue)
            r1 = tumor[root].position

            setindex!(colors, :blue, root)

            for n in shuffle!( findall(find_neighbors(cellbox, root; s=2)) )
                r2 = tumor[n].position
                d = r2.-r1

                setindex!(colors, :orange, n)
                colorsNode[] = colors
                sleep(0.01)
                # recordframe!(io)

                if norm(d) < 2.

                    setindex!(colors, :red, n)
                    colorsNode[] = colors
                    sleep(0.01)
                    # recordframe!(io)

                    r2 .+= d*(2.2/norm(d) - 1)#/2
                    #r1 .-= d*(2.2/norm(d) - 1)/2
                    cellbox[n] = pos2box(r2, dimv)
                    push!(queue, n)

                    setindex!(points, Point2f0(r2...), n)
                    pointsNode[] = points
                    sleep(0.01)
                    setindex!(colors, :grey, n)
                    # recordframe!(io)
                else
                    setindex!(colors, :lightgrey, n)
                end
            end
            setindex!(colors, :lightgrey, root)
        end
        # recordframe!(io)
    #end
end

function pushing_recursive!_animate(tumor, root, cellbox, dimv::Val)
    points = [Point2f0(cell.position...) for cell in tumor]
    pointsNode = Node(points)
    colors = fill(:lightgrey, length(tumor))
    colorsNode = Node(colors)

    scene = meshscatter(pointsNode, markersize = 1.0, color = colorsNode, scale_plot = false, shading=false)
    xlims!(-30,30)
    ylims!(-30,30)
    scene[Axis][:ticks][:ranges] = (-20:4:20, -20:4:20)
    display(scene)
    queue = Int[]
    record(scene, "test_rec.mp4") do io
        pushing_recursive!_animate_exec(points, pointsNode, colors, colorsNode, io, tumor, root, cellbox, dimv, queue)
        recordframe!(io)
    end
end

function pushing_recursive!_animate_exec(points, pointsNode, colors, colorsNode, io, tumor, root, cellbox, dimv, queue)
    r1 = tumor[root].position

    colors[root] = :blue
    push!(queue, root)

    N=length(tumor)
    for n in shuffle!( findall(find_neighbors(cellbox, root; s=2)) )

        setindex!(colors, :orange, n)
        colorsNode[] = colors
        sleep(0.01)
        recordframe!(io)

        r2 = tumor[n].position
        d = r2.-r1
        if norm(d) < 2.

            setindex!(colors, :red, n)
            colorsNode[] = colors
            sleep(0.01)
            recordframe!(io)

            r2 .+= d*(2.2/norm(d) - 1)#/2
            #r1 .-= d*(2.2/norm(d) - 1)/2
            cellbox[n] = pos2box(r2, dimv)

            setindex!(points, Point2f0(r2...), n)
            pointsNode[] = points
            sleep(0.01)
            setindex!(colors, :grey, root)
            recordframe!(io)

            pushing_recursive!_animate_exec(points, pointsNode, colors, colorsNode, io, tumor, n, cellbox, dimv, queue)
        else
            colors[n] = n in queue ? :grey : :lightgrey
        end
    end
    pop!(queue)
    colors[root] = :lightgrey
end

max_stars_repo_name: aangaji/TumorGrowth

max_issues_repo_name: aangaji/TumorGrowth

max_forks_repo_name: aangaji/TumorGrowth


text:
# add AWS   # in pkg mode
# add AWSS3, Serialization

using AWS, AWSS3, Serialization
struct MyData
  a::Int
  b::String
end

d = MyData(1,"xyz")
aws = global_aws_config(; region="us-west-2")

s3_create_bucket(aws, "my.bucket")

b = IOBuffer()

serialize(b, d)
s3_put(aws, "your-s3-bucket-name","myfile.bin", b.data)

ddat = s3_get(aws, "your-s3-bucket-name","myfile.bin")
d2 = deserialize(IOBuffer(ddat))

@assert d == d2

p = S3Path("s3://my.bucket/test1.txt")  # provides an filesystem-like interface
write(p, "some data")
read(p, byte_range=1:4)  # returns b"some"

s3_delete_bucket("my.bucket")


text:
module TestFoldsCUDATests
const HAS_FOLDSCUDATESTS = try
    using FoldsCUDATests
    true
catch err
    @info "Failed to import `FoldsCUDATests`" exception = (err, catch_backtrace())
    false
end
if HAS_FOLDSCUDATESTS
    FoldsCUDATests.runtests_unionarrays()
end
end  # module


text:
Base.Broadcast.broadcastable(a::Angle) = Ref(a)

deg(ang::Real) = Angle(deg2rad(ang))
rad(ang::Real) = Angle(ang)
deg(ang::Angle) = rad2deg(ang.rad)
rad(ang::Angle) = ang.rad

Base.:+(a1::Angle, a2::Angle) = Angle(a1.rad + a2.rad)
Base.:-(a1::Angle, a2::Angle) = Angle(a1.rad - a2.rad)
Base.:-(a::Angle) = Angle(-a.rad)
Base.:*(a::Angle, r::Real) = Angle(a.rad * r)
Base.:*(r::Real, a::Angle) = a * r
Base.:/(a::Angle, r::Real) = Angle(a.rad / r)
Base.:/(a1::Angle, a2::Angle) = a1.rad / a2.rad
Base.cos(a::Angle) = cos(a.rad)
Base.sin(a::Angle) = sin(a.rad)
Base.tan(a::Angle) = tan(a.rad)
Base.isless(a1::Angle, a2::Angle) = a1.rad < a2.rad
Base.isgreater(a1::Angle, a2::Angle) = a1.rad > a2.rad
Base.isequal(a1::Angle, a2::Angle) = a1.rad == a2.rad


text:
using LatLongDF
using Test

@testset "LatLongDF.jl" begin
    # Write your tests here.
end


text:
function shift3(a, b, c, d)
    # a = b
    # b = c
    # c = d
    return b, c, d
end

function shift2(a, b, c)
    a = b
    b = c
    return a, b
end

middle(x1, x2) = sqrt(x1 * x2)

"""
Golden Section Search for the maximum in a one dimensional interval
"""
# function findCandidate(basis, proj, gmin::Float, gmax::Float)
#     @assert gmax > gmin
#     if gmin < Float(0.01)
#         gmin = Float(0.01)
#     end

#     R::Float = (sqrt(Float(5)) - 1) / 2
#     C::Float = 1 - R
#     # tol = 0.01
#     # R = 0.4
#     # C = 0.6

#     ax, cx = log(gmin), log(gmax)
#     bx = (ax + cx) / 2
#     # bx = R * ax + C * cx
#     # bx = sqrt(ax * cx)  # if ax<1, then replace it with one

#     # x0, x3 = gmin, gmax
#     # x2 = middle(gmin, gmax)
#     # x1 = middle(x0, x2)

#     x0, x3 = ax, cx
#     # at any given time, we will keep track of four points x0, x1, x2, x3
#     # make x0 to x1 the smaller segment and fill in the new point to be tried
#     if abs(cx - bx) < abs(bx - ax)
#         x1, x2 = bx, bx + C * (cx - bx)
#     else
#         x2, x1 = bx, bx - C * (bx - ax)
#     end
#     @assert x0 <= x1 <= x2 <= x3 "inital $(exp(x0))  $(exp(x1))  $(exp(x2))  $(exp(x3))"
#     # println("$x0, $x1, $x2, $x3")
#     # the initial function evaluations. Note that we never need to evaluate the function at the original endpoints.
#     f1, f2 = -Residual(basis, proj, exp(x1)), -Residual(basis, proj, exp(x2))
#     N = 2
#     # println("$f1  vs $f2")
#     # while abs(x3 - x0) > max(Float(1), ax) / Float(100)
#     while abs(x3 - x0) > 0.1
#         println("old $(exp(x0))  $(exp(x1))  $(exp(x2))  $(exp(x3))")
#         println("old $f1 vs $f2")
#         if f2 < f1
#             # println("this one: $x0, $x1, $x2, $(R * x1 + C * x3)")
#             println("one")
#             x0, x1, x2 = shift3(x0, x1, x2, R * x1 + C * x3)
#             # x0, x1, x2 = shift3(x0, x1, x2, x1^R * x3^C)
#             # println("$x1 x $x3 --> $(sqrt(x1 * x3))")
#             # x0, x1, x2 = shift3(x0, x1, x2, sqrt(x1 * x3))
#             f1, f2 = shift2(f1, f2, -Residual(basis, proj, exp(x2)))
#             # println("after one: $x0, $x1, $x2: $(x2 - x1)")
#         else
#             println("two")
#             x3, x2, x1 = shift3(x3, x2, x1, R * x2 + C * x0)
#             # x3, x2, x1 = shift3(x3, x2, x1, sqrt(x2 * x0))
#             f2, f1 = shift2(f2, f1, -Residual(basis, proj, exp(x1)))
#         end
#         println("new $(exp(x0))  $(exp(x1))  $(exp(x2))  $(exp(x3))")
#         println("new $f1 vs $f2")
#         @assert x0 <= x1 <= x2 <= x3 "$(exp(x0))  $(exp(x1))  $(exp(x2))  $(exp(x3))"
#         N += 1
#     end
#     println("$N eval")
#     if f1 < f2
#         return exp(x1)
#     else
#         return exp(x2)
#     end
# end

function findCandidate(basis, proj, gmin::Float, gmax::Float)
    @assert gmax > gmin
    if gmin < Float(0.01)
        gmin = Float(0.01)
    end

    R::Float = (sqrt(Float(5)) - 1) / 2
    C::Float = 1 - R
    # tol = 0.01
    # R = 0.4
    # C = 0.6

    a, b = log(gmin), log(gmax)

    h = b - a
    # if h <= tol:
    #     return (a, b)

    c = a + C * h
    d = a + R * h
    yc = -Residual(basis, proj, exp(c))
    yd = -Residual(basis, proj, exp(d))
    N = 2

    # for k in range(n-1):
    # while abs(b - a) > 0.1 && abs(yc - yd) > Float(0) * 10
    while abs(b - a) > 0.1
        # println("old $(exp(a))  $(exp(c))  $(exp(d))  $(exp(b))")
        # println("$yc  vs  $yd")
        if yc < yd
            b = d
            d = c
            yd = yc
            h = R * h
            c = a + C * h
            yc = -Residual(basis, proj, exp(c))
        else
            a = c
            c = d
            yc = yd
            h = R * h
            d = a + R * h
            yd = -Residual(basis, proj, exp(d))
        end
        # println("new $(exp(a))  $(exp(c))  $(exp(d))  $(exp(b))")
        # println("$yc  vs  $yd")
        N += 1
        @assert a <= c <= d <= b "$(exp(a))  $(exp(c))  $(exp(d))  $(exp(b))"
    end

    # println("max eval = $N")
    if yc < yd
        return exp(c)
    else
        return exp(d)
    end
end

# function findCandidate(basis, proj, gmin::Float, gmax::Float)
#     @assert gmax > gmin

#     if abs(gmin) < 100 * eps(Float(0)) && gmax > 100
#         gmax = 100 # if the first grid is 0, then the maximum should be between (0, 100)
#     end
#     if abs(gmin) > 100 * eps(Float(0)) && gmax / gmin > 100
#         gmax = gmin * 100  # the maximum won't be larger than 100*gmin
#     end

#     N = 32
#     dg = (gmax - gmin) / N

#     ###################   if gmin/gmax are at the boundary 0/Λ, the maximum could be at the edge ##################
#     if abs(gmin) < 100 * eps(Float(0)) && Residual(basis, proj, gmin) > Residual(basis, proj, gmin + dg)
#         return gmin
#     end
#     if abs(gmax - basis.Λ) < 100 * eps(Float(gmax)) && Residual(basis, proj, gmax) > Residual(basis, proj, gmax - dg)
#         return gmax
#     end

#     ###################  the maximum must be between (gmin, gmax) for the remaining cases ##################
#     # check https://www.geeksforgeeks.org/find-the-maximum-element-in-an-array-which-is-first-increasing-and-then-decreasing/ for detail

#     l, r = 1, N - 1 # avoid the boundary gmin and gmax
#     while l <= r
#         m = l + Int(round((r - l) / 2))
#         g = gmin + m * dg

#         r1, r2, r3 = Residual(basis, proj, g - dg), Residual(basis, proj, g), Residual(basis, proj, g + dg)
#         if r2 >= r1 && r2 >= r3
#             # plotResidual(basis, proj, gmin, gmax)
#             return g
#         end

#         if r3 < r2 < r1
#             r = m - 1
#         elseif r1 < r2 < r3
#             l = m + 1
#         else
#             if abs(r1 - r2) < 1e-17 && abs(r2 - r3) < 1e-17
#                 return g
#             end
#             println("warning: illegl! ($l, $m, $r) with ($r1, $r2, $r3)")
#             # plotResidual(basis, proj, gmin, gmax)
#             exit(0)
#         end
#     end
#     # plotResidual(basis, proj, gmin, gmax)
#     throw("failed to find maximum between ($gmin, $gmax)!")
# end


text:
"""
    longest_repeated_substrings(t,overlap_allowed)

tests the text t for the existence of longest repeated substrings (LRSs) of length 2 or more and returns:

- the length of a longest repeated substring in t
- the set of longest repeated substrings in t
- per LRS a list of its occurences (starting positions).

If no such strings exists: an empty set is returned and length nothing.

# Arguments

* t: the text (String,Array,Tuple) to find longest repeated substrings in.
* overlap_allowed: whether the occurences of a single LRS may overlap.

# Examples

    longest_repeated_substrings("ananas") == (3, Dict("ana" => [1, 3]))

    longest_repeated_substrings("ananas"; overlap_allowed=false) == (2, Dict("an" => [1, 3], "na" => [2, 4]))

#

"""
function longest_repeated_substrings(t; overlap_allowed::Bool = true)

    (nodes, root) = mk_tree(t)

    if overlap_allowed
        lenlongest, strlongest = longest_repeated_substring_T(t, nodes, root)

        return ((lenlongest, Dict(k => sort(collect(v)) for (k, v) in strlongest)))
    else

        lenlongest, strlongest = longest_maximal_repeated_no_overlap_all(t, nodes, root)

        for (s, l) in strlongest
            t = reverse(sort(l))
            last = pop!(t)
            new = [last]

            while length(t) > 0
                p = pop!(t)
                if last + lenlongest <= p
                    push!(new, p)
                    last = p
                end
            end
            strlongest[s] = new
        end

        return (lenlongest, strlongest)
    end

end


text:
@testset "constructors, IC.flat" begin
    @test IC.composite() == Never()
    c1 = Patience(1)
    c2 = InvalidValue()
    c3 = Step(1)

    @test IC.composite(c1) == c1

    d = IC.composite(c1, c2, Never(), c3, c1)
    show(d)
    println()
    @test IC.flat(d) == (c1, c2, c3)

    # codecov:
    @test IC._in(Never(), c1)
    @test IC._in(Never(), d)
end

@testset "behaviour" begin
    a = Step(4)
    b = NumberLimit(2)
    c = InvalidValue()
    d = IC.composite(a, b, c)

    # separated:
    m = SquareRooter(4)
    state_a1 = IC.update!(a, m, 0, 1)
    state_b1 = IC.update!(b, m, 0, 1)
    state_c1 = IC.update!(c, m, 0, 1)
    state_a2 = IC.update!(a, m, 0, 2, state_a1)
    state_b2 = IC.update!(b, m, 0, 2, state_b1)
    state_c2 = IC.update!(c, m, 0, 2, state_c1)
    done_a = IC.done(a, state_a2)
    done_b = IC.done(b, state_b2)
    done_c = IC.done(c, state_c2)
    report_a = IC.takedown(a, 0, state_a2)
    report_b = IC.takedown(b, 0, state_b2)
    report_c = IC.takedown(c, 0, state_c2)

    # composed:
    m = SquareRooter(4)
    state_d1 = IC.update!(d, m, 0, 1)
    @test IC.flat(state_d1) == (state_a1, state_b1, state_c1)
    state_d2 = IC.update!(d, m, 0, 2, state_d1)
    @test IC.flat(state_d2) == (state_a2, state_b2, state_c2)
    done_d = IC.done(d, state_d2)
    @test done_d == done_a || done_b || done_c
    report_d = IC.takedown(d, 0, state_d2)
    @test report_d == ((a, report_a), (b, report_b), (c, report_c))
end

@testset "traits" begin
    needs_nothing = NumberLimit(4)
    needs_loss = Threshold(0.01)
    needs_training = WithTrainingLossesDo()
    needs_both = PQ()

    @test !IC.needs_loss(IC.composite(needs_nothing))
    @test !IC.needs_loss(IC.composite(needs_nothing, needs_training))
    @test IC.needs_loss(IC.composite(needs_nothing, needs_loss))
    @test IC.needs_loss(IC.composite(needs_nothing, needs_both))
    @test IC.needs_loss(IC.composite(needs_loss))
    @test IC.needs_loss(IC.composite(needs_loss, needs_nothing))

    @test !IC.needs_training_losses(IC.composite(needs_nothing))
    @test IC.needs_training_losses(IC.composite(needs_nothing, needs_training))
    @test !IC.needs_training_losses(IC.composite(needs_nothing, needs_loss))
    @test IC.needs_training_losses(IC.composite(needs_nothing, needs_both))
    @test IC.needs_training_losses(IC.composite(needs_training))
    @test !IC.needs_training_losses(IC.composite(needs_loss, needs_nothing))
end

true


text:
"""
# Stipple.Layout

Utilities for rendering the general layout of a Stipple app, such as of a data dashboard web page or supporting themes.
"""
module Layout

import Genie
using Stipple

export layout

"""
    `function layout(output::Union{String,Vector}; partial::Bool = false, title::String = "", class::String = "", style::String = "",
                      head_content::String = "", channel::String = Genie.config.webchannels_default_route) :: String`

Utility for creating a basic web page structure, including doctype as well as <HTML>, <HEAD>, <TITLE>, <META viewport>,
  and <BODY> tags, together with the output content.

If `partial` is `true`, the page layout HTML elements are not returned.

### Examples

```julia
julia> layout([
        span("Hello", @text(:greeting))
        ])
"<!DOCTYPE html>\n<html><head><title></title><meta name=\"viewport\" content=\"width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no, minimal-ui\" /></head><body class style><link href=\"https://fonts.googleapis.com/css?family=Material+Icons\" rel=\"stylesheet\" /><link href=\"https://fonts.googleapis.com/css2?family=Lato:ital,wght@0,400;0,700;0,900;1,400&display=swap\" rel=\"stylesheet\" /><link href=\"/css/stipple/stipplecore.min.css\" rel=\"stylesheet\" /><link href=\"/css/stipple/quasar.min.css\" rel=\"stylesheet\" /><span v-text='greeting'>Hello</span><script src=\"/js/channels.js?v=1.17.1\"></script><script src=\"/js/stipple/underscore-min.js\"></script><script src=\"/js/stipple/vue.js\"></script><script src=\"/js/stipple/quasar.umd.min.js\"></script>\n<script src=\"/js/stipple/apexcharts.min.js\"></script><script src=\"/js/stipple/vue-apexcharts.min.js\"></script><script src=\"/js/stipple/stipplecore.js\" defer></script><script src=\"/js/stipple/vue_filters.js\" defer></script></body></html>"
```

```julia
julia> layout([
        span("Hello", @text(:greeting))
        ], partial = true)
"<link href=\"https://fonts.googleapis.com/css?family=Material+Icons\" rel=\"stylesheet\" /><link href=\"https://fonts.googleapis.com/css2?family=Lato:ital,wght@0,400;0,700;0,900;1,400&display=swap\" rel=\"stylesheet\" /><link href=\"/css/stipple/stipplecore.min.css\" rel=\"stylesheet\" /><link href=\"/css/stipple/quasar.min.css\" rel=\"stylesheet\" /><span v-text='greeting'>Hello</span><script src=\"/js/channels.js?v=1.17.1\"></script><script src=\"/js/stipple/underscore-min.js\"></script><script src=\"/js/stipple/vue.js\"></script><script src=\"/js/stipple/quasar.umd.min.js\"></script>\n<script src=\"/js/stipple/apexcharts.min.js\"></script><script src=\"/js/stipple/vue-apexcharts.min.js\"></script><script src=\"/js/stipple/stipplecore.js\" defer></script><script src=\"/js/stipple/vue_filters.js\" defer></script>"
```
"""
function layout(output::Union{String,Vector}; partial::Bool = false, title::String = "", class::String = "", style::String = "",
                head_content::String = "", channel::String = Genie.config.webchannels_default_route, core_theme::Bool = true) :: String

  isa(output, Vector) && (output = join(output, '\n'))

  content = string(
    theme(; core_theme = core_theme),
    output,
    Stipple.deps(channel)
  )

  partial && return content

  Genie.Renderer.Html.doc(
    Genie.Renderer.Html.html([
      Genie.Renderer.Html.head([
        Genie.Renderer.Html.title(title)
        Genie.Renderer.Html.meta(name="viewport", content="width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no, minimal-ui")
        head_content
      ])
      Genie.Renderer.Html.body(content, class=class, style=style)
    ])
  )
end

include(joinpath("layout", "page.jl"))
include(joinpath("layout", "theme.jl"))

end


text:
# Rust (1987) bus engine replacement problem
# Notation: I use underscore to denote empirical counterparts

# Import packages
using Optim
using Distributions
using Statistics
using DataFrames
using CSV
using LinearAlgebra

function compute_U(θ::Vector, s::Vector)::Matrix
    """Compute static utility"""
    u1 = - θ[1]*s - θ[2]*s.^2       # Utility of not investing
    u2 = - θ[3]*ones(size(s))       # Utility of investing
    U = [u1 u2]                     # Combine in a matrix
    return U
end;

function compute_Vbar(θ::Vector, λ::Number, β::Number, s::Vector)::Matrix
    """Compute value function by Bellman iteration"""
    k = length(s)                                 # Dimension of the state space
    U = compute_U(θ, s)                           # Static utility
    index_λ = Int[1:k [2:k; k]];                  # Mileage index
    index_A = Int[1:k ones(k,1)];                 # Investment index
    γ = Base.MathConstants.eulergamma             # Euler's gamma

    # Iterate the Bellman equation until convergence
    Vbar = zeros(k, 2);
    Vbar1 = Vbar;
    dist = 1;
    iter = 0;
    while dist>1e-8
        V = γ .+ log.(sum(exp.(Vbar), dims=2))     # Compute value
        expV = V[index_λ] * [1-λ; λ]               # Compute expected value
        Vbar1 =  U + β * expV[index_A]             # Compute v-specific
        dist = max(abs.(Vbar1 - Vbar)...);         # Check distance
        iter += 1;
        Vbar = Vbar1                               # Update value function
    end
    return Vbar
end;

function generate_data(θ::Vector, λ::Number, β::Number, s::Vector, N::Int)::Tuple
    """Generate data from primitAes"""
    Vbar = compute_Vbar(θ, λ, β, s)             # Solve model
    ε = rand(Gumbel(0,1), N, 2)                 # Draw shocks
    St = rand(s, N)                             # Draw states
    A = (((Vbar[St,:] + ε) * [-1;1]) .> 0)     # Compute investment decisions
    δ = (rand(Uniform(0,1), N) .< λ)            # Compute mileage shock
    St1 = min.(St .* (A.==0) + δ, max(s...))   # Compute nest state
    df = DataFrame(St=St, A=A, St1=St1)
    CSV.write("../data/rust.csv", df)
    return St, A, St1
end;

function logL_Rust(θ0::Vector, λ::Number, β::Number, s::Vector, St::Vector, A::BitVector)::Number
    """Compute log-likelihood functionfor Rust problem"""
    # Compute value
    Vbar = compute_Vbar(θ0, λ_, β, s)

    # Implied choice probabilities
    EP = exp.(Vbar[:,2]) ./ (exp.(Vbar[:,1]) + exp.(Vbar[:,2]))

    # Likelihood
    logL = sum(log.(EP[St[A.==1]])) + sum(log.(1 .- EP[St[A.==0]]))
    return -logL
end;

function compute_T(k::Int, λ_::Number)::Array
    """Compute transition matrix"""
    T = zeros(k, k, 2);

    # Conditional on not investing
    T[k,k,1] = 1;
    for i=1:k-1
        T[i,i,1] = 1-λ_
        T[i,i+1,1] = λ_
    end

    # Conditional on investing
    T[:,1,2] .= 1-λ_;
    T[:,2,2] .= λ_;

    return(T)
end;

function HM_inversion(CCP::Matrix, T::Array, U::Matrix, β::Number)::Vector
    """Perform HM inversion"""

    # Compute LHS (to be inverted)
    γ = Base.MathConstants.eulergamma
    LEFT = I - β .* (CCP[:,1] .* T[:,:,1] + CCP[:,2] .* T[:,:,2])

    # Compute LHS (not to be inverted)
    RIGHT = γ .+ sum(CCP .* (U .- log.(CCP)) , dims=2)

    # Compute V
    EV_ = inv(LEFT) * RIGHT
    return vec(EV_)
end;

function from_EV_to_EP(EV_::Vector, T::Array, U::Matrix, β::Number)::Vector
    """Compute expected policy from expected value"""
    E = exp.( U + β .* [(T[:,:,1] * EV_) (T[:,:,2] * EV_)] )
    EP_ = E[:,2] ./ sum(E, dims=2)
    return vec(EP_)
end;

function logL_HM(θ0::Vector, λ::Number, β::Number, s::Vector, St::Vector, A::BitVector, T::Array, CCP::Matrix)::Number
    """Compute log-likelihood function for HM problem"""
    # Compute static utility
    U = compute_U(θ0, s)

    # Espected value by inversion
    EV_ = HM_inversion(CCP, T, U, β)

    # Implies choice probabilities
    EP_ = from_EV_to_EP(EV_, T, U, β)

    # Likelihood
    logL = sum(log.(EP_[St[A.==1]])) + sum(log.(1 .- EP_[St[A.==0]]))
    return -logL
end;

function logL_AM(θ0::Vector, λ::Number, β::Number, s::Vector, St::Vector, A::BitVector, T::Array, CCP::Matrix, K::Int)::Number
    """Compute log-likelihood function for AM problem"""
    # Compute static utility
    U = compute_U(θ0, s)
    EP_ = CCP[:,2]

    # Iterate HM mapping
    for _=1:K
        EV_ = HM_inversion(CCP, T, U, β)    # Espected value by inversion
        EP_ = from_EV_to_EP(EV_, T, U, β)   # Implies choice probabilities
        CCP = [(1 .- EP_) EP_]
    end

    # Likelihood
    logL = sum(log.(EP_[St[A.==1]])) + sum(log.(1 .- EP_[St[A.==0]]))
    return -logL
end;

## Main

# Set parameters
θ = [0.13; -0.004; 3.1];
λ = 0.82;
β = 0.95;

# State space
k = 10;
s = Vector(1:k);

# Compute value function
Vbar = compute_Vbar(θ, λ, β, s);

# Generate data
N = Int(1e5);
St, A, St1 = generate_data(θ, λ, β, s, N);
print("\n\nwe observe ", sum(A), " investment decisions in ", N, " observations")

# Estimate lambda
Δ = St1 - St;
λ_ = mean(Δ[(A.==0) .& (St.<10)]);

print("\n\nEstimated lambda: $λ_ (true = $λ)")

# True likelihood value
logL_trueθ = logL_Rust(θ, λ, β, s, St, A);
print("\n\nThe likelihood at the true parameter is $logL_trueθ")

# Select starting values
θ0 = Float64[0,0,0];

# Optimize
θ_R = optimize(x -> logL_Rust(x, λ, β, s, St, A), θ0).minimizer;
print("\n\nEstimated thetas: $θ_R (true = $θ)")

# Not all initial values are equally good
θ0 = Float64[1,1,1];

# Optimize
θ_R2 = optimize(x -> logL_Rust(x, λ, β, s, St, A), θ0).minimizer;
print("\n\nEstimated thetas: $θ_R2 (true = $θ)")

# Estimate CCP
P = [mean(A[St.==i]) for i=s]
CCP = [(1 .- P) P]

# Compute T
T = compute_T(k, λ_)

# Optimize
θ0 = Float64[0,0,0];
θ_HM = optimize(x -> logL_HM(x, λ, β, s, St, A, T, CCP), θ0).minimizer;
print("\n\nEstimated thetas: $θ_HM (true = $θ)")

# Aguirregabiria & Mira (2002)
K = 2

# Optimize
θ0 = Float64[0,0,0];
θ_AM = optimize(x -> logL_AM(x, λ, β, s, St, A, T, CCP, K), θ0).minimizer;
print("\n\nEstimated thetas: $θ_AM (true = $θ)")

# Compare times
θ0 = Float64[0,0,0];
time_Rust = optimize(x -> logL_Rust(x, λ, β, s, St, A), θ0).time_run;
time_HM = optimize(x -> logL_HM(x, λ, β, s, St, A, T, CCP), θ0).time_run;
time_AM = optimize(x -> logL_AM(x, λ, β, s, St, A, T, CCP, K), θ0).time_run;
print("Time Rust: $time_Rust\nTime HM: $time_HM\nTime AM: $time_AM")