text:
include("diagonal_sbp.jl")
include("deriv_ops_new.jl")
include("split_matrix_free_dev.jl")
include("split_matrix_free.jl")

# using CUDAdrv
# CUDAdrv.CuDevice(0)

# Solving Poisson Equation
# Δ u(x,y) = f(x,y)
# Manufactured Solution: u(x,y) = sin(πx .+ πy)

# using CuArrays, CUDAnative
using LinearAlgebra
using SparseArrays
using Plots
using CUDA
using IterativeSolvers
using BenchmarkTools
using MAT

function e(i,n)
    # A = Matrix{Float64}(I,n,n)
    # return A[:,i]
    out = spzeros(n)
    out[i] = 1.0
    return out
end

function eyes(n)
    # return Matrix{Float64}(I,n,n)
    out = spzeros(n,n)
    for i in 1:n
        out[i,i] = 1.0
    end
    return out
end

function u(x,y)
    return sin.(π*x .+ π*y)
end

function Diag(A)
    # Self defined function that is similar to Matlab Diag
    return Diagonal(A[:])
end

function Operators_2d(i, j, p=2, h_list_x = ([1/2^1, 1/2^2, 1/2^3, 1/2^4, 1/2^5, 1/2^6, 1/2^7, 1/2^8, 1/2^9, 1/2^10, 1/2^11, 1/2^12, 1/2^13]),
			 h_list_y = ([1/2^1, 1/2^2, 1/2^3, 1/2^4, 1/2^5, 1/2^6, 1/2^7, 1/2^8, 1/2^9, 1/2^10, 1/2^11, 1/2^12, 1/2^13])
			 )
    hx = h_list_x[i];
    hy = h_list_y[j];

    x = range(0,step=hx,1);
    y = range(0,step=hy,1);
    m_list = 1 ./h_list_x;
    n_list = 1 ./h_list_y;

    # Matrix Size
    N_x = Integer(m_list[i]);
    N_y = Integer(n_list[j]);

    (D1x, HIx, H1x, r1x) = diagonal_sbp_D1(p,N_x,xc=(0,1));
    (D2x, S0x, SNx, HI2x, H2x, r2x) = diagonal_sbp_D2(p,N_x,xc=(0,1));

    (D1y, HIy, H1y, r1y) = diagonal_sbp_D1(p,N_y,xc=(0,1));
    (D2y, S0y, SNy, HI2y, H2y, r2y) = diagonal_sbp_D2(p,N_y,xc=(0,1));

    # BSx = sparse(SNx - S0x);
    # BSy = sparse(SNy - S0y);
    BSx = SNx - S0x
    BSy = SNy - S0y

    # Forming 2d Operators
    # e_1x = sparse(e(1,N_x+1));
    # e_Nx = sparse(e(N_x+1,N_x+1));
    # e_1y = sparse(e(1,N_y+1));
    # e_Ny = sparse(e(N_y+1,N_y+1));
    e_1x = e(1,N_x+1);
    e_Nx = e(N_x+1,N_x+1);
    e_1y = e(1,N_x+1);
    e_Ny = e(N_y+1,N_y+1);

    # I_Nx = sparse(eyes(N_x+1));
    # I_Ny = sparse(eyes(N_y+1));
    I_Nx = eyes(N_x+1);
    I_Ny = eyes(N_y+1);

    e_E = kron(e_Nx,I_Ny);
    e_W = kron(e_1x,I_Ny);
    e_S = kron(I_Nx,e_1y);
    e_N = kron(I_Nx,e_Ny);

    E_E = kron(sparse(Diag(e_Nx)),I_Ny);   # E_E = e_E * e_E'
    E_W = kron(sparse(Diag(e_1x)),I_Ny);
    E_S = kron(I_Nx,sparse(Diag(e_1y)));
    E_N = sparse(kron(I_Nx,sparse(Diag(e_Ny))));

    D1_x = kron(D1x,I_Ny);
    D1_y = kron(I_Nx,D1y);

    D2_x = kron(D2x,I_Ny);
    D2_y = kron(I_Nx,D2y);
    D2 = D2_x + D2_y

    HI_x = kron(HIx,I_Ny);
    HI_y = kron(I_Nx,HIy);

    H_x = kron(H1x,I_Ny);
    H_y = kron(I_Nx,H1y);

    BS_x = kron(BSx,I_Ny);
    BS_y = kron(I_Nx,BSy);

    HI_tilde = kron(HIx,HIx);
    H_tilde = kron(H1x,H1y);

    return (D1_x, D1_y, D2_x, D2_y, D2, HI_x, HI_y, BS_x, BS_y, HI_tilde, H_tilde, I_Nx, I_Ny, e_E, e_W, e_S, e_N, E_E, E_W, E_S, E_N)
end

h_list_x = [1/2^1, 1/2^2, 1/2^3, 1/2^4, 1/2^5, 1/2^6, 1/2^7, 1/2^8, 1/2^9, 1/2^10, 1/2^11, 1/2^12, 1/2^13]
h_list_y = [1/2^1, 1/2^2, 1/2^3, 1/2^4, 1/2^5, 1/2^6, 1/2^7, 1/2^8, 1/2^9, 1/2^10, 1/2^11, 1/2^12, 1/2^13]

rel_errs = []
iter_errs = []
# for k in 1:length(h_list_x)
println("################### BEGIN TEST #########################")
for k in 6:13

    println()
    i = j  = k
    println("##########   Starting Test for k = ", k, "   ######################")

    hx = h_list_x[i];
    hy = h_list_y[j];

    x = range(0,step=hx,1);
    y = range(0,step=hy,1);
    m_list = 1 ./h_list_x;
    n_list = 1 ./h_list_y;

    # Matrix Size
    N_x = Integer(m_list[i]);
    N_y = Integer(n_list[j]);

    Nx = N_x + 1;
    Ny = N_y + 1;

    println("Nx = $Nx, Ny=$Ny")

    # 2D operators
    (D1_x, D1_y, D2_x, D2_y, D2, HI_x, HI_y, BS_x, BS_y, HI_tilde, H_tilde, I_Nx, I_Ny, e_E, e_W, e_S, e_N, E_E, E_W, E_S, E_N) = Operators_2d(i,j);

     # Analytical Solutions
    analy_sol = u(x,y');

    # Penalty Parameters
    tau_E = -13/hx;
    tau_W = -13/hx;
    tau_N = -1;
    tau_S = -1;

    beta = 1;

    # Forming SAT terms

    ## Formulation 1
    SAT_W = tau_W*HI_x*E_W + beta*HI_x*BS_x'*E_W;
    SAT_E = tau_E*HI_x*E_E + beta*HI_x*BS_x'*E_E;

    # SAT_S = tau_S*HI_y*E_S*D1_y
    # SAT_N = tau_N*HI_y*E_N*D1_y

    SAT_S = tau_S*HI_y*E_S*BS_y;
    SAT_N = tau_N*HI_y*E_N*BS_y;

    SAT_W_r = tau_W*HI_x*E_W*e_W + beta*HI_x*BS_x'*E_W*e_W;
    SAT_E_r = tau_E*HI_x*E_E*e_E + beta*HI_x*BS_x'*E_E*e_E;
    SAT_S_r = tau_S*HI_y*E_S*e_S;
    SAT_N_r = tau_N*HI_y*E_N*e_N;

    (alpha1,alpha2,alpha3,alpha4,beta) = (tau_N,tau_S,tau_W,tau_E,beta);

    g_W = sin.(π*y);
    g_E = -sin.(π*y);
    g_S = -π*cos.(π*x);
    # g_N = -π*cos.(π*x)
    g_N = π*cos.(π*x .+ π);

    # Solving with CPU
    A = D2 + SAT_W + SAT_E + SAT_S + SAT_N;

    b = -2π^2*u(x,y')[:] + SAT_W_r*g_W + SAT_E_r*g_E + SAT_S_r*g_S + SAT_N_r*g_N;

    A = - H_tilde*A;
    b = - H_tilde*b;

    @show nnz(A) * sizeof(Float64)

    @show Base.summarysize(A)
    @show Base.summarysize(b)
    @show Nx * Ny * sizeof(Float64)

    @show Base.summarysize(A) / Base.summarysize(b)

    println()

    A_d = CUDA.CUSPARSE.CuSparseMatrixCSC(A);
    b_d = CuArray(b);

    # @show Base.summarysize(A_d)
    # @show Base.summarysize(A_d)

    # testing matrix_split method
    b_reshaped = reshape(b,Nx,Ny);
    b_reshaped_GPU = CuArray(b_reshaped);

    x_GPU = similar(b_reshaped_GPU)
    x_full_GPU = similar(b_reshaped_GPU)
    A_d * b_d
    matrix_free_A(b_reshaped_GPU,x_GPU)
    matrix_free_A_full_GPU(b_reshaped_GPU,x_full_GPU)

    odata = similar(b_reshaped_GPU);
    r_GPU = similar(b_reshaped_GPU);
    p_GPU = similar(b_reshaped_GPU);
    Ap_GPU = similar(b_reshaped_GPU);

    CG_GPU_v2(b_reshaped_GPU,x_GPU,odata,r_GPU,p_GPU,Ap_GPU)

    ## SpMV test

    println("TESTING SpMV vs Matrix_FREE")
    iter_times = Nx + Ny

    time_matrix_free = @elapsed begin
        for _ in 1:iter_times
            matrix_free_A(b_reshaped_GPU,x_GPU)
        end
    end
    time_matrix_free = time_matrix_free * 1000 / iter_times
    @show time_matrix_free

    time_matrix_free_full_GPU = @elapsed begin
        for _ in 1:iter_times
            matrix_free_A_full_GPU(b_reshaped_GPU,x_full_GPU)
        end
    end
    time_matrix_free_full_GPU = time_matrix_free_full_GPU * 1000 / iter_times
    @show time_matrix_free_full_GPU

    time_CUBLAS = @elapsed begin
        for _ in 1:iter_times
            A_d * b_d
        end
    end
    time_CUBLAS = time_CUBLAS * 1000 / iter_times
    @show time_CUBLAS
    println()

    ## End SpMV test

    x_GPU = CuArray(zeros(Nx,Ny));
    (iter_steps,rsold_GPU) = CG_GPU(b_reshaped_GPU,x_GPU);
    println("Iteration steps till convergence for x_GPU: $iter_steps")
    @show rsold_GPU
    matrix_free_GPU_time_in_CG = (time_matrix_free * iter_steps)

    x_full_GPU = CuArray(zeros(Nx,Ny));
    (iter_steps,rsold_full_GPU) = CG_full_GPU(b_reshaped_GPU,x_full_GPU);
    matrix_free_full_GPU_time_in_CG = (time_matrix_free_full_GPU * iter_steps)
    println("Iteration steps till convergence for x_full_GPU: $iter_steps")
    @show rsold_full_GPU

    _,history= cg(A_d,b_d,log=true)
    @show history
    @show history.data[:resnorm][end]
    @show history.data[:resnorm][end-1]

    # x_GPU = CuArray(zeros(Nx,Ny));
    #         # CG_GPU_dev(b_reshaped_GPU,x_GPU)
    # CG_GPU(b_reshaped_GPU,x_GPU)

    iter_times = min(5,max(13-k,1))
    println()
    println("Starting Timing, results in ms")

    println("Time for matrix_free SpMV in CG: $matrix_free_GPU_time_in_CG")
    println("Time for matrix_free_full_GPU SpMV in CG: $matrix_free_full_GPU_time_in_CG")
    CUBLAS_SpMV_time_in_CG = (time_CUBLAS * history.iters)
    println("Time for CUBLAS SpMV in CG_IterativeSolvers: $CUBLAS_SpMV_time_in_CG\n")

    println()

    t_CG_GPU = @elapsed begin
        for i in 1:iter_times
            x_GPU = CuArray(zeros(Nx,Ny))
            # CG_GPU_dev(b_reshaped_GPU,x_GPU)
            CG_GPU(b_reshaped_GPU,x_GPU)
        end
    end
    t_CG_GPU = t_CG_GPU * 1000 / iter_times
    @show t_CG_GPU

    # t_CG_GPU_v2 = @elapsed begin # testing preallocation
    #     for i in 1:iter_times
    #         x_GPU = CuArray(zeros(Nx,Ny))
    #         # CG_GPU_dev(b_reshaped_GPU,x_GPU)
    #         CG_GPU_v2(b_reshaped_GPU,x_GPU,odata,r_GPU,p_GPU,Ap_GPU)
    #     end
    # end
    # t_CG_GPU_v2 = t_CG_GPU_v2 * 1000 / iter_times
    # @show t_CG_GPU_v2

    t_CG_full_GPU = @elapsed begin
        for i in 1:iter_times
            x_full_GPU = CuArray(zeros(Nx,Ny))
            # CG_GPU_dev(b_reshaped_GPU,x_GPU)
            CG_full_GPU(b_reshaped_GPU,x_full_GPU)
        end
    end
    t_CG_full_GPU = t_CG_full_GPU * 1000 / iter_times
    @show t_CG_full_GPU

    t_CG_GPU_IterativeSolvers = @elapsed begin
        for i in 1:iter_times
            cg(A_d,b_d,log=true)
        end
    end
    t_CG_GPU_IterativeSolvers = t_CG_GPU_IterativeSolvers * 1000 / iter_times
    @show t_CG_GPU_IterativeSolvers

    println("End time comparison\n")

    ## Compare Efficiency

    overhead_CG_matrix_free =  t_CG_GPU - matrix_free_GPU_time_in_CG
    overhead_CG_matrix_free_full_GPU = t_CG_full_GPU - matrix_free_full_GPU_time_in_CG
    overhead_CG_IterativeSolvrs = t_CG_GPU_IterativeSolvers - CUBLAS_SpMV_time_in_CG

    @show overhead_CG_matrix_free
    @show overhead_CG_matrix_free_full_GPU
    @show overhead_CG_IterativeSolvrs

    println()

    efficiency_CG_matrix_free = matrix_free_GPU_time_in_CG / t_CG_GPU
    efficiency_CG_matrix_free_full_GPU = matrix_free_full_GPU_time_in_CG / t_CG_full_GPU
    efficiency_CG_IterativeSolvrs = CUBLAS_SpMV_time_in_CG / t_CG_GPU_IterativeSolvers

    @show efficiency_CG_matrix_free
    @show efficiency_CG_matrix_free_full_GPU
    @show efficiency_CG_IterativeSolvrs

    println()

    ## End Compare Efficiency
    println("##########   Ending Test for k = ", k, "   ######################\n")

end

function check_richardson(A,ω)
    return opnorm(Matrix(I,size(A)) - ω*A)
end

function modified_richardson(A,b;ω=0.15)
    x = zeros(length(b))
    # ω = -0.15
    count = 0
    norms = []
    for _ in 1:100
        iter_norm = norm(A*x-b)
        append!(norms,iter_norm)
        if iter_norm >= 1e-6 * norms[1]
            count += 1
            x = x + ω*(b- A*x)
        end
    end
    return (x,count,norms)
end

function jacobi_diy(A,b)
    x = zeros(length(b))
    D = Diagonal(A)
    LU = A - D
    # ω = -0.15
    count = 0
    norms = []
    for _ in 1:100
        iter_norm = norm(A*x-b)
        append!(norms,iter_norm)
        if iter_norm >= 1e-6 * norms[1]
            count += 1
            x = D\(b - LU*x)
        end
    end
    return (x,count,norms)
end

function test(A,b;ω=1/13)
    (x_j,count_j,norms_j) = jacobi_diy(A,b)
    (x_r,count_r,norms_r) = modified_richardson(A,b,ω=ω)
    plot(norms_j)
    plot!(norms_r)
end


text:
#=
parametric_types:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-09
=#

module Parametric
    # struct
    mutable struct Point{Type}
        x::Type
        y::Type
    end

    # methods
    function get_x_type(p::Point{Int})
        return typeof(p.x)
    end

    function get_x_type(p::Point{Float64})
        return typeof(p.x)
    end

    function main()
        p_int = Point(2, 3)
        p_float = Point(2.5, 3.8)
        println("Type of p_int.x = ", get_x_type(p_int))
        println("Type of p_float.x = ", get_x_type(p_float))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .Parametric
    Parametric.main()
end

ext: jl

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
function func()
    x = 3
    while x < 10
        x = x + 1
    end
    println(x)
end

mutable struct test
    x
end

a = test([3])

function update!(testobj, num)
    push!(testobj.x, num)
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
abstract type Kozachenko end

abstract type Kraskov end
abstract type Kraskov1 <: Kraskov end
abstract type Kraskov2 <: Kraskov end


text:
module PolyesterWeave

using BitTwiddlingConvenienceFunctions: nextpow2
using ThreadingUtilities: _atomic_store!, _atomic_or!, _atomic_xchg!
using Static
using IfElse: ifelse
using CPUSummary: num_threads

export request_threads, free_threads!

@static if VERSION ≥ v"1.6.0-DEV.674"
  @inline function assume(b::Bool)
    Base.llvmcall(("""
      declare void @llvm.assume(i1)

      define void @entry(i8 %byte) alwaysinline {
      top:
        %bit = trunc i8 %byte to i1
        call void @llvm.assume(i1 %bit)
        ret void
      }
  """, "entry"), Cvoid, Tuple{Bool}, b)
  end
else
  @inline assume(b::Bool) = Base.llvmcall(("declare void @llvm.assume(i1)", "%b = trunc i8 %0 to i1\ncall void @llvm.assume(i1 %b)\nret void"), Cvoid, Tuple{Bool}, b)
end

const WORKERS = Ref{NTuple{8,UInt64}}(ntuple(((zero ∘ UInt64)), Val(8)))

include("unsignediterator.jl")
include("request.jl")

dynamic_thread_count() = min((Sys.CPU_THREADS)::Int, Threads.nthreads())
function worker_mask_init(x)
  if x < 0
    0x0000000000000000
  elseif x ≥ 64
    0xffffffffffffffff
  else
    ((0x0000000000000001 << x) - 0x0000000000000001)
  end
end
# function static_thread_init()
#   nt = num_threads()
#   Base.Cartesian.@ntuple 8 i -> worker_mask_init(nt - (i-1)*64)
# end
function reset_workers!()
  # workers = ntuple(((zero ∘ UInt64)), Val(8))
  nt = dynamic_thread_count() - 1
  WORKERS[] = Base.Cartesian.@ntuple 8 i -> worker_mask_init(nt - (i-1)*64)
end
function __init__()
  reset_workers!()
end

end


text:
include("dynparams.jl")
include("rne_park.jl")


text:
using Base: @kwdef
using Reexport
@reexport using KhepriBase
using Dates
using ColorTypes

# resolve conflicts
using KhepriBase:
  XYZ,
  Text

import Base:
  show

import KhepriBase:
  backend,
  backend_bounding_box,
  #backend_captured_shape,
  #backend_captured_shapes,
  #backend_changed_shape,
  #backend_create_block,
  #backend_disable_update,
  #backend_enable_update,
  #backend_extrusion,
  backend_fill_curves,
  backend_frame_at,
  #backend_generate_captured_shape,
  #backend_generate_captured_shapes,
  backend_get_family_ref,
  #backend_highlight_shape,
  #backend_highlight_shapes,
  #backend_loft_curve_point,
  #backend_loft_curves,
  #backend_loft_surface_point,
  #backend_loft_surfaces,
  #backend_map_division,
  backend_name,
  #backend_panel,
  #backend_pre_selected_shapes_from_set,
  backend_revolve_curve,
  backend_revolve_point,
  backend_revolve_surface,
  #backend_right_cuboid,
  #backend_rotation_minimizing_frames,
  #backend_set_length_unit,
  #backend_shape_from_ref,
  #backend_spotlight,
  #backend_stroke_arc,
  #backend_stroke_unite,
  #backend_surface_boundary,
  #backend_surface_domain,
  #backend_wall,
  decode,
  encode,
  intersect_ref,
  parse_signature,
  realize,
  #slice_ref,
  maybe_merged_node,
  maybe_merged_bar,
  save_shape!,
  subtract_ref,
  switch_to_layer,
  truss_bar_family_cross_section_area,
  unite_ref,
  #unite_refs,
  with_material_as_layer,
  use_material_as_layer


text:
### Logistic function in its general form

# See https://en.wikipedia.org/wiki/Logistic_function
# l scales the curve. For l>1, the curve is stretched, whereas for l<1 it is shrinked
# The curve's maximum value coincides with l
# k gives the curve's steepness. For larger k, the curve becomes more steep
# x0 is the x value of the curve's midpoint
# y0 is the y value of the curve's midpoint
# The logistic function has been defined for tuning purposes

logistic(x::Real, l::Real=1., k::Real=1., x0::Real=0., y0::Real=0.) = l/(1+exp(-k*(x-x0)))+y0


text:
#!/usr/bin/env bash
    #=
    exec julia --project="$(realpath $(dirname $(dirname $0)))" --color=yes --startup-file=no -e "include(popfirst!(ARGS))" \
    "${BASH_SOURCE[0]}" "$@"
    =#

# TAKE ONE

using Nemo, Primes

function _isprime(n)
	return Nemo.isprime(Nemo.fmpz(n))
end

function _factors(n)
	return n == 0 ? [] : Nemo.factor(Nemo.fmpz(n))
end
# function _factors(n)
# 	first(i) for i in Primes.factor(n)
# end

function primedivisors(n)
    n == 0 && return []
    _isprime(n) && return [Nemo.fmpz(n)]
    f = _factors(n)
    return sort!([p for (p, e) ∈ f])
end

function ω(n)
	nprimedivisors = 0
	if n == 0
		nprimedivisors = 0
	end
	if _isprime(n)
		nprimedivisors = length(Nemo.fmpz(n))
	end

	return length(_factors(n))

	# return fmpz(length(primedivisors(n)))
end

function Ω(n)
    n == fmpz(0) && return 0
    _isprime(n) && return fmpz(1)

    return sum(e for (_, e) ∈ _factors(n))
end

function isperfectpower(n)
	return ω(n) == 1 && Ω(n) != 1
end

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

#=
You find the smallest pair (a, n) (small with respect to n) such that a^n is your target number. Then you need to check if a is prime.
=#
function isprimepower(n::Integer)
	# return any(iszero(mod(n, p)) for p in primes(ceil(Int, sqrt(n))))

	n ≤ 1 && return false
	isprime(n) && return true
	ub = isqrt(n) + 1

	for a in primes(ub)
		for m in 2:ub # m > 1
			res = a^m
			res > n && break
			isequal(res, n) && return true
		end
	end

	return false
end

function primepower_brute_force(n::Integer)
	n ≤ 1 && return nothing
	isprime(n) && return n, 1

	for a in primes(ceil(Int, sqrt(n)))
		for m in 2:ceil(Int, sqrt(n)) # m > 1
			a^m > n && break
			isequal(a^m, n) && return a, m
		end
	end

	return nothing
end

function perfectpower_brute_force(n::Integer)
	n ≤ 1 && return nothing

	for a in 1:n
		for m in 2:n
			isequal(a^m, n) && return m, a
		end
	end

	return nothing
end

#========================================================================#

## The following function are obtained from DJB

struct DJBFloat
	a::Signed
	n::Unsigned
end

lg(x) = log2(x)

#=
From Knuths The Art of Computer Programming, semi-numerical algorithms, 2nd edition, 4.3.1, exercise 16.

Design an algorithm that divides a nonnegative n-place integer $(u_1u_2\ldots u_n)_b$ by $v$, where $v$ is a single-precision number (i.e., $0 < v < b$), producing the quotient $(w_1w_2\ldots 2_n)_b$ and remainder $r$.
=#

function Base.divrem(u::Integer, v::Real, b::Integer)
	r, j = 0, 1
	n = ndigits(u) #; base = b
	u̲ = round.(reverse(digits(u)), base = b)
	w̲ = Vector{Int}(undef, n)

	while j ≤ n
		w̲[j] = floor(Int, (r*b + u̲[j]) / v)
		r = mod(r*b + u̲[j], v)
		j += 1
	end

	# w = sum(w̲[i] * 10 ^ (length(w̲) - i) for i in eachindex(w̲))
	w = parse(Int, join(string.(w̲)))
	# w = foldr((i, j) -> 10 * j + i, reverse(w̲))
	#
	# n = 0
	# reduce(w̲) do i, j
	# 	10 * j + i
	# 	# global n = n + ndigits(w̲[i])
	# end
	# n = 0
	# w = sum(eachindex(w̲)) do i
	# 	w̲[i] * 10 ^ (length(w̲) - i + n)
	# 	global n = n + ndigits(w̲[i])
	# end
	# n = 0

	return w, r
end

Base.div(u::Integer, v::Real, b::Integer) = first(divrem(u, v, b))
Base.rem(u::Integer, v::Real, b::Integer) = last(divrem(u, v, b))

# function Base.div(a::Integer, n::Integer, k::Integer, b::Integer)
# 	f = nextpow(2, a)
# 	return a + f - ceil(Int, lg(k)) - b, floor(Int, n / (2^(f - ceil(Int, lg(k)) - b) * k))
# end

#=
In this section, I define truncation to $b$ bits, written $\text{trunc}_b$, and show that $\frac{r}{\text{trunc}_br}$ is between $1$ and $1+2^{1-b}$.  More generally, for any positive integer $k$, I define $\text{div}_b(r,k)$ as a floating point approximation to $\frac{r}{k}$, so that $\frac{r}{k\text{div}_b(r,k)}$ is between $1$ and $1+2^{1-b}$.\\ Fix $b\geq 1$.  Set $\text{div}_b(a,n,k)=(a+f-\left\lceil\text{lg} k\right\rceil - b, \left\lfloor\frac{n}{2^{f-\left\lceil\text{lg} k\right\rceil-b}k}\right\rfloor)$, where $2^{f-1}\leq n<w^f$.  (Note that $f-\left\lceil \text{lg} k\right\rceil - b$ may be negative.)  This map induces a map, also denoted $\text{div}_b$, upon positive floating-point numbers: \begin{equation}\text{div}_b(2^an,k)=2^{a+f-\left\lceil\text{lg} k\right\rceil - b}\left\lfloor\frac{n}{2^{f-\left\lceil\text{lg}k\right\rceil - b}k}\right\rfloor\quad\quad\quad\text{if }2^{f-1}\leq n<2^f.\end{equation}\\ To compute $\text{div}_b(r,k)$ I also use an algorithm designed for dividing by small integers; time spent computing $\text{div}_b$ is not $M$-time.
=#
function Base.trunc(n::Real, b::Integer)
    return n & (1 << b - 1)
	# return div(n, 1, b)
end

#=
Write a MIX program that multiplies $(u_1u_2\ldots u_n)_b$ by $v$, where $v$ is a single-precision number (i.e., $0 \leq v < b$), producing the answer $(w_0w_1\ldots w_n)_b$.  How much running time is required?

		ENTX	N		1	ADD		CARRY	N
		JOV		OFLO	1	JNOV	*+2		N
		ENTX	0		1	INCX	1		K
	2H	STX		CARRY	N	STA		W,1		N
		LDA		U,1		N	DEC1	1		N
		MUL		V		N	J1P		2B		N
		SLC		5		N	STX		W		1

The running time is 23N + K + 5 cycles, and K is rougly (1/2)N.
=#
function mul(u::Integer, v::Real, b::Integer)
	n = ndigits(u) #; base = b
	u̲ = round.(reverse(digits(u)), base = b)
	w̲ = Vector{Int}(undef, n)
end

#=
Given a positive floating-point number r and two positive integers b, k, to prime pow_b(r, k).  Let $\text{pow}_b(r,k)$ be the $b$-bit approximate $k$th power of $r$.
\begin{equation}\begin{aligned}
\text{pow}_b(r,1) &= \text{trunc}_br; \\
\text{pow}_b(r,2k) &= \text{trunc}_b(\text{pow}_b(r,k)^2); \\
\text{pow}_b(r,2k+1) &= \text{trunc}_b(\text{pow}_b(r,2k)\text{pow}_b(r,1))
\end{aligned}\end{equation}
=#
function __pow(r::Real, b::Integer, k::Integer)
	if r ≤ 0 || b ≤ 0 || k ≤ 0
		throw(error("r, b, and k must be positive numbers to compute pow_b(r, k)."))
	end
	println("pow: $r, $b, $k")
	isone(k) && return trunc(r, b)
	iszero(mod(k, 2)) && return trunc(__pow(r, b, k ÷ 2)^2, b) # use \div because k is even
	return trunc(__pow(r, b, k - 1) * trunc(r, b), b)
end

#=
I am trying to find a root $z$ of $z^ky-1$.
To compute nroot_b(y, k) for 1 \le b \le \ceil{lg8k}: In advance, find the exponent g satisfying 2^{g-  1} < y < 2^g, and set a = \floor{-g / k}, so that 2^a \le y^{-1/k} < 2^{a+1}.  Also set B = \ceil{lg(66(2k + 1))}
=#
function __nroot(y::Real, k::Real, b::Real)
	# find the exponent g satisfying 2^(g - 1) < y < 2^g
	g = nothing
	i = 0
	while true
		if 2^(float(i) - 1) < y && y ≤ float(2)^i
			g = i
			break
		end
		i += 1
	end

	a = floor(Int, -g / k)
	B = ceil(Int, lg(66 * (2*k + 1)))
	z = 2^float(a) + 2^(float(a) - 1)
	println("nroot: $z")
	j = 1

	while ! isequal(j, b)
		println("nroot: $y, $k, $j, $B")
		z = __nroot(y, k, j)
		r = trunc(__pow(z, k, B) * trunc(y, B), B)
		z = ifelse(r ≤ 993//1024, z + 2^(a - j + 1), z)
		z = ifelse(r > 1, z - 2^(a - j - 1), z)
		j = j + 1
	end

	return z
end

#=
To compute nroot_b(y, k) for b \ge 4 + \ceil{lg 8k} + 1: In advance set b' = \ceil{lg 2k} + \ceil{(b - \ceil{lg 2k}) / 2} and B = 2b' + 4 - \ceil{lg k}.  Note that b; < b.
=#
function __newton_nroot(y::Real, k::Real, b::Real)
	b′ = ceil(Int, (b - ceil(Int, lg(2*k))) / 2)
	B = 2*b′ + 4 - ceil(Int, lg(l))
	z = ifelse(b′ ≤ ceil(Int, lg(8*k)), __nroot(__nroot(y, k, b′)), __newton_nroot(y, k, b′)) # if else, then b′ ≥ ceil(Int, lg(8*k)) + 1
	r₂ = trunc(z, B) * (k + 1)
	r₃ = trunc(__pow(z, k + 1, B) * trunc(y, B), B)

	return r₄ = div(r₂ - r₃, k, B)
end

#=
Given positive integers n, x, k, to compute the sign of n - x^k: In advance, set f = \ceil{lg 2n}
=#
function __sign(n::Integer, x::Integer, k::Integer)
	if n ≤ 0 || x ≤ 0 || k ≤ 0
		throw(error("n, x, and k must be positive numbers to compute the sign of n - x^k using this algorithm."))
	end

	f = ceil(Int, lg(2n))
	b = 1
	while true
		r = __pow(x, k, b + ceil(Int, lg(8k)))

		n < r && return -1
		r*(1 + 2^(-b)) ≤ n && return 1
		b ≥ f && return 0

		b = min(2*b, f)
	end
end

#=
Given integers n \ge 2 and k \ge 2, and a positive floating-point number y, to see if n is a kth power: in advance, set f = \floor{lg 2n} and b = 3 + \ceil{f / k}
=#
function __iskthpower(n::Integer, k::Integer, y::Real)
	if n < 2 || k < 2 || y ≤ 0
		thow(error("n, k, and y must be positive (and n and k must be more than one) to see if n is a kth power using this algorithm."))
	end

	f = floor(Int, lg(2*n))
	b = 3 + ceil(Int, f / k)
	r = __nroot(y, k, b)

	x = nothing
	for i in 1:r
		if abs(r - x) ≤ 5//8
			x = i
		end
	end

	(iszero(x) || abs(r - x) ≥ 1//4) && return 0

	sign = __sign(n, x, k) # why do we do this step?!
	isequal(n, x^k) && return x

	return 0
end

#=
Given an integer n \ge 2, to decompose n as a perfect power if possible: in advance, set f = \floor{lg 2n}
=#
function __find_prime_power(n::Integer)
	if n < 2
		throw(error("This algorithm to compute the prime power if possible only works for integers more than one."))
	end

	# isequal(n, 2) && return n, 1

	f = floor(Int, lg(2*n))
	y = __nroot(n, 1, 3 + ceil(Int, f / 2))

	for p in primes(f)
		x = __iskthpower(n, p, y)
		x > 0 && return x, p
	end

	return x, 1
end


text:
"""
    DiabaticModels

Models defined within this module subtype the [`DiabaticModel`](@ref) and provide
potentials as Hermitian matrices and derivatives as arrays of Hermitian matrices.
"""
module DiabaticModels

using ..NQCModels: NQCModels
using Unitful: @u_str
using UnitfulAtomic: austrip

using Parameters: Parameters
using LinearAlgebra: Hermitian
using StaticArrays: SMatrix, SVector

"""
    DiabaticModel <: Model

`DiabaticModel`s are used when a system has multiple electronic states that are
presented in the diabatic representation. This is the case for the majority
of model systems.

# Implementation

`DiabaticModel`s should implement:
* `potential(model, R)`
* `derivative!(model, D, R)`
* `nstates(model)`
* `ndofs(model)`

# Example

In this example we create a simple 2 state, 1 dimensional diabatic model `MyModel`.
As noted above, we implement the 4 relevant functions then evaluate the potential.
Since this is a 1D model the argument `R` accepts a `Real` value.

```jldoctest
using StaticArrays: SMatrix
using LinearAlgebra: Hermitian

struct MyModel <: NQCModels.DiabaticModels.DiabaticModel end

NQCModels.nstates(::MyModel) = 2
NQCModels.ndofs(::MyModel) = 1

function NQCModels.potential(::MyModel, R::Real)
    V11 = R
    V22 = -R
    V12 = 1
    return Hermitian(SMatrix{2,2}(V11, V12, V12, V22))
end

function NQCModels.derivative!(::MyModel, D, R::Real)
    return Hermitian(SMatrix{2,2}(1, 0, 0, 1))
end

model = MyModel()
NQCModels.potential(model, 10)

# output

2×2 Hermitian{Int64, SMatrix{2, 2, Int64, 4}}:
 10    1
  1  -10
```
"""
abstract type DiabaticModel <: NQCModels.Model end

"""
    LargeDiabaticModel <: DiabaticModel

Same as the `DiabaticModels` but uses normal Julia arrays instead of StaticArrays and must
implement the inplace `potential!` rather than `potential`.
This is useful when `nstates` is very large and StaticArrays are no longer efficient.
"""
abstract type LargeDiabaticModel <: DiabaticModel end

function NQCModels.derivative!(model::LargeDiabaticModel, D, R::AbstractMatrix)
    if NQCModels.ndofs(model) == 1
        if size(R, 2) == 1
            NQCModels.derivative!(model, D[1], R[1])
            return D
        else
            NQCModels.derivative!(model, view(D, 1, :), view(R, 1, :))
            return D
        end
    elseif size(R, 2) == 1
        NQCModels.derivative!(model, view(D, :, 1), view(R, :, 1))
        return D
    else
        throw(MethodError(NQCModels.derivative!, (model, D, R)))
    end
end

"""
    DiabaticFrictionModel <: LargeDiabaticModel

These models are defined identically to the `LargeDiabaticModel` but
allocate extra temporary arrays when used with `NQCDynamics.jl`.

This allows for the calculation of electronic friction
internally from the diabatic potential after diagonalisation
and calculation of nonadiabatic couplings.
"""
abstract type DiabaticFrictionModel <: LargeDiabaticModel end

function matrix_template(model::DiabaticModel, eltype)
    n = NQCModels.nstates(model)
    return SMatrix{n,n}(zeros(eltype, n, n))
end
function matrix_template(model::LargeDiabaticModel, eltype)
    zeros(eltype, NQCModels.nstates(model), NQCModels.nstates(model))
end

function vector_template(model::DiabaticModel, eltype)
    n = NQCModels.nstates(model)
    return SVector{n}(zeros(eltype, n))
end
function vector_template(model::LargeDiabaticModel, eltype)
    zeros(eltype, NQCModels.nstates(model))
end

function NQCModels.zero_derivative(model::DiabaticModel, R::AbstractMatrix)
    [Hermitian(matrix_template(model, eltype(R))) for _ in CartesianIndices(R)]
end

function NQCModels.potential(model::LargeDiabaticModel, R::AbstractMatrix)
    V = Hermitian(matrix_template(model, eltype(R)))
    NQCModels.potential!(model, V, R)
    return V
end

include("double_well.jl")
export DoubleWell
include("tully_models.jl")
export TullyModelOne
export TullyModelTwo
export TullyModelThree
include("three_state_morse.jl")
export ThreeStateMorse
include("spin_boson.jl")
export OhmicSpectralDensity
export DebyeSpectralDensity
export SpinBoson
export BosonBath
include("1D_scattering.jl")
export Scattering1D
include("ouyang_models.jl")
export OuyangModelOne
include("gates_holloway_elbow.jl")
export GatesHollowayElbow
include("subotnik.jl")
export MiaoSubotnik
include("ananth_models.jl")
export AnanthModelOne
export AnanthModelTwo

end # module


text:
module TestConversions
    using Base.Test
    using DataFrames

    df = DataFrame()
    df[:A] = 1:5
    df[:B] = [:A, :B, :C, :D, :E]
    @test isa(convert(Array, df), Matrix{Any})
    @test convert(Array, df) == convert(Array, convert(DataArray, df))
    @test isa(convert(Array{Any}, df), Matrix{Any})

    df = DataFrame()
    df[:A] = 1:5
    df[:B] = 1.0:5.0
    @test isa(convert(Array, df), Matrix{Real})
    @test convert(Array, df) == convert(Array, convert(DataArray, df))
    @test isa(convert(Array{Any}, df), Matrix{Any})
    @test isa(convert(Array{Float64}, df), Matrix{Float64})

    df = DataFrame()
    df[:A] = 1.0:5.0
    df[:B] = 1.0:5.0
    a = convert(Array, df)
    aa = convert(Array{Any}, df)
    ai = convert(Array{Int}, df)
    @test isa(a, Matrix{Float64})
    @test a == convert(Array, convert(DataArray, df))
    @test a == convert(Matrix, df)
    @test isa(aa, Matrix{Any})
    @test aa == convert(Matrix{Any}, df)
    @test isa(ai, Matrix{Int})
    @test ai == convert(Matrix{Int}, df)

    df[1,1] = NA
    @test_throws ErrorException convert(Array, df)
    da = convert(DataArray, df)
    daa = convert(DataArray{Any}, df)
    dai = convert(DataArray{Int}, df)
    @test isa(da, DataMatrix{Float64})
    @test isequal(da, convert(DataMatrix, df))
    @test isa(daa, DataMatrix{Any})
    @test isequal(daa, convert(DataMatrix{Any}, df))
    @test isa(dai, DataMatrix{Int})
    @test isequal(dai, convert(DataMatrix{Int}, df))

    a = [1.0,2.0]
    b = [-0.1,3]
    c = [-3.1,7]
    di = Dict([("a", a), ("b", b), ("c", c)])

    df = convert(DataFrame,di)
    @test isa(df,DataFrame)
    @test names(df) == Symbol[x for x in sort(collect(keys(di)))]
    @test df[:a] == a
    @test df[:b] == b
    @test df[:c] == c

    a = [1.0]
    di = Dict([("a", a), ("b", b), ("c", c)])
    @test_throws ArgumentError convert(DataFrame,di)

end


text:
struct PosixPath <: AbstractPath
    parts::Tuple
end

PosixPath() = PosixPath(tuple())

function PosixPath(str::AbstractString)
    str = String(str)

    if isempty(str)
        return PosixPath(tuple("."))
    end

    tokenized = split(str, POSIX_PATH_SEPARATOR)
    if isempty(tokenized[1])
        tokenized[1] = POSIX_PATH_SEPARATOR
    end
    return PosixPath(tuple(map(String, tokenized)...))
end

# The following should be implemented in the concrete types
==(a::PosixPath, b::PosixPath) = parts(a) == parts(b)
Base.String(path::PosixPath) = joinpath(parts(path)...)
parts(path::PosixPath) = path.parts

Base.show(io::IO, path::PosixPath) = print(io, "p\"$(join(parts(path), '/'))\"")

function isabs(path::PosixPath)
    if parts(path)[1] == POSIX_PATH_SEPARATOR
        return true
    else
        return false
    end
end

drive(path::PosixPath) = ""

function root(path::PosixPath)
    if parts(path)[1] == POSIX_PATH_SEPARATOR
        return POSIX_PATH_SEPARATOR
    else
        return ""
    end
end

function expanduser(path::PosixPath)
    p = parts(path)

    if p[1] == "~"
        if length(p) > 1
            return PosixPath(tuple(homedir(), p[2:end]...))
        else
            return PosixPath(tuple(homedir()))
        end
    end

    return path
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to string processing.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

"""
    _crop_str(str::AbstractString, cols_beg::Int, cols_end::Int = -1)

Crop the string `str` between the columns `cols_beg` and `cols_end`. If any of
the values are negative, then it will be selected as the string limits.

"""
function _crop_str(str::AbstractString, cols_beg::Int, cols_end::Int = -1)
    str_width = textwidth(str)

    cols_beg < 1 && (cols_beg = 1)

    cropped_str = ""
    num_cols = 1
    init = false

    for c in str
        cw = textwidth(c)

        if !init
            Δ = num_cols - cols_beg

            if Δ < 0
                if Δ + cw > 0
                    cropped_str *= " "^(Δ + cw)
                    num_cols += Δ + cw
                    init = true
                    continue
                end

                num_cols += cw
                continue
            end
        end

        num_cols += cw
        (cols_end > 0) && (num_cols > cols_end + 1) && break
        cropped_str *= c
    end

    return cropped_str
end


text:
#=
   1 RM Calculations (strength training)

   See https://en.wikipedia.org/wiki/One-repetition_maximum

   Here we calculate two versions of 1 RM (1 repetition),
   adjusted for reps and also adjusted for age (which in my
   case is necessary :-))

   For age:64 weight:102 reps:5 we get

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

        onerm_epley   119.0048    0.5035     0.0050    0.0049   10111.3967    0.9999     5230.9347
    onerm_age_epley   193.8589    0.8203     0.0082    0.0080   10111.3967    0.9999     5230.9347
      onerm_brzycki   114.7478    0.5052     0.0051    0.0057    9583.0426    1.0007     4957.6009
  onerm_age_brzycki   186.9242    0.8230     0.0082    0.0092    9583.0426    1.0007     4957.6009

  The age adjustments should be interpreted as the weight a person at age 30 would have to lift.

  And for my deadlift PR at 1x117kg: age:64 weight:117 reps:1

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

        onerm_epley   120.9058    0.5066     0.0051    0.0050   9629.8855    1.0000    10031.1308
    onerm_age_epley   196.9556    0.8252     0.0083    0.0082   9629.8855    1.0000    10031.1308
      onerm_brzycki   116.9990    0.4978     0.0050    0.0061   9253.4812    1.0002     9639.0429
  onerm_age_brzycki   190.5914    0.8110     0.0081    0.0099   9253.4812    1.0002     9639.0429

  So, for my 1 rep 117kg lift, a person at age would have to lift about 190kg to be as strong as I am
  now (at age 64).

   Cf ~/blog/OneRM.blog
     ~/webppl/OneRM.wppl

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

@model function one_rm(age=63,weight=102,reps=5)
    age_coeff = Dict(
         30=> 1.000, 31=> 1.016, 32=> 1.031, 33=> 1.046, 34=> 1.059,
         35=> 1.072, 36=> 1.083, 37=> 1.096, 38=> 1.109, 39=> 1.122,
         40=> 1.135, 41=> 1.149, 42=> 1.162, 43=> 1.176, 44=> 1.189,
         45=> 1.203, 46=> 1.218, 47=> 1.233, 48=> 1.248, 49=> 1.263,
         50=> 1.279, 51=> 1.297, 52=> 1.316, 53=> 1.338, 54=> 1.361,
         55=> 1.385, 56=> 1.411, 57=> 1.437, 58=> 1.462, 59=> 1.488,
         60=> 1.514, 61=> 1.541, 62=> 1.568, 63=> 1.598, 64=> 1.629,
         65=> 1.663, 66=> 1.699, 67=> 1.738, 68=> 1.779, 69=> 1.823,
         70=> 1.867, 71=> 1.910, 72=> 1.953, 73=> 2.004, 74=> 2.060,
         75=> 2.117, 76=> 2.181, 77=> 2.255, 78=> 2.336, 79=> 2.419,
         80=> 2.540, 81=> 2.597, 82=> 2.702, 83=> 2.831, 84=> 2.981,
         85=> 3.153, 86=> 3.352, 87=> 3.580, 88=> 3.843, 89=> 4.145,
         90=> 4.493
        )

    weight ~ Uniform(1,200)
    reps   ~ DiscreteUniform(1,30)
    age    ~ Uniform(30,90)

    # OneRM adjusted for reps
    onerm_epley ~ Normal(weight*(1.0 + reps/30.0),0.5)

    # OneRM adjusted for age
    onerm_age_epley ~ Dirac(onerm_epley * age_coeff[round(Int,age)])

    # OneRM adjusted for reps
    onerm_brzycki ~ Normal(weight * (36.0/(37.0 - reps)),0.5)

    # OneRM adjusted for age
    onerm_age_brzycki ~ Dirac(onerm_brzycki * age_coeff[round(Int,age)])

end

function run_one_rm(age,weight,reps)
    println("age:$age weight:$weight reps:$reps")
    model = one_rm(age,weight,reps)
    # chns = sample(model, Prior(), 10_000)
    # chns = sample(model, MH(), 10_000)
    # chns = sample(model, PG(5), 10_000)
    chns = sample(model, SMC(), 10_000)
    # chns = sample(model, IS(), 10_000)

    display(chns)
    # display(plot(chns))
    println()
end

age = 64
weight=102
reps=5
run_one_rm(age,weight,reps)

age = 64
weight=117
reps=1
run_one_rm(age,weight,reps)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
immutable ZeroTo{T<:Signed} <: AbstractUnitRange{T}
    stop::T
    function ZeroTo(stop)
        Base.depwarn("ZeroTo is deprecated, use ZeroRange instead", :ZeroTo)
        new(max(T(-1), stop))
    end
end
ZeroTo{T<:Signed}(stop::T) = ZeroTo{T}(stop)

Base.length(r::ZeroTo) = r.stop+1

Base.length{T<:Union{Int,Int64}}(r::ZeroTo{T}) = T(r.stop+1)

let smallint = (Int === Int64 ?
                Union{Int8,UInt8,Int16,UInt16,Int32,UInt32} :
                Union{Int8,UInt8,Int16,UInt16})
    Base.length{T <: smallint}(r::ZeroTo{T}) = Int(r.stop)+1
    Base.start{T<:smallint}(r::ZeroTo{T}) = 0
end

Base.first{T}(r::ZeroTo{T}) = zero(T)

Base.start{T}(r::ZeroTo{T}) = zero(T)

@inline function Base.getindex{T}(v::ZeroTo{T}, i::Integer)
    @boundscheck ((i > 0) & (i <= length(v))) || Base.throw_boundserror(v, i)
    convert(T, i-1)
end

@inline function Base.getindex{T}(r::ZeroTo{T}, s::ZeroTo)
    @boundscheck checkbounds(r, s)
    ZeroTo(T(s.stop))
end

@inline function Base.getindex{T}(r::ZeroTo{T}, s::AbstractUnitRange)
    @boundscheck checkbounds(r, s)
    T(first(s)-1):T(last(s)-1)
end

Base.intersect(r::ZeroTo, s::ZeroTo) = ZeroTo(min(r.stop,s.stop))

Base.promote_rule{T1,T2}(::Type{ZeroTo{T1}},::Type{ZeroTo{T2}}) =
    ZeroTo{promote_type(T1,T2)}
Base.convert{T<:Real}(::Type{ZeroTo{T}}, r::ZeroTo{T}) = r
Base.convert{T<:Real}(::Type{ZeroTo{T}}, r::ZeroTo) = ZeroTo{T}(r.stop)

Base.show(io::IO, r::ZeroTo) = print(io, typeof(r).name, "(", r.stop, ")")


text:
"""
outparser: CAMEOS output parser
"""

using JLD # Should be placed here too to avoid later crash loading JLD file

module main

include("optimize.jl")
include("types.jl")

using ArgParse, JLD

function parse_commandline()
	s = ArgParseSettings()
	@add_arg_table s begin
        "mark_gene"
		    help = "1st positional argument, name of gene 'mark'"
            arg_type = String
			default = "cysJ"
        "deg_gene"
		    help = "2nd positional argument, name of gene 'deg'"
            arg_type = String
			default = "infA"
		"runid"
			help = "3rd positional argument, run id code created by CAMEOS"
            arg_type = String
            default = "LgolUs7k"
        "--frame"
			help = "optional argument, CAMEOS frame"
            arg_type = String
		    default = "p1"
        "--fasta"
            help = "generate a FASTA file too"
            action = :store_true
        "--just-fullseq"
            help = "the FASTA will only contain the full sequence"
            action = :store_true
	end
	return parse_args(s)
end

function outparse_cameos()
    println("=-= CAMEOS output parser =-= v0.4 - Mar 2021 =-= by LLNL =-=")

    # Parse arguments
   	parsed_args = parse_commandline()
	mark_gene = parsed_args["mark_gene"]
	deg_gene = parsed_args["deg_gene"]
	runid = parsed_args["runid"]
	frame = parsed_args["frame"]
    fasta = parsed_args["fasta"]
    justfullseq = parsed_args["just-fullseq"]

    # Get complete path for input and output files
    subdir = string(mark_gene, "_", deg_gene, "_", frame)
    jld_file = string("saved_pop_", runid, ".jld")
    csv_file = string("summary_", runid, ".csv")
    fa_file = string("variants_", runid, ".fasta")
    in_path = string("output/", subdir, "/", jld_file)
    out_path = string("output/", subdir, "/", csv_file)
    fa_path = string("output/", subdir, "/", fa_file)

    # Load JLD file
    print("Loading data from JLD file ", in_path, " ... ")
    variants = load(in_path)["variants"];
    println("OK!")

    # Save CSV file
    print("Saving data to CSV file ", out_path, " : ")
    parsed = 0
    open(out_path, "w") do io
        println(io, string("full_seq,",
                           mark_gene,"_psls,",
                           mark_gene,"_seq,",
                           deg_gene,"_psls,",
                           deg_gene,"_seq"))
        for var in variants
            println(io, var.full_sequence,",",
                    var.mark_prob,",",
                    var.mark_seq,",", #var.mark_nuc,",",
                    var.deg_prob,",",
                    var.deg_seq) #,",",var.deg_nuc)
            parsed += 1
            if parsed % 200 == 0
                print(".")
            end
        end
    end
    println(" OK!")
    println(parsed, " variants parsed for ", runid)

    # Save fasta file
    if fasta
        print("Saving data to FASTA file ", fa_path, " : ")
        variant = 1
        open(fa_path, "w") do io

            function print2fasta(var, field::String, nucs::Bool=false)
                # Prints header and sequence to fasta file
                println(io, ">CAMEOS run:", runid,
                        " mark_gene:", mark_gene,
                        " deg_gene:", deg_gene,
                        " variant:", variant, "/", parsed,
                        " ", string(field))
                if nucs
                    println(io, getfield(getfield(var, Symbol(field)), :nucs))
                else
                    println(io, getfield(var, Symbol(field)))
                end
                return nothing
            end

            for var in variants
                if justfullseq
                    print2fasta(var, "full_sequence")  # equal to mark_nuc.nucs
                else
                    print2fasta(var, "mark_seq")
                    print2fasta(var, "mark_nuc", true)
                    print2fasta(var, "deg_seq")
                    print2fasta(var, "deg_nuc", true)
                end
                variant += 1
                if variant % 200 == 0
                    print(".")
                end
            end
        end
        println(" OK!")
        println(variant-1, " variants in FASTA file for ", runid)
    end
end

@time outparse_cameos()
end

max_stars_repo_name: khyox/CAMEOS

max_issues_repo_name: BiosecSFA/CAMEOS

max_forks_repo_name: BiosecSFA/CAMEOS


text:
import Colors
export lsimplot, stepplot, impulseplot, bodeplot, nyquistplot, sigmaplot, marginplot, setPlotScale, gangoffour, gangoffourplot, gangofseven, pzmap, pzmap!, nicholsplot

getColorSys(i,Nsys)   = convert(Colors.RGB,Colors.HSV(360*((i-1)/Nsys)^1.5,0.9,0.8))
function getStyleSys(i,Nsys)
    Ncmax = 5
    if Nsys <= Ncmax
        return Dict(:c => convert(Colors.RGB,Colors.HSV(360*((i-1)/Nsys)^1.5,0.9,0.8)) , :l => :solid)
    end
    styles = [:solid, :dash, :dot, :dashdot]
    Nstyles = min(ceil(Int,Nsys/Ncmax), length(styles))
    Nc = ceil(Int,Nsys / Nstyles)
    istyle = min(floor(Int,i/Nc)+1 , length(styles))
    c = convert(Colors.RGB,Colors.HSV(360*((mod(i-1,Nc)/Nc))^1.5,0.9,0.8))

    return Dict(:c => c, :l => styles[istyle])
end

_PlotScale = "dB"
_PlotScaleFunc = :identity
_PlotScaleStr = "(dB)"

@doc """`setPlotScale(str)`

Set the default scale of magnitude in `bodeplot` and `sigmaplot`.
`str` should be either `"dB"` or `"log10"`.""" ->
function setPlotScale(str::AbstractString)
    if str == "dB"
        plotSettings = (str, :identity, "(dB)")
    elseif str == "log10"
        plotSettings = (str, :log10, "")
    else
        error("Scale must be set to either \"dB\" or \"log10\"")
    end
    global _PlotScale, _PlotScaleFunc, _PlotScaleStr
    _PlotScale, _PlotScaleFunc, _PlotScaleStr = plotSettings
end

@doc """`fig = lsimplot(sys::StateSpace, u, t[, x0, method]; kwargs...), lsimplot(sys::TransferFunction,u,t[,method]; kwargs...)`

`lsimplot(StateSpace[sys1, sys2...], u, t[, x0, method]; kwargs...), lsimplot(TransferFunction[sys1, sys2...], u, t[, method]; kwargs...)`

Calculate the time response of the `LTISystem`(s) to input `u`. If `x0` is
ommitted, a zero vector is used.

Continuous time systems are discretized before simulation. By default, the
method is chosen based on the smoothness of the input signal. Optionally, the
`method` parameter can be specified as either `:zoh` or `:foh`.

`kwargs` is sent as argument to Plots.plot.""" ->
function lsimplot{T<:StateSpace}(systems::Vector{T}, u::Union{AbstractVecOrMat,Function},
    t::AbstractVector, x0::VecOrMat=zeros(systems[1].nx, 1),
    method::Symbol=_issmooth(u) ? :foh : :zoh; kwargs...)
    if !_same_io_dims(systems...)
        error("All systems must have the same input/output dimensions")
    end
    ny, nu = size(systems[1])
    fig = Plots.plot(layout=(ny,1))
    s2i(i,j) = sub2ind((ny,1),j,i)
    for (si, s) in enumerate(systems)
        y = lsim(s, u, t, x0, method)[1]
        for i=1:ny
            ydata = reshape(y[:, i], size(t, 1))
            style = iscontinuous(s) ? :path : :steppost
            ytext = (ny > 1) ? "Amplitude to: y($i)": "Amplitude"
            Plots.plot!(fig, t, ydata, l=style, xlabel="Time (s)", ylabel=ytext, title="System Response", subplot=s2i(1,i), lab="\$G_\{$(si)\}\$"; getStyleSys(si,length(systems))..., kwargs...)
        end
    end
    return fig
end
lsimplot(sys::LTISystem, u::Union{AbstractVecOrMat,Function}, t::AbstractVector, args...; kwargs...) =
    lsimplot(StateSpace[sys], u, t, args...; kwargs...)
lsimplot{T<:LTISystem}(sys::Vector{T}, u::Union{AbstractVecOrMat,Function}, t::AbstractVector, args...; kwargs...) =
    lsimplot(StateSpace[s for s in sys], u, t, args...; kwargs...)

for (func, title) = ((:step, "Step Response"), (:impulse, "Impulse Response"))
    funcname = Symbol("$(func)plot")
    @eval begin
        function $funcname{T<:LTISystem}(systems::Vector{T}, Ts_list::Vector, Tf::Real; kwargs...)
            if !_same_io_dims(systems...)
                error("All systems must have the same input/output dimensions")
            end
            ny, nu = size(systems[1])
            fig = Plots.plot(layout=(ny,nu))
            titles = fill("", 1, ny*nu)
            s2i(i,j) = sub2ind((ny,nu),j,i)
            for (si,(s, Ts)) in enumerate(zip(systems, Ts_list))
                t = 0:Ts:Tf
                y = ($func)(s, t)[1]
                for i=1:ny
                    for j=1:nu
                        ydata = reshape(y[:, i, j], size(t, 1))
                        style = iscontinuous(s) ? :path : :steppost
                        ttext = (nu > 1 && i==1) ? $title*" from: u($j) " : $title
                        titles[s2i(i,j)] = ttext
                        ytext = (ny > 1 && j==1) ? "Amplitude to: y($i)": "Amplitude"
                        Plots.plot!(fig, t, ydata, l=style, xlabel="Time (s)", ylabel=ytext, subplot=s2i(i,j), lab="\$G_\{$(si)\}\$"; getStyleSys(si,length(systems))..., kwargs...)
                    end
                end
            end
            Plots.plot!(fig, title=titles)
            return fig
        end
        $funcname{T<:LTISystem}(systems::Vector{T}, Tf::Real; kwargs...) =
            $funcname(systems, map(_default_Ts, systems), Tf; kwargs...)
        $funcname{T<:LTISystem}(systems::Vector{T}; kwargs...) =
            $funcname(systems, _default_time_data(systems)...; kwargs...)
        $funcname{T<:LTISystem}(systems::Vector{T}, t::AbstractVector; kwargs...) =
            $funcname(systems, repmat([t[2] - t[1]], length(systems)), t[end]; kwargs...)
        $funcname(sys::LTISystem, args...; kwargs...) = $funcname(LTISystem[sys], args...; kwargs...)
    end
end

@doc """`fig = stepplot(sys, args...)`, `stepplot(LTISystem[sys1, sys2...], args...)`

Plot the `step` response of the `LTISystem`(s). A final time `Tf` or a
time vector `t` can be optionally provided.""" -> stepplot

@doc """`fig = impulseplot(sys, args...)`, `impulseplot(LTISystem[sys1, sys2...], args...)`

Plot the `impulse` response of the `LTISystem`(s). A final time `Tf` or a
time vector `t` can be optionally provided.""" -> impulseplot

## FREQUENCY PLOTS ##
@doc """`fig = bodeplot(sys, args...)`, `bodeplot(LTISystem[sys1, sys2...], args...; plotphase=true, kwargs...)`

Create a Bode plot of the `LTISystem`(s). A frequency vector `w` can be
optionally provided.

`kwargs` is sent as argument to Plots.plot.""" ->
function bodeplot{T<:LTISystem}(systems::Vector{T}, w::AbstractVector; plotphase=true, kwargs...)
    if !_same_io_dims(systems...)
        error("All systems must have the same input/output dimensions")
    end
    ny, nu = size(systems[1])
    fig = Plots.plot(layout=((plotphase?2:1)*ny,nu))
    s2i(i,j) = sub2ind((nu,(plotphase?2:1)*ny),j,i)
    nw = length(w)
    for (si,s) = enumerate(systems)
        mag, phase = bode(s, w)[1:2]
        if _PlotScale == "dB"
            mag = 20*log10(mag)
        end
        xlab = plotphase ? "" : "Frequency (rad/s)"
        for j=1:nu
            for i=1:ny
                magdata = vec(mag[:, i, j])
                if all(magdata .== -Inf)
                    # 0 system, don't plot anything
                    continue
                end
                phasedata = vec(phase[:, i, j])
                Plots.plot!(fig, w, magdata, grid=true, yscale=_PlotScaleFunc, xscale=:log10, xlabel=xlab, subplot=s2i((plotphase?(2i-1):i),j), title="Bode plot from: u($j)", ylabel="Magnitude $_PlotScaleStr", lab="\$G_\{$(si)\}\$"; getStyleSys(si,length(systems))..., kwargs...)
                plotphase && Plots.plot!(fig, w, phasedata, grid=true, xscale=:log10, ylabel="Phase (deg)", subplot=s2i(2i,j), xlabel="Frequency (rad/s)", lab="\$G_\{$(si)\}\$"; getStyleSys(si,length(systems))..., kwargs...)
            end
        end
    end

    return fig
end
bodeplot{T<:LTISystem}(systems::Vector{T}; plotphase=true, kwargs...) =
    bodeplot(systems, _default_freq_vector(systems, :bode); plotphase=plotphase, kwargs...)
bodeplot(sys::LTISystem, args...; plotphase=true, kwargs...) = bodeplot([sys], args...; plotphase=plotphase, kwargs...)

@doc """`fig = nyquistplot(sys; kwargs...)`, `nyquistplot(LTISystem[sys1, sys2...]; kwargs...)`

Create a Nyquist plot of the `LTISystem`(s). A frequency vector `w` can be
optionally provided.

`kwargs` is sent as argument to Plots.plot.""" ->
function nyquistplot{T<:LTISystem}(systems::Vector{T}, w::AbstractVector; neg=false, kwargs...)
    if !_same_io_dims(systems...)
        error("All systems must have the same input/output dimensions")
    end
    ny, nu = size(systems[1])
    nw = length(w)
    fig = Plots.plot(layout=(ny,nu))
    s2i(i,j) = sub2ind((ny,nu),j,i)
    # Ensure that `axes` is always a matrix of handles
    for (si,s) = enumerate(systems)
        re_resp, im_resp = nyquist(s, w)[1:2]
        for j=1:nu
            for i=1:ny
                redata = re_resp[:, i, j]
                imdata = im_resp[:, i, j]
                ylim = (min(max(-20,minimum(imdata)),-1), max(min(20,maximum(imdata)),1))
                xlim = (min(max(-20,minimum(redata)),-1), max(min(20,maximum(redata)),1))
                Plots.plot!(fig,redata, imdata, title="Nyquist plot from: u($j)", ylabel="To: y($i)", ylims=ylim, xlims=xlim, subplot=s2i(i,j), lab="\$G_\{$(si)\}\$"; getStyleSys(si,length(systems))..., kwargs...)

                if si == length(systems)
                    v = linspace(0,2π,100)
                    S,C = sin(v),cos(v)
                    Plots.plot!(fig,C,S,l=:dash,c=:black, lab="", subplot=s2i(i,j))
                    Plots.plot!(fig,C-1,S,l=:dash,c=:red, grid=true, lab="", subplot=s2i(i,j))
                    # neg && Plots.plot!(fig[i, j],redata, -imdata, args...)
                end
            end
        end
    end
    return fig
end

nyquistplot{T<:LTISystem}(systems::Vector{T}; kwargs...) =
    nyquistplot(systems, _default_freq_vector(systems, :nyquist); kwargs...)
nyquistplot(sys::LTISystem, args...; kwargs...) = nyquistplot([sys], args...; kwargs...)

@doc """
fig = `nicholsplot{T<:LTISystem}(systems::Vector{T}, w::AbstractVector; kwargs...)`

Create a Nichols plot of the `LTISystem`(s). A frequency vector `w` can be
optionally provided.

Keyword arguments:
```
text = true
Gains = [12, 6, 3, 1, 0.5, -0.5, -1, -3, -6, -10, -20, -40, -60]
pInc = 30
sat = 0.4
val = 0.85
fontsize = 10
```

`pInc` determines the increment in degrees between phase lines.

`sat` ∈ [0,1] determines the saturation of the gain lines

`val` ∈ [0,1] determines the brightness of the gain lines

Additional keyword arguments are sent to the function plotting the systems and can be
used to specify colors, line styles etc. using regular Plots.jl syntax

This function is based on code subject to the two-clause BSD licence
Copyright 2011 Will Robertson
Copyright 2011 Philipp Allgeuer

""" ->
function nicholsplot{T<:LTISystem}(systems::Vector{T}, w::AbstractVector;
    text=true,
    Gains = [12, 6, 3, 1, 0.5, -0.5, -1, -3, -6, -10, -20, -40, -60],
    pInc = 30,
    sat = 0.4,
    val = 0.85,
    fontsize = 10,
    kwargs...)
    if !_same_io_dims(systems...)
        error("All systems must have the same input/output dimensions")
    end
    ny, nu = size(systems[1])

    if !iscontinuous(systems[1])
        w_nyquist = 2π/systems[1].Ts
        w = w[w.<= w_nyquist]
    end
    nw = length(w)

    # Gain circle functions
    angle(x)        = unwrap(atan2(imag(x),real(x)))
    RadM(m)         = abs(m/(m^2-1))
    CentreM(m)      = m^2/(1-m^2)
    Ny(mdb,t)       = CentreM(10^(mdb/20))+RadM(10^(mdb/20)).*(cosd(t)+im.*sind(t))
    Niϕ(mdb,t)      = rad2deg((angle(Ny(mdb,t))))
    Ni_Ga(mdb,t)    = 20.*log10(abs(Ny(mdb,t)))

    # Phase circle functions
    Radϕ(ϕ)         = 1./(2.*abs(sind(ϕ)))
    Nyℜ(ϕ,t)        = -0.5+Radϕ(ϕ).*cosd(t+mod(ϕ,180)-90)
    Nyℑ(ϕ,t)        = 1./(2.*tand(ϕ))+Radϕ(ϕ).*sind(t+mod(ϕ,180)-90)
    Niϕϕ(ϕ,t)       = rad2deg((angle(Nyℜ(ϕ,t)+im*Nyℑ(ϕ,t))))+360*floor(ϕ/360)
    Ni_Gaϕ(ϕ,t)     = 20.*log10(abs(Nyℜ(ϕ,t)+im*Nyℑ(ϕ,t)))
    Ni_La(ϕ)        = 0.090*10^(ϕ/60)
    getColor(mdb)   = convert(Colors.RGB,Colors.HSV(360*((mdb-minimum(Gains))/(maximum(Gains)-minimum(Gains)))^1.5,sat,val))

    fig             = Plots.plot()
    megaangles      = vcat(map(s -> 180/π*angle(squeeze(freqresp(s, w)[1],(2,3))), systems)...)
    filter!(x-> !isnan(x), megaangles)
    PCyc            = Set{Int}(floor(Int,megaangles/360))
    PCyc            = sort(collect(PCyc))

    #  Gain circles
    for k=Gains
        ϕVals   =Niϕ(k,0:0.1:360)
        GVals   =Ni_Ga(k,0:0.1:360)
        for l in PCyc
            Plots.plot!(fig,ϕVals+l*360,GVals,c=getColor(k), grid=false)
            if text
                offset  = (l+1)*360
                TextX   = Niϕ(k,210)+offset
                TextY   = Ni_Ga(k,210)
                Plots.plot!(fig,ann=(TextX,TextY,Plots.text("$(string(k)) dB",fontsize)))
            end
        end
    end

    #  Phase circles
    Phi=PCyc[1]*360:pInc:PCyc[end]*360
    T1=logspace(-4,log10(180),300)
    T2=[T1; 360-flipdim(T1,1)]

    for k=Phi
        if abs(sind(k))<1e-3
            Plots.plot!(fig,[k,k],[-110,25],c=Colors.RGB(0.75*[1, 1, 1]...))
            if cosd(5)>0
                TextX=k
                TextY=1
            else
                TextX=k
                TextY=-46.5
            end
        else
            Plots.plot!(fig,Niϕϕ(k,T2),Ni_Gaϕ(k,T2),c=Colors.RGB(0.75*[1,1,1]...))
            Offset=k-180*floor(Int,k/180);
            if sign(sind(k))==1
                TextX=Niϕϕ(k,Ni_La(180-Offset))
                TextY=Ni_Gaϕ(k,Ni_La(180-Offset))
            else
                TextX=Niϕϕ(k,-Ni_La(Offset))+360
                TextY=Ni_Gaϕ(k,-Ni_La(Offset))
            end
        end
        if text
            Plots.plot!(fig,ann=(TextX,TextY,Plots.text("$(string(k))°",fontsize)))
        end

        Plots.plot!(fig, title="Nichols chart", grid=false, legend=false)

    end
    dKwargs = Dict(kwargs)
    LW = "linewidth" ∈ keys(dKwargs) ? pop!(dKwargs,"linewidth") : 2

    # colors = [:blue, :cyan, :green, :yellow, :orange, :red, :magenta]
    for (sysi,s) = enumerate(systems)
        ℜresp, ℑresp        = nyquist(s, w)[1:2]
        ℜdata               = squeeze(ℜresp, (2,3))
        ℑdata               = squeeze(ℑresp, (2,3))
        mag                 = 20*log10(sqrt(ℜdata.^2 + ℑdata.^2))
        angles              = 180/π*angle(im*ℑdata.+ℜdata)
        Plots.plot!(fig,angles, mag; linewidth = LW, getStyleSys(sysi,length(systems))..., kwargs...)
    end

    return fig
end

nicholsplot{T<:LTISystem}(systems::Vector{T};kwargs...) =
    nicholsplot(systems, _default_freq_vector(systems, :nyquist);kwargs...)
nicholsplot(sys::LTISystem, args...; kwargs...) = nicholsplot([sys],args...; kwargs...)

@doc """`sigmaplot(sys, args...)`, `sigmaplot(LTISystem[sys1, sys2...], args...)`

Plot the singular values of the frequency response of the `LTISystem`(s). A
frequency vector `w` can be optionally provided.

`kwargs` is sent as argument to Plots.plot.""" ->
function sigmaplot{T<:LTISystem}(systems::Vector{T}, w::AbstractVector; kwargs...)
    if !_same_io_dims(systems...)
        error("All systems must have the same input/output dimensions")
    end
    ny, nu = size(systems[1])
    nw = length(w)
    fig = Plots.plot()
    for (si, s) in enumerate(systems)
        sv = sigma(s, w)[1]
        if _PlotScale == "dB"
            sv = 20*log10(sv)
        end
        for i in 1:size(sv, 2)
            Plots.plot!(fig, w, sv[:, i], xscale=:log10, yscale=_PlotScaleFunc; getStyleSys(si,length(systems))..., kwargs...)
        end
    end
    Plots.plot!(fig, title="Sigma Plot", xlabel="Frequency (rad/s)",
        ylabel="Singular Values $_PlotScaleStr")
    return fig
end
sigmaplot{T<:LTISystem}(systems::Vector{T}; kwargs...) =
    sigmaplot(systems, _default_freq_vector(systems, :sigma); kwargs...)
sigmaplot(sys::LTISystem, args...; kwargs...) = sigmaplot([sys], args...; kwargs...)

@doc """`fig = marginplot(sys::LTISystem [,w::AbstractVector];  kwargs...)`, `marginplot(sys::Vector{LTISystem}, w::AbstractVector;  kwargs...)`

Plot all the amplitude and phase margins of the system(s) `sys`.
A frequency vector `w` can be optionally provided.

`kwargs` is sent as argument to Plots.plot.""" ->
function marginplot{T<:LTISystem}(systems::Vector{T}, w::AbstractVector; kwargs...)
    if !_same_io_dims(systems...)
        error("All systems must have the same input/output dimensions")
    end
    ny, nu = size(systems[1])
    fig = bodeplot(systems, w, kwargs...)
    s2i(i,j) = sub2ind((ny,nu),j,i)
    titles = Array(AbstractString,nu,ny,2,2)
    titles[:,:,1,1] = "Gm: "
    titles[:,:,2,1] = "Pm: "
    titles[:,:,1,2] = "Wgm: "
    titles[:,:,2,2] = "Wpm: "
    for (si, s) in enumerate(systems)
        for j=1:nu
            for i=1:ny
                wgm, gm, wpm, pm, fullPhase = sisomargin(s[i,j],w, full=true, allMargins=true)
                if _PlotScale == "dB"
                    mag = 20*log10(1./gm)
                    oneLine = 0
                else
                    mag = 1./gm
                    oneLine = 1
                end
                for k=1:length(wgm)
                    #Plot gain margins
                    Plots.plot!(fig, [wgm[k];wgm[k]], [1;mag[k]], lab="", subplot=s2i(2i-1,j); getStyleSys(si,length(systems))...)
                end
                #Plot gain line at 1
                Plots.plot!(fig, [w[1],w[end]], [oneLine,oneLine], l=:dash, c=:gray, lab="", subplot=s2i(2i-1,j))
                titles[j,i,1,1] *= "["*join([@sprintf("%2.2f",v) for v in gm],", ")*"] "
                titles[j,i,1,2] *= "["*join([@sprintf("%2.2f",v) for v in wgm],", ")*"] "
                for k=1:length(wpm)
                    #Plot the phase margins
                    Plots.plot!(fig, [wpm[k];wpm[k]],[fullPhase[k];fullPhase[k]-pm[k]], lab="", subplot=s2i(2i,j); getStyleSys(si,length(systems))...)
                    #Plot the line at 360*k
                    Plots.plot!(fig, [w[1],w[end]],(fullPhase[k]-pm[k])*ones(2), l=:dash, c=:gray, lab="", subplot=s2i(2i,j))
                end
                titles[j,i,2,1] *=  "["*join([@sprintf("%2.2f",v) for v in pm],", ")*"] "
                titles[j,i,2,2] *=  "["*join([@sprintf("%2.2f",v) for v in wpm],", ")*"] "
            end
        end
    end
    for j = 1:nu
        for i = 1:ny
            Plots.title!(fig, titles[j,i,1,1]*" "*titles[j,i,1,2], subplot=s2i(2i-1,j))
            Plots.title!(fig, titles[j,i,2,1]*" "*titles[j,i,2,2], subplot=s2i(2i,j))
        end
    end
    return fig
end
marginplot{T<:LTISystem}(systems::Vector{T}; kwargs...) =
    marginplot(systems, _default_freq_vector(systems, :bode); kwargs...)
marginplot(sys::LTISystem, args...; kwargs...) = marginplot([sys], args...; kwargs...)

# HELPERS:

function _same_io_dims(systems::LTISystem...)
    sizes = map(size, systems)
    return all(s -> s == sizes[1], sizes)
end

function _default_time_data{T<:LTISystem}(systems::Vector{T})
    sample_times = [_default_Ts(i) for i in systems]
    Tf = 100*maximum(sample_times)
    return sample_times, Tf
end
_default_time_data(sys::LTISystem) = _default_time_data(LTISystem[sys])

@doc """`fig = pzmap!(fig, system, args...; kwargs...)`

Create a pole-zero map of the `LTISystem`(s) in figure `fig`, `args` and `kwargs` will be sent to the `scatter` plot command.""" ->
function pzmap!(fig, system::LTISystem, args...; kwargs...)
    if system.nu + system.ny > 2
        warn("pzmap currently only supports SISO systems. Only transfer function from u₁ to y₁ will be shown")
    end

    z,p,k = zpkdata(system)
    !isempty(z[1]) && Plots.scatter!(fig,real(z[1]),imag(z[1]),m=:c,markersize=15., markeralpha=0.5, args...; kwargs...)
    !isempty(p[1]) && Plots.scatter!(fig,real(p[1]),imag(p[1]),m=:x,markersize=15., args...; kwargs...)
    Plots.title!("Pole-zero map")

    if system.Ts > 0
        v = linspace(0,2π,100)
        S,C = sin(v),cos(v)
        Plots.plot!(fig,C,S,l=:dash,c=:black, grid=true)
    end
    Plots.plot!(fig,legend=false)

    return fig
end

@doc """`fig = pzmap(system, args...; kwargs...)`

Create a pole-zero map of the `LTISystem`(s), `args` and `kwargs` will be sent to the `scatter` plot command.""" ->
pzmap(system::LTISystem, args...; kwargs...) = pzmap!(Plots.plot(), system, args...; kwargs...)

@doc """`fig = gangoffourplot(P::LTISystem, C::LTISystem)`, `gangoffourplot(P::Union{Vector, LTISystem}, C::Vector; plotphase=false)`

Gang-of-Four plot.

`kwargs` is sent as argument to Plots.plot.""" ->
function gangoffourplot(P::Union{Vector, LTISystem}, C::Vector, args...; plotphase=false, kwargs...)
    S,D,N,T = gangoffour(P,C)
    fig = bodeplot(LTISystem[[S[i] D[i]; N[i] T[i]] for i = 1:length(C)], args..., plotphase=plotphase; kwargs...)
    lower = plotphase ? 3 : 2
    s2i(i,j) = sub2ind((2,(plotphase?4:2)),j,i)
    Plots.plot!(fig,title="\$S = 1/(1+PC)\$", subplot=s2i(1,1))
    Plots.plot!(fig,title="\$D = P/(1+PC)\$", subplot=s2i(1,2))
    Plots.plot!(fig,title="\$N = C/(1+PC)\$", subplot=s2i(lower,1))
    Plots.plot!(fig,title="\$T = PC/(1+PC\$)", subplot=s2i(lower,2))
    return fig
end

function gangoffourplot(P::LTISystem,C::LTISystem, args...; plotphase=false, kwargs...)
    gangoffourplot(P,[C], args...; plotphase=plotphase, kwargs...)
end


text:
lattice = load_lattice(joinpath(@__DIR__, "CrFe2Ga_qe.in"))
atoms   = load_atoms(  joinpath(@__DIR__, "CrFe2Ga_qe.in"))
atoms   = attach_pseudos(atoms, Cr="hgh/pbe/cr-q14.hgh",
                                Fe="hgh/pbe/fe-q16.hgh",
                                Ga="hgh/pbe/ga-q3.hgh")
magnetic_moments = build_magnetic_moments(atoms; Ga=0.0, Fe=5.0, Cr=5.0)

cases = Dict(
    "standard_02" => (damping=DFTK.FixedDamping(0.3), ),
)

config[@__DIR__] = (
    directory        =  @__DIR__,
    lattice          = lattice,
    atoms            = atoms,
    magnetic_moments = magnetic_moments,
    #
    kwargs_model     = (smearing=Smearing.Gaussian(), temperature=0.01),
    Ecut             = 45,
    kgrid            = [13, 13, 13],
    supersampling    = 2.0,
    potential        = nothing,
    kwargs_scf       = (acceleration=DFTK.AndersonAcceleration(; maxcond=1e2), ),
    mixing           = KerkerMixing(),
    cases            = cases,
)

max_stars_repo_name: mfherbst/supporting-adaptive-damping

max_issues_repo_name: mfherbst/supporting-adaptive-damping

max_forks_repo_name: mfherbst/supporting-adaptive-damping


text:
using Media
import REPL

using Logging: with_logger
using .Progress: JunoProgressLogger

ends_with_semicolon(x) = REPL.ends_with_semicolon(split(x,'\n',keepempty = false)[end])

withpath(f, path) =
  CodeTools.withpath(f, path == nothing || isuntitled(path) ? nothing : path)

const evallock = ReentrantLock()

handle("evalshow") do data
  fixjunodisplays()
  @dynamic let Media.input = Editor()
    @destruct [text, line, path, mod] = data
    mod = getmodule(mod)

    lock(evallock)
    result = hideprompt() do
      with_logger(JunoProgressLogger()) do
        withpath(path) do
          res = @errs include_string(mod, text, path, line)

          if res isa EvalError
              Base.showerror(IOContext(stderr, :limit => true), res)
          elseif res !== nothing && !ends_with_semicolon(text)
            display(res)
          end
        end
      end
    end
    unlock(evallock)

    Base.invokelatest() do
      display = Media.getdisplay(typeof(result), Media.pool(Editor()), default = Editor())
      !isa(result, EvalError) && ends_with_semicolon(text) && (result = nothing)
      display ≠ Editor() && result !== nothing && @ierrs render(display, result)
    end

    nothing
  end
end

handle("eval") do data
  fixjunodisplays()
  @dynamic let Media.input = Editor()
    @destruct [text, line, path, mod, errorInRepl || false] = data
    mod = getmodule(mod)

    lock(evallock)
    result = hideprompt() do
      with_logger(JunoProgressLogger()) do
        withpath(path) do
          res = @errs include_string(mod, text, path, line)
          if errorInRepl && res isa EvalError
            try
              Base.showerror(IOContext(stderr, :limit => true), res)
            catch err
              show(stderr, err)
            end
          end
          return res
        end
      end
    end
    unlock(evallock)

    Base.invokelatest() do
      !isa(result, EvalError) && ends_with_semicolon(text) && (result = nothing)
      @ierrs displayandrender(result)
    end
  end
end

# dummy handler for Revise compat
handle("evalrepl") do data
  @warn "Juno's evalrepl handler is deprecated."
end

handle("evalall") do data
  fixjunodisplays()
  @dynamic let Media.input = Editor()
    @destruct [setmod = :module || nothing, path || "untitled", code] = data
    mod = if setmod !== nothing
       getmodule(setmod)
    elseif isabspath(path)
      getmodule(CodeTools.filemodule(path))
    else
      Main
    end

    lock(evallock)
    hideprompt() do
      with_logger(JunoProgressLogger()) do
        withpath(path) do
          result = nothing
          try
            result = include_string(mod, code, path)
          catch err
            ee = EvalError(err, stacktrace(catch_backtrace()))

            # show error in REPL:
            Base.showerror(IOContext(stderr, :limit => true), ee)
            # show notification (if enabled in Atom):
            @msg error(d(:msg => "Error evaluating $(basename(path))",
                         :detail => string(ee),
                         :dismissable => true))
          end

          Base.invokelatest() do
            @ierrs displayandrender(result)
          end
        end
      end
    end
    unlock(evallock)
  end
  return
end

handle("docs") do data
  @destruct [
    word,
    mod || "Main"
  ] = data
  docs(word, mod)
end

function docs(word, mod = "Main")
  docstring = @errs getdocs(mod, word)
  docstring isa EvalError && return Dict(:error => true)

  mtable = try getmethods(mod, word)
    catch e
      []
    end

  Dict(:error    => false,
       :type     => :dom,
       :tag      => :div,
       :contents =>  map(x -> render(Inline(), x), [docstring; mtable]))
end


text:
rand(d::UnivariateDistribution) = quantile(d, rand())
rand(rng::AbstractRNG, d::UnivariateDistribution) = quantile(d, rand(rng))

##### specific distributions #####

const discrete_distributions = [
    "categorical"
]

for dname in discrete_distributions
    include(joinpath("univariate", "discrete", "$(dname).jl"))
end


text:
# # An equation of second order in two dimmensions
#
# Let \\(\Omega = (0, 1)^2\\) the unit square and denote \\( \Gamma = \partial \Omega\\) its boundary.
# Consider the following problem. Given \\(f \in C(\Omega)\\), find a function \\(u\\) satisfying
# -\Delta u + u = f in \Omega
# u = 0 on \Gamma

using SymFEL
using SymPy
using LinearAlgebra
using SparseArrays
using WriteVTK
using JLD2

## discretization parameters
# we use the mesh square.msh (in gmsh format)
# obtained from square.geo using gmsh
# this mesh is formed by quad elements
# we use gmsh module for read the mesh
# there are some problems using Threads and gmsh (in linux)
# the file is square.jld2 is prepared separetely
@load "square.jld2"

# elementary matrices - P2 x P2
elem_Mxy = SymFEL.get_square_lagrange_em((2, 2), (0, 0), (0, 0))
elem_Kxy = SymFEL.get_square_lagrange_em((2, 2), (1, 0), (1, 0)) +
    SymFEL.get_square_lagrange_em((2, 2), (0, 1), (0, 1))

dx = norm(nodes_coordinate[:, elements_bound[1,1]] - nodes_coordinate[:, elements_bound[2,1]])

elem_Kxy_dx = convert(Matrix{Float64}, elem_Kxy.subs(h, dx))
elem_Mxy_dx = convert(Matrix{Float64}, elem_Mxy.subs(h, dx));

# global matrices
K = SymFEL.assemble_squaremesh_FE_matrix(elem_Kxy_dx, elements_int, order1=2, order2=2)
M = SymFEL.assemble_squaremesh_FE_matrix(elem_Mxy_dx, elements_int, order1=2, order2=2)

f = (2*pi^2 + 1) * sin.(pi * nodes_coordinate[1,:]) .* sin.(pi * nodes_coordinate[2,:])

F = M * f

A = K + M

# boundary condition
tgv = 1e30
A[nodes_boundary, nodes_boundary] += tgv * sparse(Matrix{Float64}(I, nodes_boundary_N, nodes_boundary_N))

u = A \ F
u_exact = sin.(pi*nodes_coordinate[1,:]) .* sin.(pi*nodes_coordinate[2,:])
err = u - u_exact

println("L2 error : ", sqrt(err' * M * err))
println("H1 error : ", sqrt(err' * K * err))

# export to vtk
cells = [MeshCell(VTKCellTypes.VTK_QUADRATIC_QUAD, elements_int[1:8, i]) for i = 1:elements_int_N]

points_x = nodes_coordinate[1, :]
points_y = nodes_coordinate[2, :]
vtkfile = vtk_grid("ex3-output", points_x, points_y, cells)

vtkfile["my_point_data", VTKPointData()] = u
outfiles = vtk_save(vtkfile)


text:
using Random

export pecuzal_embedding

"""
    pecuzal_embedding(s; kwargs...) → Y, τ_vals, ts_vals, ΔLs, ⟨ε★⟩
A unified approach to properly embed a timeseries or a set of timeseries
(`Dataset`) based on the recent PECUZAL algorithm due to Kraemer et al.[^Kraemer2021].

For more details, see the description below.

## Keyword arguments
* `τs = 0:50`: Possible delay values `τs` (in sampling time units). For each of
  the `τs`'s the continuity statistic ⟨ε★⟩ gets computed and further processed
  in order to find optimal delays `τᵢ` for each embedding cycle `i`.
* `w::Int = 0`: Theiler window (neighbors in time with index `w` close to the point,
  that are excluded from being true neighbors). `w=0` means to exclude only the
  point itself, and no temporal neighbors.
* `samplesize::Real = 1.0`: Fraction of state space points
  to be considered (fiducial points v) to average ε★ over, in order to
  produce `⟨ε★⟩`. Lower fraction value decreases accuracy as well as computation time.
* `K::Int = 13`: the amount of nearest neighbors in the δ-ball (read algorithm
  description). Must be at least 8 (in order to gurantee a valid statistic).
  `⟨ε★⟩` is computed taking the minimum result over all `k ∈ K`.
* `KNN::Int = 3`: the amount of nearest neighbors considered, in order to compute
  σ_k^2 (read algorithm description [`uzal_cost`]@ref). If given a vector, the
  minimum result over all `knn ∈ KNN` is returned.
* `L_threshold::Real = 0`: The algorithm breaks, when this threshold is exceeded
  by `ΔL` in an embedding cycle (set as a positive number, i.e. an absolute value
  of `ΔL`).
* `α::Real = 0.05`: The significance level for obtaining the continuity statistic
* `p::Real = 0.5`: The p-parameter for the binomial distribution used for the
  computation of the continuity statistic ⟨ε★⟩.
* `max_cycles = 50`: The algorithm will stop after that many cycles no matter what.
* `econ::Bool = false`: Economy-mode for L-statistic computation. Instead of
  computing L-statistics for time horizons `2:Tw`, here we only compute them for
  `2:2:Tw`, see description for further details.
* `verbose = true`: Print information about the process.

## Description
The method works iteratively and gradually builds the final embedding vectors
`Y`. Based on the `⟨ε★⟩`-statistic (of [`pecora`](@ref)) the algorithm picks an
optimal delay value `τᵢ` for each embedding cycle `i`.
For achieving that, we take the inpute time series `s`, denoted as the actual
phase space trajectory `Y_actual` and compute the continuity statistic `⟨ε★⟩`.

1. Each local maxima in `⟨ε★⟩` is used for constructing a
   candidate embedding trajectory `Y_trial` with a delay corresponding to that
   specific peak in `⟨ε★⟩`.
2. We then compute the `L`-statistic (of [`uzal_cost`](@ref))
   for `Y_trial` (`L-trial`) and `Y_actual` (`L_actual`) for increasing prediction
   time horizons (free parameter in the `L`-statistic) and save the maximum
   difference `max(L-trial - L_actual)` as `ΔL` (Note that this is a
   negative number, since the `L`-statistic decreases with better reconstructions).
3. We pick the `τ`-value, for which `ΔL` is minimal (=maximum decrease of
   the overall `L`-value) and construct the actual embedding trajectory
   `Y_actual` (steps 1.-3. correspond to an embedding cycle).
4. We repeat steps 1.-3. with `Y_actual` as input and stop the algorithm when `ΔL` is > 0,
   i.e. when and additional embedding component would not lead to a lower overall
   L-value. `Y_actual` -> `Y`.

In case of multivariate embedding, i.e. when embedding a set of M time series
(`s::Dataset`), in each embedding cycle the continuity statistic `⟨ε★⟩` gets
computed for all M time series available. The optimal delay value `τ` in each
embedding cycle is chosen as the peak/`τ`-value for which `ΔL` is
minimal under all available peaks and under all M `⟨ε★⟩`'s. In the first
embedding cycle there will be M² different `⟨ε★⟩`'s to consider, since it is
not clear a priori which time series of the input should consitute the first
component of the embedding vector and form `Y_actual`.

The range of considered delay values is determined in `τs` and for the
nearest neighbor search we respect the Theiler window `w`. The final embedding
vector is stored in `Y` (`Dataset`). The chosen delay values for each embedding
cycle are stored in `τ_vals` and the according time series numbers chosen for
each delay value in `τ_vals` are stored in `ts_vals`. For univariate embedding
(`s::Vector`) `ts_vals` is a vector of ones of length `τ_vals`, because there is
simply just one timeseries to choose from. The function also returns the
`ΔLs`-values for each embedding cycle and the continuity statistic `⟨ε★⟩`
as an `Array` of `Vector`s.

For distance computations the Euclidean norm is used.

[^Kraemer2021]: Kraemer, K.H., Datseris, G., Kurths, J., Kiss, I.Z., Ocampo-Espindola, Marwan, N. (2021). [A unified and automated approach to attractor reconstruction. New Journal of Physics 23(3), 033017](https://iopscience.iop.org/article/10.1088/1367-2630/abe336).
"""
function pecuzal_embedding(s::Vector{T}; τs = 0:50 , w::Int = 0,
    samplesize::Real = 1, K::Int = 13, KNN::Int = 3, L_threshold::Real = 0,
    α::Real = 0.05, p::Real = 0.5, max_cycles::Int = 50, econ::Bool = false,
    verbose = true,
    ) where {T<:Real}

    verbose && println("Initializing PECUZAL algorithm for univariate input...")
    @assert 0 < samplesize ≤ 1
    @assert all(x -> x ≥ 0, τs)
    @assert L_threshold ≥ 0
    threshold = -L_threshold # due to the negativity of L-decrease
    metric = Euclidean()

    s_orig = s
    s = standardize(s) # especially important for comparative L-statistics
    # define actual phase space trajectory
    Y_act = Dataset(s)

    # set a flag, in order to tell the while loop when to stop. Each loop
    # stands for encountering a new embedding dimension
    flag, counter = true, 1

    # preallocate output variables
    τ_vals = Int64[0]
    ts_vals = Int64[1]
    Ls = Float64[]
    ε★s = Array{T}(undef, length(τs), max_cycles)

    # loop over increasing embedding dimensions until some break criterion will
    # tell the loop to stop/break
    while flag
        verbose && println("Starting $counter-th embedding cycle...")
        Y_act = pecuzal_embedding_cycle!(
                Y_act, flag, s, τs, w, counter, ε★s, τ_vals, metric,
                Ls, ts_vals, samplesize, K, α, p, KNN, econ)
        flag = pecuzal_break_criterion(Ls, counter, max_cycles, threshold)
        counter += 1
    end
    # construct final reconstruction vector
    NN = (length(s)-sum(τ_vals[1:end-1]))
    Y_final = s_orig
    for i = 2:length(τ_vals[1:end-1])
        Y_final = DelayEmbeddings.hcat_lagged_values(Y_final, s_orig, τ_vals[i])
    end
    return Y_final, τ_vals[1:end-1], ts_vals[1:end-1], Ls[1:end-1], ε★s[:,1:counter-1]

end

function pecuzal_embedding(Y::Dataset{D, T}; τs = 0:50 , w::Int = 1,
    samplesize::Real = 1, K::Int = 13, KNN::Int = 3, L_threshold::Real = 0,
    α::Real = 0.05, p::Real = 0.5, max_cycles::Int = 50, econ::Bool = false,
    verbose = true,
    ) where {D, T<:Real}

    verbose && println("Initializing PECUZAL algorithm for multivariate input...")
    @assert 0 < samplesize ≤ 1
    @assert all(x -> x ≥ 0, τs)
    @assert L_threshold ≥ 0
    threshold = -L_threshold # due to the negativity of L-decrease
    metric = Euclidean()

    Y_orig = Y
    Y = standardize(Y) # especially important for comparative L-statistics

    # define actual phase space trajectory
    Y_act = []

    # set a flag, in order to tell the while loop when to stop. Each loop
    # stands for encountering a new embedding dimension
    flag, counter = true, 1

    # preallocate output variables
    τ_vals = Int64[0]
    ts_vals = Int64[]
    Ls = Float64[]
    ε★s = fill(zeros(T, length(τs), size(Y,2)), 1, max_cycles)

    # loop over increasing embedding dimensions until some break criterion will
    # tell the loop to stop/break
    while flag
        verbose && println("Starting $counter-th embedding cycle...")
        Y_act = pecuzal_multivariate_embedding_cycle!(
                Y_act, flag, Y, τs, w, counter, ε★s, τ_vals, metric,
                Ls, ts_vals, samplesize, K, α, p, KNN, econ)

        flag = pecuzal_break_criterion(Ls, counter, max_cycles, threshold)
        counter += 1
    end
    # construct final reconstruction vector
    Y_final = Y_orig[:,ts_vals[1]]
    for i = 2:length(τ_vals[1:counter-1])
        Y_final = DelayEmbeddings.hcat_lagged_values(Y_final,Y_orig[:,ts_vals[i]],τ_vals[i])
    end

    return Y_final, τ_vals[1:end-1], ts_vals[1:end-1], Ls[1:end-1], ε★s[:,1:counter-1]

end

"""
Perform one univariate embedding cycle on `Y_act`. Return the new `Y_act`
"""
function pecuzal_embedding_cycle!(Y_act::Dataset{D, T}, flag::Bool, s::Vector,
        τs, w::Int, counter::Int, ε★s::AbstractArray, τ_vals::Vector{Int}, metric,
        Ls::Vector{T}, ts_vals::Vector{Int}, samplesize::Real, K::Int, α::Real,
        p::Real, KNN::Int, econ::Bool) where {D, T}

    ε★, _ = pecora(s, Tuple(τ_vals), Tuple(ts_vals); delays = τs, w = w,
                samplesize = samplesize, K = K, metric = metric, α = α,
                p = p, undersampling = false)
    ε★s[:,counter] = ε★

    # zero-padding of ⟨ε★⟩ in order to also cover τ=0 (important for the multivariate case)
    ε★ = vec([0; ε★])
    # get the L-statistic-decrease for each peak in ⟨ε★⟩ and take the maximum one
    L_trials, max_idx = local_L_statistics(ε★, Y_act, s, τs, KNN, w, samplesize, metric, econ)
    L_min, min_idx = findmin(L_trials)

    push!(τ_vals, τs[max_idx[min_idx]-1])
    push!(ts_vals, 1)
    push!(Ls, L_min)

    # create phase space vector for this embedding cycle
    Y_act = DelayEmbeddings.hcat_lagged_values(Y_act, s, τ_vals[counter+1])

    return Y_act
end

"""
Perform one embedding cycle on `Y_act` with a multivariate set Ys
"""
function pecuzal_multivariate_embedding_cycle!(Y_act, flag::Bool,
        Ys::Dataset{DT, T}, τs, w::Int, counter::Int, ε★s::AbstractMatrix,
        τ_vals::Vector{Int}, metric, Ls::Vector{T}, ts_vals::Vector{Int},
        samplesize::Real, K::Int, α::Real, p::Real, KNN::Int, econ::Bool
        ) where {D, DT, T}

    M = size(Ys,2)
    # in the 1st cycle we have to check all (size(Y,2)^2 combinations and pick
    # the tau according to maximum L-statistic decrease)
    if counter == 1
        Y_act = first_embedding_cycle_pecuzal!(Ys, M, τs, w, samplesize, K,
                                metric, α, p, KNN, τ_vals, ts_vals, Ls, ε★s, econ)
    # in all other cycles we just have to check (size(Y,2)) combinations and pick
    # the tau according to maximum L-statistic decrease
    else
        Y_act = embedding_cycle_pecuzal!(Y_act, Ys, counter, M, τs, w, samplesize,
                            K, metric, α, p, KNN, τ_vals, ts_vals, Ls, ε★s, econ)
    end
    return Y_act
end

"""
Perform the first embedding cycle of the multivariate embedding. Return the
actual reconstruction vector `Y_act`.
"""
function first_embedding_cycle_pecuzal!(Ys::Dataset{D, T}, M::Int, τs, w::Int,
            samplesize::Real, K::Int, metric, α::Real, p::Real, KNN::Int,
            τ_vals::Vector{Int}, ts_vals::Vector{Int}, Ls::Vector{T},
            ε★s::AbstractMatrix, econ::Bool) where {D, T}
    counter = 1
    L_min = zeros(M)
    L_act = zeros(M)
    L_min_idx = zeros(Int, M)
    ε★ = zeros(length(τs), M*M)
    idx = zeros(Int, M)
    ξ_min = zeros(M)
    for ts = 1:M
        ε★[:,1+(M*(ts-1)):M*ts], _ = pecora(Ys, (0,), (ts,); delays = τs,
                    w = w, samplesize = samplesize, K = K, metric = metric,
                    α = α, p = p, undersampling = false)
        L_min[ts], L_min_idx[ts], idx[ts]  = choose_right_embedding_params(
                                        ε★[:,1+(M*(ts-1)):M*ts], Ys[:,ts],
                                        Ys, τs, KNN, w, samplesize,
                                        metric, econ)
    end
    L_mini, min_idx = findmin(L_min)
    # update τ_vals, ts_vals, Ls, ε★s
    push!(τ_vals, τs[L_min_idx[min_idx]])
    push!(ts_vals, min_idx)             # time series to start with
    push!(ts_vals, idx[min_idx])        # result of 1st embedding cycle
    push!(Ls, L_mini)                   # L-value of 1st embedding cycle
    ε★s[counter] = ε★[:,1+(M*(ts_vals[1]-1)):M*ts_vals[1]]

    # create phase space vector for this embedding cycle
    Y_act = DelayEmbeddings.hcat_lagged_values(Ys[:,ts_vals[counter]],
                                 Ys[:,ts_vals[counter+1]],τ_vals[counter+1])

    return Y_act
end

"""
Perform an embedding cycle of the multivariate embedding, but the first one.
Return the actual reconstruction vector `Y_act`.
"""
function embedding_cycle_pecuzal!(Y_act::Dataset{D, T}, Ys::Dataset{DT, T},
            counter::Int, M::Int, τs, w::Int, samplesize::Real, K::Int, metric,
            α::Real, p::Real, KNN::Int, τ_vals::Vector{Int}, ts_vals::Vector{Int},
            Ls::Vector{T}, ε★s::AbstractMatrix, econ::Bool) where {D, DT, T}

    ε★, _ = pecora(Ys, Tuple(τ_vals), Tuple(ts_vals); delays = τs, w = w,
            samplesize = samplesize, K = K, metric = metric, α = α,
            p = p, undersampling = false)
    # update τ_vals, ts_vals, Ls, ε★s
    choose_right_embedding_params!(ε★, Y_act, Ys, τ_vals, ts_vals, Ls, ε★s,
                                counter, τs, KNN, w, samplesize, metric, econ)
    # create phase space vector for this embedding cycle
    Y_act = DelayEmbeddings.hcat_lagged_values(Y_act, Ys[:, ts_vals[counter+1]],
                                                        τ_vals[counter+1])
    return Y_act
end

"""
    Choose the maximum L-decrease and corresponding τ for each ε★-statistic,
based on picking the peak in ε★, which corresponds to the minimal `L`-statistic.
"""
function choose_right_embedding_params!(ε★::AbstractMatrix, Y_act,
            Ys::Dataset{D, T}, τ_vals::Vector{Int}, ts_vals::Vector{Int},
            Ls::Vector{T}, ε★s::AbstractMatrix, counter::Int, τs, KNN::Int,
            w::Int, samplesize::Real, metric, econ::Bool) where {D, T}

    L_min_ = zeros(size(Ys,2))
    τ_idx = zeros(Int,size(Ys,2))
    for ts = 1:size(Ys,2)
        # zero-padding of ⟨ε★⟩ in order to also cover τ=0 (important for the multivariate case)
        # get the L-statistic for each peak in ⟨ε★⟩ and take the one according to L_min
        L_trials_, max_idx_ = local_L_statistics(vec([0; ε★[:,ts]]), Y_act, Ys[:,ts],
                                        τs, KNN, w, samplesize, metric, econ)
        L_min_[ts], min_idx_ = findmin(L_trials_)
        τ_idx[ts] = max_idx_[min_idx_]-1
    end
    idx = sortperm(L_min_)
    L_mini, min_idx = findmin(L_min_)
    push!(τ_vals, τs[τ_idx[min_idx]])
    push!(ts_vals, min_idx)
    push!(Ls, L_mini)

    ε★s[counter] = ε★
end

"""
    Choose the right embedding parameters of the ε★-statistic in the first
embedding cycle. Return the `L`-decrease-value, the corresponding index value of
the chosen peak `τ_idx` and the number of the chosen time series to start with
`idx`.
"""
function choose_right_embedding_params(ε★::AbstractMatrix, Y_act,
            Ys::Dataset{D, T}, τs, KNN::Int, w::Int, samplesize::Real, metric,
            econ::Bool) where {D, T}
    L_min_ = zeros(size(Ys,2))
    τ_idx = zeros(Int,size(Ys,2))
    for ts = 1:size(Ys,2)
        # zero-padding of ⟨ε★⟩ in order to also cover τ=0 (important for the multivariate case)
        # get the L-statistic for each peak in ⟨ε★⟩ and take the one according to L_min
        L_trials_, max_idx_ = local_L_statistics(vec([0; ε★[:,ts]]), Dataset(Y_act), Ys[:,ts],
                                        τs, KNN, w, samplesize, metric, econ)
        L_min_[ts], min_idx_ = findmin(L_trials_)
        τ_idx[ts] = max_idx_[min_idx_]-1
    end
    idx = sortperm(L_min_)
    return L_min_[idx[1]], τ_idx[idx[1]], idx[1]
end

"""
    Return the maximum decrease of the L-statistic `L_decrease` and corresponding
delay-indices `max_idx` for all local maxima in ε★
"""
function local_L_statistics(ε★::Vector{T}, Y_act::Dataset{D, T}, s::Vector{T},
        τs, KNN::Int, w::Int, samplesize::Real, metric, econ::Bool
        ) where {D, T}
    _, max_idx = get_maxima(ε★) # determine local maxima in ⟨ε★⟩
    L_decrease = zeros(Float64, length(max_idx))
    for (i,τ_idx) in enumerate(max_idx)
        # create candidate phase space vector for this peak/τ-value
        Y_trial = DelayEmbeddings.hcat_lagged_values(Y_act, s, τs[τ_idx-1])
        # compute L-statistic for Y_act and Y_trial and get the maximum decrease
        L_decrease[i] = uzal_cost_pecuzal(Y_act, Y_trial, τs[end]; K = KNN,
                    w , metric, econ, samplesize)
    end
    return L_decrease, max_idx
end

function pecuzal_break_criterion(Ls::Vector{T}, counter::Int,
        max_num_of_cycles::Int, threshold::Real) where {T}
    flag = true
    if counter == 1 && Ls[end] > threshold
        println("Algorithm stopped due to increasing L-values in the first embedding cycle. "*
                "Valid embedding NOT achieved ⨉.")
        flag = false
    end
    if counter > 1 && Ls[end] > threshold
        println("Algorithm stopped due to increasing L-values. "*
                "VALID embedding achieved ✓.")
        flag = false
    end
    if max_num_of_cycles == counter
        println("Algorithm stopped due to hitting max cycle number. "*
                "Valid embedding NOT achieved ⨉.")
        flag = false
    end
    return flag
end

"""
Return the maxima of the given time series s and its indices
"""
function get_maxima(s::Vector{T}) where {T}
    maximas = T[]
    maximas_idx = Int[]
    N = length(s)
    flag = false
    first_point = 0
    for i = 2:N-1
        if s[i-1] < s[i] && s[i+1] < s[i]
            flag = false
            push!(maximas, s[i])
            push!(maximas_idx, i)
        end
        # handling constant values
        if flag
            if s[i+1] < s[first_point]
                flag = false
                push!(maximas, s[first_point])
                push!(maximas_idx, first_point)
            elseif s[i+1] > s[first_point]
                flag = false
            end
        end
        if s[i-1] < s[i] && s[i+1] == s[i]
            flag = true
            first_point = i
        end
    end
    # make sure there is no empty vector returned
    if isempty(maximas)
        maximas, maximas_idx = findmax(s)
    end
    return maximas, maximas_idx
end

"""
    uzal_cost_pecuzal(Y1::Dataset, Y2::Dataset, Tw; kwargs...) → L_decrease
This function is based on the functionality of [`uzal_cost`](@ref), here
specifically tailored for the needs in the PECUZAL algorithm.
Compute the L-statistics `L1` and `L2` for the input datasets `Y1` and `Y2` for
increasing time horizons `T = 1:Tw`. For each `T`, compute `L1` and `L2` and
decrease `L_decrease = L2 - L1`. If `L_decrease` is a negative value, then `Y2`
can be regarded as a "better" reconstruction that `Y1`. Break, when `L_decrease`
reaches the 1st local minima, since this will typically also be the global
minimum. Return the according minimum `L_decrease`-value.

## Keyword arguments

* `K = 3`: the amount of nearest neighbors considered, in order to compute σ_k^2
  (read algorithm description).
  If given a vector, minimum result over all `k ∈ K` is returned.
* `metric = Euclidean()`: metric used for finding nearest neigbhors in the input
  state space trajectory `Y.
* `w = 1`: Theiler window (neighbors in time with index `w` close to the point,
  that are excluded from being true neighbors). `w=0` means to exclude only the
  point itself, and no temporal neighbors.
* `econ::Bool = false`: Economy-mode for L-statistic computation. Instead of
  computing L-statistics for time horizons `2:Tw`, here we only compute them for
  `2:2:Tw`.
* `samplesize::Real = 1`: determine the fraction of all phase space points
  to be considered (fiducial points v) for computing the L-statistic
"""
function uzal_cost_pecuzal(Y::Dataset{D, ET}, Y_trial::Dataset{DT, ET}, Tw::Int;
        K::Int = 3, w::Int = 1, econ::Bool = false, metric = Euclidean(),
        samplesize::Real = 1) where {D, DT, ET}

    @assert DT == D+1
    @assert Tw ≥ 0
    @assert 0 < samplesize ≤ 1.0 "`samplesize` must be ∈ (0,1]"

    if econ
        tws = 2:2:Tw # start at 2 will eliminate bad results for noise
    else
        tws = 2:Tw # start at 2 will eliminate bad results for noise
    end

    NNN = length(Y_trial)-1
    # preallocation for 1st dataset
    ϵ² = zeros(NNN)             # neighborhood size
    E² = zeros(NNN, length(tws))         # conditional variance
    ϵ_ball = zeros(ET, K+1, D)  # epsilon neighbourhood
    u_k = zeros(ET, D)          # center of mass
    # preallocation for 2nd dataset
    ϵ²_trial = zeros(NNN)             # neighborhood size
    E²_trial = zeros(NNN, length(tws))         # conditional variance
    ϵ_ball_trial = zeros(ET, K+1, DT) # epsilon neighbourhood
    u_k_trial = zeros(ET, DT)         # center of mass

    dist_former = 9999999 # intial L-decrease

    NN = length(Y_trial)-tws[end]
    if NN < 1
        error("Time series too short for given possible delays and Theiler window to find enough nearest neighbours")
    end
    if samplesize==1
        ns = 1:NN
        Nfp = length(ns)
    else
        Nfp = Int(floor(samplesize*NN)) # number of considered fiducial points
        ns = sample(vec(1:NN), Nfp, replace = false)  # indices of fiducial points
    end

    vs = Y[ns] # the fiducial points in the data set
    vs_trial = Y_trial[ns] # the fiducial points in the data set

    vtree = KDTree(Y[1:NN], metric)
    allNNidxs, allNNdist = DelayEmbeddings.all_neighbors(vtree, vs, ns, K, w)
    vtree_trial = KDTree(Y_trial[1:NN], metric)
    allNNidxs_trial, allNNdist_trial = DelayEmbeddings.all_neighbors(vtree_trial, vs_trial, ns, K, w)

    # loop over each time horizon
    cnt = 1
    for T in tws
        # compute conditional variances and neighborhood-sizes
        compute_conditional_variances!(ns, vs, vs_trial, allNNidxs,
            allNNidxs_trial, Y, Y_trial, ϵ_ball, ϵ_ball_trial, u_k, u_k_trial,
            T, K, metric, ϵ², ϵ²_trial, E², E²_trial, cnt)

        # compute distance of L-values and check whether that distance can be
        # increased
        dist = compute_L_decrease(E², E²_trial, ϵ², ϵ²_trial, cnt, Nfp)
        if isnan(dist)
            error("Computed 0-distances, due to duplicate datapoints in your data. Try to add minimal additive noise to the signal you wish to embed and try again.")
        end
        if dist > dist_former && dist_former<0
            break
        else
            dist_former = dist
        end
        cnt += 1
    end
    return dist_former
end

function compute_conditional_variances!(ns, vs, vs_trial, allNNidxs::Vector{Array{Int64,1}},
        allNNidxs_trial::Vector{Array{Int64,1}}, Y::Dataset{D, P},
        Y_trial::Dataset{DT, P}, ϵ_ball::Array{P, 2}, ϵ_ball_trial::Array{P, 2},
        u_k::Vector{P}, u_k_trial::Vector{P}, T::Int, K::Int, metric, ϵ²::Vector,
        ϵ²_trial::Vector, E²::Array{P, 2}, E²_trial::Array{P, 2}, cnt::Int) where {P, D, DT}

    # loop over each point on the trajectories
    for (i,v) in enumerate(vs)
        NNidxs = allNNidxs[i] # indices of k nearest neighbors to v of 1st dataset
        NNidxs_trial = allNNidxs_trial[i] # indices of k nearest neighbors to v of 2nd dataset
        # pairwise distance of fiducial points and all neighbours
        pdsqrd = DelayEmbeddings.fiducial_pairwise_dist_sqrd(view(Y.data, NNidxs), v, metric)
        ϵ²[i] = (2/(K*(K+1))) * pdsqrd  # Eq. 16
        pdsqrd_trial = DelayEmbeddings.fiducial_pairwise_dist_sqrd(view(Y_trial.data, NNidxs_trial), vs_trial[i], metric)
        ϵ²_trial[i] = (2/(K*(K+1))) * pdsqrd_trial  # Eq. 16

        E²[i, cnt] = comp_Ek2!(ϵ_ball, u_k, Y, ns[i], NNidxs, T, K, metric) # Eqs. 13 & 14
        E²_trial[i, cnt] = comp_Ek2!(ϵ_ball_trial, u_k_trial, Y_trial, ns[i], NNidxs_trial, T, K, metric) # Eqs. 13 & 14
    end
end

function compute_L_decrease(E²::Array{P, 2}, E²_trial::Array{P, 2}, ϵ²::Vector{P},
        ϵ²_trial::Vector{P}, T::Int, NN::Int) where {P}
    # 1st dataset
    # Average E²[T] over all prediction horizons
    E²_avrg = mean(E²[1:NN,1:T], dims=2)                   # Eq. 15
    σ² = E²_avrg ./ ϵ²[1:NN] # noise amplification σ², Eq. 17
    σ²_avrg = mean(σ²) # averaged value of the noise amplification, Eq. 18
    α² = 1 / mean(ϵ²[1:NN].^(-1)) # for normalization, Eq. 21
    L = log10(sqrt(σ²_avrg)*sqrt(α²))
    # 2nd dataset
    # Average E²[T] over all prediction horizons
    E²_avrg_trial = mean(E²_trial[1:NN,1:T], dims=2)                   # Eq. 15
    σ²_trial = E²_avrg_trial ./ ϵ²_trial[1:NN] # noise amplification σ², Eq. 17
    σ²_avrg_trial = mean(σ²_trial) # averaged value of the noise amplification, Eq. 18
    α²_trial = 1 / mean(ϵ²_trial[1:NN].^(-1)) # for normalization, Eq. 21
    L_trial = log10(sqrt(σ²_avrg_trial)*sqrt(α²_trial))

    return L_trial - L
end

"""
    comp_Ek2!(ϵ_ball,u_k,Y,v,NNidxs,T,K,metric) → E²(T)
Returns the approximated conditional variance for a specific point in state space
`ns` (index value) with its `K`-nearest neighbors, which indices are stored in
`NNidxs`, for a time horizon `T`. This corresponds to Eqs. 13 & 14 in [^Uzal2011].
The norm specified in input parameter `metric` is used for distance computations.
"""
function comp_Ek2!(ϵ_ball::Array{P, 2}, u_k::Vector{P}, Y::Dataset{D, P},
        ns::Int, NNidxs::Vector, T::Int, K::Int, metric
        ) where {D, P}
    # determine neighborhood `T` time steps ahead
    ϵ_ball[1, :] .= Y[ns+T]
    @inbounds for (i, j) in enumerate(NNidxs)
        ϵ_ball[i+1, :] .= Y[j + T]
    end

    # compute center of mass
    @inbounds for i in 1:size(Y)[2]; u_k[i] = sum(view(ϵ_ball, :, i))/(K+1); end # Eq. 14

    E²_sum = 0
    @inbounds for j = 1:K+1
        E²_sum += (evaluate(metric, view(ϵ_ball, j, :), u_k))^2
    end
    E² = E²_sum / (K+1)         # Eq. 13
end


text:
"""
The base type used by all nodes within a Koza tree.
"""
abstract KozaNode <: KozaParent

"""
Returns a randomly selected node from a given sub-tree.
"""
sample(n::KozaNode) = search(n, rand(1:num_nodes(n)))

"""
Calculates the depth of a node within the tree which it belongs to.
"""
function depth(n::KozaNode)
  i = 1
  while !isroot(n)
    n = n.parent
    i += 1
  end
  return i
end

"""
Performs a breadth-first search of a sub-tree rooted at a given node, n, for
the node at location p.
"""
function search(n::KozaNode, p::Int)
  i = 1
  q = KozaNode[n]

  # We don't actually have to check this...
  # We could just keep popping nodes?
  while !isempty(q)
    n = pop!(q) # pop vs. shift
    if i == p
      return n
    end
    if isa(n, KozaNonTerminal)
      append!(q, n.children)
    end
    i += 1
  end
end


text:
function operation(s, n, op)
	if s == undef
		s = n
	else
		if op == '+'
			s += n
		elseif op == '*'
			s *= n
		end
	end
	s
end

function eval1(expr, start = 1)
	ops = split(expr, " ")
	s = undef
	lastop = undef
	i = start
	while i <= length(ops)
		try
			n = parse(Int, ops[i])
			s = operation(s, n, lastop)
		catch
			op = ops[i]
			if op[1] in ['+','*']
				lastop = op[1]
			elseif op[1] == '('
				(x,i) = eval1(expr, i+1)
				s = operation(s, x, lastop)
			elseif op[1] == ')'
				return (s, i)
			end
		end
		i += 1
	end
	(s, length(ops))
end

function replaceplusmul(e)
	for i=1:length(e.args)
		if e.args[i] == :+
			e.args[i] = :*
		elseif e.args[i] == :*
			e.args[i] = :+
		elseif typeof(e.args[i]) == Expr
			e.args[i] = replaceplusmul(e.args[i])
		end
	end
	e
end

function eval2(expr)
	expr = replace(expr, "+" => "#")
	expr = replace(expr, "*" => "+")
	expr = replace(expr, "#" => "*")
	e = Meta.parse(expr)
	e = replaceplusmul(e)
	x = eval(e)
end

#prints solution of both parts
function x()
	f = open("18.txt")
	s1 = 0
	s2 = 0
	for i=1:377
		#println("line {$i}")
		e = readline(f)
		s1 += eval1(e)[1]
		s2 += eval2(e)[1]
	end
	println(s1)
	println(s2)
end


text:
module Slomo

include("utils.jl")
include("constants.jl")
include("cosmology.jl")
include("integrate.jl")
include("models.jl")
include("halos/halos.jl")
include("tracers.jl")
include("anisotropy.jl")
include("jeans.jl")

using Slomo.Models
using Slomo.Halos
using Slomo.Tracers
using Slomo.Anisotropy
using Slomo.Jeans

export Halos, SersicModel, ConstantBetaModel, RSBetaModel, JeansModel
export mass, density, density2d, beta, update
export sigma_los, sigma_los_parallel, kappa_los

end