text:
# utilities for piece-wise linear functions

"slope and offset for line through two points"
function linecoefs(x1, y1, x2, y2)
    x1 ≈ x2 && error("line too steep!")

    # line equation:             y = slope*x + offset
    # from Taylor of 1st order:      f'(x0)*(x - x0) + f(x0)
    # use x1 for x0:             y = slope*(x - x1) + f(x1)
    slope = (y2 - y1)/(x2 - x1)
    offset = y1 - slope*x1
    slope, offset
end

"""
Return coefficient matrix [a b c] of overestimating inequalities of the form
  a*x + b*y ≥ c
for a piece-wise linear function given by
"""
function pwl_ineqs(xs, ys)
    @assert length(xs) == length(ys)
    m = length(xs) - 1

    ineqs = fill(0.0, (m, 3))
    for i = 1:m
        slope, offset = linecoefs(xs[i], ys[i], xs[i+1], ys[i+1])
        ineqs[i, :] = [-slope, 1.0, offset]
    end
    ineqs
end

"Compute the inverse of a monotone PWL function at a point"
function pwl_inverse(xs, ys, y)
    @assert length(xs) == length(ys)
    @assert issorted(xs)
    increasing, decreasing = issorted(ys), issorted(ys, rev=true)
    @assert increasing || decreasing
    if decreasing # make it increasing
        y, ys = -y, -ys
    end
    @assert isapproxin(y, ys[1], ys[end], atol=1e-4)

    right = findfirst(yi -> y ≤ yi, ys)
    if right == 1 # on the left boundary
        return xs[1]
    end

    left = right - 1
    @assert isapproxin(y, ys[left], ys[right])
    if ys[left] ≈ ys[right]
        return xs[left]
    end
    # y = (1-λ)*ys[left] + λ*ys[right]
    λ = (y - ys[left])/(ys[right] - ys[left])
    x = (1-λ)*xs[left] + λ*xs[right]
    @assert isapproxin(x, xs[left], xs[right])
    x
end


text:
using Test
using Receive: receive, Break

@testset "Both channels imidiatelly block at `take!`" begin let

    quit = Channel()
    results = Channel()

    t3 = @async begin
        sleep(0.5)
        for i in 1:10000
            put!(results, i)
        end
    end

    t4 = @async begin
        sleep(1)
        put!(quit, nothing)
    end

    state = 0

    receive(results, quit) do ch, value

        if ch==quit
            close(quit)
            return Break()
        end

        state += 1
    end

    @test state == 10000
    @test istaskdone(t3) == true

end
end

@testset "One channel blocks at `put!` other at `take!`" begin let

    quit = Channel()
    results = Channel()

    t3 = @async begin
        for i in 1:10000
            put!(results, i)
        end
    end

    t4 = @async begin
        sleep(1)
        put!(quit, nothing)
    end

    #take!(results)
    sleep(0.5)

    state = 0

    receive(results, quit) do ch, value

        if ch==quit
            close(quit)
            return Break()
        end

        state += 1

    end

    @test state == 10000
    @test istaskdone(t3) == true

end
end

@testset "Interupting while being flooded" begin let

    quit = Channel()
    results = Channel()

    t3 = @async begin
        for i in 1:1000000
            put!(results, i)
        end
    end

    t4 = @async begin
        sleep(1)
        put!(quit, nothing)
    end

    state = 0
    curvalue = 0

    receive(results, quit) do ch, value

        if ch==quit
            close(quit)
            return Break()
        end

        state += 1
        curvalue = value

    end

    @test state == curvalue
    @test state > 10000
    @test istaskdone(t3) == false
    @test take!(results) == state + 1

end
end

@testset "One channel blocks at `take!` frequently and other blocks" begin let

    quit = Channel()
    results = Channel()

    t3 = @async begin
        for i in 1:10
            sleep(0.1)
            put!(results, i)
        end
    end

    t4 = @async begin
        sleep(0.5)
        put!(quit, nothing)
    end

    state = 0

    receive(results, quit) do ch, value

        if ch==quit
            close(quit)
            return Break()
        end

        state = value
    end

    @test state > 3
    @test istaskdone(t3) == false
    @test take!(results) == state + 1

end
end

ext: jl


text:
__precompile__(true)

module PPQI

using PyCall

export load_config, load_model, model_forward

const inference = PyNULL()

function __init__()
    copy!(inference, pyimport("paddle.inference"))
end

struct InferenceModel
    predictor::PyObject
    input_handles::Vector{PyObject}
    output_handles::Vector{PyObject}
end

function load_config(modelpath::String, use_gpu::Bool=false, gpu_id::Int=0, use_mkldnn::Bool=false, cpu_treads::Int=1)::PyObject
    if isdir(modelpath)

        if isfile(joinpath(modelpath, "__params__"))
            model = joinpath(modelpath, "__model__")
            params = joinpath(modelpath, "__params__")
            config = inference.Config(model, params)

        elseif isfile(joinpath(modelpath, "params"))
            model = joinpath(modelpath, "model")
            params = joinpath(modelpath, "params")
            config = inference.Config(model, params)

        elseif isfile(joinpath(modelpath, "__model__"))
            config = inference.Config(modelpath)

        end

    elseif isfile(modelpath + ".pdmodel")
        model = modelpath + ".pdmodel"
        params = modelpath + ".pdiparams"
        config = inference.Config(model, params)

    end

    if use_gpu
        config.enable_use_gpu(100, gpu_id)

    else
        config.disable_gpu()
        config.set_cpu_math_library_num_threads(cpu_treads)

        if use_mkldnn
            config.enable_mkldnn()
        end

    end

    return config

end

function load_model(config::PyObject)::InferenceModel
    predictor = inference.create_predictor(config)
    input_names = predictor.get_input_names()
    output_names = predictor.get_output_names()
    input_handles = PyObject[]

    for input_name in input_names
        input_handle = predictor.get_input_handle(input_name)
        push!(input_handles, input_handle)
    end

    output_handles = PyObject[]

    for output_name in output_names
        output_handle = predictor.get_output_handle(output_name)
        push!(output_handles, output_handle)
    end

    model = InferenceModel(predictor, input_handles, output_handles)

    return model
end

function model_forward(model::InferenceModel, input_datas::Vector)::Vector
    for input_handle in model.input_handles, data in input_datas
        input_handle.copy_from_cpu(data)
    end

    model.predictor.run()
    output_datas = Array[]

    for output_handle in model.output_handles
        data = output_handle.copy_to_cpu()
        push!(output_datas, data)
    end

    return output_datas

end

end


text:
using Flux, DiffEqFlux, DifferentialEquations, Plots

## Setup ODE to optimize
function lotka_volterra(du,u,p,t)
  x, y = u
  α, β, δ, γ = p
  du[1] = dx = α*x - β*x*y
  du[2] = dy = -δ*y + γ*x*y
end
u0 = [1.0,1.0]
tspan = (0.0,10.0)
p = [1.5,1.0,3.0,1.0] #True Params
prob = ODEProblem(lotka_volterra,u0,tspan,p)

# Verify ODE solution
sol = solve(prob,Tsit5())
plot(sol)

# Generate data from the ODE
data_sol = solve(prob,Tsit5(),saveat=0.1)
A1 = data_sol[1,:] # length 101 vector
A2 = data_sol[2,:] # length 101 vector
t = 0:0.1:10.0
scatter!(t,A1,color=[1],label = "rabbits")
scatter!(t,A2,color=[2],label = "wolves")

# Build a neural network that sets the cost as the difference from the
# generated data and true data

p = param([4., 1.0, 2.0, 0.4]) # Initial Parameter Vector
function predict_rd() # Our 1-layer neural network
  diffeq_rd(p,prob,Tsit5(),saveat=0.1)
end
loss_rd() = sum(abs2,predict_rd()-data_sol) # loss function

# Optimize the parameters so the ODE's solution stays near 1

data = Iterators.repeated((), 100)
opt = ADAM(0.1)
cb = function () #callback function to observe training
  #= display(loss_rd()) =#
  # using `remake` to re-create our `prob` with current parameters `p`
  scatter(t,A1,color=[1],label = "rabbit data")
  scatter!(t,A2,color=[2],label = "wolves data")
  display(plot!(solve(remake(prob,p=Flux.data(p)),Tsit5(),saveat=0.1),ylim=(0,6),labels = ["rabbit model","wolf model"],color=[1 2]))
end
# Display the ODE with the initial parameter values.
cb()
Flux.train!(loss_rd, [p], data, opt, cb = cb)
p

# ----- Testing with proposed objective function -----

p = param([4., 1.0, 2.0, 0.4]) # Initial Parameter Vector
data_sim = reduce(hcat, data_sol.u)
function loss_rd() # Our 1-layer neural network
  sum(abs2.(adams_moulton_estimator(p, data_sim, t, lotka_volterra)))
end

# Optimize the parameters so the ODE's solution stays near 1

data = Iterators.repeated((), 100)
opt = ADAM(0.1)
cb = function () #callback function to observe training
  #= display(loss_rd()) =#
  # using `remake` to re-create our `prob` with current parameters `p`
  scatter(t,A1,color=[1],label = "rabbit data")
  scatter!(t,A2,color=[2],label = "wolves data")
  display(plot!(solve(remake(prob,p=Flux.data(p)),Tsit5(),saveat=0.1),ylim=(0,6),labels = ["rabbit model","wolf model"],color=[1 2]))
end
# Display the ODE with the initial parameter values.
cb()
Flux.train!(loss_rd, [p], data, opt, cb = cb)

# ----- Testing with Floudas' dataset -----

i = 1
prob = problem_set[i]
f = prob.fun
exp_data = prob.data
lb = prob.bounds[1]
ub = prob.bounds[end]
true_par = prob.phi
t = prob.t
shared_rand = rand(1:1234)
rng = MersenneTwister(shared_rand)
p_am = param([rand(rng, Uniform(lb[i], ub[i])) for i in 1:length(true_par)])
rng = MersenneTwister(shared_rand)
p_rd = param([rand(rng, Uniform(lb[i], ub[i])) for i in 1:length(true_par)])
println(p_am)
println(p_rd)

function loss_am() # Our 1-layer neural network
  sum(abs2.(adams_moulton_estimator(p_am, exp_data, t, f)))
end

# Optimize the parameters so the ODE's solution stays near 1

data = Iterators.repeated((), 1000)
opt_method = ADAM(2.0)
u0 = exp_data[:,1]
tspan = (t[1], t[end])
prob_instance = ODEProblem(f, u0, tspan)

function predict_rd() # Our 1-layer neural network
  diffeq_rd(p_rd,prob_instance,Tsit5(),saveat=t)
end
loss_rd() = sum(soft_l1.(abs2.((predict_rd()-exp_data)))) # loss function

cb_am = function () #callback function to observe training
  # @show(loss_am())
  # using `remake` to re-create our `prob` with current parameters `p`
  plot_canvas = plot()
  for i in 1:length(u0)
    scatter!(plot_canvas,t,transpose(exp_data[i,:]), color=[i], label="[A_$i]")
  end
  display(plot!(plot_canvas,solve(remake(prob_instance,p=Flux.data(p_am)),Tsit5(),saveat=t)))
end
# Display the ODE with the initial parameter values.
cb_am()
Flux.train!(loss_am, [p_am], data, opt_method, cb = cb_am)
println("True:\n$true_par\nEst:\n$p_am\nDistance:\n$(euclidean(true_par,p_am))")

cb_rd = function () #callback function to observe training
  # @show(loss_rd())
  # using `remake` to re-create our `prob` with current parameters `p`
  plot_canvas = plot()
  for i in 1:length(u0)
    scatter!(plot_canvas,t,transpose(exp_data[i,:]), color=[i], label="[A_$i]")
  end
  display(plot!(plot_canvas,solve(remake(prob_instance,p=Flux.data(p_rd)),Tsit5(),saveat=t)))
end
# Display the ODE with the initial parameter values.
cb_rd()
Flux.train!(loss_rd, [p_rd], data, opt_method, cb = cb_rd)

println("True:\n$true_par\nEst:\n$p_rd\nDistance:\n$(euclidean(true_par,p_rd))")
m = Chain(
          p_rd -> diffeq_rd(p_rd, prob_instance, Tsit5(), saveat=t),
          )
ps = Flux.params(m)
Flux.train!(loss_rd, ps, data, opt_method)

true_par
p_rd

# ----- Second example -----
u0 = Float32[2.; 0.]
datasize = 30
tspan = (0.0f0,1.5f0)

function trueODEfunc(du,u,p,t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end

function ODEfunc(du,u,p,t)
    du .= ((u.^3)'*p)'
end

t = range(tspan[1],tspan[2],length=datasize)
prob = ODEProblem(trueODEfunc,u0,tspan)
ode_data = Array(solve(prob,Tsit5(),saveat=t))
plot(transpose(ode_data))

par = [param(1.0) param(1.0); param(1.0) param(1.0)]
#par = param([1. 1.; 1. 1.])
function f1(x)
  data_shooting_estimator_node(par, x, t, ODEfunc; steps=1)
end
function f2(x)
  data_shooting_estimator_node(par, x, t, ODEfunc; steps=2)
end

function fe(x)
    sum(abs2.(data_shooting_estimator_endnode(par, x, t, ODEfunc; steps=3)))
end

dudt = Chain(x -> f1(x), x -> f2(x))
ps = Flux.params(dudt)

pred = n_ode(ode_data) # Get the prediction using the correct initial condition
scatter(t,ode_data[1,:],label="data")
scatter!(t,Flux.data(pred[1,:]),label="prediction")
function predict_n_ode()
  n_ode(u0)
end
loss_n_ode() = sum(abs2,ode_data .- predict_n_ode())

data = Iterators.repeated((), 100)
opt_method = ADAM(0.1)
cb = function () #callback function to observe training
  display(loss_n_ode())
  # plot current prediction against data
  cur_pred = Flux.data(predict_n_ode())
  pl = scatter(t,ode_data[1,:],label="data")
  scatter!(pl,t,cur_pred[1,:],label="prediction")
  display(plot(pl))
end

# Display the ODE with the initial parameter values.
cb()

dudt(ode_data)
loss_testing() = sum(abs2, dudt(ode_data).-ode_data)
Flux.train!(loss_testing, par, data, opt_method)
params(dudt)
par

# ----- Try simples example -----

# Parameter that I want to learn
learn_p = param([4.01;1.1])
# Function used for generating output using learn_p
function f_square_cube(x)
    x.^learn_p
end
# Single layer neural network
model_chain = Chain(x -> f_square_cube(x),
                    x -> f_square_cube(x)
                    )
tracking = params(learn_p)
loss(x,y) = sum(abs2, model_chain(x) - y)
rand_data = rand((1.0:0.1:100.0),2,300)
square_cube_rand = rand_data.^[4;9]
data = Iterators.repeated((rand_data, square_cube_rand), 2000)
Flux.train!(loss, tracking, data, ADAM(0.1))
tracking

#----- Now to ODEs -----
i = 1
p_solve = problem_set[i]
ode_fun = p_solve.fun
phi = p_solve.phi
t = p_solve.t
ode_data = p_solve.data
learn_p = param([rand(1:100) for i in 1:length(phi)])
learn_p
t = range(t[1], stop=t[end], length=10)
tspan = (t[1], t[end])
ini_cond = p_solve.data[:,1]
ode_prob = ODEProblem(p_solve.fun, ini_cond, tspan, p_solve.phi)
ode_sol  = solve(ode_prob, lsoda(), saveat=reduce(vcat, t))
data_original = reduce(hcat, ode_sol.u)
ode_data = copy(data_original)

function run_for_steps(steps, learn_p)
    function f_learn(x)
        data_shooting_estimator_node(learn_p,x,t,ode_fun; steps=steps)
    end
    # Single layer neural network
    model_chain = Chain(x -> f_learn(x))

    tracking = params(learn_p)
    function loss(x,y)
        res = Flux.mse(model_chain(x), y)
        if res == NaN
            res = 10^5
        end
        return res
    end
    data = Iterators.repeated((ode_data, ode_data), 500)
    Flux.train!(loss, tracking, data, ADAM(0.1))
    tracking
end

@time run_for_steps(1, learn_p)
run_for_steps(2, learn_p)
run_for_steps(3, learn_p)
run_for_steps(40, learn_p)

run_for_steps(40, learn_p)

# One step only
# Function used for generating output using learn_p
learn_p = param([7.71, 8.81])
function f_learn(x)
    data_shooting_estimator_node(learn_p,x,t,ode_fun; steps=1)
end
# Single layer neural network
model_chain = Chain(x -> f_learn(x))

tracking = params(learn_p)
loss(x,y) = sum(abs2, model_chain(x) - y)
data = Iterators.repeated((ode_data, ode_data), 10000)
Flux.train!(loss, tracking, data, ADAM(0.1))
tracking

# ----- Testing Adams-Moulton -----
i = 6
p_solve = problem_set[i]
ode_fun = p_solve.fun
phi = p_solve.phi
t = p_solve.t
ode_data = p_solve.data
t = p_solve.t
t = range(t[1], stop=t[end], length=50)
tspan = (t[1], t[end])
ini_cond = p_solve.data[:,1]
ode_prob = ODEProblem(p_solve.fun, ini_cond, tspan, p_solve.phi)
ode_sol  = solve(ode_prob, lsoda(), saveat=reduce(vcat, t))
data_original = reduce(hcat, ode_sol.u)
ode_data = copy(data_original)
plot(transpose(ode_data))
learn_p = param([rand(0.0001:0.0001:100) for i in 1:length(phi)])
tracking = params(learn_p)
function f_learn(x)
    adams_moulton_fourth_estimator(learn_p,x,t,ode_fun; return_estimated=true)
end
# Single layer neural network
model_chain = Chain(x -> f_learn(x))
function loss_mse(x,y)
    res = sum(abs2, model_chain(x) - y) + 10^10
    return res
end
data = Iterators.repeated((ode_data, ode_data), 1000)
function_error = []
function cb()
    err = Tracker.data(loss_mse(ode_data, ode_data))
    if err < 100
        push!(function_error, err)
    end
end
@time Flux.train!(loss_mse, tracking, data, ADAM(0.1), cb = cb)
plot(function_error[2:end])
tracking
euclidean(learn_p, phi)

# ----- Adams Moulton One on Each Neuron -----
i = 7
p_solve = problem_set[i]
ode_fun = p_solve.fun
phi = p_solve.phi
t = p_solve.t
ode_data = p_solve.data
t = p_solve.t
t = range(t[1], stop=t[end], length=50)
tspan = (t[1], t[end])
ini_cond = p_solve.data[:,1]
ode_prob = ODEProblem(p_solve.fun, ini_cond, tspan, p_solve.phi)
ode_sol  = solve(ode_prob, lsoda(), saveat=reduce(vcat, t))
data_original = reduce(hcat, ode_sol.u)
ode_data = copy(data_original)
plot(transpose(ode_data))
learn_p = param([rand(0.0001:0.0001:100) for i in 1:length(phi)])
tracking = params(learn_p)
function f_learn(x)
    adams_moulton_estimator(learn_p,x,t,ode_fun; return_estimated=true)
end
# Single layer neural network
function concat_neurons(args...)
    reestructure = zeros(eltype(args[1]), length(args[1][:,end]), length(args))
    for i in 1:length(args)
        reestructure[:,i] = args[i][:,end]
    end
    reestructure
end
a1 = adams_moulton_estimator(phi, ode_data[:,1:2], t[1:2], ode_fun; return_estimated=true)
a2 = adams_moulton_estimator(phi, ode_data[:,2:3], t[2:3], ode_fun; return_estimated=true)
a3 = adams_moulton_estimator(phi, ode_data[:,3:4], t[3:4], ode_fun; return_estimated=true)
con = concat_neurons(a1, a2, a3)
ode_data[:,2:4]
con - ode_data[:,2:4]

model_chain = Chain(x -> f_learn(x), x -> concat_neurons(x))
function loss_mse(x,y)
    res = Flux.mse(model_chain(x), y)
end

function loss_mse()
    proc_data = [ode_data[:,i:i+1] for i in 1:size(ode_data)[2]-1]
    segment = [f_learn(data_piece) for data_piece in proc_data]
    estimated_data = concat_neurons(segment...)
    res = Flux.mse(estimated_data, ode_data[:,2:end])
end

#data = Iterators.repeated([(ode_data[:,i:i+1], ode_data[:,2:end]) for i in 1:size(ode_data)[2]-1], 10000)
#data = [(ode_data[:,i:i+1], ode_data[:,2:end]) for i in 1:size(ode_data)[2]-1]
data = Iterators.repeated((), 1000)
function_error = []
function cb()
    err = Tracker.data(loss(ode_data, ode_data))
    if err < 100
        push!(function_error, err)
    end
end
@time Flux.train!(loss_mse, tracking, data, ADAM(0.1))
tracking
euclidean(learn_p, phi)
phi

# ----- Another test, adding layers  -----
i = 7
p_solve = problem_set[i]
ode_fun = p_solve.fun
phi = p_solve.phi
t = p_solve.t
ode_data = p_solve.data
t = p_solve.t
t = range(t[1], stop=t[end], length=10)
tspan = (t[1], t[end])
ini_cond = p_solve.data[:,1]
ode_prob = ODEProblem(p_solve.fun, ini_cond, tspan, p_solve.phi)
ode_sol  = solve(ode_prob, lsoda(), saveat=reduce(vcat, t))
data_original = reduce(hcat, ode_sol.u)
ode_data = copy(data_original)
plot(transpose(ode_data))
rand_ini = [rand(0.0001:0.0001:100) for i in 1:length(phi)]
learn_p = param(rand_ini)
#learn_p = param(phi)
#tracking = params(learn_p)
function f_learn(p)
    vec(adams_moulton_estimator(p, ode_data, t, ode_fun; return_estimated=true))
end
# 10 layer neural network
size(ode_data)
model_chain = Chain(Dense(length(ode_data),length(phi), σ),
                    Dense(length(phi),100, σ),
                    Dense(100,100, σ),
                    Dense(100,100, σ),
                    Dense(100,100, σ),
                    Dense(100,100, σ),
                    Dense(100,100, σ),
                    Dense(100,length(phi), abs),
                    x -> f_learn(x))
tracking = params(model_chain)
function loss(x,y)
    res = Flux.mse(model_chain(x), y)
    if res == NaN
        res = 10^10
    end
    return res
end
data = Iterators.repeated((vec(ode_data), vec(ode_data)), 1000000)
function_error = []
function cb()
    println(Flux.data(loss(vec(ode_data), vec(ode_data))))
end
@time Flux.train!(loss, tracking, data, ADAM(0.1))
print(params(model_chain))
# ----- Simpler model -----
learn_p = param([7.71]) |> gpu
function f(x)
    x.*learn_p
end
# Single layer neural network
model_chain = Chain(x -> f(x)) |> gpu
tracking = params(learn_p)
loss(x,y) = sum(abs2, model_chain(x) - y)
rand_data = rand((1.0:0.1:100.0),500,500) |> gpu
label = rand_data.*2 |> gpu
data = Iterators.repeated((rand_data, label), 5000)
@time Flux.train!(loss, tracking, data, ADAM(0.1))
tracking

# ----- Other one -----
model_chain = Chain(Dense(100, 50, σ),
                    Dense(50,100), softmax) |> gpu
tracking = params(model_chain) |> gpu
loss(x,y) = sum(abs2, model_chain(x) - y)
rand_data = rand((1.0:0.1:100.0),100,100)
label = rand_data.*[2]
data = Iterators.repeated((rand_data, label), 500)
@time Flux.train!(loss, tracking, data, ADAM(0.1))
loss(rand_data, label)

# ----- CuArrays are fast as they get big -----
a = rand((1:0.001:1000),10000,10000)
b = cu(a)
@time a*a
@time b*b

max_stars_repo_name: atgmello/DEParamEst

max_issues_repo_name: atgmello/DEParamEst

max_forks_repo_name: atgmello/DEParamEst


text:
using DynamicHMCModels, LinearAlgebra, StatsFuns

ProjDir = @__DIR__
cd(ProjDir)

delim = ';'
df = CSV.read(joinpath("..", "..", "data", "chimpanzees.csv"), DataFrame; delim)
df.pulled_left = convert(Array{Int64}, df.pulled_left)
df.prosoc_left = convert(Array{Int64}, df.prosoc_left)
first(df, 5)

Base.@kwdef mutable struct Chimpanzees{Ty <: AbstractVector,
  Tx <: AbstractMatrix}
    "Observations."
    y::Ty
    "Covariates"
    x::Tx
    "Number of observations"
    N::Int
end

# Write a function to return properly dimensioned transformation.

function make_transformation(model::Chimpanzees)
  as( (β = as(Array, size(model.x, 2)), ) )
end

# Instantiate the model with data and inits.

N = size(df, 1)
x = hcat(ones(Int64, N), df[:, :prosoc_left]);
y = df[:, :pulled_left]
model = Chimpanzees(;y=y, x=x, N=N);

# Make the model callable with a single argument.

function (model::Chimpanzees)(θ)
    @unpack y, x, N = model   # extract the data
    @unpack β = θ  # works on the named tuple too
    ll = 0.0
    ll += sum(logpdf.(Normal(0, 10), β)) # a & bp
    ll += sum([loglikelihood(Binomial(1, logistic(dot(x[i, :], β))), [y[i]]) for i in 1:N])
    ll
end

println()
θ = (β = [1.0, 2.0],)
model(θ)
println()

# Wrap the problem with a transformation, then use Flux for the gradient.

P = TransformedLogDensity(make_transformation(model), model)
∇P = ADgradient(:ForwardDiff, P)
results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000)
posterior = P.transformation.(results.chain)

p = as_particles(posterior)
display(p)

stan_result = "
Iterations = 1:1000
Thinning interval = 1
Chains = 1,2,3,4
Samples per chain = 1000

Empirical Posterior Estimates:
      Mean        SD       Naive SE       MCSE      ESS
 a 0.05103234 0.12579086 0.0019889282 0.0035186307 1000
bp 0.55711212 0.18074275 0.0028577937 0.0040160451 1000

Quantiles:
       2.5%        25.0%       50.0%      75.0%      97.5%
 a -0.19755400 -0.029431425 0.05024655 0.12978825 0.30087758
bp  0.20803447  0.433720250 0.55340400 0.67960975 0.91466915
";

# End of `10/m10.2d.jl`


text:
"""
    SupportedMarksDistributions

Union type for the marks distributions supported by the package for the SEMPP models.
"""
SupportedMarksDistributions = Union{Type{Distributions.GeneralizedPareto}, Type{EGPD.EGPpower}}

"""
SEMPPExpKern(data, μ, ϕ, γ, markdens, ξ, α, β, κ)

Self-exciting marked point process with exponential kernel model type such that for all t, and all m there is an event in t with mark m with probability λ_g(t)f(m|t)) where :
```
λ_g(t) = μ + ϕ * ν(t)
```

with `ν(t) = \\sum_{t_k < t} (1 + δ*m_k) (\\exp(-γ(t-t_k))`

f(m|t) is markdens with scale `σ_t = β + α * ν(t)`

`t_k` are the timestamps of all events and `m_k` their marks.
"""
mutable struct SEMPPExpKern <: SEPP
    data::Union{MarkedTimeSeries, Nothing}
    μ::Real
    ϕ::Real
    γ::Real
    δ::Real
    impact_function::Function
    markdens::SupportedMarksDistributions
    ξ::Real
    α::Real
    β::Real
    κ::Real
    cov_mat::Union{Array{Float64, 2}, Nothing}

    function SEMPPExpKern(data::Union{MarkedTimeSeries, Nothing} = nothing; μ::Real = rand(), ϕ::Real = rand(), γ::Real = rand(), δ::Real = rand(), impact_function::Function = (d -> 1 + tanh(d)), markdens::SupportedMarksDistributions = Distributions.GeneralizedPareto, ξ::Real = rand(), α::Real = rand(), β::Real = rand(), κ::Real = rand(), cov_mat::Union{Array{Float64, 2}, Nothing} = nothing)
        if any((μ, ϕ, γ, α, β, δ, κ) .< 0)
            error("paramaters except for ξ must be positive or zero")
        else
            new(data, μ, ϕ, γ, δ, impact_function, markdens, ξ, α, β, κ, cov_mat)
        end
    end

end

params(sepp::SEMPPExpKern) = Dict(:μ => sepp.μ, :ϕ => sepp.ϕ, :γ => sepp.γ, :impact_function => sepp.impact_function, :δ => sepp.δ, :markdens => sepp.markdens, :ξ => sepp.ξ, :α => sepp.α, :β => sepp.β, :κ => sepp.κ)

"""
    lin_coeff_impact(sepp::SEMPPExpKern)::Real

Get the linear coefficient impact of the marks on the rate.
"""
lin_coeff_impact(sepp::SEMPPExpKern)::Real = sepp.δ

"""
    marks_scale_params(sepp::SEMPPExpKern)::Real

Get the paramaters weighing the volatility function into the scale parameter of the mark density: `\\sigma_t = \\beta + \\alpha \\nu(t)`.
"""
marks_scale_params(sepp::SEMPPExpKern)::Real = (sepp.α, sepp.β)

shape(sepp::SEMPPExpKern) = sepp.κ

decay(sepp::SEMPPExpKern) = sepp.ξ


text:
using PRLeg

PRLeg.update()


text:
function fcsc_fastmulti{T}(X::Array{T, 3}, dictSize::Int, filterSize::Int, optParams::FCSCOpts{T})
# Input:
#   X: HxWxN stack of images (padded already)
#
# Alg. parameters:
#   mu_max: 1e5
#   tau: [1.01 1.5] larger means faster convergence but sacrifices primal feasibility,
#        smaller means slower convergence

  const mu_max = 1e5;
  const tau = 1.01;

  const Dinit    = optParams.Dinit;
  const beta::T  = optParams.beta;
        mu_s::T  = optParams.mu_s;
        mu_t::T  = optParams.mu_t;
  const max_iter = optParams.max_iter;
  const tol      = optParams.tol;
  const DEBUG    = optParams.DEBUG;
  const FIXED_D  = optParams.FIXED_D;

  (convH, convW, numImgs) = size(X);
  convLen = convH*convW;
  assert(!FIXED_D || !isempty(Dinit), "Dinit must be provided when FIXED_D is enabled.");

  # Initialization
  if !FIXED_D;
    mask = ones(T, convH, convW, dictSize);
    mask[1:filterSize, 1:filterSize, :] = 0;
    mask = find(mask);

    if isempty(Dinit);
      Dinit = reshape(mod([convert(T, 1):convH*convW*dictSize], 13), (convH, convW, dictSize));
    else;
      Dinit = padarray(Dinit, (convH-size(Dinit, 1), convW-size(Dinit, 2)), convert(T, 0), "post");
    end
    d = Dinit;
    s = copy(d);
    s[mask] = 0;
    normalize!(s);
    Zt_fft = SharedArray(Complex{T}, (dictSize, convLen));
    ZtZ_fft = SharedArray(Complex{T}, (dictSize, dictSize, convLen));
  else;
    d = padarray(Dinit, (convH-size(Dinit, 1), convW-size(Dinit, 2)), convert(T, 0), "post");
    s = d;
  end
  lambda_s = zeros(T, convH, convW, dictSize);

  # Transform images to the frequency domain
  x_fft = reshape(fft2(X), convLen, numImgs);
  z = zeros(T, convLen*dictSize, numImgs);
  t = zeros(T, convLen*dictSize, numImgs);
  lambda_t = zeros(T, convLen*dictSize, numImgs);
  Dt_fft = zeros(Complex{T}, dictSize, convLen);

  workerIds = workers();
  numWorkers = min(convLen, length(workerIds));
  numPerWorker = int(convLen/numWorkers);
  inds = [{ (i-1)*numPerWorker+(1:numPerWorker) for i=1:numWorkers-1 }; { (numWorkers-1)*numPerWorker+1:convLen }];

  @printf("iter=000");
  if DEBUG;
    debugData = {"x"=>X,
                 "d"=>d, "s"=>s,
                 "lambda_s"=>lambda_s,
                 "mu_s"=>mu_s,
                 "z"=>z, "t"=>t,
                 "lambda_t"=>lambda_t,
                 "mu_t"=>mu_t,
                 "beta"=>beta,
                 "function"=>compute_lagrangian_dictonly};
    (prev_obj,rerr,spar) = compute_lagrangian(debugData);
    @printf(": obj=%.10f (rerr=%.10f, spar=%.10f)\n", prev_obj, rerr, spar);
  end

  iter = 0;
  converged = false;
  if !FIXED_D;
    learnData = LearningData(zeros(Complex{T}, dictSize, convLen), zeros(Complex{T}, dictSize, dictSize, convLen),
      convH, convW, filterSize);
    learnParams = LearningOpts(T); learnParams.mu_s = mu_s; learnParams.max_iter = 1; learnParams.tol = tol;
    if DEBUG;
      learnParams.DEBUG = DEBUG;
      learnParams.debugData = debugData;
    end
  end
  codeData = FastCodingData(zeros(Complex{T}, dictSize, numImgs, convLen), zeros(Complex{T},
    dictSize, dictSize, convLen), convH, convW);
  codeParams = FastCodingOpts(T); codeParams.beta = beta; codeParams.max_iter = 1; codeParams.tol = tol;
    codeParams.t_fft = zeros(Complex{T}, dictSize, numImgs, convLen);
    codeParams.lambda_t_fft = zeros(Complex{T}, dictSize, numImgs, convLen);
  if DEBUG;
    codeParams.DEBUG = DEBUG;
    codeParams.debugData = debugData;
  end
  while ~converged;
    iter += 1;

    @printf("\b\b\b\b\b\b\b\biter=%03d", iter);
    if DEBUG;
      debugData["z"] = z; debugData["t"] = t;
      debugData["lambda_t"] = lambda_t;
      debugData["mu_t"] = mu_t;
      debugData["d"] = d; debugData["s"] = s;
      debugData["lambda_s"] = lambda_s;
      debugData["mu_s"] = mu_s;
      (prev_obj,rerr,spar) = compute_lagrangian(debugData);
      @printf(": obj=%.10f (rerr=%.10f, spar=%.10f)\n", prev_obj, rerr, spar);
    end

    #
    # CODING
    #
    copy!(Dt_fft, reshape(fft2(d), convLen, dictSize)');
    for i=1:convLen;
      # codeData.DtD_fft[:, :, i] = Dt_fft[:, i]*Dt_fft[:, i]';
      codeData.DtD_fft[:, :, i] = broadcast(*, Dt_fft[:, i], Dt_fft[:, i]');
    end
    # NOTE: if numImgs gets too large, we can do this in mini-batches
    for j=1:numImgs;
      # codeData.Dtx_fft[:, j, :] = Dt_fft.*repmat(x_fft[:, j].', dictSize, 1);
      codeData.Dtx_fft[:, j, :] = broadcast(.*, Dt_fft, x_fft[:, j].');
      codeParams.t_fft[:, j, :] = reshape(fft2(reshape(t[:, j], convH, convW, dictSize)), convLen, dictSize).';
      codeParams.lambda_t_fft[:, j, :] = reshape(fft2(reshape(lambda_t[:, j], convH, convW, dictSize)), convLen, dictSize).';
    end
    codeParams.z = z;
    codeParams.t = t;
    codeParams.lambda_t = lambda_t;
    codeParams.mu_t = mu_t;
    (z, t) = fcsc_coding_fastmulti(codeData, codeParams);

    #
    # LEARNING
    #
    if !FIXED_D;
      fill!(learnData.ZtZ_fft, 0);
      fill!(learnData.Ztx_fft, 0);
      for j=1:numImgs;
        copy!(Zt_fft, reshape(fft2(reshape(z[:, j], convH, convW, dictSize)), convLen, dictSize)');
#        for i=1:convLen;
#          # learnData.ZtZ_fft[:, :, i] += Zt_fft[:, i]*Zt_fft[:, i]';
#          learnData.ZtZ_fft[:, :, i] += broadcast(*, Zt_fft[:, i], Zt_fft[:, i]');
#        end
        @sync begin
          for p=1:numWorkers;
            @async begin
              remotecall_wait(workerIds[p], outerprod!, ZtZ_fft, Zt_fft, inds[p]);
            end
          end
        end
        broadcast!(+, learnData.ZtZ_fft, learnData.ZtZ_fft, ZtZ_fft);
        # learnData.Ztx_fft += Zt_fft.*repmat(x_fft[:, j].', dictSize, 1);
        learnData.Ztx_fft += broadcast(.*, Zt_fft, x_fft[:, j].');
      end
      learnParams.Dinit = d;
      learnParams.s = s;
      learnParams.lambda_s = lambda_s;
      learnParams.mu_s = mu_s;
      (d, s) = fcsc_learning_multi(learnData, learnParams);
    end

    # Update the Lagrangian multiplier
    for i=1:convLen*dictSize*numImgs; lambda_t[i] += mu_t*(z[i]-t[i]); end
    if !FIXED_D;
      for i=1:convLen*dictSize; lambda_s[i] += mu_s*(d[i]-s[i]); end
    end
#    debugData["lambda_t"] = lambda_t;
#    debugData["lambda_s"] = lambda_s;
#    (obj,rerr,spar) = compute_lagrangian(debugData);
#    is_mono = (prev_obj>=obj) ? 'O' : 'X';
#    @printf("          L: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
#    prev_obj = obj;

    # Update penalty term
    mu_t = min(tau*mu_t, mu_max);
    mu_s = min(tau*mu_s, mu_max);
#    debugData["mu_t"] = mu_t;
#    debugData["mu_s"] = mu_s;
#    (obj,rerr,spar) = compute_lagrangian(debugData);
#    is_mono = (prev_obj>=obj) ? 'O' : 'X';
#    @printf("          M: %c obj=%.10f (rerr=%.10f, spar=%.10f)\n", is_mono, obj, rerr, spar);
#    prev_obj = obj;

    learnConverged = FIXED_D || mean((d-s).^2)<=tol;
    codeConverged = true;
    for j=1:numImgs;
      if mu_t>0 && mean((z[:, j]-t[:, j]).^2)>tol;
        codeConverged = false;
        break
      end
    end
    if iter>=max_iter || (learnConverged && codeConverged);
      converged = true;
    end
  end
  s = s[1:filterSize, 1:filterSize, :];

  return (s, z, t);
end

max_stars_repo_name: bkong/FCSC


text:
J = j[1].RB2.J;
β2 = j[1].RB2.x[4:7]
# β2 = []
function M_β(β::Vector{T}) where T<:Real
    E = genE(β)
    M = 4*transpose(E)*J*E
    if typeof(β[1]) == Float64
        Ms = real(sqrt(M))
        # C = cholesky(M)
        # println("L = ", norm(Ms - C.L))
        # println(norm(Ms - C.U))
        # Ms = (sqrt.(diag(M))).*Matrix{T}(I,size(M))
    else # ForwardDiff.Dual variable
        # C = cholesky(M)
        # Ms = C.U

        # Ms = (sqrt.(diag(M))).*Matrix{T}(I,size(M))

        # diag_Mfinal = diag(Mfinal);
        # vals_sq = zeros(T,size(Mfinal,1))
        # for i=1:size(Mfinal,1)
        #     vals_sq[i] = sqrt(diag_Mfinal[i])
        # end
        # Ms = diagm(vals_sq);

        # Ms = (Mfinal^(0.5))

        vals, vecs = eigen(M);
        vals_sq = zeros(T,size(vals))
        for i=1:length(vals)
            vals_sq[i] = sqrt(vals[i])
        end
        Ms = vecs*diagm(vals_sq)*(vecs\I(size(vecs,1)))
    end
    return M, Ms
end
# M_β(β2)
m = central_fdm(10,1)
temp1 = FiniteDifferences.jacobian(m, z->M_β(z)[2], β2)[1]
temp2 = ForwardDiff.jacobian(z->M_β(z)[1], β2)
temp1_nz = findall(x -> abs(x)>1e-6, temp1)
temp2_nz = findall(x -> abs(x)>1e-6, temp2)
norm(temp1 - temp2)

# function kronSum(A::Array{T,2},B::Array{T,2}) where T<:Real
#     out = kron(A, Matrix{T}(I,(size(B)))) + kron(Matrix{T}(I,(size(A))),B)
#     return out
# end

M_check, Mhalf_check = M_β(β2)
dM_half = inv(kronSum(permutedims(Mhalf_check), Mhalf_check))
norm(dM_half*temp2 - temp1)

## trying svd for square root
# clearconsole()
# using  FiniteDifferences
# using LinearAlgebra
# using GenericLinearAlgebra
# using ForwardDiff
# using Revise
# function sqrtM(x::Array{T}) where T<:Real
#     # M = assembleM(x,j)
#     M = genMatM(x,j[1].RB2)
#     M1 = M[4:7,4:7]#[5:7,5:7]
#     # E = genE(q)
#     # M = 4*transpose(E)*J*E
#     if typeof(M[1]) == Float64
#         Ms = real(sqrt(M1))
#         # println("Ms = ", Ms)
#         # F = GenericLinearAlgebra.svd!(M)
#         # Ms2 = permutedims(F.Vt)*diagm(sqrt.(F.S))*F.Vt
#         # println("Mserr = ", norm(Ms - Ms2))
#     else
#         F = GenericLinearAlgebra.svd!(M1)
#         Ms = permutedims(F.Vt)*diagm(sqrt.(F.S))*F.Vt
#     end
#     return M, Ms
# end
# function sqrtMq(q::Array{T}) where T<:Real
#     E = genE(q)
#     M = 4*transpose(E)*j[1].RB2.J*E
#     if typeof(M[1]) == Float64
#         out = real(sqrt(M))
#         # F = GenericLinearAlgebra.svd!(M)
#         # Ms = permutedims(F.Vt)*diagm(sqrt.(F.S))*F.Vt
#         # println("Mserr = ", norm(out - Ms))
#     else
#         F = GenericLinearAlgebra.svd!(M)
#         out = permutedims(F.Vt)*diagm(sqrt.(F.S))*F.Vt
#     end
#     return M, out
# end
# # x0 = [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.0, 0.0, 0.24999999999999997, 0.0, -0.43301270189221935, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.294735414289422];
# # q2 = [sqrt(2)/2; sqrt(2)/2;zeros(2)]
# x0 = [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.0; 0.0; 0.24999999999999997; 0.0; -0.43301270189221935; [1;zeros(3)]; zeros(7)];
# # x0 = [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.0; 0.0; 0.24999999999999997; 0.0; -0.43301270189221935; randomQuat(); zeros(7)];
# updateRigidBody!(j[1].RB2,x0)
# # x0 = randomQuat();
# _,tempMs = sqrtM(x0)
# # temp = rand(10)
# m = central_fdm(10,1)
# finJ = FiniteDifferences.jacobian(m, z->sqrtM(z)[2], x0)[1]
# fdJ = ForwardDiff.jacobian(z->sqrtM(z)[2],x0)
# finJ_nz = findall(x -> abs(x)>1e-6, finJ)
# fdJ_nz = findall( x-> abs(x)>0, fdJ)
# println()
# println(norm(fdJ - finJ))
#
# # q2 = randomQuat();
# q2 = (x0[18:21])
# sqrtMq(q2)
# finJ2 = FiniteDifferences.jacobian(m, z->sqrtMq(z)[2], q2)[1]
# fdJ2 = ForwardDiff.jacobian(z->sqrtMq(z)[2],q2)
# println(norm(fdJ2 - finJ2))
##
# J = j[1].RB2.J;
# function pinvM(q::Array{T}) where T<:Real
#     E = genE(q)
#     M = 4*transpose(E)*J*E
#     Ms = sqrt(diagm(M))
#     F = GenericLinearAlgebra.svd!(diagm(M))
#     out = pinv(diagm(M))
#     return out
# end
# ForwardDiff.jacobian(pinvM, temp)
## symbolic
# using SymPy
#
# u0, u1, u2, u3 = symbols("u0 u1 u2 u3")
#
# E = [u0 u1 u2 u3;
#     -u1 u0 u3 -u2;
#     -u2 -u3 u0 u1;
#     -u3 u2 -u1 u0];
#
# J = j[1].RB2.J;
# M = 4*permutedims(E)*J*E
# F = GenericLinearAlgebra.svd!(M)
# # Ms = (M)^(2)
# diffM = [M.diff((u0,1))[:] M.diff((u1,1))[:] M.diff((u2,1))[:] M.diff((u3,1))[:]]
#
# diffM_eval = zeros(size(diffM))
# for i=1:length(diffM)
#     diffM_eval[i] = diffM[i](1.0, 0.0, 0.0, 0.0)
# end
# diffM_eval
# M[1,1](1.0, 0.0, 0.0, 0.0)
##
# function invM(x)
#     M = diagm(x)
#     return inv(M)
# end
# r=rand(3)
# finJ_M = FiniteDifferences.jacobian(m,z->invM(z), r)[1]
# fdJ_M = ForwardDiff.jacobian(z->invM(z), r)
# println("err= ", norm(finJ_M - fdJ_M))
## Mdiff function
clearconsole()
function kronSum(A::Array{T,2},B::Array{T,2}) where T<:Real
    out = kron(A, Matrix{T}(I,(size(B)))) + kron(Matrix{T}(I,(size(A))),B)
    return out
end
# x0 = [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.0; 0.0; 0.24999999999999997; 0.0; -0.43301270189221935; [1;zeros(3)]; zeros(7)];
# x0 = [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.0; 0.0; 0.24999999999999997; 0.0; -0.43301270189221935; randomQuat(); zeros(7)];

function Mdiff(x::Vector{T}) where T<:Real
    M = assembleM(x,j)
    Ms = real(sqrt(M))
    M_inv = inv(M)
    Ms_inv = real(sqrt(M_inv))

    dM = ForwardDiff.jacobian(z -> assembleM(z,j),x)
    dM_inv = -(kron(permutedims(inv(M)),inv(M)))*ForwardDiff.jacobian(z->assembleM(z,j), x)
    # dM_inv = ForwardDiff.jacobian(z -> inv(assembleM(z,j)),x)
    dMs = inv(kronSum(permutedims(Ms), Ms))*dM
    dMs_inv = -kron(permutedims(Ms_inv),Ms_inv)*(dMs);
    return M, Ms, Ms_inv, dM, dMs, dMs_inv, M_inv, dM_inv
end
m = central_fdm(5,1);
# finJ_dMs = FiniteDifferences.jacobian(m, z->Mdiff(z)[2], x0)[1]
# finJ_dMs_inv = FiniteDifferences.jacobian(m, z->Mdiff(z)[3], x0)[1]
# println()
# println(norm(finJ_dMs - Mdiff(x0)[5]))
# println(norm(finJ_dMs_inv - Mdiff(x0)[6]))
#
# _,_,Ms_inv,_,dMs, fdJ_dMs_inv = Mdiff(x0);
# # fdJ_dMs_inv = -kron(permutedims(Ms_inv),Ms_inv)*(dMs);
# println(norm(finJ_dMs_inv - fdJ_dMs_inv))
##
clearconsole()
function t1Fn(x)
    A, b = Ab_VecOfMat(x,j)
    M = assembleM(x,j)
    M, _, _, _, _, _, M_inv, dM_inv = Mdiff(x);
    # M = genMatM(x,j[1].RB2)
    # M1 = M#[4:7,4:7]
    # F = svd(M);
    # M_inv = F.V*diagm(1./(F.S))*(permutedims(F.U))
    # println("err = ", norm(inv(M) - pinv(M)))
    # sleep(1000)
    # M, _, Ms_inv, _, _, dMs_inv = Mdiff(x)
    sz_A1 = size(A,1)
    # t1 = permutedims(A)*((A*(M\(permutedims(A))))\I(sz_A1))
    t1 = permutedims(A)*inv(A*inv(M)*permutedims(A))

    dA = ForwardDiff.jacobian(z->Ab_VecOfMat(z,j)[1],x);
    dA_tr = ForwardDiff.jacobian(z->permutedims(Ab_VecOfMat(z,j)[1]),x);

    Y = A*M_inv*permutedims(A); invY = inv(Y);
    dY1 = kron(A*M_inv,1.0I(sz_A1))*dA;
    dY2 = kron(A,A)*dM_inv;
    dY3 = kron(1.0I(sz_A1), A*M_inv)*dA_tr;
    dY = dY1 + dY2 + dY3;
    dY_inv = -kron(permutedims(invY),invY)*dY;

    dt1_1 = kron(permutedims(invY),1.0I(size(A,2)))*dA_tr;
    dt1_2 = kron(1.0I(sz_A1),permutedims(A))*dY_inv;
    dt1 = dt1_1 + dt1_2;
    return t1, dt1
end
function Gdiff(x::Vector{T}) where T<:Real
    A, b = Ab_VecOfMat(x,j)
    M, _, Ms_inv, _, _, dMs_inv, _, _ = Mdiff(x)

    G = A*Ms_inv;
    Gp = permutedims(G)*((G*permutedims(G))\I(size(G,1)))

    t1, dt1 = t1Fn(x)
    dGp1 = kron(permutedims(t1),Matrix{T}(I,size(M)))*dMs_inv

    sz_t2 =size(t1,2)
    t2 = kron(Matrix{T}(I,(sz_t2,sz_t2)), Ms_inv)
    dGp2 = t2*dt1;

    dGp = dGp1 + dGp2
    return Gp, dGp
end
# Gdiff(x0)
finJ_G = FiniteDifferences.jacobian(m,z->Gdiff(z)[1],x0)[1]
fdJ_G = Gdiff(x0)[2]#ForwardDiff.jacobian(z->Gdiff(z)[2],x0)
println("\nerr_Gp", norm(finJ_G - fdJ_G))
#
# finJ_t1 = FiniteDifferences.jacobian(m,z->t1Fn(z)[1], x0)[1]
# fdJ_t1 = t1Fn(x0)[2]
# println("\nerr_t1= ", norm(finJ_t1 - fdJ_t1))

##
x0, u = getXU_0(j);

function hdiff(x)
    Fu = assembleF(x,u,j,GravityInInertial)
    A, b = Ab_VecOfMat(x,j)
    M = assembleM(x,j)

    a = inv(M)*Fu
    h = b-A*a
    return h
end
function Fcdiff(x::Vector{T}) where T<:Real
    Fu = assembleF(x,u,j,GravityInInertial)
    A, b = Ab_VecOfMat(x,j)
    M, Ms, _, _, dMs, _, _, _ = Mdiff(x)

    a = M\Fu
    h = b-A*a

    Gp, dGp = Gdiff(x)

    Fc = Ms*Gp*h

    dFc1 = kron(permutedims(Gp*h),I(size(M,1)))*dMs
    dFc2 = kron(permutedims(h),Ms)*dGp
    dFc3 = Ms*Gp*ForwardDiff.jacobian(z->hdiff(z),x)

    dFc = dFc1 + dFc2 + dFc3;

    dFu = ForwardDiff.jacobian(z->assembleF(z,u,j,GravityInInertial),x)
    return Fc, dFc, Fu, dFu
end
finJ_Fc = FiniteDifferences.jacobian(m,z->Fcdiff(z)[1],x0)[1]
fdJ_Fc = Fcdiff(x0)[2]
println("\nerr_Fc = ", norm(finJ_Fc - fdJ_Fc))

finJ_Fu = FiniteDifferences.jacobian(m,z->Fcdiff(z)[3],x0)[1]
fdJ_Fu = Fcdiff(x0)[4]
println("\nerr_Fu = ", norm(finJ_Fu - fdJ_Fu))

## fxdotdiff
invJ = j[1].RB2.invJ;
nJ = length(j); # also equal to number of bodies (excluding inertial frame)
function accdiff(x::Vector{T}) where T<:Real
    len_x = length(x)-14
    # dacc = zeros(T,(len_x,len_x))

    M = assembleM(x,j)
    Fc, dFc, Fu, dFu = Fcdiff(x)

    # dM_inv = ForwardDiff.jacobian(z ->inv(assembleM(z,j)), x)

    # dacc1 = kron(permutedims(Fc+Fu),1.0I(size(M,1)))*dM_inv
    # dacc2 = inv(M)*(dFc + dFu)
    #
    # dacc = dacc1 + dacc2
    dacc = zeros(T,(7*nJ,length(x)))
    for k=1:nJ
        ct1 = 7*(k-1); ct2 = 14*k;
        dacc[ct1+1:ct1+3,:] = 1/j[k].RB2.m*(dFc[ct1+1:ct1+3,:] +   dFu[ct1+1:ct1+3,:])

        E = genE(x[ct2+4:ct2+7])
        fdJ_E = ForwardDiff.jacobian(z->genE(z[ct2+4:ct2+7]),x)
        fdJ_trE = ForwardDiff.jacobian(z->transpose(genE(z[ct2+4:ct2+7])),x)

        dacc2_1 = kron(permutedims(invJ*E*((Fc+Fu)[ct1+4:ct1+7])),I(size(E,1)))*fdJ_trE
        dacc2_2 = kron(permutedims((Fc+Fu)[ct1+4:ct1+7]),permutedims(E)*invJ)*fdJ_E
        dacc2_3 = permutedims(E)*invJ*E*(dFc[ct1+4:ct1+7,:] + dFu[ct1+4:ct1+7,:])

        # println("\ndacc2_1 = ", dacc2_1)
        # println("\ndacc2_2 = ", dacc2_2)
        # println("\ndacc2_3 = ", dacc2_3)
        # sleep(1000);
        dacc[ct1+4:ct1+7,:] = 1/4*(dacc2_1 + dacc2_2 + dacc2_3)
    end
    return dacc
end
# accdiff(x0)
##
function fxdotdiff(x::Vector{T}) where T<:Real
    len_x = length(x)
    dFx = zeros(T,(len_x,len_x))
    # dFx = Matrix{T}(undef,(len_x,len_x))

    dacc = accdiff(x);
    # println("dacc = ", dacc)

    for k=1:length(j)
        dFx[14*k+1:14*k+7,14*k+8:14*k+14] = 1.0*I(7)
        dFx[14*k+8:14*k+14,:] = dacc[7*(k-1)+1:7*k,:]
    end
    return dFx
end

finJ_fx = FiniteDifferences.jacobian(m,z->fxdot(z,u,j,GravityInInertial),x0)[1];

fdJ_fx = fxdotdiff(x0);
# println("\nerr_fx = ", norm(finJ_fx[25:end,:] - fdJ_fx[25:end,:]))
println("\nerr_fx = ", norm(finJ_fx- fdJ_fx))
## M for β = [1.0;zeros(3)]
q = deepcopy(x0); q[18:21] = [1.0;zeros(3)]
function Mq(x::AbstractVector{T}) where T<:Real
    M = assembleM(x,j)
    M1 = M[4:7,4:7]
    out = (M1)#\I(size(M1,1))
    return out
end
finJ_M = FiniteDifferences.jacobian(m,z->inv(Mq(z)), q)[1]
M1 = Mq(q);
fdJ_M = -(kron(permutedims(inv(M1)),inv(M1)))*ForwardDiff.jacobian(Mq, q)
println("\nerr_invM = ", norm(finJ_M - fdJ_M))
## symbolic trying in login node (hprc)
# using SymPy
# using LinearAlgebra
# u0, u1, u2, u3 = symbols("u0 u1 u2 u3")
# uSym = [u0;u1;u2;u3];
# E = [u0 u1 u2 u3;
#     -u1 u0 u3 -u2;
#     -u2 -u3 u0 u1;
#     -u3 u2 -u1 u0];
#
# J = j[1].RB2.J;
# M = 4*permutedims(E)*J*E
# invM = (M)^(-1)s
## forwarddiff benchmarks
# function tempFn(z)
#     out = z.^2;
#     return out
# end
# y = rand(100);
# fdJ_fn = x->ForwardDiff.jacobian(tempFn,x)
# @btime fdJ_fn(y);
# @btime ForwardDiff.jacobian(tempFn,y);

jb1_nz = findall(x -> abs(x) > 1e-6, jb1)
jb2_nz = findall(x -> abs(x) > 0, jb2)

max_stars_repo_name: isrlab/Multibody-Dynamics

max_issues_repo_name: isrlab/Multibody-Dynamics

max_forks_repo_name: isrlab/Multibody-Dynamics


text:
#--------------------------------------------------------------
# MultilevelTT package
# example code for https://arxiv.org/abs/1802.09062
#--------------------------------------------------------------
# Auxiliary routines for evaluating representation condition numbers
# (or respective bounds for matrix representations)
module TTCondition

using TT

export repampvals, repamp, repampbound, repcondbound

function repampvals(t::Tensor, l::Int64)
  if l == 1
    T2 = unfold(t[2:end])
    return 1., norm(t[1][:]),
      maximum(svdvals(reshape(T2, (size(T2,1), size(T2,2)))))
  elseif l == length(t)
    T1 = unfold(t[1:(l-1)])
    return maximum(svdvals(reshape(T1, (size(T1,2), size(T1,3))))),
      norm(t[l][:]), 1.
  else
    T1 = unfold(t[1:(l-1)])
    T2 = unfold(t[(l+1):end])
    return maximum(svdvals(reshape(T1, (size(T1,2), size(T1,3))))),
      norm(t[l][:]),
      maximum(svdvals(reshape(T2, (size(T2,1), size(T2,2)))))
  end
end

# computes the exact value of the sum of ramp_l(t) over l
function repamp(t::Tensor)
  s = 0.;
  for l = 1:length(t)
    a, b, c = repampvals(t, l)
    s += a*b*c
  end
  return s
end

# computes the exact value of ramp_l(t)
function repamp(t::Tensor, l::Integer)
  a, b, c = repampvals(t, l)
  return a*b*c
end

function leftUnfold(T::TensorMatrix, l::Int64, i::Int64)
  s = [sizes(T)[i][j] for j = 1:2, i in eachindex(sizes(T))]
  s = s[:,1:(l-1)]
  t = Tensor(T)[1:(l-1)]
  t[l-1] = reshape(t[l-1][:,:,i], (size(t[l-1],1),size(t[l-1],2),1))
  M = unfold(t)[1,:,1]
  L = l-1
  dim = zeros(Int64,2L)
  dim[1:L] = (2L-1):-2:1
  dim[(L+1):(2L)] = (2L):-2:2
  return reshape(dimpermute(reshape(M, (s[:]...)), dim),
                (prod(s[1,:]), prod(s[2,:])))
end

function rightUnfold(T::TensorMatrix, l::Int64, i::Int64)
  s = [sizes(T)[i][j] for j = 1:2, i in eachindex(sizes(T))]
  s = s[:,(l+1):length(T)]
  t = Tensor(T)[(l+1):length(T)]
  t[1] = reshape(t[1][i,:,:], (1,size(t[1],2),size(t[1],3)))
  M = unfold(t)[1,:,1]
  L = length(T) - l
  dim = zeros(Int64,2L);
  dim[1:L] = (2L-1):-2:1;
  dim[(L+1):(2L)] = (2L):-2:2;
  return reshape(dimpermute(reshape(M, (s[:]...)), dim),
            (prod(s[1,:]), prod(s[2,:])))
end

function repampvals(T::TensorMatrix, l::Int64)
  cnorms = [norm(T[l][i,:,:,j]) for i=1:size(T[l],1), j=1:size(T[l],4)]
  if l == 1
    rnorms = [norm(rightUnfold(T, l, i)) for i = 1:size(T[l],4)]
    return [1.], cnorms, rnorms
  elseif l == length(T)
    lnorms = [norm(leftUnfold(T, l, i)) for i = 1:size(T[l],1) ]
    return lnorms, cnorms, [1.]
  else
    lnorms = [norm(leftUnfold(T, l, i)) for i = 1:size(T[l],1) ]
    rnorms = [norm(rightUnfold(T, l, i)) for i = 1:size(T[l],4) ]
    return lnorms, cnorms, rnorms
  end
end

function repampbound(T::TensorMatrix, l::Int64)
  a, b, c = repampvals(T, l)
  return norm(a)*vecnorm(b)*norm(c)
end

function repampbound(T::TensorMatrix)
  return sum(repampbound(T,l) for l = 1:length(T))
end

function repcondbound(T::TensorMatrix)
  return repampbound(T) / minimum(svdvals(unfold(T)))
end

end


text:
#=
    The File contains `Densik<T>` and related funcs implementation

    `Densik` implements `AbstractSparseObject` interface
=#

#------------------------------------------------------------------------------

include("sparsik.jl")

#------------------------------------------------------------------------------

"""
Densik
"""
mutable struct Densik{T<:Field} <: AbstractSparseVector{T}
    field::T
    data::Vector{<:FieldElem}
end

#------------------------------------------------------------------------------

density(v::Densik) = error("O(n)")
nnz(v::Densik) = error("O(n)")

base_ring(v::Densik) = v.field

#------------------------------------------------------------------------------

function Base.deepcopy_internal(x::Densik, stackdict::IdDict)
    y = Densik(x.field, Base.deepcopy_internal(x.data, stackdict))
    stackdict[x] = y
    return y
end

#------------------------------------------------------------------------------

first_nonzero(v::Densik) = error("O(n)")

#------------------------------------------------------------------------------

scale!(v::Densik, c) = error("O(n)")
scale(v::Densik, c) = error("O(n)")

#------------------------------------------------------------------------------

# returns v + c * u
# mutates v
# if vector `u` has k nonzeroes
# O(k)
function reduce!(v::Densik{T}, u::AbstractSparseObject{T}, c) where {T}
    for (idx, val) in u
        v[idx] += val * c
    end
    return v
end

# returns v + c * u
# if vector `u` has k nonzeroes
# O(k)
function Base.reduce(v::Densik{T}, u::AbstractSparseObject{T}, c) where {T}
    return reduce!(deepcopy(v), u, c)
end

#------------------------------------------------------------------------------

# returns standard inner product of v and u
# if vector `u` has k nonzeroes
# O(k)
function inner(v::Densik{T}, u) where {T}
    ans = zero(base_ring(v))

    for (idx, val) in u
        ans += val * v[idx]
    end

    return ans
end

inner(v, u::Densik{T}) where {T} = inner(u, v)

#------------------------------------------------------------------------------

Base.empty!(v::Densik) = error("O(n)")

#------------------------------------------------------------------------------

#------------------------------------------------------------------------------

Base.zero(v::Densik) = error("O(n)")
Base.iszero(v::Densik{T}) where {T} = error("O(n)")

Base.get(v::Densik, i::Int) = v.data[i]
Base.getindex(v::Densik, i::Int) = v.data[i]

Base.setindex!(v::Densik, val, i::Int) = v.data[i] = val

Base.size(v::Densik) = (length(v.data), )
Base.size(v::Densik, i::Int) = size(v)[i]
Base.length(v::Densik) = length(v.data)
dim(v::Densik) = length(v)

Base.isempty(v::Densik) =  error("O(n)")

# -----------------------------------------------------------------------------

==(v::Densik{T}, u::Densik{T}) where {T} = (v.data = u.data)
!=(v::Densik{T}, u::Densik{T}) where {T} = !(v == u)

# -----------------------------------------------------------------------------

Base.repr(::MIME"text/plain", A::Densik) = "($(join(map(idx -> A[idx], 1 : dim(A)), ", ")))"
Base.show(io::IO, A::Densik) = print(io, repr(MIME("text/plain"), A))

# -----------------------------------------------------------------------------

function Base.iterate(v::Densik, state=1)
    if state <= length(v)
        state += 1
        return ((state, v[state]), state)
    else
        return nothing
    end
end

# -----------------------------------------------------------------------------

Base.valtype(v::Densik) = elem_type(base_ring(v))
Base.eltype(v::Densik) = (Int, valtype(v))

#-----------------------------------------------------------------------------

# Returns a new Densik object consisting of elements
# from the given `v` each converted to the `field`
# Throws if any of the denominators reduces to zero
#
# O(n)
function modular_reduction(v::Densik, field)
    ans = Densik(zeros(dim(v)), field)

    for (i, x) in v.data
        y = modular_reduction(x, field)
        ans.data[i] = y
    end

    return ans
end

#-----------------------------------------------------------------------------

# returns a new Sparsik object consisting of elements
# from the `v` each reconstructed from v.field to QQ
#
# O†(k) if k is the number of nonzero in `v`
# † reconstruction cost not included
function rational_reconstruction(v::Densik)
    ans = Densik(zeros(dim(v)), QQ)

    tmp_type = BigInt
    if isa(base_ring(v), GaloisField)
        tmp_type = Int
    end
    ch = convert(tmp_type, characteristic(v.field))

    for (i, x) in v.data
        y = rational_reconstruction(convert(tmp_type, x), ch)
        ans.data[i] = y
    end

    return ans
end

#-----------------------------------------------------------------------------

#-----------------------------------------------------------------------------

function random_densik(sz::Int, field::T) where {T<:Field}
    if T <: FracField
        # todo!  ??????
        dense_repr = broadcast(QQ, rand(1 : 2^31-1, sz))
        return  Densik(field, dense_repr)
    elseif T <: Union{GaloisField, GaloisFmpzField}
        dense_repr = rand(field, sz)
        return Densik(field, dense_repr)
    end

    error("random sparsik over $field is not implemented!")
end

function random_densik(sz::Tuple{Int, Int}, field::T) where {T<:Field}
    random_densik(prod(sz), field)
end

#-----------------------------------------------------------------------------

# v[idx] += c
# mutates v
# O(1)
function reduce!(v::Densik{T}, idx::Int, c) where {T}
    v[idx] += c
    return v
end

#-----------------------------------------------------------------------------

# returns a new Densik object consisting of elements
# from the `v` each reconstructed from v.field to QQ
#
# O(n)
function rational_reconstruction(v::Densik)
    ans = Densik(QQ, zeros(QQ, size(v)...))
    tmp_type = BigInt
    if isa(base_ring(v), GaloisField)
        tmp_type = Int
    end
    ch = convert(tmp_type, characteristic(base_ring(v)))

    ans.data = [
        rational_reconstruction(convert(tmp_type, x), ch)
        for x in v.data
    ]
    return ans
end


text:
include("Schollwock.jl")
using LinearAlgebra
using Arpack
using TensorOperations

W = Dict{Int,AbstractArray}()
function W_constructor(l::Int;j=-1.0,λ = 2)
    I = [1. 0.
        0. 1.]
    Z = [1. 0.
       0. -1.]
   for k in 2:l-1
       W[k] = zeros(Float64,(2,2,3,3))
       W[k][:,:,1,1] = I
       W[k][:,:,3,3] = I
       W[k][:,:,2,1] = Z
       W[k][:,:,3,1] = λ*Z
       W[k][:,:,3,2] = j*Z
   end
   W[1] = zeros(Float64,(2,2,1,3))
   W[1][:,:,1,1] = λ*Z
   W[1][:,:,1,2] = j*Z
   W[1][:,:,1,3] = I
   W[l] = zeros(Float64,(2,2,3,1))
   W[l][:,:,1,1] = I
   W[l][:,:,2,1] = Z
   W[l][:,:,3,1] = λ*Z
end
Le = 12
W_constructor(Le)

function W_constructor2(l::Int;j=-1.0,λ=2.0)
    Z = [1. 0.
       0. -1.]
    Hss = λ*Z
    Hnn = j*kron(Z,Z)
    return NNHandSSHtoMPO(Hnn,Hss,l)
end

function W_constructor3(l::Int;j=-1.0,λ=2.0)
    II = [1. 0.
        0. 1.]
    Z = [1. 0.
       0. -1.]
    Hss = λ*Z
    Hnn = j*kron(Z,Z)
    return NNHtoMPO(Hnn+kron(Hss,II)+kron(II,Hss),l)
end

# display(W[1])
# display(W_constructor2(Le)[1])
# println(norm(W[4][:,:,3,2]- W_constructor2(Le)[4][:,:,3,2]))
# display(W[4][:,:,3,2]*W[4][:,:,2,1])
# display(W_constructor2(Le)[4][:,:,3,2]*W_constructor2(Le)[4][:,:,2,1])
# gr1 = DMRG(W,Le,maxiteration=20,D=5)#,ψ0 = ψ₀)
# res = deepcopy((gr.tensors[1])[1,:,:])
# for i in 2:Le
#     global res
#     res = res*(gr.tensors[i])[1,:,:]
# end
# println(res*res'/inner(gr,gr))
# res = deepcopy((gr.tensors[1])[2,:,:])
# for i in 2:Le
#     global res
#     res = res*(gr.tensors[i])[2,:,:]
# end
# println(res*res'/inner(gr,gr))

gr2 = DMRG(W_constructor2(Le,j=1),Le,ϵ=1e-20,maxiteration=100,D=4)
gr2p = DMRG(W_constructor2(Le,j=1),Le,ϵ=1e-20,maxiteration=100,D=4)
println("First:")
mydump(gr2)
println("Second")
mydump(gr2p)
println(abs2(inner(gr2,gr2p)))
# gr3 = DMRG(W_constructor3(Le),Le,maxiteration=20,D=5)
# println("correct answer")
# mydump(gr1)
# println("NNHandSSHtoMPO")
# mydump(gr2)
# println("This should do nothing:")
# gr2p = CutEntanglement(gr2,Le-3,after=false)
# mydump(gr2p)
# println("NNHtoMPO")
# mydump(gr3)

#GHZ state MPS:
# function GHZ_constructor(l::Int)::MPS2
#     A₀ = [1.0 0.0; 0.0 0.0]
#     A₁ = [0.0 0.0; 0.0 1.0]
#     t = Dict{Int, Array{Float64,3}}()
#     for i in 2:l-1
#         t[i] = zeros(Float64,(2,2,2))
#         t[i][1,:,:] = A₀
#         t[i][2,:,:] = A₁
#     end
#     t[1] = zeros(Float64,(2,1,2))
#     t[1][1,:,:] = [1.0 0.0]
#     t[1][2,:,:] = [0.0 1.0]
#     t[l] = zeros(Float64,(2,2,1))
#     t[l][1,:,:] = [1.0; 0.0]
#     t[l][2,:,:] = [0.0; 1.0]
#     return MPS2(t,l,2,l,2)
# end
#
# println("GHZ before cut")
# mydump(GHZ_constructor(13))
# cut = CutEntanglement(GHZ_constructor(13),10,after=true)
# println("GHZ after cut")
# mydump(cut)
#
# println(inner(GHZ_constructor(13),cut))

max_stars_repo_path: example/Ising model test.jl

max_stars_repo_name: CNOT/DMRG

max_issues_repo_path: example/Ising model test.jl

max_issues_repo_name: CNOT/DMRG

max_forks_repo_path: example/Ising model test.jl

max_forks_repo_name: CNOT/DMRG


text:
function gencodeelem2(filename::String, f, xdg, udg, odg, wdg, uinf, param, time)

foldername = "app";
opufile = "opu" * filename;
cpufile = "cpu" * filename;
gpufile = "gpu" * filename;

stropu = "template <typename T> void " * opufile;
strgpu = "template <typename T>  __global__  void kernel" * gpufile;

tmp = "(T *f, T *xdg, T *udg, T *odg, T *wdg, T *uinf, T *param, T time, int modelnumber, int ng, int nc, int ncu, int nd, int ncx, int nco, int ncw, int nce, int npe, int ne)\n";

stropu = stropu * tmp * "{\n";
stropu = stropu * "\tfor (int i = 0; i <ng; i++) {\n";

strgpu = strgpu * tmp * "{\n";
strgpu = strgpu * "\tint i = threadIdx.x + blockIdx.x * blockDim.x;\n";
strgpu = strgpu * "\twhile (i<ng) {\n";

str = "\t\tint j = i%npe;\n";
str = str * "\t\tint k = (i-j)/npe;\n";

fstr = string(f[:]);
str = varsassign(str, "param", length(param), 0, fstr);
str = varsassign(str, "uinf", length(uinf), 0, fstr);

varname = "xdg";
for i = 1:length(xdg)
    str1 = varname * string(i);
    str2 = varname * "[j+npe*" * string(i-1) * "+npe*ncx*k" * "]";
    str = str * "\t\tT " * str1 * " = " * str2 * ";\n";
end
varname = "udg";
for i = 1:length(udg)
    str1 = varname * string(i);
    str2 = varname * "[j+npe*" * string(i-1) * "+npe*nc*k" * "]";
    str = str * "\t\tT " * str1 * " = " * str2 * ";\n";
end
varname = "odg";
for i = 1:length(odg)
    str1 = varname * string(i);
    str2 = varname * "[j+npe*" * string(i-1) * "+npe*nco*k" * "]";
    str = str * "\t\tT " * str1 * " = " * str2 * ";\n";
end
varname = "wdg";
for i = 1:length(wdg)
    str1 = varname * string(i);
    str2 = varname * "[j+npe*" * string(i-1) * "+npe*ncw*k" * "]";
    str = str * "\t\tT " * str1 * " = " * str2 * ";\n";
end

n = length(f);
#for i=1:n
#    f[i] = sympy.simplify(f[i]);
#end
ts, fs = sympy.cse(f);
m = length(ts);
for i = 1:m
    str1 = sympy.ccode(ts[i][1]);
    str2 = sympy.ccode(ts[i][2]);
    str = str * "\t\tT " * str1 * " = " * str2 * ";\n";
end
for i = 1:n
    str1 = "\t\tf[j+npe*" * string(i-1) * "+npe*nce*k" * "]";
    str2 = sympy.ccode(fs[1][i]);
    str = str * str1 * " = " * str2 * ";\n";
end

stropu = stropu * str * "\t}\n" * "}\n\n";
tmp = "template void " * opufile;
tmp = tmp * "(double *, double *, double *, double *, double *, double *, double *, double, int, int, int, int, int, int, int, int, int, int, int);\n";
tmp = tmp * "template void " * opufile;
tmp = tmp * "(float *, float *, float *, float *, float *, float *, float *, float, int, int, int, int, int, int, int, int, int, int, int);\n";
stropu = stropu * tmp;

ioopu = open(foldername * "/" * opufile * ".cpp", "w");
write(ioopu, stropu);
close(ioopu);

strgpu = strgpu * str * "\t\ti += blockDim.x * gridDim.x;\n";
strgpu = strgpu * "\t}\n" * "}\n\n";
tmp = "template <typename T> void " * gpufile;
tmp = tmp * "(T *f, T *xdg, T *udg, T *odg, T *wdg, T *uinf, T *param, T time, int modelnumber, int ng, int nc, int ncu, int nd, int ncx, int nco, int ncw, int nce, int npe, int ne)\n";
tmp = tmp * "{\n";
tmp = tmp * "\tint blockDim = 256;\n";
tmp = tmp * "\tint gridDim = (ng + blockDim - 1) / blockDim;\n";
tmp = tmp * "\tgridDim = (gridDim>1024)? 1024 : gridDim;\n";
tmp = tmp * "\tkernel" * gpufile * "<<<gridDim, blockDim>>>(f, xdg, udg, odg, wdg, uinf, param, time, modelnumber, ng, nc, ncu, nd, ncx, nco, ncw, nce, npe, ne);\n";
tmp = tmp * "}\n\n";
tmp = tmp * "template void " * gpufile;
tmp = tmp * "(double *, double *, double *, double *, double *, double *, double *, double, int, int, int, int, int, int, int, int, int, int, int);\n";
tmp = tmp * "template void " * gpufile;
tmp = tmp * "(float *, float *, float *, float *, float *, float *, float *, float, int, int, int, int, int, int, int, int, int, int, int);\n";
strgpu = strgpu * tmp;

iogpu = open(foldername * "/" * gpufile * ".cu", "w");
write(iogpu, strgpu);
close(iogpu);

strcpu = replace(stropu, "opu" => "cpu");
strcpu = replace(strcpu, "for (int i = 0; i <ng; i++) {" => "#pragma omp parallel for\n\tfor (int i = 0; i <ng; i++) {");

iocpu = open(foldername * "/" * cpufile * ".cpp", "w");
write(iocpu, strcpu);
close(iocpu);

end

max_stars_repo_name: rloekvh/Exasim

max_issues_repo_name: rloekvh/Exasim

max_forks_repo_name: rloekvh/Exasim


text:
module TestBase

include("preamble.jl")
using Kaleido: prefer_singleton_callable, BijectionLens

lenses_as_shown = include("lenses_as_shown.jl") :: Array
desired_show = filter(x -> !occursin(x, "[]"),
                      readlines(joinpath(@__DIR__, "lenses_as_shown.jl")))
desired_show = strip.(desired_show)

struct ShowTestCase
    lens::Lens
    source::String
    desired_show::Bool
end

Base.show(io::IO, case::ShowTestCase) = print(io, case.source)

showtestcases = ShowTestCase.(lenses_as_shown, desired_show, Ref(true))

@testset "repr" begin
    @testset for case in showtestcases
        l0 = case.lens
        code = repr(l0)
        l1 = Base.include_string(@__MODULE__, code)
        @test l1 == l0
    end
end

@testset "show" begin
    @testset for case in showtestcases
        l = case.lens

        str1 = sprint(show, l)
        @debug """
        Show of $(case.source):
        $str1
        """
        @test occursin("Kaleido.", str1)

        str2 = sprint(show, l; context=:limit => true)
        @debug """
        Show of $(case.source):
        $str2
        """
        @test !occursin("Kaleido.", str2)
        if case.desired_show
            @test str2 == case.source
        end
    end
end

@testset "prefer_singleton_callable" begin
    @test sizeof((Int,)) > 0
    @test sizeof((prefer_singleton_callable(Int),)) == 0
    @test sizeof((prefer_singleton_callable(identity),)) == 0
end

end  # module


text:
function crra_test()
    @testset "CRRA" begin
        println("CRRA_test")
        for pSigma in [1.0, 2.0]
            crraS = CRRA(pSigma, true);
            cM = collect(range(2.0, 3.0, length = 4)) * collect(range(1.5, 2.5, length = 5))';
            uM = utility(crraS, cM);
            u5 = utility(crraS, cM[5]);
            @test u5 ≈ uM[5]
            muM = marginal_utility(crraS, cM);
            @test size(cM) == size(uM)
            @test all(muM .> 0)

            dc = 1e-5;
            u2M = utility(crraS, cM .+ dc);
            mu2M = (u2M .- uM) ./ dc;
            @test muM ≈ mu2M  atol = 1e-5
        end
    end
end

function pv_test()
    @testset "PV consumption" begin
        for pSigma in [1.0, 2.0]
            crraS = CRRA(pSigma, true)
            R = 1.07;
            beta = 0.97;
            betaR = beta * R;
            T = 3;
            pv = pv_consumption(crraS, beta, R, T);
            g = cons_growth(crraS, betaR);
            # println("Consumption growth $g")
            pvTrue = 1 + g/R + (g/R)^2;
            @test pv ≈ pvTrue

            ltIncome = 17.3;
            cV = cons_path(crraS, beta, R, T, ltIncome);
            @test length(cV) == T  &&  all(cV .> 0.0)
            cGrowthV = cV[2 : T] ./ cV[1 : (T-1)];
            @test all(abs.(cGrowthV .- g) .< 1e-4)
            discV = (1/R) .^ (0 : (T-1));
            pv = sum(cV .* discV);
            @test pv ≈ ltIncome

            devV = euler_deviation(crraS, beta, R, cV);
            @test all(abs.(devV) .< 1e-4)

            utilV = utility(crraS, cV);
            # uGrowth = util_growth(crraS, g);
            # uGrowthV = utilV[2 : T] ./ utilV[1 : (T-1)];
            # @test all(abs.(uGrowthV .- uGrowth) .< 1e-4)
            betaV = beta .^ (0 : (T-1));
            ltUtil = sum(utilV .* betaV);
            @test ltUtil ≈ lifetime_utility(crraS, beta, R, T, ltIncome)

            # Two ways of computing lifetime utility from utility
            util0 = -1.23;
            u1 = lifetime_utility(util0, beta, T);
            u2 = lifetime_utility(fill(util0, T), beta, T);
            @test u1 ≈ u2

            # Marginal utility of wealth
            ltUtil ≈ lifetime_utility(crraS, beta, R, T, ltIncome);
            dy = 1e-6;
            ltUtil2 = lifetime_utility(crraS, beta, R, T, ltIncome + dy);
            mu = (ltUtil2 - ltUtil) / dy;
            mu2 = UtilityFunctionsLH.mu_wealth(crraS, beta, R, T, ltIncome);
            @test isapprox(mu, mu2, atol = 1e-5)
        end
    end
end

function euler_test()
    @testset "Euler" begin
        for pSigma in [1.0, 2.0]
            crraS = CRRA(pSigma, true)
            R = 1.07;
            beta = 0.97;
            T = 5;
            g = cons_growth(crraS, beta .* R);
            devV = euler_deviation(crraS, beta, R, g .^ (1:T));
            @test all(abs.(devV) .< 1e-4)
        end
    end
end

@testset "CRRA" begin
    crra_test()
    pv_test()
    euler_test()
end

# --------------


text:
using FastAI
using Flux: mse, @functor, Dense, Descent, train!, Params, gradient, update!
using Flux.Data: DataLoader

"A simple DataBunch where `x` is random and `y = a*x + b` plus some noise."
function test_dbunch(a=2, b=3, bs=16, n_train=10, n_valid=2)

    function get_data(n)
        xy = []
        for i in 1:bs*n
            x = rand()
            y = a*x .+ b + 0.1*rand()
            push!(xy,([x],[y]))
        end
        return DataLoader(xy, batchsize=bs, shuffle=true)
    end

    train_dl = get_data(n_train)
    valid_dl = get_data(n_valid)
    return DataBunch(train_dl, valid_dl)
end

struct TestCallback <: AbstractCallback end

function test_learner(n_train=10, n_valid=2, cuda=false, lr=0.01)
    data = test_dbunch() #n_train=n_train,n_valid=n_valid)
    return Learner(data, Dense(1,1), loss=mse, opt=Descent(0.001))
end

function one_batch(dl::DataLoader)
    for d in dl
        return d
    end
end

function is_close(a,b;eps=1e-5)
    #@show a
    #@show b
    return (a-b)^2 < eps^2
end

function run_learner(learn)

    xys = learn |> data_bunch |> train |> one_batch
    lf = loss(learn)
    mf = model(learn)

    function lff(xy)
        x = xy[1]
        p = mf(x)[1]
        y = xy[2]
        return lf(p,y)
    end

    init_loss = sum(lff.(xys))
    fit!(learn,16)
    final_loss = sum(lff.(xys))
    @test final_loss < init_loss
end


text:
"""
A categorical distribution, P(x|parents(x)) where all parents are discrete integers 1:N.

The ordering of `distributions` array follows the convention in Decision Making Under Uncertainty.
Suppose a variable has three discrete parents. The first parental instantiation
assigns all parents to their first bin. The second will assign the first
parent (as defined in `parents`) to its second bin and the other parents
to their first bin. The sequence continues until all parents are instantiated
to their last bins.

This is equivalent to:

X,Y,Z
1,1,1
2,1,1
1,2,1
2,2,1
1,1,2
...
"""
struct CategoricalCPD{D} <: CPD{D}
    target::NodeName
    parents::NodeNames
    # list of instantiation counts for each parent, in same order as parents
    parental_ncategories::Vector{Int}
    # a vector of distributions in DMU order
    distributions::Vector{D}
end
CategoricalCPD(target::NodeName, d::D) where {D<:Distribution} =
        CategoricalCPD(target, NodeName[], Int[], D[d])

name(cpd::CategoricalCPD) = cpd.target
parents(cpd::CategoricalCPD) = cpd.parents
nparams(cpd::CategoricalCPD) = sum(d->paramcount(params(d)), cpd.distributions)

function (cpd::CategoricalCPD)(a::Assignment=Assignment())
    if isempty(cpd.parents)
        return first(cpd.distributions)
    else
        sub = [a[p] for p in cpd.parents]
        shape = ntuple(i -> cpd.parental_ncategories[i],
                length(cpd.parental_ncategories))
        ind = LinearIndices(shape)[sub...]
        return cpd.distributions[ind]
    end
end
(cpd::CategoricalCPD)(pair::Pair{NodeName}...) =
        (cpd)(Assignment(pair)) # cpd(:A=>1)

"""
    Distributions.ncategories(cpd::CategoricalCPD)

Return the number of categories for a cpd.
"""
Distributions.ncategories(cpd::CategoricalCPD) =
    ncategories(first(cpd.distributions))

function Distributions.fit(::Type{CategoricalCPD{D}},
                           data::DataFrame,
                           target::NodeName,
                           ) where {D}

    # no parents

    d = convert(D, fit(D, data[!,target]))
    CategoricalCPD(target, NodeName[], Int[], D[d])
end
function Distributions.fit(::Type{CategoricalCPD{D}},
                           data::DataFrame,
                           target::NodeName,
                           parents::NodeNames,
                           ) where {D}

    # with parents

    if isempty(parents)
        return fit(CategoricalCPD{D}, data, target)
    end

    # ---------------------
    # pull discrete dataset
    # 1st row is all of the data for the 1st parent
    # 2nd row is all of the data for the 2nd parent, etc.
    # calc parent_instantiation_counts

    nparents = length(parents)
    parental_ncategories = map!(p->infer_number_of_instantiations(data[!,p]), Array{Int}(undef, length(parents)), parents)
    dims = [1:parental_ncategories[i] for i in 1:nparents]
    distributions = Array{D}(undef, prod(parental_ncategories))
    for (q, parent_instantiation) in enumerate(Iterators.product(dims...))
        arr = Array{eltype(data[!,target])}(undef, 0)
        for i in 1 : nrow(data)
            if all(j->data[i,parents[j]]==parent_instantiation[j], 1:nparents) # parental instantiation matches
                push!(arr, data[i, target])
            end
        end
        distributions[q] = fit(D, arr)
    end

    CategoricalCPD(target, parents, parental_ncategories, distributions)
end

#####

const DiscreteCPD = CategoricalCPD{Categorical{Float64,Vector{Float64}}}

DiscreteCPD(target::NodeName, prob::AbstractVector{T}) where {T<:Real} = CategoricalCPD(target, Categorical(prob ./ sum(prob)))

function Distributions.fit(::Type{DiscreteCPD},
    data::DataFrame,
    target::NodeName;
    ncategories::Int = infer_number_of_instantiations(data[!,target]),
    )

    d = convert(Categorical{Float64,Vector{Float64}}, fit_mle(Categorical, ncategories, data[!,target]))
    CategoricalCPD(target, NodeName[], Int[], Categorical{Float64,Vector{Float64}}[d])
end
function Distributions.fit(::Type{DiscreteCPD},
    data::DataFrame,
    target::NodeName,
    parents::NodeNames;
    parental_ncategories::Vector{Int} = map!(p->infer_number_of_instantiations(data[!,p]), Array{Int}(undef, length(parents)), parents),
    target_ncategories::Int = infer_number_of_instantiations(data[!,target]),
    )

    # with parents

    if isempty(parents)
        return fit(DiscreteCPD, data, target, ncategories=target_ncategories)
    end

    nparents = length(parents)
    dims = [1:parental_ncategories[i] for i in 1:nparents]
    distributions = Array{Categorical{Float64,Vector{Float64}}}(undef, prod(parental_ncategories))
    arr = Array{eltype(data[!,target])}(undef, 0)
    for (q, parent_instantiation) in enumerate(Iterators.product(dims...))
        empty!(arr)
        for i in 1 : nrow(data)
            if all(j->data[i,parents[j]]==parent_instantiation[j], 1:nparents) # parental instantiation matches
                push!(arr, data[i, target])
            end
        end
        if !isempty(arr)
            distributions[q] = fit_mle(Categorical, target_ncategories, arr)
        else
            distributions[q] = Categorical(target_ncategories)
        end
    end

    CategoricalCPD(target, parents, parental_ncategories, distributions)
end


text:
println("reminder to Russell")
println("to install packages use")
println("launch julia using sudo")
println("try to import module first using, 'using'")
println("include('ml.jl')")
#println("Pkg.add('PyCall')")
println("\n\n\n\n\n\n")

function trypy()
	try
        ENV["PYTHON"] = "/anaconda3/bin/python"

        Pkg.test("PyCall")
        plt = pyimport("matplotlib.pyplot")
        x = range(0;stop=2*pi,length=1000); y = sin.(3*x + 4*cos.(2*x));
        plt.plot(x, y, color="red", linewidth=2.0, linestyle="--")
        plt.show()

		#return unreliable_connect()
	catch ex
        ENV["PYTHON"] = "/opt/conda/bin/python"
        Pkg.test("PyCall")
        math = pyimport("math")
        math.sin(math.pi / 4) # returns ≈ 1/√2 = 0.70710678...
	   end
end
#trypy()
#try
try
	using Knet
	using PyCall
    using Plots;
    using StatsPlots;
    using DataFrames
    using IJulia
	using PyPlot
	using Plots
	using GR

catch
    using Pkg
	#Pkg.rm("StatPlots")
	Pkg.add("StatsPlots")
	#Pkg.add("Gallium")
	Pkg.add(“GR”)
    Pkg.add("IJulia");
    Pkg.add("Plots");
    Pkg.add("StatsPlots"); #to install the StatPlots package.
    Pkg.add("DataFrames");
	Pkg.add("Seaborn")
	Pkg.add("PyPlot")
    using DataFrames
    using IJulia
	#using Gallium
	using Knet

    using Plots;
    using StatsPlots;
    using PyCall
end

#=

	include(Knet.dir("data","housing.jl"));
	x,y = housing();

	predict(w, x) = w[1] * x .+ w[2];
	loss(w, x, y) = mean(abs2, predict(w, x)-y);
	lossgradient = grad(loss);

	function train(w, data; lr=0.01)
	    for (x,y) in data
	        dw = lossgradient(w, x, y)
	        for i in 1:length(w)
	            w[i] -= dw[i]*lr
	        end
	    end
	    return w
	end;

	#plotly();
	gr();
	scatter(x', y[1,:], layout=(3,5), reg=true, size=(950,500))
	savefig("dnn_scatter_plot.png")
	w = Any[ 0.1*randn(1,13), 0.0 ];
	errdf = DataFrame(Epoch=1:20, Error=0.0);
	#BigInt cntr;

	println("try to import module first using, 'using'")
	#println("Pkg.add('PyCall')")
	println("\n\n\n\n\n\n")

	for i=1:200
		cntr = 1
		#println("cntr",cntr)

	    println(i,cntr)
	    train(w, [(x,y)])
	    if mod(i, 10) == 0
	        println("Epoch $i: $(round(loss(w,x,y)))")
	        errdf[cntr, :Epoch]=i
	        errdf[cntr, :Error]=loss(w,x,y)
	        cntr+=1
	    end
	end;
	print(errdf)

=#

@pyimport pickle
f = pybuiltin("open")("scraped_new.p","rb")
p = pickle.Unpickler(f)
scraped_new = p[:load]()
f[:close]()
scraped_new[1]["wcount"]
sn = [ sn["wcount"] for sn in scraped_new ]
#varinfo(sn)
histogram(sn)
#using PyPlot
#h = PyPlot.plt.hist(sn)
savefig("document_length_distribution.png")
varinfo()

#exit();

max_stars_repo_name: mcgurrgurr/ScienceAccessibility

max_issues_repo_name: mcgurrgurr/ScienceAccessibility

max_forks_repo_name: mcgurrgurr/ScienceAccessibility


text:
struct FaceQuadratures
    quads::Any
    facetoquad::Any
    ncells::Any
    nfaces::Any
    function FaceQuadratures(quads, facetoquad)
        nphase, nfaces, ncells = size(facetoquad)
        @assert nfaces == 4
        @assert nphase == 2
        @assert all(facetoquad .>= 0)
        @assert all(facetoquad .<= length(quads))
        new(quads, facetoquad, ncells, nfaces)
    end
end

function FaceQuadratures(cutmesh, levelset, numqp)

    ncells = number_of_cells(cutmesh)

    quads = face_quadratures(numqp)
    quad1d = ReferenceQuadratureRule(numqp)
    facetoquad = zeros(Int, 2, 4, ncells)

    for cellid = 1:ncells
        s = cell_sign(cutmesh,cellid)
        if s == +1
            facetoquad[1, :, cellid] .= 1:4
        elseif s == -1
            facetoquad[2, :, cellid] .= 1:4
        elseif s == 0
            nodeids = nodal_connectivity(background_mesh(cutmesh),cellid)
            load_coefficients!(levelset,cellid)

            pquad = face_quadratures(interpolater(levelset), +1, quad1d)
            idxstart = length(quads) + 1
            append!(quads, pquad)
            facetoquad[1, :, cellid] .= idxstart:(idxstart+3)

            nquad = face_quadratures(interpolater(levelset), -1, quad1d)
            idxstart = length(quads) + 1
            append!(quads, nquad)
            facetoquad[2, :, cellid] .= idxstart:(idxstart+3)
        else
            error("Expected cellsign ∈ {-1,0,+1}, got cellsign = $s")
        end
    end
    return FaceQuadratures(quads, facetoquad)
end

function Base.getindex(facequads::FaceQuadratures, s, faceid, cellid)
    phaseid = cell_sign_to_row(s)
    return facequads.quads[facequads.facetoquad[phaseid, faceid, cellid]]
end

function uniform_face_quadratures(facequads::FaceQuadratures)
    return facequads.quads[1:facequads.nfaces]
end

function number_of_faces_per_cell(facequads::FaceQuadratures)
    return facequads.nfaces
end

function Base.show(io::IO, facequads::FaceQuadratures)
    ncells = facequads.ncells
    nuniquefacequads = length(facequads.quads)
    str = "FaceQuadratures\n\tNum. Cells: $ncells\n\tNum. Unique Quadratures: $nuniquefacequads"
    print(io, str)
end

function face_quadratures(numqp)
    quad1d = tensor_product_quadrature(1, numqp)
    facequads = [extend_to_face(quad1d, faceid) for faceid = 1:4]
    return facequads
end

function face_quadrature(faceid, levelset, signcondition, quad1d)
    dir, coordval = reference_face(faceid)
    quad = QuadratureRule(ImplicitDomainQuadrature.one_dimensional_quadrature(
        [x -> levelset(extend(x, dir, coordval))],
        [signcondition],
        -1.0,
        +1.0,
        quad1d,
    ))
    return extend_to_face(quad, faceid)
end

function face_quadratures(levelset, signcondition, quad1d)
    quads = [face_quadrature(faceid, levelset, signcondition, quad1d) for faceid = 1:4]
    return quads
end

function extend_to_face(quad::QuadratureRule, faceid)
    extp = extend_to_face(quad.points, faceid)
    return QuadratureRule(extp, quad.weights)
end

function extend_to_face(points, faceid)
    dir, coordval = reference_face(faceid)
    @assert dir == 1 || dir == 2
    flipdir = dir == 1 ? 2 : 1
    return extend([coordval], flipdir, points)
end

function reference_face(faceid)
    if faceid == 1
        return (2, -1.0)
    elseif faceid == 2
        return (1, +1.0)
    elseif faceid == 3
        return (2, +1.0)
    elseif faceid == 4
        return (1, -1.0)
    else
        error("Expected faceid ∈ {1,2,3,4}, got faceid = $faceid")
    end
end

function update_face_quadrature!(facequads::FaceQuadratures, s, faceid, cellid, quad)
    row = cell_sign_to_row(s)
    idx = facequads.facetoquad[row, faceid, cellid]
    @assert idx > 4 "Attempting to modify the uniform cell face quadrature rules"
    facequads.quads[idx] = quad
end


text:
#export check_vector_input
export colspace, eye, ewprod, khatrirao, krontkron, kron, krontv, krtv, tkrtv, perfect_shuffle, lanczos, lanczos_tridiag, randsvd
export VectorCell, MatrixCell, TensorCell, CoreCell

"""
    VectorCell(undef,N)

Cell of vectors of length N.
"""
const VectorCell = Array{AbstractVector{<:Number},1}
"""
    MatrixCell(undef,N)

Cell of matrices of length N.
"""
const MatrixCell = Array{AbstractMatrix{<:Number},1}
"""
    TensorCell(undef,N)

Cell of multidimensional arrays of length N.
"""
const TensorCell = Array{AbstractArray{<:Number},1}
export CoreCell

"""
    CoreCell(undef,N)

Cell of 3D tensors of length N. Suitable for the cores of TTtensor.
"""
const CoreCell = Array{AbstractArray{<:Number,3},1}

# """
#     ewprod(I1::CartesianIndex,I2::CartesianIndex)
#
# Element-wise product of two inputs of type CartesianIndex.
# """
function ewprod(I1::CartesianIndex,I2::CartesianIndex)
  N=length(I1)
  @assert(length(I2)==N,"Dimension mismatch.")
  prod=zeros(Int,N)
  [prod[n]=I1[n]*I2[n] for n=1:N]
  CartesianIndex{N}(tuple(prod...))
end

"""
    check_vector_input(input,dim,default_value)

Check whether the input vector is of appropriate size or if input is number create vector out of it.
"""
function check_vector_input(input,dim::Integer,default_value::Number)
  if length(input)>0
     if isa(input,Number)
        input=repeat([input],dim)
    else
      @assert(dim==length(input),"Dimensions mismatch.")
    end
  else #fixed-precision problem
    input=repeat([default_value],dim)
  end
  input
end

"""

    colspace(X; <keyword arguments>)

Column space basis.
## Arguments:
- `X`: Matrix.
- `method` ∈ {"svd","lanczos","randsvd"} Method for SVD. Default: "svd".
- `maxrank::Integer`: Maximal rank. Optional.
- `atol::Number`: Drop singular values below atol.  Default: 1e-8.
- `rtol::Number`: Drop singular values below rtol*sigma_1. Optional.
- `p::Integer`: Oversampling parameter used by lanczos and randsvd methods. Defaul p=10.
"""
function colspace(X::AbstractMatrix{<:Number};method="svd",maxrank=0,atol=1e-8,rtol=0,p=10)
  if method == "lanczos"
    U,S=lanczos(X,tol=atol,reqrank=maxrank,p=p)
  elseif method == "randsvd"
    U,S=randsvd(X,tol=atol,reqrank=maxrank,p=p)
  else
    U,S,V=svd(X)
  end
  if maxrank!=0 && size(U,2)>maxrank
    U=U[:,1:maxrank]
    if maxrank<length(S)
      S=S[1:maxrank]
    end
  end
  rtol != 0 ? tol=rtol*S[1] : tol=atol
  K=findall(x-> x>tol ? true : false,S)
  U[:,K]
end

"""
    eye(n::Integer[,full=1])

Identity matrix of size nxn. If full=0, returns type Diagonal.
"""
function eye(n::Integer,full=1)
  full==1 ? Matrix(1.0I, n, n) : Diagonal(ones(n))
end

"""
    khatrirao(M,t='n')
    khatrirao(M,n,t='n')

Khatri-Rao product of matrices from M:  M₁⊙ ⋯ ⊙ Mₙ. Optionally skip nth matrix.
If t='t', compute transpose Khatri-Rao product: M₁⊙ᵀ ⋯ ⊙ᵀ Mₙ.
"""
function khatrirao(M::MatrixCell,t='n')
    N=length(M)
    if t== 't'
      sz=[size(M[n],1) for n=1:N]
      @assert(any(map(Bool,sz .-size(M[1],1)))==0,"Matrices must have the same number of rows.")
      [M[n]=M[n]' for n=1:N]
      khatrirao(M)'
    else
      K=[size(M[n],2) for n=1:N]
      @assert(any(map(Bool,K .-size(M[1],2)))==0,"Matrices must have the same number of columns.")
      J,K=size(M[end])
      X=reshape(M[end],J, 1, K)
      for n=N-1:-1:1
        sz=size(M[n],1)
        Y=reshape(M[n],1,sz,K)
        X=reshape(Y.*X,sz*J,1,K)
        J=sz*J
      end
      reshape(X,size(X,1),K)
    end
end
khatrirao(M::Array{Matrix{<:Number}},t='n')=khatrirao(MatrixCell(M),t)
function khatrirao(M1::AbstractMatrix{<:Number}, M2::AbstractMatrix{<:Number},t='n')
  M=MatrixCell(undef,2)
  M[1]=M1;M[2]=M2
  khatrirao(M,t)
end
function khatrirao(M1::MatrixCell, M2::MatrixCell,t='n')
    @assert(length(M1)==length(M2),"Matrix cells must be of the same length.")
    N=length(M1)
    M=MatrixCell(undef,N)
    for n=1:N
      M[n]=khatrirao(M1[n],M2[n],t)
    end
    M
end
#Skip nth matrix from M
function khatrirao(M::MatrixCell,n::Integer,t='n')
    khatrirao(deleteat!(M,n),t)
end
khatrirao(M::Array{Matrix{<:Number}},n::Integer,t='n') =khatrirao(MatrixCell(M),n,t)

#Extension of Base.kron to work with MatrixCell
function kron(M::MatrixCell,t='n')
  N=length(M)
  t=='n' ? K=M[1] : K=M[1]'
  for n=2:N
      t=='n' ? K=kron(K,M[n]) : K=kron(K,M[n]')
  end
  K
end
#If array of matrices isn't defined as MatrixCell, but as M=[M1,M2,...,Mn]:
kron(M::Array{AbstractMatrix{<:Number}},t='n')=kron(MatrixCell(M),t)

"""
    krontkron(A,v,t='n')

Kronecker product of matrices multiplied by Kronecker product of vectors.
"""
function krontkron(A::MatrixCell,v::VectorCell,t='n')
  if t=='t'
      [A[n]=A[n]' for n=1:length(A)]
  end
  @assert(length(A)==length(v),"There must be the same number of matrices and vectors.")
  N=length(A)
  @assert(!any(map(Bool,[size(A[n],2)-length(v[n]) for n=1:N])),"Dimension mismatch.")
  w=A[1]*v[1]
  [w=kron(w,A[n]*v[n]) for n=2:N]
  w
end
#If arrays are not defined as VectorCell/MatrixCell, but as [M1,M2,...,Mn]:
krontkron(A::Array{Matrix{<:Number}},v::Array{Vector{<:Number}},t='n')=krontkron(MatrixCell(A),VectorCell(v),t)
krontkron(A::Array{Matrix{<:Number}},v::VectorCell,t='n')=krontkron(MatrixCell(A),v,t)
krontkron(A::MatrixCell,v::Array{Vector{<:Number}},t='n')=krontkron(A,VectorCell(v),t)

"""
    krontv(A,B,v)

Kronecker product times vector: (A ⊗ B)v.
If v is a matrix, multiply column by column.
"""
function krontv(A::AbstractMatrix{<:Number},B::AbstractMatrix{<:Number},v::AbstractVector{<:Number})
  m,n=size(A);
  p,q=size(B);
  @assert(length(v)==q*n, "Dimensions mismatch.")
  if n*q*p+m*n*p <= m*q*p+q*n*m
    vec(B*reshape(v,q,n)*A');
  else
    vec(B*(reshape(v,q,n)*A'));
  end
end
function krontv(A::AbstractMatrix{<:Number},B::AbstractMatrix{<:Number},M::AbstractMatrix{<:Number})
  if sort(collect(size(vec(M))))[1]==1
        return krontv(A,B,vec(M))
  end
  m,n=size(A);
  p,q=size(B);
  @assert(size(M,1)==q*n, "Dimensions mismatch.");
  Mprod=zeros(m*p,size(M,2))
  for j=1:size(M,2)
    Mprod[:,j]=krontv(A,B,M[:,j])
  end
  Mprod
end

"""
    krtv(A,B,v)

Khatri-Rao product times vector: (A ⊙ B)v.
If v is a matrix, multiply column by column.
"""
function krtv(A::AbstractMatrix{<:Number},B::AbstractMatrix{<:Number},v::AbstractVector{<:Number})
  @assert(size(A,2)==size(B,2),"Dimension mismatch.")
  m,n=size(A);
  p=size(B,1);
  @assert(length(v)==n, "Dimensions mismatch.")
  if n<=p
    vec(B*(v.*A'));
  else
    vec((B.*v')*A');
  end
end
function krtv(A::AbstractMatrix{<:Number},B::AbstractMatrix{<:Number},M::AbstractMatrix{<:Number})
  @assert(size(A,2)==size(B,2),"Dimension mismatch.")
  if sort(collect(size(vec(M))))[1]==1
    return krtv(A,B,vec(M));
  end
  m,n=size(A);
  p=size(B,1);
  @assert(size(M,1)==n, "Dimensions mismatch.");
  Mprod=zeros(m*p,size(M,2))
  for j=1:size(M,2)
    Mprod[:,j]=krtv(A,B,M[:,j])
  end
  Mprod
end

"""
   tkrtv(A,B,v)

Transpose Khatri-Rao product times vector: (A ⊙ᵀ B)v.
If v is a matrix, multiply column by column.
"""
function tkrtv(A::AbstractMatrix{<:Number},B::AbstractMatrix{<:Number},v::AbstractVector{<:Number})
  @assert(size(A,1)==size(B,1),"Dimension mismatch.")
  m,n=size(A)
  p=size(B,2)
  @assert(length(v)==n*p, "Dimensions mismatch.")
  if n<=p
    sum((B*reshape(v,p,n)).*A,dims=2)
  else
    sum(B.*(A*reshape(v,p,n)'),dims=2)
  end
end
function tkrtv(A::AbstractMatrix{<:Number},B::AbstractMatrix{<:Number},M::AbstractMatrix{<:Number})
  @assert(size(A,1)==size(B,1),"Dimension mismatch.")
  if sort(collect(size(vec(M))))[1]==1
    return tkrtv(A,B,vec(M))
  end
  m,n=size(A);
  p=size(B,2);
  @assert(size(M,1)==n*p, "Dimensions mismatch.");
  Mprod=zeros(m,size(M,2))
  for j=1:size(M,2)
    Mprod[:,j]=tkrtv(A,B,M[:,j])
  end
  Mprod
end

"""
    perfect_shuffle(piles_nmbr,piles_len,piles_block=0)

Perfect shuffle vector.
"""
function perfect_shuffle(piles_nmbr::Integer,piles_len::Integer,piles_block=[])
    if piles_block==[]
        p=zeros(Int,piles_nmbr*piles_len)
        let
            i=1
            for l=1:piles_len, n=1:piles_nmbr
                p[i]=(n-1)*piles_len+l
                i+=1
            end
        end
    else
        @assert(length(piles_block)==piles_nmbr,"Dimension mismatch.")
        #p=zeros(Int,piles_nmbr*piles_len*piles_block)
        p=zeros(Int,1,0)
        for l=1:piles_len, n=1:piles_nmbr
            p=[p sum(piles_block[1:n-1])*piles_len .+ collect((l-1)*piles_block[n]+1:l*piles_block[n])...]
        end
    end
    p[:]
end

"""
    lanczos(A; <keyword arguments>)

Lanczos based SVD - computes left singular vectors and singular values of a matrix.

## Arguments:
- `A::Matrix`
- `tol`: Tolerance - discard singular values below tol. Default: 1e-8.
- `maxit`: Maximal number of iterations. Default: 1000.
- `reqrank`: Number of singular values and singular vectors to compute. Optional.
- `p`: Oversampling parameter. Default: p=10.
"""
function lanczos(A::AbstractMatrix{<:Number};tol=1e-8,maxit=1000,reqrank=0,p=10)
  Q,T=lanczos_tridiag(A,tol=tol,maxit=maxit,reqrank=reqrank,p=p)
  E=eigen(T,tol,Inf);
  U=E.vectors[:,end:-1:1];
  S=sqrt.(abs.(E.values[end:-1:1]));
  if reqrank!=0
    U=Q*U[:,1:reqrank];
    S=S[1:reqrank];
  else
    K=findall(x-> x>tol ? true : false,S)
    U=Q*U[:,K];
    S=S[K];
  end
  U,S
end

"""
    lanczos_tridiag(A; <keyword arguments>)

Lanczos tridiagonalization algorithm - returns orthonormal Q and symmetric tridiagonal T such that A≈Q*T*Q'.

## Arguments:
- `A::Matrix`
- `tol`: Tolerance. Default: 1e-8.
- `maxit`: Maximal number of iterations. Default: 1000.
- `reqrank`: Number of singular values and singular vectors to compute. Optional.
- `p`: Oversampling parameter. Default: p=10.
"""
function lanczos_tridiag(A::AbstractMatrix{<:Number};tol=1e-8,maxit=1000,reqrank=0,p=10)
  m,n=size(A)
  K=min(m,maxit);
  if reqrank != 0
      K=min(reqrank+p,K);
  end
  isreal(A) ? α=zeros(K) : α=zeros(Complex,K)
  β=zeros(K)
  v=randn(m);
  q=v/norm(v);
  Q=zeros(m,1)
  Q[:,1]=q;
  k=0; #needed if stopping criterion is met
  for k=1:K
    r=A*(A'*Q[:,k])
    α[k]=dot(r,Q[:,k])
    r=r-α[k]*Q[:,k] #without orthogonalization: r=r-α[k]*Q[:,k]-β[k-1]*Q[:,k-1]
    [r=r-Q*(Q'*r) for i=1:3]
    β[k]=norm(r)
    if β[k] < tol
      K=k
      break
    end
    if k!=K
      Q=[Q r/β[k]]
    end
  end
  isreal(α) ? T=SymTridiagonal(α[1:K], β[1:K-1]) : T=SymTridiagonal(real(α[1:K]), β[1:K-1]) +im*Diagonal(imag(α[1:K])) #T=Hermitian(Matrix(SymTridiagonal(real(α[1:K]), β[1:K-1]) +im*Diagonal(imag(α[1:K]))))
  Q,T
  end

function randrange(A::AbstractMatrix{<:Number},gram=true,t='t';tol=1e-8,maxit=1000,reqrank=0,r=10,p=10)
  m,n=size(A)
  if reqrank!=0
    if t=='t'
      Y=A'*randn(m,reqrank+p)
      if gram
        Y=A*Y
      end
    elseif t=='n'
      Y=A*randn(n,reqrank+p)
      if gram
        Y=A'*Y
      end
    else
      error("Wrong input.")
    end
    Q=Matrix(qr(Y).Q)
  else
    maxit=min(m,n,maxit);
    rangetol=tol*sqrt.(pi/2)/10; #assures ||A-Q*Q'*A||<=tol
    if t=='t'
      Y=A'*randn(m,r)
      if gram
        Y=A*Y
        Q=zeros(m,0)
      else
        Q=zeros(n,0)
      end
    elseif t=='n'
      Y=A*randn(n,r)
      if gram
        Y=A'*Y
        Q=zeros(n,0)
      else
        Q=zeros(m,0)
      end
    end

    j=0
    maxcolnorm=maximum([norm(Y[:,i]) for i=1:r]);
    while maxcolnorm > rangetol && j<maxit
      j+=1;
      v=Q'*Y[:,j];
      Y[:,j]-=Q*v;
      q=Y[:,j]/norm(Y[:,j]);
      Q=[Q q];
      if t=='t'
        w=A'*randn(m)
        if gram
          w=A*w
        end
      elseif t=='n'
        w=A*randn(n)
        if gram
          w=A'*w
        end
      end
      v=Q'*w;
      Y=[Y w-Q*v]; #Y[:,j+r]=w-Q*v;
      v=q'*Y[:,j+1:j+r-1]
      Y[:,j+1:j+r-1]-=q*v;
      maxcolnorm=maximum([norm(Y[:,i]) for i=j+1:j+r]);
    end
  end
  Q
end

"""
   randsvd(A; <keyword arguments>)

Randomized SVD algorithm - returns left singular vectors and singular values of a matrix.

## Arguments:
- `A::Matrix`
- `tol`: Tolerance - discard singular values below tol. Default: 1e-8.
- `maxit`: Maximal number of iterations. Default: 1000.
- `reqrank`: Number of singular values and singular vectors to compute. Optional.
- `r`: Number of samples for stopping criterion. Default: r=10.
- `p`: Oversampling parameter. Default: p=10.
"""
function randsvd(A::AbstractMatrix{<:Number},svdvecs="left";tol=1e-8,maxit=1000,reqrank=0,r=10,p=10)
  m,n=size(A)
  if svdvecs=="left" || svdvecs=="right"
    svdvecs=="left" ? t ='t' : t='n'
    Q=randrange(A,true,t,tol=tol,maxit=maxit,reqrank=reqrank,r=r,p=p)
    svdvecs=="left" ? B=A'*Q : B=A*Q
    if isreal(B)
      B=Symmetric(B'*B)
    else
      B=Hermitian(B'*B)
    end
    E=eigen(B,tol,Inf)
    U=E.vectors[:,end:-1:1];
    S=sqrt.(abs.(E.values[end:-1:1]))
    if reqrank != 0
      if size(U,2)<reqrank
	      @warn "Requested rank exceeds the actual rank. Try changing tolerance."
        U=Q*U;
      else
	      U=Q*U[:,1:reqrank]
      	S=S[1:reqrank]
      end
    else
      K=findall(x-> x>tol ? true : false,S)
      U=Q*U[:,K]
      S=S[K]
    end
    U,S
  elseif svdvecs=="both"
    Q=randrange(A,false,'n',tol=tol,maxit=maxit,reqrank=reqrank,r=r,p=p)
    B=transpose(Q)*A
    U,S,V=svd(B)
    if reqrank != 0
      if size(U,2)<reqrank
	      @warn "Requested rank exceeds the actual rank. Try changing tolerance."
        U=Q*U;
      else
	      U=Q*U[:,1:reqrank]
      	S=S[1:reqrank]
        V=V[:,1:reqrank]
      end
    else
      K=findall(x-> x>tol ? true : false,S)
      U=Q*U[:,K]
      S=S[K]
      V=V[:,K]
    end
    U,S,V
  end
end