text:
RealOrComplex = Union{Real,Complex}

## User should just use TransferFunction
immutable SisoZpk <: SisoTf
    z::Vector{Complex{Float64}}
    p::Vector{Complex{Float64}}
    k::Float64
    function SisoZpk(z::Vector{Complex{Float64}}, p::Vector{Complex{Float64}}, k::Float64)
        if k == zero(k)
            p = []
            z = []
        end
        new(z, p, k)
    end
end

SisoZpk{T<:RealOrComplex,S<:RealOrComplex}(z::AbstractArray{T}, p::AbstractArray{S}, k::Real) = SisoZpk(Complex128[z...], Complex128[p...], Float64(k))

# Taking care of empty vectors being of type Any
function SisoZpk(z::AbstractArray, p::AbstractArray, k::Real)
    if length(z) > 0
        if !(eltype(z) <: RealOrComplex)
            error("Zeros must be real or complex")
        end
    else
        z = Array(Float64,0)
    end
    if length(p) > 0
        if !(eltype(p) <: RealOrComplex)
            error("poles must be real or complex")
        end
    else
        p = Array(Float64,0)
    end
    SisoZpk(z, p, k)
end

function minreal(sys::SisoZpk, eps::Real)
    newZ = copy(sys.z)
    newP = Vector{Complex{Float64}}()
    doubles = Vector{Int64}()
    newZ = copy(sys.z)
    for  p in sys.p
        if !isempty(newZ)
            val, zi = findmin(abs(p-newZ))
        else
            val = Inf #Keep looping over p, adding poles
        end
        if val < eps
            deleteat!(newZ, zi)
            continue
        else
            push!(newP, p)
        end
    end
    SisoZpk(newZ, newP, sys.k)
end

function Base.num(t::SisoZpk)
    return copy(t.z)
end

function Base.den(t::SisoZpk)
    return copy(t.p)
end

tzero(sys::SisoZpk) = num(sys)
pole(sys::SisoZpk) = den(sys)

function zp2polys(vec)
    polys = Array{Poly{Float64},1}(0)
    polesiLeft = Set(1:length(vec))
    while length(polesiLeft) > 0
        p = vec[pop!(polesiLeft)]
        if abs(imag(p)) < sqrt(eps())
            push!(polys,Poly(float([1, -real(p)])))
        else
            polesiLeftVec = [i for i in polesiLeft]
            polesTest = Complex128[vec[polesiLeftVec]...]
            val, i = findmin(abs(polesTest-conj(p)))
            val > 2*sqrt(eps()) && error("Could not find conjugate to pole")
            push!(polys,Poly(float([1, -2*real(p), real(p)^2+imag(p)^2])))
            pop!(polesiLeft,polesiLeftVec[i])
        end
    end
    polys
end

function numpoly(G::SisoZpk)
    zpolys = zp2polys(G.z)
end

function denpoly(G::SisoZpk)
    ppolys = zp2polys(G.p)
end

function evalfr(sys::SisoZpk, s::Number)
    S = promote_type(typeof(s), Float64)
    den = reduce(*, (poly) -> polyval(poly, s), zp2polys(sys.p))
    if den == zero(S)
        return convert(S, Inf)
    else
        num = reduce(*, (poly) -> polyval(poly, s), zp2polys(sys.z))
        return sys.k*num/den
    end
end

function print_siso(io::IO, t::SisoZpk, var=:s)
    zpolys = zp2polys(t.z)
    ppolys = zp2polys(t.p)
    # Convert the numerator and denominator to strings
    if length(zpolys) < 2
        numstr = ( length(zpolys) == 0 ) ? "1.0" : sprint(print_poly, zpolys[1], var)
    else
        numstr = reduce(*,"",["("*sprint(print_poly, z, var)*")" for z in zpolys])
    end
    if length(ppolys) < 2
        denstr = ( length(ppolys) == 0 ) ? "1.0" : sprint(print_poly, ppolys[1], var)
    else
        denstr = reduce(*,"",["("*sprint(print_poly, p, var)*")" for p in ppolys])
    end
    # Figure out the length of the separating line
    len_num = length(numstr)
    len_den = length(denstr)
    dashcount = max(len_num, len_den)

    # Center the numerator or denominator
    if len_num < dashcount
        numstr = "$(repeat(" ", div(dashcount - len_num, 2)))$numstr"
    else
        denstr = "$(repeat(" ", div(dashcount - len_den, 2)))$denstr"
    end

    gainstr = string(t.k)
    #Add spaces to account for gain string
    numstr = " "^(length(gainstr))*numstr
    denstr = " "^(length(gainstr))*denstr
    println(io, numstr)
    println(io, gainstr*repeat("-", dashcount))
    println(io, denstr)
end

Base.promote_rule{T<:Real}(::Type{SisoZpk}, ::Type{T}) = SisoZpk
Base.convert(::Type{SisoZpk}, b::Real) = SisoZpk([], [], b)

Base.zero(::Type{SisoZpk}) = SisoZpk([],[],0.0)
Base.zero(::SisoZpk) = Base.zero(SisoZpk)

Base.length(t::SisoZpk) = max(length(t.z), length(t.p))

==(t1::SisoZpk, t2::SisoZpk) = (t1-t2).k == 0.0
function isapprox(t1::SisoZpk, t2::SisoZpk; rtol::Real=sqrt(eps()), atol::Real=sqrt(eps()))
    tdiff = t1 - t2
    isapprox(tdiff.k, 0, atol=atol, rtol=rtol)
end

function +(t1::SisoZpk, t2::SisoZpk)
  numPoly = t1.k*prod(zp2polys(t1.z))*prod(zp2polys(t2.p))+t2.k*prod(zp2polys(t2.z))*prod(zp2polys(t1.p))
  z = roots(numPoly)
  if length(numPoly) > 0
      k = numPoly[1]
      p = [t1.p;t2.p]
  else
      k = 0
      p = []
  end
  SisoZpk(z,p,k)
end

+(t::SisoZpk, n::Real) = t + SisoZpk([],[],n)
+(n::Real, t::SisoZpk) = SisoZpk([],[],n) + t
.+(t::SisoZpk, n::Real) = t + n
.+(n::Real, t::SisoZpk) = t + n

-(t1::SisoZpk, t2::SisoZpk) = +(t1,-t2)
-(n::Real, t::SisoZpk) = SisoZpk([],[],n) - t
-(t::SisoZpk, n::Real) = +(t, -n)
.-(t::SisoZpk, n::Real) = -(t, n)
.-(n::Real, t::SisoZpk) = -(n, t)

-(t::SisoZpk) = SisoZpk(t.z, t.p, -t.k)

*(t1::SisoZpk, t2::SisoZpk) = SisoZpk([t1.z;t2.z], [t1.p;t2.p], t1.k*t2.k)
*(t::SisoZpk, n::Real) = SisoZpk(t.z, t.p, t.k*n)
*(n::Real, t::SisoZpk) = *(t, n)
.*(t1::SisoZpk, t2::SisoZpk) = *(t1, t2)
.*(t::SisoZpk, n::Real) = *(t, n)
.*(n::Real, t::SisoZpk) = *(t, n)

/(n::Real, t::SisoZpk) = SisoZpk(t.p, t.z, n/t.k)
/(t::SisoZpk, n::Real) = SisoZpk(t.z, t.p, t.k/n)
/(t1::SisoZpk, t2::SisoZpk) = t1*(1/t2)


text:
#=

  Euler #11 in Julia.

  Problem 11
  """
  In the 20x20 grid below, four numbers along a diagonal line have
  been marked in red.

  ...

  The product of these numbers is 26 x 63 x 78 x 14 = 1788696.

  What is the greatest product of four adjacent numbers in any direction
  (up, down, left, right, or diagonally) in the 20 x 20 grid?
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http:www.hakank.org/julia/

=#

include("Euler.jl")

m = [
     8  2 22 97 38 15  0 40  0 75  4  5  7 78 52 12 50 77 91  8;
    49 49 99 40 17 81 18 57 60 87 17 40 98 43 69 48  4 56 62  0;
    81 49 31 73 55 79 14 29 93 71 40 67 53 88 30  3 49 13 36 65;
    52 70 95 23  4 60 11 42 69 24 68 56  1 32 56 71 37  2 36 91;
    22 31 16 71 51 67 63 89 41 92 36 54 22 40 40 28 66 33 13 80;
    24 47 32 60 99  3 45  2 44 75 33 53 78 36 84 20 35 17 12 50;
    32 98 81 28 64 23 67 10 26 38 40 67 59 54 70 66 18 38 64 70;
    67 26 20 68  2 62 12 20 95 63 94 39 63  8 40 91 66 49 94 21;
    24 55 58  5 66 73 99 26 97 17 78 78 96 83 14 88 34 89 63 72;
    21 36 23  9 75  0 76 44 20 45 35 14  0 61 33 97 34 31 33 95;
    78 17 53 28 22 75 31 67 15 94  3 80  4 62 16 14  9 53 56 92;
    16 39  5 42 96 35 31 47 55 58 88 24  0 17 54 24 36 29 85 57;
    86 56  0 48 35 71 89  7  5 44 44 37 44 60 21 58 51 54 17 58;
    19 80 81 68  5 94 47 69 28 73 92 13 86 52 17 77  4 89 55 40;
     4 52  8 83 97 35 99 16  7 97 57 32 16 26 26 79 33 27 98 66;
    88 36 68 87 57 62 20 72  3 46 33 67 46 55 12 32 63 93 53 69;
     4 42 16 73 38 25 39 11 24 94 72 18  8 46 29 32 40 62 76 36;
    20 69 36 41 72 30 23 88 34 62 99 69 82 67 59 85 74  4 36 16;
    20 73 35 29 78 31 90  1 74 31 49 71 48 86 81 16 23 57  5 54;
     1 70 54 71 83 51 54 69 16 92 33 48 61 43 52  1 89 19 67 48
]

# 0.07859598s
function euler11a()
    len = size(m)[1]
    mt = m'
    n = 4
    maxp = 0

    # Rows and columns
    for row = 1:len-1
        for col in 1:len-n
            p1 = prod(m[row,col:col+n-1])
            if p1 > maxp
                maxp = p1
            end
            p2 = prod(mt[row,col:col+n-1])
            if p2 > maxp
                maxp = p2
            end
        end
    end

    # diag_down
    for j in 1:17
        for i in 1:17
            pa = []
            for a in 0:n-1
                push!(pa,m[a+i,a+j])
            end
            p = prod(pa)
            if p > maxp
                maxp = p
            end
        end
    end

    # diag_up
    for j in 2:17
        for i = 5:20
            pa = []
            for a = 0:3
                push!(pa,m[i-a,j+a])
            end
            p = prod(pa)
            if p > maxp
                maxp = p
            end
        end
    end

    return maxp;
end

run_euler(euler11a)

ext: jl

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# # Convex Optimization in Julia
#
# ## Madeleine Udell | ISMP 2015

#-

# ## Convex.jl team
#
# * [Convex.jl](https://github.com/cvxgrp/Convex.jl): Madeleine Udell, Karanveer Mohan, David Zeng, Jenny Hong

#-

# ## Collaborators/Inspiration:
#
# * [CVX](http://www.cvxr.com): Michael Grant, Stephen Boyd
# * [CVXPY](https://github.com/cvxgrp/cvxpy): Steven Diamond, Eric Chu, Stephen Boyd
# * [JuliaOpt](https://github.com/JuliaOpt): Miles Lubin, Iain Dunning, Joey Huchette

## initial package installation

#-

## Make the Convex.jl module available
using Convex, SparseArrays, LinearAlgebra
using SCS # first order splitting conic solver [O'Donoghue et al., 2014]

## Generate random problem data
m = 50;  n = 100
A = randn(m, n)
x♮ = sprand(n, 1, .5) # true (sparse nonnegative) parameter vector
noise = .1*randn(m)    # gaussian noise
b = A*x♮ + noise      # noisy linear observations

## Create a (column vector) variable of size n.
x = Variable(n)

## nonnegative elastic net with regularization
λ = 1
μ = 1
problem = minimize(square(norm(A * x - b)) + λ*square(norm(x)) + μ*norm(x, 1),
                   x >= 0)

## Solve the problem by calling solve!
solve!(problem, SCS.Optimizer(verbose=0))

println("problem status is ", problem.status) # :Optimal, :Infeasible, :Unbounded etc.
println("optimal value is ", problem.optval)

#-

using Interact, Plots
## Interact.WebIO.install_jupyter_nbextension() # might be helpful if you see `WebIO` warnings in Jupyter
@manipulate throttle=.1 for λ=0:.1:5, μ=0:.1:5
    global A
    problem = minimize(square(norm(A * x - b)) + λ*square(norm(x)) + μ*norm(x, 1),
                   x >= 0)
    solve!(problem, SCS.Optimizer(verbose=0))
    histogram(evaluate(x), xlims=(0,3.5), label="x")
end

# # Quick convex prototyping

#-

# ## Variables

## Scalar variable
x = Variable()

#-

## (Column) vector variable
y = Variable(4)

#-

## Matrix variable
Z = Variable(4, 4)

# # Expressions

#-

# Convex.jl allows you to use a [wide variety of functions](http://convexjl.readthedocs.org/en/latest/operations.html) on variables and on expressions to form new expressions.

x + 2x

#-

e = y[1] + logdet(Z) + sqrt(x) + minimum(y)

# ### Examine the expression tree

e.children[2]

# # Constraints
#
# A constraint is convex if convex combinations of feasible points are also feasible. Equivalently, feasible sets are convex sets.
#
# In other words, convex constraints are of the form
#
# * `convexExpr <= 0`
# * `concaveExpr >= 0`
# * `affineExpr == 0`

x <= 0

#-

square(x) <= sum(y)

#-

M = Z
for i = 1:length(y)
    global M += rand(size(Z)...)*y[i]
end
M ⪰ 0

# # Problems

x = Variable()
y = Variable(4)
objective = 2*x + 1 - sqrt(sum(y))
constraint = x >= maximum(y)
p = minimize(objective, constraint)

#-

## solve the problem
solve!(p, SCS.Optimizer(verbose=0))
p.status

#-

x.value

#-

## can evaluate expressions directly
evaluate(objective)

# ## Pass to solver
#
# call a `MathProgBase` solver suited for your problem class
#
# * see the [list of Convex.jl operations](http://convexjl.readthedocs.org/en/latest/operations.html) to find which cones you're using
# * see the [list of solvers](http://www.juliaopt.org/) for an up-to-date list of solvers and which cones they support

#-

# to solve problem using a different solver, just import the solver package and pass the solver to the `solve!` method: eg
#
#     using Mosek
#     solve!(p, Mosek.Optimizer())

#-

# ## Warmstart

## Generate random problem data
m = 50;  n = 100
A = randn(m, n)
x♮ = sprand(n, 1, .5) # true (sparse nonnegative) parameter vector
noise = .1*randn(m)    # gaussian noise
b = A*x♮ + noise      # noisy linear observations

## Create a (column vector) variable of size n.
x = Variable(n)

## nonnegative elastic net with regularization
λ = 1
μ = 1
problem = minimize(square(norm(A * x - b)) + λ*square(norm(x)) + μ*norm(x, 1),
                   x >= 0)
@time solve!(problem, SCS.Optimizer(verbose=0))
λ = 1.5
@time solve!(problem, SCS.Optimizer(verbose=0), warmstart = true)

# # DCP examples

## affine
x = Variable(4)
y = Variable(2)
sum(x) + y[2]

#-

2*maximum(x) + 4*sum(y) - sqrt(y[1] + x[1]) - 7 * minimum(x[2:4])

#-

## not dcp compliant
log(x) + square(x)

#-

## $f$ is convex increasing and $g$ is convex
square(pos(x))

#-

## $f$ is convex decreasing and $g$ is concave
invpos(sqrt(x))

#-

## $f$ is concave increasing and $g$ is concave
sqrt(sqrt(x))


text:
using MyPackage
using Documenter

makedocs(;
    modules=[MyPackage],
    authors="nobuta05 <nobuta05@gmail.com> and contributors",
    repo="https://github.com/nobuta05/MyPackage.jl/blob/{commit}{path}#L{line}",
    sitename="MyPackage.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://nobuta05.github.io/MyPackage.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/nobuta05/MyPackage.jl",
)


text:
###########################################################
# Generate discretization matrices for FINITE DIFFERENCE
"""
    generate_fd_interior_mat

Generate the discretization matrices for the interior for Finite Difference.
"""
function generate_fd_interior_mat( nx, nz, hx, hz)
    ex = ones(nx)
    ez = ones(nz)

    # DISCRETIZATION OF THE SECOND DERIVATIVE
    Dxx = spdiagm(-1 => ex[1:end-1], 0 => -2*ex, 1 => ex[1:end-1])
    Dzz = spdiagm(-1 => ez[1:end-1], 0 => -2*ez, 1 => ez[1:end-1])
    #IMPOSE PERIODICITY IN Z-DIRECTION
    Dzz[1, end] = 1;
    Dzz[end, 1] = 1;

    Dxx = Dxx/(hx^2);
    Dzz = Dzz/(hz^2);

    # DISCRETIZATION OF THE FIRST DERIVATIVE
    Dz  = spdiagm(-1 => -ez[1:end-1], 1 => ez[1:end-1])

    #IMPOSE PERIODICITY
    Dz[1, end] = -1;
    Dz[end, 1] = 1;

    Dz = Dz/(2*hz);

    return Dxx, Dzz, Dz
end

"""
    generate_fd_boundary_mat

Generate the discretization matrices for the boundary for Finite Difference.
"""
function generate_fd_boundary_mat( nx, nz, hx, hz)
    # BUILD THE SECOND BLOCK C1
    e1 = spzeros(nx,1)
    e1[1] = 1
    en = spzeros(nx,1)
    en[end] = 1
    Iz = sparse(1.0I, nz, nz)
    C1 = [kron(e1,Iz) kron(en,Iz)]/(hx^2);

    # BUILD THE THIRD BLOCK C2^T
    d1 = 2/hx;
    d2 = -1/(2*hx);
    vm = spzeros(1,nx);
    vm[1] = d1;
    vm[2] = d2;
    vp = spzeros(1,nx);
    vp[end] = d1;
    vp[end-1] = d2;
    C2T = [kron(vm,Iz); kron(vp,Iz)];

    return C1, C2T
end

###########################################################
# Generate Wavenumber FINITE DIFFERENCE
"""
    generate_wavenumber_fd

Generate a wavenumber for Finite Difference.
"""
function generate_wavenumber_fd( nx::Integer, nz::Integer, wg::String, delta::Number)
    if wg == "TAUSCH"
        return generate_wavenumber_fd_tausch( nx, nz, delta)
    elseif wg == "JARLEBRING"
        return generate_wavenumber_fd_jarlebring( nx, nz, delta)
    end
    # Use early-bailout principle. If a supported waveguide is found, compute and return. Otherwise end up here and throw and error
    error("No wavenumber loaded: The given Waveguide '", wg ,"' is not supported in 'FD' discretization.")
end

"""
    generate_wavenumber_fd_tausch

Generate the wavenumber in FD discretization for the waveguide
described by TAUSCH.
"""
function generate_wavenumber_fd_tausch( nx::Integer, nz::Integer, delta::Number)
    xm = 0;
    xp = (2/pi) + 0.4;
    zm = 0;
    zp = 1;

    xm = xm - delta;
    xp = xp + delta;

    # Domain (First generate including the boundary)
    X = range(xm, stop = xp, length = nx+2)
    hx = step(X);
    X = collect(X);
    Z = range(zm, stop = zp, length = nz+1)
    hz = step(Z);
    Z = collect(Z);
    # Removing the boundary
    X = X[2:end-1];
    Z = Z[2:end];

    # The actual wavenumber
    k1 = sqrt(2.3)*pi;
    k2 = sqrt(3)*pi;
    k3 = pi;
    k = function(x,z)
            z_ones = ones(size(z)) #Placeholder of ones to expand x-vector
            k1*(x .<= 0) .* z_ones +
            k2*(x.>0) .* (x.<=2/pi) .* z_ones +
            k2*(x.>2/pi) .* (x.<=(2/pi+0.4)) .* (z.>0.5) +
            k3*(x.>2/pi) .* (z.<=0.5) .* (x.<=(2/pi+0.4)) +
            k3*(x.>(2/pi+0.4)) .* z_ones;
        end

    K = k(X', Z).^2;
    Km = k(-Inf, 1/2)[1];
    Kp = k(Inf, 1/2)[1];
    return K, hx, hz, Km, Kp
end

"""
    generate_wavenumber_fd_jarlebring

Generate the wavenumber in FD discretization for the waveguide
described by JARLEBRING.
"""
function generate_wavenumber_fd_jarlebring( nx::Integer, nz::Integer, delta::Number)
    xm = -1;
    xp = 1;
    zm = 0;
    zp = 1;

    xm = xm - delta;
    xp = xp + delta;

    # Domain (First generate including the boundary)
    X = range(xm, stop = xp, length = nx+2)
    hx = step(X);
    X = collect(X);
    Z = range(zm, stop = zp, length = nz+1)
    hz = step(Z);
    Z = collect(Z);
    # Removing the boundary
    X = X[2:end-1];
    Z = Z[2:end];

    # The actual wavenumber
    k1 = sqrt(2.3)*pi;
    k2 = 2*sqrt(3)*pi;
    k3 = 4*sqrt(3)*pi;
    k4 = pi;
    k = function(x,z)
            z_ones = ones(size(z)) #Placeholder of ones to expand x-vector
            x_ones = ones(size(x)) #Placeholder of ones to expand z-vector
            k1 *(x.<=-1) .* z_ones  +
            k4 *(x.>1) .* z_ones  +
            k4 *(x.>(-1+1.5)) .* (x.<=1) .* (z.<=0.4) +
            k3 *(x.>(-1+1)) .* (x.<=(-1+1.5)) .* z_ones +
            k3 *(x.>(-1+1.5)) .* (x.<=1) .* (z.>0.4) +
            k3 *(x.>-1) .* (x.<=(-1+1)) .* (z.>0.5) .* (z.*x_ones-(x.*z_ones)/2 .<= 1) +
            k2 *(x.>-1) .* (x.<=(-1+1)) .* (z.>0.5) .* (z.*x_ones-(x.*z_ones)/2 .> 1) +
            k3 *(x.>-1) .* (x.<=(-1+1)) .* (z.<=0.5) .* (z.*x_ones+(x.*z_ones)/2 .> 0) +
            k2 *(x.>-1) .* (x.<=(-1+1)) .* (z.<=0.5) .* (z.*x_ones+(x.*z_ones)/2 .<= 0);
        end

    K = k(X', Z).^2;
    Km = k(-Inf, 1/2)[1];
    Kp = k(Inf, 1/2)[1];
    return K, hx, hz, Km, Kp
end


text:
include("../roblib.jl")

function mark_all_cells(r)
    num_steps = through_rectangles_into_angle(r, (Down, Left))

    side = Right
    while (isborder(r, Up) == false) || (isborder(r, Right) == false)
        putmarkers!(r, side)
        if (isborder(r, Up) == false)
            move!(r, Up)
        end
        side = inverse(side)
    end

    through_rectangles_into_angle(r, (Down, Left))
    movements!(r, (Up, Right), reverse(num_steps))
end


text:
"""
Given an integer n = 4, the output should be

   #
  ##
 ###
####

ref: https://www.hackerrank.com/challenges/staircase/problem
"""

function staircase(size = nothing,builder = nothing, pad_char = nothing)
    """
    The method takes the size, builder and padding character of a staircase
    and prints the result to console
    if size is not provided, size is defaulted to 6
    if builder is not provided, defaulted to #
    if pad_char is not provided, defaulted to ' ', white space
    """

    # check if values are not provided and assign appropriate defaults
    size = (size === nothing) ? 6 : size
    builder = (builder === nothing) ? "#" : builder
    pad_char = (pad_char === nothing) ? " " : pad_char

    # Julia is 1-index based, this means from 1 to size
    # for example, for size = 6, we get 1,2,3,4,5,6
    for loop in 1:size
        # print line to console, pad white spaces and then builder
        println(join([pad_char for lp in 1:(size-loop)]),join([builder for rl in 1:loop]))
    end
end

# printing an empty line
println() # printing an empty line at the end, to look nice
staircase(9)
println() # printing an empty line at the end, to look nice

max_stars_repo_name: csendranshi/Hackerrank-Codes

max_issues_repo_name: csendranshi/Hackerrank-Codes

max_forks_repo_name: csendranshi/Hackerrank-Codes


text:
module TestLTSA
    using ManifoldLearning
    using Test

    @testset "LTSA" begin

	k = 12
	d = 2
	X, L = swiss_roll()
	Y = transform(LTSA, X; k=k, d=d)

	@test outdim(Y) == d
	@test size(projection(Y), 2) == size(X, 2)
	@test neighbors(Y) == k
	@test length(eigvals(Y)) == d
    end
end


text:
module ImportTweetsTask

using Genie, Twitter, SearchLight, App, Toolbox, Logger, SentimentAnalysis, Channels

function description()
  """
  Imports a new batch of tips from the Twitter API
  """
end

function run_task!() :: TaskResult
  while true
    for search in SearchLight.find(Search)
      if haskey(Channels.SUBSCRIPTIONS, "tweets-" * search.search)
        import_tweets("#" * search.search)
      end
    end

    sleep(.5)
  end

  VoidTaskResult()
end

function import_tweets(search_hash)
  latest_tweet = SearchLight.find(Tweet, SQLQuery(order = "id DESC", limit = 1))
  latest_tweet_id = isempty(latest_tweet) ? "1" : first(latest_tweet).tweet_id

  for tweet in tweets(search_hash, latest_tweet_id)
    t = App.Tweet(tweet_id = tweet.id_str,
                  text = tweet.full_text,
                  user_name = tweet.user["name"],
                  screen_name = tweet.user["screen_name"],
                  avatar_url = tweet.user["profile_image_url"],
                  favorite_count = tweet.favorite_count,
                  retweet_count = tweet.retweet_count)

    sentiment = SentimentAnalysis.analyze(t.text)
    t.polarity = sentiment.polarity
    t.subjective = sentiment.subjective
    t.search_hash = search_hash[2:end]

    SearchLight.create_or_update_by!!(t, :tweet_id)

    sleep(.5)
  end
end

function tweets(hash, latest_tweet_id)
  get_search_tweets(hash, options = Dict("tweet_mode" => "extended", "result_type" => "recent", "lang" => "en", "since_id" => latest_tweet_id))
end

end

max_stars_repo_name: essenciary/genie-demo-tweet-stats

max_issues_repo_name: essenciary/genie-demo-tweet-stats

max_forks_repo_name: essenciary/genie-demo-tweet-stats


text:
#
# Correctness Tests
#

using Base.Test
using DataFrames

my_tests = ["extras.jl",
            "data.jl",
            "index.jl",
            "dataframe.jl",
            "operators.jl",
            "io.jl",
            # "formula.jl",
            "datastream.jl",
            "constructors.jl",
            "indexing.jl",
            # "indexedvector.jl",
            "RDA.jl",
            "sort.jl",
            "grouping.jl",
            "iteration.jl",
            #"test/colfuncs.jl",
            "duplicates.jl"]

println("Running tests:")

for my_test in my_tests
    println(" * $(my_test)")
    include(my_test)
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.4
#   kernelspec:
#     display_name: Julia 1.3.1
#     language: julia
#     name: julia-1.3
# ---

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# # Transport Streamfunction And Divergent Component
#
# Here we apply such [Helmholtz Decomposition](https://en.wikipedia.org/wiki/Helmholtz_decomposition) to vertically integrated transports defined over a Global Ocean model (`C-grid`).
#
# 1. read vertically integrated transport from file
# 2. compute its convergence and apply land mask
# 3. decompose into rotational and divergent components
# 4. derive streamfunction from the rotational component

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Read Model Grid & Transports From File
#
# 1. pre-requisites
# 2. read variables

# +
#]add MITgcmTools#master

# + {"slideshow": {"slide_type": "-"}, "cell_style": "center"}
using MeshArrays, Plots, Statistics, MITgcmTools
include("helper_functions.jl")
get_grid_if_needed()

Γ=GridLoad(GridSpec("LatLonCap","../inputs/GRID_LLC90/"))
(Tx,Ty,τx,τy,η)=trsp_read("LatLonCap","../inputs/GRID_LLC90/")
SPM,lon,lat=read_SPM("../inputs/GRID_LLC90/")
msk=1.0 .+ 0.0 * mask(view(Γ["hFacC"],:,1),NaN,0.0);

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Helmholtz Decomposition
#
# 1. convergence & land mask
# 2. compute scalar potential
# 3. subtract divergent component
# 4. compute vector potential / streamfunction
#
# _note: masking avoids singularities related to isolated canyons_

# + {"slideshow": {"slide_type": "subslide"}}
#convergence & land mask
TrspCon=msk.*convergence(Tx,Ty);

#scalar potential
TrspPot=ScalarPotential(TrspCon)

#Divergent transport component
(TxD,TyD)=gradient(TrspPot,Γ)
TxD=TxD.*Γ["DXC"]
TyD=TyD.*Γ["DYC"]

#Rotational transport component
TxR = Tx-TyD
TyR = Ty-TyD

#vector Potential
TrspPsi=VectorPotential(TxR,TyR,Γ);

# + {"slideshow": {"slide_type": "subslide"}, "cell_type": "markdown"}
# ### Verify The Results Consistency
#
# _TrspCon-tmpCon should be negligible compared with TrspCon_

# + {"slideshow": {"slide_type": "-"}}
tmpCon=convergence(TxD,TyD)
tmp1=TrspCon[3]
tmp2=tmp1[findall(isfinite.(tmp1))]
errCon=1/sqrt(mean(tmp2.^2)).*(tmpCon[3]-TrspCon[3])
heatmap(errCon)

# + {"slideshow": {"slide_type": "slide"}, "cell_type": "markdown"}
# ### Map Out Streamfunction And Scalar Potential
#
# _Interpolation is used to create global maps_

# + {"slideshow": {"slide_type": "-"}, "cell_style": "split"}
TrspPsiI=MatrixInterp(write(1e-6*msk*TrspPsi),SPM,size(lon))
contourf(vec(lon[:,1]),vec(lat[1,:]),transpose(TrspPsiI),
    title="Streamfunction",clims=(-50,50))

# + {"slideshow": {"slide_type": "-"}, "cell_style": "split"}
TrspPotI=MatrixInterp(write(1e-6*msk*TrspPot),SPM,size(lon))
contourf(vec(lon[:,1]),vec(lat[1,:]),transpose(TrspPotI),
    title="Scalar Potential",clims=(-0.4,0.4))

max_stars_repo_name: abelolivelanna/GlobalOceanNotebooks

max_issues_repo_name: abelolivelanna/GlobalOceanNotebooks

max_forks_repo_name: abelolivelanna/GlobalOceanNotebooks


text:
using Random, AutoGrad, Knet, KnetLayers
"""
    intro(args)

    Create configuration dictionary given program arguments.
"""
function intro(args)
    s = ArgParseSettings()
    @add_arg_table s begin
        ("--dataSet"; arg_type=String; default="UDDataSet"; help="TRDataSet or UDDataSet")
        ("--langcode"; arg_type=String; default="fi"; help="lang code for UDDataSet")
        ("--version"; arg_type=Int; default=2018; help="2006 | 2016 | 2018 for TRDataSet")
        ("--bestModel"; default=Morse.dir("checkpoints/bestModel"); help="path to best model")
        ("--modelFile"; arg_type=String; help="load a pretrained model")
        ("--sourceModel"; arg_type=String)
        ("--trainSize"; arg_type=Int; default=typemax(Int))
        ("--logFile"; arg_type=String; default=Morse.dir("logs/morse"); help="log file")
        ("--genFile"; arg_type=String; default=Morse.dir("generations/morse"); help="file for generations")
        ("--hiddenSizes"; nargs='+';arg_type=Int; default=[512,512,512,512];
            help="char-encoder, context-encoder, decoder-l1, decoder-l2")
        ("--embedSizes"; nargs='+';arg_type=Int; default=[64,256])
            help="char and output embeddings"
        ("--batchSize"; arg_type=Int; default=1; help="batch size for training")
        ("--atype"; default=(gpu()>=0 ? "KnetArray{Float32}" : "Array{Float32}"))
        ("--seed"; arg_type=Int; default=31; help="random number seed.")
        ("--previous"; arg_type=Int; default=2; help="output encoder window")
        ("--dropouts"; arg_type=Float64; default= 0.3; help="dropout probabilities")
        ("--decayRate"; arg_type=Float64; default=0.8; help="learning rate decay")
        ("--optimizer"; arg_type=String; default="SGD(;lr=1.6,gclip=60.0)")
        ("--patience"; arg_type=Int; default=10; help="number of validation to wait")
        ("--mode"; arg_type=Int;default=1; help="1=training, 2=evaluate, 3=transfer")
        ("--epochs"; arg_type=Int; default=1; help="# of epochs for training.")
        ("--wordLimit"; arg_type=Int; default=100; help="maximum # of words in a sentence.")
        ("--threshold"; arg_type=Float64; default=0.0; help="threshold to be considered")
        ("--modelType"; arg_type=String; default="MorseModel";
            help="MorseModel | S2S | S2SContext | Classifier")
        ("--lemma"; action = :store_true;
            help = "predict lemma+morph feats (only for sequential models)")
    end
    config = parse_args(args,s;as_symbols=true)

    lang   = config[:dataSet]   == "TRDataSet"  ? string("TR-tr",config[:version]) : string("UD-",config[:langcode])
    tsize  = config[:trainSize] == typemax(Int) ? "full"  : string(config[:trainSize],"sent")

    postfix = string(".",config[:modelType],"_lemma_",config[:lemma],"_lang_",lang,"_size_",tsize)
    config[:bestModel] = config[:bestModel]*postfix
    config[:logFile] = config[:logFile]*postfix
    config[:genFile] = config[:genFile]*postfix
    config[:bestTestAcc] = config[:bestDevAcc] = 0.0

    if isfile(config[:bestModel]*".jld2") && config[:modelFile] === nothing
        config[:modelFile] = config[:bestModel]*".jld2"
    end

    return config
end

"""
    main(ARGS=[]; config = intro(ARGS))

    Main procedure according to config. see `intro`
"""
function main(ARGS=[]; config = intro(ARGS))

    logFile = config[:logFile]!==nothing ? open(config[:logFile]*".log","a+") : nothing

    DataSet  = eval(Meta.parse(config[:dataSet]))
    download(DataSet) # downloads if dataset is not available in the data/

    ModelType = eval(Meta.parse(config[:modelType]))

    if DataSet === UDDataSet
        code = config[:langcode]
        lang = CODE_TO_LANG[code]
        files = Morse.dir("data","ud-treebanks-v2.1/UD_$(lang)") *
                        "/$(code)-ud-" .* ["train","dev","test"] .* ".conllu"
    else # DataSet === TRDataSet
        vers = config[:version]
        splits = vers==2018 ? [".train"] : [".train",".test"]
        files = Morse.dir("data","TrMor2018")*"/TrMor$(vers)/trmor$(vers)".*splits
    end

    data, vocab, parser = prepareData(files, DataSet; vers=config[:version],
                                                      tsize=config[:trainSize],
                                                      seed=config[:seed],
                                                      withLemma=config[:lemma],
                                                      parseAll=(ModelType<:Disambiguator))

    data = miniBatch(data, vocab)

    ArrayType = eval(Meta.parse(config[:atype]))
    KnetLayers.settype!(ArrayType)
    KnetLayers.Knet.seed!(config[:seed])

    printConfig(logFile,config)
    printLog(logFile,"dataLengths (train,dev,test): ",length.(data))
    printLog(logFile,"vocabLengths", (chars=length(vocab.chars),
                                      tags=length(vocab.tags),
                                      comptags=length(vocab.comptags),
                                      words=length(vocab.words)))

    if config[:mode] == 1 # train
        if  config[:modelFile] !== nothing
            model,o,_,_ = loadModel(config[:modelFile])
            config[:bestDevAcc]  = o[:bestDevAcc]
            config[:bestTestAcc] = o[:bestTestAcc]
        else
            model = ModelType(config, vocab)
            setoptim!(model, eval(Meta.parse(config[:optimizer])))
        end
    elseif config[:mode] == 2 # generate
	source,_,sv,_ = loadModel(config[:modelFile])
        if isa(source,ModelType)
            model = source
        else
            model = ModelType(config, vocab)
            transfer!(model,vocab,source,sv)
        end
	source=nothing; sv=nothing;
    elseif config[:mode] == 3
        # transfer learning
        model = ModelType(config, vocab)
        setoptim!(model, eval(Meta.parse(config[:optimizer])))
        if config[:sourceModel] !== nothing
	    !isfile(config[:sourceModel]) && error("couldn't found sourceModel in its location!")
            source,_,sv,_ = loadModel(config[:sourceModel])
            transfer!(model,vocab,source,sv)
            source=nothing; sv=nothing;
        else
            error("source model is not given")
        end
    end
    Knet.gc();
    #Train
    if config[:mode] != 2
        trainmodel!(model, data, config, vocab, parser; logFile=logFile)
        model=nothing; Knet.gc();
        if isfile(config[:bestModel]*".jld2")
            model,_,_,_ = loadModel(config[:bestModel]*".jld2")
        end
    end

    #Final Generation for Test
    open(config[:genFile]*".gen","w") do f
        scores = evaluate(model, last(data), vocab, parser, file=f)
        printLog(logFile, "epoch=final | set=Test | scores: ", scores)
    end

    #Final Evaluation for Logging
    for (set,name) in zip(data,(:Train, :Dev))
        scores = evaluate(model, set, vocab, parser, file=nothing)
        printLog(logFile, "epoch=final | set=$name | scores: ", scores)
    end

    return model, vocab, parser, config, data
end

"""
    trainmodel!(M::Model, data::Vector, o::Dict{Symbol,Any}, v::Vocabulary,
                           p::Parser; logFile::IO=stdout, patience=o[:patience])

    Trains a `Model` on given data.
    Make sure that you `setoptim!` for model parameters before calling `trainepoch!`

    It applies learning rate decay and early stopping according to patience argument
    It also does evaluations and logs them after each epoch.
"""

function trainmodel!(M::Model, data::Vector, o::Dict{Symbol,Any}, v::Vocabulary,
                            p::Parser; logFile::IO=stdout, patience=o[:patience])
    for i=1:o[:epochs]
        trnloss = trainepoch!(M, data[1], v; wordLimit=o[:wordLimit])
        printLog(logFile, "epoch=$i | set=Train | scores: ", (loss=trnloss[1]/trnloss[2],))
        for (k,set) in enumerate((:Dev,:Test))
            scores = evaluate(M, data[k+1], v, p, file=nothing)
            printLog(logFile, "epoch=$i | set=$set | scores: ", scores)
            acc = scores.accuracies.complete
            if acc > o[Symbol(:best,set,:Acc)]
                o[Symbol(:best,set,:Acc)] = acc;
                if set == :Dev
                    patience = o[:patience]
                    saveModel(o[:bestModel]*".jld2", M, o, v, p)
                end
                printLog(logFile,"Current Best $set Accuracies: ", scores.accuracies)
            end
        end

        i%10==0 && saveModel(o[:bestModel]*"_epoch$i.jld2", M, o, v, p)

        if (patience -= 1) == o[:patience] ÷ 2
            lrdecay!(M,o[:decayRate])
        elseif patience < 0
            printLog(logFile,"patience < 0, training stops")
            break
        end

    end
end

"""
    trainepoch!(M::Model, data::Vector{SentenceBatch}, v::Vocabulary; wordLimit=100)

    Trains a `Model` one epoch on given data.
    Make sure that you `setoptim!` for model parameters before calling `trainepoch!``
"""
function trainepoch!(M::Model, data::Vector{SentenceBatch}, v::Vocabulary; wordLimit=100)
    parameters = params(M)
    trainloss = zeros(2)
    for (i,d) in progress(enumerate(shuffle(data)))
        length(d.encodedIOs) > wordLimit && continue
        J = @diff loss(M, d; v=v)
        isnan(value(J)) && continue
        if isa(M,Sequential)
            trainloss .+= [value(J),1] .* sum(d.masks)
        else
            trainloss .+= [value(J),1] .* length(d.encodedIOs)
        end
        for p in parameters
            update!(value(p), grad(J,p), p.opt)
        end
    end
    return trainloss
end

# Knet ambiguity fix (TODO: open Knet issue)
# ERROR: MethodError: update!(::Knet.KnetArray{Float32,3}, ::Nothing, ::Knet.SGD) is ambiguous. Candidates:
#   update!(w::Knet.KnetArray{Float32,N} where N, g, p::Knet.SGD) in Knet at /kuacc/users/dyuret/.julia/dev/Knet/src/update.jl:475
#   update!(w::Knet.KnetArray{Float32,N} where N, g::Nothing, p) in Knet at /kuacc/users/dyuret/.julia/dev/Knet/src/update.jl:547
# Possible fix, define
#   update!(::Knet.KnetArray{Float32,N} where N, ::Nothing, ::Knet.SGD)
Knet.update!(w::KnetArray{Float32}, ::Nothing, ::Knet.SGD) = w


text:
"""
$(TYPEDSIGNATURES)

High-level place-transition nets (HL-PTNet) have a toolspecific structure
defined for token graphics. Contains <tokenposition> tags.
"""
function parse_tokengraphics(node, pntd; kw...)
    nn = nodename(node)
    nn == "tokengraphics" || error("element name wrong: $nn")
    positions = allchildren("tokenposition", node)
    if isempty(positions)
        TokenGraphics() # Empty is legal.
    else
        #TODO: Enforce type sameness of position coordinates? How?
        TokenGraphics(parse_tokenposition.(positions, Ref(pntd); kw...))
    end
end

"""
$(TYPEDSIGNATURES)

Cartesian coordinate relative to containing element.
"""
function parse_tokenposition(node, pntd; kw...)
    nn = nodename(node)
    nn == "tokenposition" || error("element name wrong: $nn")

    parse_graphics_coordinate(node, pntd; kw...)
end

"""
$(TYPEDSIGNATURES)

Arcs, Annotations and Nodes have different graphics semantics.
Return a [`Graphics`](@ref) holding the union of possibilities.
"""
function parse_graphics(node, pntd; kw...)
    @debug node
    nn = nodename(node)
    nn == "graphics" || error("element name wrong: $nn")

    d = PnmlDict(:tag => Symbol(nn),
                 :line => nothing, :positions => Coordinate[], :dimension => nothing,
                 :fill => nothing, :font => nothing, :offset => nothing)
    foreach(elements(node)) do child
        @match nodename(child) begin
            "dimension" => (d[:dimension] = parse_graphics_coordinate(child, pntd; kw...))
            "fill"      => (d[:fill] = parse_graphics_fill(child, pntd; kw...))
            "font"      => (d[:font] = parse_graphics_font(child, pntd; kw...))
            "line"      => (d[:line] = parse_graphics_line(child, pntd; kw...))
            "offset"    => (d[:offset] = parse_graphics_coordinate(child, pntd; kw...))
            "position"  => (push!(d[:positions], parse_graphics_coordinate(child, pntd; kw...)))
            _ => @warn "ignoring <graphics> child '$(child)'"
        end
    end
    Graphics(;
             dim=d[:dimension],
             fill=d[:fill],
             font=d[:font],
             line=d[:line],
             offset=d[:offset],
             position=d[:positions])
end

"""
$(TYPEDSIGNATURES)

Return [`Line`](@ref).
"""
function parse_graphics_line(node, pntd; kw...)
    nn = nodename(node)
    (nn == "line") || error("element name wrong: $nn")

    color = EzXML.haskey(node, "color") ? node["color"] : nothing
    shape = EzXML.haskey(node, "shape") ? node["shape"] : nothing
    style = EzXML.haskey(node, "style") ? node["style"] : nothing
    width = EzXML.haskey(node, "width") ? node["width"] : nothing

    Line(; shape, color, width, style)
end

"""
$(TYPEDSIGNATURES)

Return [`Coordinate`](@ref).
Specification seems to only use integers, we also allow real numbers.
"""
function parse_graphics_coordinate(node, pntd; kw...)
    nn = nodename(node)
    (nn=="position" || nn=="dimension" ||
     nn=="offset" || nn=="tokenposition") || error("element name wrong: $nn")

    EzXML.haskey(node, "x") || throw(MalformedException("$nn missing x", node))
    EzXML.haskey(node, "y") || throw(MalformedException("$nn missing y", node))

    Coordinate(number_value(node["x"]), number_value(node["y"]))
end

"""
$(TYPEDSIGNATURES)

Return [`Fill`](@ref)
"""
function parse_graphics_fill(node, pntd; kw...)
    nn = nodename(node)
    (nn == "fill") || error("element name wrong: $nn")

    clr  = EzXML.haskey(node, "color") ? node["color"] : nothing
    img  = EzXML.haskey(node, "image") ? node["image"] : nothing
    gclr = EzXML.haskey(node, "gradient-color")    ? node["gradient-color"] : nothing
    grot = EzXML.haskey(node, "gradient-rotation") ? node["gradient-rotation"] : nothing

    Fill(color=clr, image=img, gradient_color=gclr, gradient_rotation=grot)
end

"""
$(TYPEDSIGNATURES)

Return [`Font`](@ref).
"""
function parse_graphics_font(node, pntd; kw...)
    nn = nodename(node)
    (nn == "font") || error("element name wrong: $nn")

    align  = EzXML.haskey(node, "align")      ? node["align"] : nothing
    deco   = EzXML.haskey(node, "decoration") ? node["decoration"] : nothing
    family = EzXML.haskey(node, "family")     ? node["family"] : nothing
    rot    = EzXML.haskey(node, "rotation")   ? node["rotation"] : nothing
    size   = EzXML.haskey(node, "size")       ? node["size"] : nothing
    style  = EzXML.haskey(node, "style")      ? node["style"] : nothing
    weight = EzXML.haskey(node, "weight")     ? node["weight"] : nothing

    Font(; family, style, weight, size, decoration=deco, align, rotation=rot)
end


text:
default_beam_family(
    beam_family_instance(
        backend_family(
            revit => RevitFamily(
                "C:\\ProgramData\\Autodesk\\RVT 2017\\Libraries\\US Metric\\Structural Framing\\Wood\\M_Timber.rfa",
                Dict(:width => "b", :height => "d")),
            autocad => )),
        width=1,
        height=2)


text:
# Types
abstract type AbstractCollectiveInferer end

"""
Relaxation labeling object. Stores the parameters necessary for
the algorithm.
"""
struct RelaxationLabelingInferer <: AbstractCollectiveInferer
	maxiter::Int
	tol::Float64
	tf::Function
	κ::Float64
	α::Float64
end

"""
Iterative classification object. Stores the parameters necessary for
the algorithm.
"""
struct IterativeClassificationInferer <: AbstractCollectiveInferer
	maxiter::Int
	tol::Float64
	tf::Function
end

"""
Gibbs sapmpling object. Stores the parameters necessary for
the algorithm.
"""
struct GibbsSamplingInferer <: AbstractCollectiveInferer
	maxiter::Int
	tol::Float64
	tf::Function
	burniter::Int
end

# Show methods
Base.show(io::IO, ci::RelaxationLabelingInferer) = print(io, "Relaxation labeling, maxiter=$(ci.maxiter), tol=$(ci.tol), κ=$(ci.κ), α=$(ci.α)")
Base.show(io::IO, ci::IterativeClassificationInferer) = print(io, "Iterative classification, maxiter=$(ci.maxiter), tol=$(ci.tol)")
Base.show(io::IO, ci::GibbsSamplingInferer) = print(io, "Gibbs sampling, maxiter=$(ci.maxiter), tol=$(ci.tol), burniter=$(ci.burniter)")
Base.show(io::IO, vci::T) where T<:AbstractVector{S} where S<:AbstractCollectiveInferer =
	print(io, "$(length(vci))-element Vector{$S} ...")

# Transform methods
function transform!(Xo::T, Ci::RelaxationLabelingInferer, obsdim::OD, Mr::M, fr_exec::E, RL::R, Adj::A, offset::Int, Xr::S,
		    update::BitVector=trues(nobs(Xo))) where {
		M, E,
		OD<:LearnBase.ObsDimension,
		T<:AbstractMatrix, R<:Vector{<:AbstractRelationalLearner},
		A<:Vector{<:AbstractAdjacency}, S<:AbstractMatrix}

	# Initializations
	n = nobs(Xo, obsdim)									# number of observations
	κ = Ci.κ										# a constant between 0 and 1
	α = Ci.α										# decay
	β = κ											# weight of current iteration estimates
	maxiter = Ci.maxiter									# maximum number of iterations
	tol = Ci.tol										# maximum error
	f_targets = Ci.tf									# function used to obtain targets
	size_out = nvars(Xo,obsdim)								# ouput size (corresponds to the number of classes)
	Xl = copy(Xo)										# local estimates
	ŷ = f_targets(Xo)									# Obtain first the labels corresponding to the local model
	ŷₒ = similar(ŷ)										#   and the 'previous' iteration estimates
	AV = adjacency_matrix.(Adj)								# Pre-calculate adjacency matrices
	Xrᵢ = matrix_prealloc(n, size_out, obsdim, 0.0)						# Initialize temporary storage

	# Iterate
	_Xo = datasubset(Xo, update, obsdim)
	_Xr = datasubset(Xr, update, obsdim)
	for it in 1:maxiter
		β = β * α									# Update learning rate
		copy!(ŷₒ, ŷ);									# Update 'previous iteration' estimates

		# Obtain relational dataset for the current iteration
		@inbounds for (i,(RLᵢ,Aᵢ)) in enumerate(zip(RL,AV))

			# Apply relational learner
			transform!(Xrᵢ, RLᵢ, Aᵢ, Xo, ŷ)

			# Update relational data output
			_Xrᵢ= datasubset(Xr, offset+(i-1)*size_out+1 : offset+i*size_out, oppdim(obsdim))
			_Xrᵢ[:] = Xrᵢ
		end

		# Update estimates

		_Xo[:] = β.*fr_exec(Mr, _Xr) + (1.0-β).*_Xo
		ŷ = f_targets(Xo)

		# Convergence check
		if isequal(ŷ,ŷₒ) || mean(abs.(ŷ-ŷₒ))<=tol
			@print_verbose 1 "Convergence reached at iteration $it."
			break
		else
			@print_verbose 2 "\tIteration $it: $(sum(ŷ.!= ŷₒ)) estimates changed"
	   	end

		# Replace non-converging estimates with local estimates
		if (it == maxiter) && (maxiter != 1)
			_nc = ŷ.!=ŷₒ 		# positions of non-converging estimates
			datasubset(Xo, _nc, obsdim)[:] = datasubset(Xl, _nc, obsdim)[:]
			@print_verbose 1 "Maximum level of iterations reached, $(sum(_nc)) estimates did not converge."
		end
	end

	return Xo
end

function transform!(Xo::T, Ci::IterativeClassificationInferer, obsdim::OD, Mr::M, fr_exec::E, RL::R, Adj::A, offset::Int, Xr::S,
		    update::BitVector=trues(nobs(Xo))) where {
		M, E,
		OD<:LearnBase.ObsDimension,
		T<:AbstractMatrix, R<:Vector{<:AbstractRelationalLearner},
		A<:Vector{<:AbstractAdjacency}, S<:AbstractMatrix}

	# Initializations
	n = nobs(Xr,obsdim)									# number of observations
	ordering = [i:i for i in findall(update)]						# observation estimation order
	maxiter = Ci.maxiter									# maximum number of iterations
	tol = Ci.tol										# maximum error
	f_targets = Ci.tf									# function used to obtain targets
	size_out = nvars(Xo,obsdim)								# ouput size (corresponds to the number of classes)
	Xl = copy(Xo)										# local estimates
	ŷ = f_targets(Xo)									# Obtain first the labels corresponding to the local model
	ŷₒ = similar(ŷ)										#   and the 'previous' iteration estimates
	AV = adjacency_matrix.(Adj)								# Pre-calculate adjacency matrices
	Xrᵢⱼ = matrix_prealloc(1, size_out, obsdim, 0.0)					# Initialize temporary storage
	second_dim = ObsDim.Constant{2}()

	# Iterate
	for it in 1:maxiter
		shuffle!(ordering)								# Randomize observation order
		copy!(ŷₒ, ŷ);									# Update 'previous iteration' estimates

		# Loop over observations and obtain individual estimates
		for rⱼ in ordering

			# Get data subsets pertinent to the current observation
			Xrⱼ = datasubset(Xr, rⱼ, obsdim)
			Xoⱼ = datasubset(Xo, rⱼ, obsdim)
			ŷⱼ = datasubset(ŷ, rⱼ)

			# Obtain relational data for the current observation
			@inbounds for (i,(RLᵢ,Aᵢ)) in enumerate(zip(RL,AV))

				# Apply relational learner
				Aᵢⱼ = adjacency_obs(Aᵢ, rⱼ, second_dim)
				transform!(Xrᵢⱼ, RLᵢ, Aᵢⱼ, Xo, ŷ)

				# Update relational data output for the current sample
				_Xrⱼ = datasubset(Xrⱼ, offset+(i-1)*size_out+1 : offset+i*size_out, oppdim(obsdim))
				_Xrⱼ[:] = Xrᵢⱼ
			end

			# Update estimates
			Xoⱼ[:] = fr_exec(Mr, Xrⱼ)
			ŷⱼ[:] = f_targets(Xoⱼ)
		end

		# Convergence check
		if isequal(ŷ,ŷₒ) || mean(abs.(ŷ-ŷₒ))<=tol
			@print_verbose 1 "Convergence reached at iteration $it."
			break
		else
			@print_verbose 2 "\tIteration $it: $(sum(ŷ.!= ŷₒ)) estimates changed"
	   	end

		# Replace non-converging estimates with local estimates
		if (it == maxiter) && (maxiter != 1)
			_nc = ŷ.!=ŷₒ 		# positions of non-converging estimates
			datasubset(Xo, _nc, obsdim)[:] = datasubset(Xl, _nc, obsdim)[:]
			@print_verbose 1 "Maximum level of iterations reached, $(sum(_nc)) estimates did not converge."
		end
	end

	return Xo
end

# Version of Gibbs sampling (experimental) similar to iterative classification (i.e. no sampling) from:
# P. Sen, G. Namata, M. Bilgic, L. Getoor, B. Gallagher, T. Eliassi-Rad "Collective classification in network data", AI Magazine 29(3), 2008

# Another (slower) alternative would be to assign to each observation a class sampled
# in accordance to the class-wise probabilities of its neighbourhood. This however
# implies 1. sampling (slow), 2. this would work only for relational learners that
# make use of the neighbourhood class estimates (i.e. :rn and :bayesrn only)
function transform!(Xo::T, Ci::GibbsSamplingInferer, obsdim::OD, Mr::M, fr_exec::E, RL::R, Adj::A, offset::Int, Xr::S,
		    update::BitVector=trues(nobs(Xo))) where {
		M, E,
		OD<:LearnBase.ObsDimension,
		T<:AbstractMatrix, R<:Vector{<:AbstractRelationalLearner},
		A<:Vector{<:AbstractAdjacency}, S<:AbstractMatrix}

	# Initializations
	n = nobs(Xr, obsdim)									# number of observations
	updateable = findall(update)
	ordering = [i:i for i in updateable]							# observation estimation order
	maxiter = Ci.maxiter									# maximum number of iterations
	burniter = Ci.burniter									# number of burn-in iterations
	tol = Ci.tol										# maximum error
	f_targets = Ci.tf									# function used to obtain targets
	size_out = nvars(Xo,obsdim)								# ouput size (corresponds to the number of classes)
	Xl = copy(Xo)										# local estimates
	ŷ = f_targets(Xo)									# Obtain first the labels corresponding to the local model
	ŷₒ = similar(ŷ)										#   and the 'previous' iteration estimates
	AV = adjacency_matrix.(Adj)								# Pre-calculate adjacency matrices
	Xrᵢⱼ = matrix_prealloc(1, size_out, obsdim, 0.0)					# Initialize temporary storage
	second_dim = ObsDim.Constant{2}()

	# Burn-in
	@print_verbose 2 "\tRunning $burniter burn-in iterations ..."
	for it in 1:burniter
		shuffle!(ordering)								# Randomize observation order

		# Loop over observations and obtain individual estimates
		for rⱼ in ordering

			# Get data subsets pertinent to the current observation
			Xrⱼ = datasubset(Xr, rⱼ, obsdim)
			Xoⱼ = datasubset(Xo, rⱼ, obsdim)
			ŷⱼ = datasubset(ŷ, rⱼ)

			# Obtain relational data for the current observation
			@inbounds for (i,(RLᵢ,Aᵢ)) in enumerate(zip(RL,AV))

				# Apply relational learner
				Aᵢⱼ = adjacency_obs(Aᵢ, rⱼ, second_dim)
				transform!(Xrᵢⱼ, RLᵢ, Aᵢⱼ, Xo, ŷ)

				# Update relational data output for the current sample
				_Xrⱼ = datasubset(Xrⱼ, offset+(i-1)*size_out+1 : offset+i*size_out, oppdim(obsdim))
				_Xrⱼ[:] = Xrᵢⱼ
			end

			# Update estimates
			Xoⱼ[:] = fr_exec(Mr, Xrⱼ)
			ŷⱼ[:] = f_targets(Xoⱼ)
		end
	end

	# Initialize class-counting structure
	class_counts = matrix_prealloc(n, size_out, obsdim, 0.0)

	# Small function that makes the class count work (even though it does not make sense)
	# for cases outside classification (i.e. input labels are floats)
	_idx_(x::AbstractVector{Int}) = x[1]
	_idx_(x::AbstractVector) = 1
	# Iterate
	@print_verbose 2 "\tRunning $maxiter iterations ..."
	for it in 1:maxiter
		shuffle!(ordering)								# Randomize observation order
		copy!(ŷₒ, ŷ);									# Update 'previous iteration' estimates

		# Loop over observations and obtain individual estimates
		for rⱼ in ordering

			# Get data subsets pertinent to the current observation
			Xrⱼ = datasubset(Xr, rⱼ, obsdim)
			Xoⱼ = datasubset(Xo, rⱼ, obsdim)
			ŷⱼ = datasubset(ŷ, rⱼ)
			cⱼ = datasubset(class_counts, rⱼ, obsdim)

			# Obtain relational data for the current observation
			@inbounds for (i,(RLᵢ,Aᵢ)) in enumerate(zip(RL,AV))

				# Apply relational learner
				Aᵢⱼ = adjacency_obs(Aᵢ, rⱼ, second_dim)
				transform!(Xrᵢⱼ, RLᵢ, Aᵢⱼ, Xo, ŷ) 				# TODO: Find a better compromise for adjacency access; views - slow for sparse matrices
												#	slicing - increases the number of allocations.
				# Update relational data output for the current sample
				_Xrⱼ = datasubset(Xrⱼ, offset+(i-1)*size_out+1 : offset+i*size_out, oppdim(obsdim))
				_Xrⱼ[:] = Xrᵢⱼ
			end

			# Update estimates
			Xoⱼ[:] = fr_exec(Mr, Xrⱼ)
			ŷⱼ[:] = f_targets(Xoⱼ)
			cⱼ[_idx_(ŷⱼ)]+=1.0
		end

		# Convergence check
		if isequal(ŷ,ŷₒ) || mean(abs.(ŷ-ŷₒ))<=tol
			@print_verbose 1 "Convergence reached at iteration $it."
			break
		else
			@print_verbose 2 "\tIteration $it: $(sum(ŷ.!= ŷₒ)) estimates changed"
	   	end
	end

	# Assign new estimates
	_Xo = datasubset(Xo, updateable, obsdim)
	_cc = datasubset(class_counts, updateable, obsdim)
	_Xo[:] = _cc ./ sum(_cc, dims=intdim(oppdim(obsdim)))

	return Xo
end


text:
timet = LinRange(0,1,101);
M = length(timet);
N = 5;
lam = 0.001;
center = [.3; .35];
center2 = [4; 4];
sd1 = .1;
gam_sd = 1;
num_comp = 5;
f_orig = zeros(M,N*length(center));
omega = 2*pi;
cnt = 1;
for ii = 1:length(center)
    global cnt
    tmp = Normal(center[ii], .1);
    tmp2 = Normal(center2[ii], sd1);
    for jj = 1:N
        f_orig[:, cnt] = rand(tmp2) * pdf.(tmp, timet);
        cnt+= 1;
    end
end

q_orig = f_to_srsf(f_orig,collect(timet));
y_orig = ones(size(q_orig,2));
y_orig[N+1:end] .= -1;

f = zeros(M, size(f_orig,2));
q = zeros(M, size(f_orig,2));
gam_orig = rgam(M, gam_sd, 2*N);
cnt = 1;
for ii = 1:length(center)
    global cnt
    for jj = 1:N
        f[:,cnt] = warp_f_gamma(collect(timet), f_orig[:,cnt], gam_orig[:,jj]);
        q[:,cnt] = f_to_srsf(f[:,cnt], collect(timet));
        cnt += 1;
    end
end


text:
using Dates, Random, Test
using MultiStartLH

mdl = MultiStartLH;

function test_optimizer(guessV)
    solnV = guessV .- 0.1;
    fVal = sum(solnV);
    exitFlag = :success;
    fValStart = fVal + 1.0;
    nEvals = round(Int, sum(guessV));
    p = Point(guessV, solnV, fValStart, fVal, :success, Second(4), nEvals);
    return p
end

function opt_test()
    @testset "MultiStart" begin
        rng = MersenneTwister(434);
        historyFn = mdl.test_history_path();
        maxHours = 0.7;
        fTol = 0.05;
        fTolN = 3;
        nPoints = 20;
        nParams = 4;
        guessesV = [1.0 .+ rand(rng, nParams)  for j = 1 : nPoints];
        ttr = TikTakRule(; nPoints);

        m = MultiStart(test_optimizer, ttr, guessesV, historyFn, maxHours,
            fTol, fTolN);
        solnV, fVal, exitFlag = optimize(m);
        @test solnV isa Vector{Float64};
        @test fVal > 0.0;
        @test exitFlag isa Symbol;
    end
end

@testset "Optimization" begin
    opt_test();
end

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


text:
using LESbrary.Utils

@testset "fit_cubic" begin
    c = fit_cubic((0, 0), (1, 1), 0, 0)

    # We should get back 3x^2 - 2x^3
    @test c[1] ≈ 0 atol=1e-12
    @test c[2] ≈ 0 atol=1e-12
    @test c[3] ≈ 3
    @test c[4] ≈ -2
end


text:
######### function for a sigmoid curve
### dose input
# x is a vector of concentrations
### parameters for a sigmoid curve
# slope
# upper asymptote (fixed)
# lower asymptote (varying)
# inflection point (xmid)

function log_logistic(upper = 1.0; dose, slope::Float64, lower::Float64, xmid::Float64)
    num = upper - lower
    denom = 1 .+ ((dose ./ xmid) .^ slope)
    return lower .+ num ./ denom
end

doses = [0:1:100.0;]

# example to test
responses = log_logistic(dose = doses, slope = 5.0, lower = 0.0, xmid = 50.0)
plot(doses, responses, ylim = 1.0, label = "b=5, c=0, d=1, e=50")

### Vary slope ###
plt1 = plot([NaN], [NaN], label = "", title = "Varying slopes", xlabel = "dose", ylabel = "response")

# range of slope values
slopes = [4.0, 8.0, 16.0, 32.0, 64.0]
for sl in slopes
    responses = log_logistic(dose = doses, slope = sl, lower = 0.0, xmid = 50.0)
    plot!(doses, responses, label = "$sl")
end
plt1

### Vary asymptote
plt2 = plot([NaN], [NaN], label = "", title = "Varying asymptotes"
, xlabel = "dose", ylabel = "response")

# range of asymptote values
asymptotes = [0.0, 0.1, 0.2, 0.3, 0.4]
for as in asymptotes
    responses = log_logistic(dose = doses, slope = 8.0, lower = as, xmid = 50.0)
    plot!(doses, responses, label = "$as")
end
plt2

### Vary inflection point
plt3 = plot([NaN], [NaN], label = "", title = "Varying inflection points"
, xlabel = "dose", ylabel = "response")

# range of xmid values
inflections = [20.0, 30.0, 40.0, 50.0, 60.0]
for infl in inflections
    responses = log_logistic(dose = doses, slope = 8.0, lower = 0.0, xmid = infl)
    plot!(doses, responses, label = "$infl")
end
plt3

plot(plt1, plt2, plt3)

### loop for high low combos
plt4 = plot([NaN], [NaN], label = "", title = "High/low parameters")

slopes2 = [5.0, 50.0]
asymps2 = [0.0, 0.5]
inflex2 = [10.0, 50.0]

for b in slopes2
    for d in asymps2
        for e in inflex2
            responses = log_logistic(dose = doses, slope = b, lower = d, xmid = e)
            plot!(doses, responses, label = "b = $b, d = $d, e = $e")
        end
    end
end
plt4

# function Hill_equation(; dose::Vector{Float64}, slope::Float64, asym::Float64, xmid::Float64)
#     xb = dose .^ (-slope)
#     num = (1-asym) .* xb
#     denom = xb .+ (xmid .^ (-slope))
#     sigmoid = num ./ denom
#     return sigmoid
# end

# responses = Hill_equation(x = doses, slope = 10.0, asym = 1.0, xmid = 30.0)
# plot(doses, responses, ylim = 1.0)

max_stars_repo_name: hanamayall/simple-feeding-DR

max_issues_repo_name: hanamayall/simple-feeding-DR

max_forks_repo_name: hanamayall/simple-feeding-DR


text:
#=
   algorithms from
   Mioara Joldes, Jean-Michel Muller, Valentina Popescu.
   Tight and rigourous error bounds for basic building blocks of double-word arithmetic.
   ACM Transactions on Mathematical Software, Association for Computing Machinery,
   2017, 44 (2), pp.1 - 27. <10.1145/3121432>. <hal-01351529v3>
   they recommend
   > For adding two double-word numbers, use Algorithm 6.
     Use Algorithm 5 only if both operands have the same sign.
  > For multiplying a double-word number by a floating-point number,
    where FMA is available, use Algorithm 9.
  > For multiplying two double-word numbers, if an FMA instruction
    is available, then Algorithm 12 is to be favored.
  > For dividing a double-word number by a floating-point number,
    use Algorithm 15.
  > For dividing two double-word numbers: Algorithm 17 is suggested.
    If an FMA instruction is available, and accuracy is important,
    prefer Algorithm 18 (it is slower).
=#

# Algorithm 1 in ref: error-free transformation

@inline function Fast2Sum(a::T, b::T) where {T<:AbstractFloat}
    s = a + b
    z = s - a
    t = b - z
    return s, t
end

# Algorithm 2 in ref: error-free transformation

function TwoSum(a::T, b::T) where {T<:AbstractFloat}
    s  = a + b
    a1 = s - b
    b1 = s - a1
    da = a - a1
    db = b - b1
    t  = da + db
    return s, t
end

function TwoDiff(a::T, b::T) where {T<:AbstractFloat}
    s  = a - b
    a1 = s + b
    b1 = s - a1
    da = a - a1
    db = b - b1
    t  = da + db
    return s, t
end

# Algorithm 3 in ref: error-free transformation

@inline function Fast2Mult(a::T, b::T) where {T<:AbstractFloat}
    s = a * b
    t = fma(a, b, -s)
    return s, t
end

# Algorithm 4 in ref: relerr 2u²  [reltime 26]

@inline function DWPlusFP(xₕᵢ::T, xₗₒ::T, y::T) where {T<:AbstractFloat}
    sₕᵢ, sₗₒ = TwoSum(xₕᵢ, y)
    v = xₗₒ + sₗₒ
    zₕᵢ, zₗₒ = TwoSum(sₕᵢ, v)
    return zₕᵢ, zₗₒ
end

@inline function DWMinusFP(xₕᵢ::T, xₗₒ::T, y::T) where {T<:AbstractFloat}
    sₕᵢ, sₗₒ = TwoDiff(xₕᵢ, y)
    v = xₗₒ + sₗₒ
    zₕᵢ, zₗₒ = TwoSum(sₕᵢ, v)
    return zₕᵢ, zₗₒ
end

# Algorithm 6 in ref: relerr  3u² + 13u³  [reltime 35]

function AccurateDWPlusDW(xₕᵢ::T, xₗₒ::T, yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
   sₕᵢ, sₗₒ = TwoSum(xₕᵢ, yₕᵢ)
   tₕᵢ, tₗₒ = TwoSum(xₗₒ, yₗₒ)
   c = sₗₒ + tₕᵢ
   vₕᵢ, vₗₒ = Fast2Sum(sₕᵢ, c)
   w = tₗₒ + vₗₒ
   zₕᵢ, zₗₒ = Fast2Sum(vₕᵢ, w)
   return zₕᵢ, zₗₒ
end

function AccurateDWMinusDW(xₕᵢ::T, xₗₒ::T, yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
   sₕᵢ, sₗₒ = TwoDiff(xₕᵢ, yₕᵢ)
   tₕᵢ, tₗₒ = TwoDiff(xₗₒ, yₗₒ)
   c = sₗₒ + tₕᵢ
   vₕᵢ, vₗₒ = Fast2Sum(sₕᵢ, c)
   w = tₗₒ + vₗₒ
   zₕᵢ, zₗₒ = Fast2Sum(vₕᵢ, w)
   return zₕᵢ, zₗₒ
end

# Algorithm 7 in ref: relerr (³/₂)u² + 4u³  [reltime 18]

@inline function DWTimesFP1(xₕᵢ::T, xₗₒ::T, y::T) where {T<:AbstractFloat}
    cₕᵢ, c1 = Fast2Mult(xₕᵢ, y)
    c2 = xₗₒ * y
    tₕᵢ, t1 = Fast2Sum(cₕᵢ, c2)
    t2 = t1 + c1
    zₕᵢ, zₗₒ = Fast2Sum(tₕᵢ, t2)
    return zₕᵢ, zₗₒ
end

# Algorithm 9 in ref: relerr 2u²  [reltime 15]

@inline function DWTimesFP3(xₕᵢ::T, xₗₒ::T, y::T) where {T<:AbstractFloat}
    cₕᵢ, cₗₒ = Fast2Mult(xₕᵢ, y)
    c = fma(xₗₒ, y, cₗₒ)
    zₕᵢ, zₗₒ = Fast2Sum(cₕᵢ, c)
    return zₕᵢ, zₗₒ
end

# Algorithm 11 in ref: relerr 6u²  [reltime 16]

function DWTimesDW2(xₕᵢ::T, xₗₒ::T, yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
   cₕᵢ, c1 = Fast2Mult(xₕᵢ, yₕᵢ)
   t0 = xₗₒ * yₗₒ
   c2 = fma(xₗₒ, yₕᵢ, t0)
   c3 = c1 + c2
   zₕᵢ, zₗₒ = Fast2Sum(cₕᵢ, c3)
   return zₕᵢ, zₗₒ
end

# Algorithm 12 in ref: relerr 5u²  [reltime 17]

function DWTimesDW3(xₕᵢ::T, xₗₒ::T, yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
   cₕᵢ, c1 = Fast2Mult(xₕᵢ, yₕᵢ)
   t0 = xₗₒ * yₗₒ
   t1 = fma(xₕᵢ, yₗₒ, t0)
   c2 = fma(xₗₒ, yₕᵢ, t1)
   c3 = c1 + c2
   zₕᵢ, zₗₒ = Fast2Sum(cₕᵢ, c3)
   return zₕᵢ, zₗₒ
end

# Algorithm 15 in ref: relerr 3u²  [reltime 27]

function DWDivFP3(xₕᵢ::T, xₗₒ::T, y::T) where {T<:AbstractFloat}
   tₕᵢ = xₕᵢ / y
   pₕᵢ, pₗₒ = Fast2Mult(tₕᵢ, y)
   dₕᵢ = xₕᵢ - pₕᵢ
   dₗₒ = dₕᵢ - pₗₒ
   d = dₗₒ + xₗₒ
   tₗₒ = d / y
   zₕᵢ, zₗₒ = Fast2Sum(tₕᵢ, tₗₒ)
   return zₕᵢ, zₗₒ
end

# Algorithm 17 in ref: relerr 15u² + 56u³ [reltime 50]

function DWDivDW2(xₕᵢ::T, xₗₒ::T, yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
   tₕᵢ = xₕᵢ / yₕᵢ
   rₕᵢ, rₗₒ = DWTimesFP1(yₕᵢ, yₗₒ, tₕᵢ)
   dₕᵢ = xₕᵢ - rₕᵢ
   dₗₒ = xₗₒ - rₗₒ
   d = dₕᵢ + dₗₒ
   tₗₒ = d / yₕᵢ
   zₕᵢ, zₗₒ = Fast2Sum(tₕᵢ, tₗₒ)
   return zₕᵢ, zₗₒ
end

# Algorithm 18 in ref: relerr < 10u² (6u² seen) [reltime 107]
# (note DWTimesDW3 replaces DWTimesDW2 per ref)

function DWDivDW3(xₕᵢ::T, xₗₒ::T, yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
   tₕᵢ = inv(yₕᵢ)
   rₕᵢ = fma(yₕᵢ, -tₕᵢ, one(T))
   rₗₒ = -(yₗₒ * tₕᵢ)
   eₕᵢ, eₗₒ = Fast2Sum(rₕᵢ, rₗₒ)
   dₕᵢ, dₗₒ = DWTimesFP3(eₕᵢ, eₗₒ, tₕᵢ)
   mₕᵢ, mₗₒ = DWPlusFP(dₕᵢ, dₗₒ, tₕᵢ)
   zₕᵢ, zₗₒ = DWTimesDW3(xₕᵢ, xₗₒ, mₕᵢ, mₗₒ)
   return zₕᵢ, zₗₒ
end

# inv(...) using Algorithms 17 and 18

# Algorithm 17 in ref: relerr 15u² + 56u³  [reltime 48]

function DWInvDW2(yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
   tₕᵢ = one(T) / yₕᵢ
   rₕᵢ, rₗₒ = DWTimesFP1(yₕᵢ, yₗₒ, tₕᵢ)
   dₕᵢ = one(T) - rₕᵢ
   dₗₒ = -rₗₒ
   d = dₕᵢ + dₗₒ
   tₗₒ = d / yₕᵢ
   zₕᵢ, zₗₒ = Fast2Sum(tₕᵢ, tₗₒ)
   return zₕᵢ, zₗₒ
end

# Algorithm 18 in ref: relerr < 10u² (6u² seen)  [reltime 72]
# (note DWTimesDW3 replaces DWTimesDW2 per ref)

function DWInvDW3(yₕᵢ::T, yₗₒ::T) where {T<:AbstractFloat}
   tₕᵢ = inv(yₕᵢ)
   rₕᵢ = fma(yₕᵢ, -tₕᵢ, one(T))
   rₗₒ = -(yₗₒ * tₕᵢ)
   eₕᵢ, eₗₒ = Fast2Sum(rₕᵢ, rₗₒ)
   dₕᵢ, dₗₒ = DWTimesFP3(eₕᵢ, eₗₒ, tₕᵢ)
   zₕᵢ, zₗₒ = DWPlusFP(dₕᵢ, dₗₒ, tₕᵢ)
   return zₕᵢ, zₗₒ
end


text:
using Test

#--- tests on SGN cnoidal wave
@testset "SGN cnoidal wave" begin
    # Build in a first way a cnoidal wave
    param = (ϵ  = 0.1, μ = 0.1, h₀ = 0.5, h₁ = 0.8, h₂ = 1.5, N=2^9)
    (η,u,v,mesh,para)=CnoidalWaveSerreGreenNaghdi(param;P=1)
    # Build the cnoidal wave with two periods
    (η2,u2,v2,mesh2,para2)=CnoidalWaveSerreGreenNaghdi(param;P=2)
    # Build in a second way the cnoidal wave
    init = CnoidalSGN(param)
    # Check they correspond
    @test init.η(mesh.x)≈η
    @test init.v(mesh.x)≈v
    @test init.η(mesh2.x)≈η2
    @test init.v(mesh2.x)≈v2
    @test η2[end÷2+1:3*end÷4]≈η[end÷2+1:2:end]
    @test v2[end÷2+1:3*end÷4]≈v[end÷2+1:2:end]

    # Solve the corresponding initial-value problem
    param=merge(param,(L=para.λ,T=1e-2,dt=1e-3))
    model = SerreGreenNaghdi(param)
    problem = Problem(model, init, param)
    solve!(problem;verbose=false)
    # Check the solution is translated
    (η,v,x,t)=solution(problem)
    @test η≈init.η(x.-para.c*param.T)
    @test v≈init.v(x.-para.c*param.T)

end

#--- tests on SGN solitary wave
@testset "SGN solitary wave" begin
    # Build in a first way a solitary wave
    param = (ϵ  = 0.1, μ = 0.1, c=1.1, L = 10, N=2^8)
    (η,u,v,mesh)=SolitaryWaveSerreGreenNaghdi(param)
    # Build the translated solitary wave on a different mesh
    param2 = (ϵ  = 0.1, μ = 0.1, c=1.1, L = 20, N=2^9);x₀=-2param.L/param.N;
    (η2,u2,v2,mesh2)=SolitaryWaveSerreGreenNaghdi(param2;x₀=x₀)
    # Build in a second way the cnoidal wave
    init = SolitarySGN(param2)
    init2 = SolitarySGN(param2.c,ϵ=param2.ϵ,μ=param2.μ,N=2^10)

    # Check they correspond
    @test init.η(mesh.x)≈η
    @test init.v(mesh.x)≈v
    @test init.η(mesh.x)≈init2.η(mesh.x)
    @test init.v(mesh.x)≈init2.v(mesh.x)
    @test init.η(mesh2.x .-x₀)≈η2
    @test init.v(mesh2.x .-x₀)≈v2
    @test η2[end÷4:3*end÷4-1]≈η[1:end]
    @test u2[end÷4:3*end÷4-1]≈u[1:end]
    @test v2[end÷4:3*end÷4-1]≈v[1:end]

    # Solve the corresponding initial-value problem
    param=merge(param2,(T=1e-2,dt=1e-3))
    model = SerreGreenNaghdi(param)
    problem = Problem(model, init, param)
    solve!(problem;verbose=false)
    # Check the solution is translated
    (η,v,x,t)=solution(problem)
    @test η≈init.η(x.-param.c*param.T)
    @test v≈init.v(x.-param.c*param.T)

end

#--- tests on Whitham-Boussinesq solitary wave
@testset "Whitham-Boussinesq solitary wave" begin
    # Build a solitary wave with different parameters
    param = (ϵ  = 0.1, μ = 0.1, c=1.1, L = 10, N=2^8)
    waves = Tuple[]
    push!(waves,SolitaryWaveWhithamBoussinesq(param))
    push!(waves,SolitaryWaveWhithamBoussinesq(param;iterative=true))
    push!(waves,SolitaryWaveWhithamBoussinesq(param;
            iterative=true,max_iter = 15, tol = 1e-9, gtol = 1e-9, ktol = 1e-9))
    push!(waves,SolitaryWaveWhithamBoussinesq(param;dealias=1))
    push!(waves,SolitaryWaveWhithamBoussinesq(param;q=0.9,β=1/2,max_iter=30))
    push!(waves,SolitaryWaveWhithamBoussinesq(param;β=1/2))

    # Check they are consistent
    (η,v,mesh)=waves[1]
    for wave in waves
        (η2,v2,mesh2)=wave
        @test η2≈η
        @test v2≈v
    end

    # Build the translated solitary wave on a different mesh
    param2 = (ϵ  = 0.1, μ = 0.1, c=1.1, L = 20, N=2^9);x₀=-2param.L/param.N;
    (η2,v2,mesh2)=SolitaryWaveWhithamBoussinesq(param2;x₀=x₀)

    # Build in a second way the solitary wave
    init = SolitaryWB(param2)
    init2 = SolitaryWB(param2.c,ϵ=param2.ϵ,μ=param2.μ,N=2^10)

    # Check they correspond
    @test init.η(mesh.x)≈init2.η(mesh.x)
    @test init.v(mesh.x)≈init2.v(mesh.x)
    @test init.η(mesh.x)≈η
    @test init.v(mesh.x)≈v
    @test init.η(mesh2.x .-x₀)≈η2
    @test init.v(mesh2.x .-x₀)≈v2
    @test η2[end÷4:3*end÷4-1]≈η[1:end]
    @test v2[end÷4:3*end÷4-1]≈v[1:end]

    # Solve the corresponding initial-value problem
    param=merge(param2,(T=1e-2,dt=1e-3))
    model = WhithamBoussinesq(param)
    problem = Problem(model, init, param)
    solve!(problem;verbose=false)
    # Check the solution is translated
    (η,v,x,t)=solution(problem)
    @test η≈init.η(x.-param.c*param.T)
    @test v≈init.v(x.-param.c*param.T)

    # Do it again with other Whitham-Boussinesq model
    init = SolitaryWB(param;α=1/2)
    model = WhithamBoussinesq(param;α=1/2)
    problem = Problem(model, init, param)
    solve!(problem;verbose=false)
    (η,v,x,t)=solution(problem)
    @test η≈init.η(x.-param.c*param.T)
    @test v≈init.v(x.-param.c*param.T)

    # Do it again with regular Boussinesq
    init = SolitaryWB(param;Boussinesq=true,a=-1/4,b=1/2)
    model = WhithamBoussinesq(param;Boussinesq=true,a=-1/4,b=1/2)
    problem = Problem(model, init, param)
    solve!(problem;verbose=false)
    (η,v,x,t)=solution(problem)
    @test η≈init.η(x.-param.c*param.T)
    @test v≈init.v(x.-param.c*param.T)

end

#--- tests on Whitham-Green-Naghdi solitary wave
@testset "Whitham-Green-Naghdi solitary wave" begin
    # Build a solitary wave with different parameters
    param = (ϵ  = 0.1, μ = 0.1, c=1.1, L = 10, N=2^8)
    waves = Tuple[]
    push!(waves,SolitaryWaveWhithamGreenNaghdi(param))
    push!(waves,SolitaryWaveWhithamGreenNaghdi(param;iterative=true))
    push!(waves,SolitaryWaveWhithamGreenNaghdi(param;
            iterative=true,max_iter = 15, tol = 1e-9, gtol = 1e-9, ktol = 1e-9))
    push!(waves,SolitaryWaveWhithamGreenNaghdi(param;dealias=1))
    push!(waves,SolitaryWaveWhithamGreenNaghdi(param;q=0.9,α=1/2,max_iter=30))
    push!(waves,SolitaryWaveWhithamGreenNaghdi(param;α=1/2))
    push!(waves,SolitaryWaveWhithamGreenNaghdi(param;method=1))
    push!(waves,SolitaryWaveWhithamGreenNaghdi(param;method=2))
    push!(waves,SolitaryWaveWhithamGreenNaghdi(param;method=3))
    push!(waves,SolitaryWaveWhithamGreenNaghdi(param;method=4))

    # Check they are consistent
    (η,u,v,mesh)=waves[1]
    for wave in waves
        (η2,u2,v2,mesh2)=wave
        @test η2≈η
        @test v2≈v
        @test u2≈u
    end

    # Build the translated solitary wave on a different mesh
    param2 = (ϵ  = 0.1, μ = 0.1, c=1.1, L = 20, N=2^9);x₀=-2param.L/param.N;
    (η2,u2,v2,mesh2)=SolitaryWaveWhithamGreenNaghdi(param2;x₀=x₀)

    # Build in a second way the solitary wave
    init = SolitaryWGN(param2)
    init2 = SolitaryWGN(param2.c,ϵ=param2.ϵ,μ=param2.μ,N=2^10)

    # Check they correspond
    @test init.η(mesh.x)≈init2.η(mesh.x)
    @test init.v(mesh.x)≈init2.v(mesh.x)
    @test init.η(mesh.x)≈η
    @test init.v(mesh.x)≈v
    @test init.η(mesh2.x .-x₀)≈η2
    @test init.v(mesh2.x .-x₀)≈v2
    @test η2[end÷4:3*end÷4-1]≈η[1:end]
    @test u2[end÷4:3*end÷4-1]≈u[1:end]
    @test v2[end÷4:3*end÷4-1]≈v[1:end]

    # Solve the corresponding initial-value problem
    param=merge(param2,(T=1e-2,dt=1e-3))
    model = WhithamGreenNaghdi(param)
    problem = Problem(model, init, param)
    solve!(problem;verbose=false)
    # Check the solution is translated
    (η,v,x,t)=solution(problem)
    @test η≈init.η(x.-param.c*param.T)
    @test v≈init.v(x.-param.c*param.T)

    # Do it again with regular Green-Naghdi
    init = SolitaryWGN(param;SGN=true)
    init2 = SolitarySGN(param)
    @test init.η(x)≈init2.η(x)
    @test init.v(x)≈init2.v(x)
    model = WhithamGreenNaghdi(param;SGN=true)
    problem = Problem(model, init, param)
    solve!(problem;verbose=false)
    (η,v,x,t)=solution(problem)
    @test η≈init.η(x.-param.c*param.T)
    @test v≈init.v(x.-param.c*param.T)

end

#--- tests on Whitham solitary wave
@testset "Whitham solitary wave" begin
    # Build a solitary wave with different parameters
    param = (ϵ  = 0.1, μ = 0.1, c=1.1, L = 10, N=2^9)
    waves = Tuple[]
    push!(waves,SolitaryWaveWhitham(param;α=1))
    push!(waves,SolitaryWaveWhitham(param;iterative=true))
    push!(waves,SolitaryWaveWhitham(param;
            iterative=true,max_iter = 15, tol = 1e-9, gtol = 1e-9, ktol = 1e-9))
    push!(waves,SolitaryWaveWhitham(param;dealias=1))
    push!(waves,SolitaryWaveWhitham(param;q=0.9,α=1/2,max_iter=30))
    push!(waves,SolitaryWaveWhitham(param;α=1/2))

    # Check they are consistent
    (u,mesh)=waves[1]
    for wave in waves
        (u2,mesh2)=wave
        @test u2≈u
    end

    # Build the translated solitary wave on a different mesh
    param2 = (ϵ  = 0.1, μ = 0.1, c=1.1, L = 20, N=2^10);x₀=-2param.L/param.N;
    (u2,mesh2)=SolitaryWaveWhitham(param2;x₀=x₀,α=1)

    # Build in a second way the solitary wave
    init = SolitaryWhitham(param2)
    init2 = SolitaryWhitham(param2.c,ϵ=param2.ϵ,μ=param2.μ,N=2^10)

    # Check they correspond
    @test init.η(mesh.x)≈init2.η(mesh.x)
    @test init.v(mesh.x)≈init2.v(mesh.x)
    @test init.v(mesh.x)≈u
    @test init.v(mesh2.x .-x₀)≈u2
    @test u2[end÷4:3*end÷4-1]≈u[1:end]

end

#--- tests on randomly generated wave
@testset "Random initial data" begin
    param=(N=20,L=2);mesh=Mesh(param);x=mesh.x;
    # Build initial data
    inits=InitialData[]
    push!(inits,Random(param;L=1,s=Inf,λ=Inf,a=(1,1)))
    push!(inits,Random(param;L=1,s=2,λ=Inf,a=(1,1)))
    push!(inits,Random(param;L=1,s=Inf,λ=1,a=(1,1)))
    push!(inits,Random(param;L=1,s=Inf,λ=Inf,a=(-1,0)))

    # Check the data have correct size and none of them yields NaN
    for init in inits
        @test !any(isnan,init.η(x))
        @test !any(isnan,init.v(x))
        @test length(init.η(x))==param.N
        @test length(init.v(x))==param.N
    end
end


text:
using TableScraper
using DataFrames
using Test

@testset "TableScraper.jl" begin
    table = scrape_tables("https://www.agenas.gov.it/covid19/web/index.php?r=site%2Fprovvedimento&q=010")[1] |> DataFrame;

    @test nrow(table) > 0
end

@testset "TableScraper.jl goratings" begin
    table = scrape_tables("https://goratings.org")[2] |> DataFrame;
    @test nrow(table) > 0
end


text:
for (i, testfile) in enumerate((
    "generic.jl",
    "wrapped_rng.jl",
    "PCG/runtests.jl",
    "MersenneTwisters/runtests.jl",
    "Xorshifts/runtests.jl"))
    @eval module $(Symbol("T$i"))
        include("common.jl")
        include($testfile)
    end
end


text:
# function to allocate incidents from the current queue of each department
# to a second queue that checks whether other departments could help out
function fill_exchange_queue!(sim_data::DataFrame,
                                incident_queue::Array{Union{Missing,Int64},3},
                                exchange_queue::Array{Union{Missing,Int64},1},
                                max_queue::Int64,
                                card_districts::Int64,
                                priority::Int64)
    place_in_queue = 1
    exchange_queue .= missing
    for i = 1:card_districts
        for j = 1:max_queue
            # allocate the incident from the main queues to the exchange queue
            if ismissing(incident_queue[j,i,priority]) == false
                current_case = incident_queue[j,i,priority]
                if place_in_queue > length(exchange_queue)
                    error("w_aus is too small!")
                    break
                else
                    if sim_data[current_case,:cars_missing] > 0
                        exchange_queue[place_in_queue] = current_case
                        place_in_queue += 1
                    end
                end
            else
                break
            end
        end
    end
end


text:
# # Construct a skyscraper

using PiCraft

# With the help of programming we can build entire cities with ease. Writing a program to build structures is a great way to minimize repetitive work while experimenting with different parameters. If you are making a small house odds are that writing a program is inefficient and manually building it will get the job done much faster. However, if you are trying to build a skyscraper which spawns 20 stories high then manually it could take several hours depending on your speed but it will definitely be very boring. Writing code to do the same will take a dozen minutes with the added advantage that making changes will be extremely easy.

# To get there first we will start small. In this tutorial we will learn how to build a simple skyscraper.

# Our first task is to decide the dimensions and the parameters. Then think of a strategy/algorithm which will do the construction.

# ## makeSkyscraper()

# * Decide the location for the building
#     `p = getTile()`
# * Decide the dimensions and parameters for the building
#     - floors = 5
#     - length = 10
#     - width = 15
#     - height = 6
# * Write the pseudo-code/algorithm
#     Pseudo-code means to write informal code rather than _actual_ code which a human can understand. It can be helpful in designing algorithms this way.
#
#         x = start_position
#         for i in 1:floors
#             buildFloor(x)
#             increment x's height by floor's height
#         end
#         Build a ceiling for the top floor
#
#         Now to build a floor
#            1) Make all the walls out of glass panes
#            2) Set the bottom level as the flooring block. (Cobblestone for instance)
#            3) Add pillars to the corners

# Let's first define a few basic parameters for our building
p = getTile()
floors = 5
# length will be in the +ve x direction and width in the +ve z direction
len = 10
width = 15
height = 6

# We're first going to write code to construct a single floor relative to a point
# Then we will use a loop to create a skyscraper
# pseudocode:
#    x = start_position
#    for i in 1:floors
#        buildFloor(x)
#        increase x's height by floor's height
#    end
#    Build a ceiling for the top floor
#
# Now to build a floor
#    1) Make all the walls out of glass panes
#    2) Set the bottom level as floor
#    3) Get Pillars on the corners

function buildFloor(p)
    #Set all walls to glass panes
    setBlocks(p .+ (1, 0, 1), p .+ (1, height, width), PiCraft.GLASS_PANE)
    setBlocks(p .+ (1, 0, 1), p .+ (len, height, 1), PiCraft.GLASS_PANE)
    setBlocks(p .+ (1, height, width), p .+(len, 0, width), PiCraft.GLASS_PANE)
    setBlocks(p .+ (len, height, width), p .+ (len, 0, 1), PiCraft.GLASS_PANE)

    #Set entire floor to Stone Brick
    setBlocks(p .+ (1, 0, 1), p .+ (len, 0, width), PiCraft.STONE_BRICK)

    #Set 4 Pillars to the corners
    setBlocks(p .+ (1, 0, 1), p .+ (1, height, 1), PiCraft.STONE_BRICK)
    setBlocks(p .+ (1, 0, width), p .+ (1, height, width), PiCraft.STONE_BRICK)
    setBlocks(p .+ (len, 0, 1), p .+ (len, height, 1), PiCraft.STONE_BRICK)
    setBlocks(p .+ (len, 0 , width), p .+ (len, height, width), PiCraft.STONE_BRICK)
end

function makeSkyscraper()
    for i in 1:floors
        buildFloor(p .+ (0, height*i - height, 0))
    end
    setBlocks(p .+ (1, height*floors, 1), p .+ (len, height*floors, width), PiCraft.IRON_BLOCK) # Set Roof to Iron Block
end

# Executing `makeSkyscraper()` once will construct a building in front of you. Try experimenting around with different dimensions.
#
# ![makeSkyscraper.gif](../assets/img/Skyscraper/makeSkyscraper.gif)
#
# ![makeSkyscraper.png](../assets/img/Skyscraper/makeSkyscraper.png)

# ***

# ## makeSkyscraper2()
# This building serves as a demonstration of how we can break down a construction task. When programming we will generally have a choice of algorithms/strategies to choose among. In our next example we will choose a different algorithm.
#
# * Decide the location, parameters and dimensions like last time.
#     -  `p = getTile()`
#     - floors = 5
#     - length = 10
#     - width = 15
#     - height = 6
# * Pseudocode
#     - Place a cuboid of blocks on the building space
#     - Hollow out the inside
#     - For each floor
#         - Add flooring
#         - Add windows
#
# This building isn't particularly nice and has too much glass. Let's have another go.
# When we use programming to solve a problem there can be multiple solutions to the same
# problem.

# This time let's make the building first and then add windows later

function makeSkyscraper2()
    # Make a filled cuboid of Blocks filling the entire building space
    setBlocks(p .+ (1, 0, 1), p .+ (len, height*floors + 1, width),Block(251, 0))

    # Hollow out the inside of the building
    setBlocks(p .+ (2, 1, 2), p .+ (len - 1, height*floors, width - 1), PiCraft.AIR)

    # Add wooden flooring
    for i in 0:floors
        setBlocks(p .+ (2, height*i, 2), p .+ (len - 1, height*i, width - 1), PiCraft.WOOD_PLANKS)
    end

    # Add window Panes, careful if you change the given building dimensions, the change in coordinates may result
    # in unintended consequences.
    for i in 0:(floors - 1)
        setBlocks(p .+ (1, height*i + 2, floor.((width + 1)/2 - 2)),
        p .+ (1, height*i + 4, floor.((width + 1)/2 + 2)), PiCraft.GLASS_PANE)
        setBlocks(p .+ (1, height*i + 2, floor.((width + 1)/2 - 2)),
        p .+ (1, height*i + 4, floor.((width + 1)/2 + 2)), PiCraft.GLASS_PANE)
    end
end

# ![makeSkyscraper2.png](../assets/img/Skyscraper/makeSkyscraper2.png)

# ---

# To get better at building you should practice and try new structures and dimensions. Make your building beautiful by experimenting with a variety of different colored blocks and parameters. You can find a list of challenges below.

# In our second iteration we made the building better, the code is easier to
# follow through and is less in comparison.
# Now, experiment with writing your own code architecturing your own designs.
# Check out the challenges listed below
# Challenges :
# 1) Make the building well lit with torches
# 2) Add a balcony to each floor.
# 3) A common ladder access for all floors.
# 4) Add randomness to your design
#    Hint: for all blocks in building if block is cobblestone and rand() < 0.1 set block to Mossy Cobblestone
# 5) Try constructing a building with the "setback" design in mind. Given below is a documentary on how the setback principle came about.

# [![How New York got its Skyline](http://img.youtube.com/vi/lGroIrQmwyw/0.jpg)](http://www.youtube.com/watch?v=lGroIrQmwyw)


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using JuLIP, ACE, JuLIP.Potentials, LinearAlgebra
using ACE.Testing: lsq, get_V0
using LinearAlgebra: qr, cond
using JuLIP: Atoms
using Plots

function rand_config(V; rattle = 0.2, nrepeat = 3)
   at = bulk(:W, cubic=true, pbc=false) * nrepeat
   return rattle!(at, rattle)
end

function get_2b_basis(species; maxdeg = 10, rcut = 8.0 )
   basis = pair_basis(species = species,
      r0 = rnn(Symbol(species)),
      maxdeg = maxdeg,
      rcut = rcut)
   return basis
end

function calc_rmse(train, Vfit)
   sqe = 0.0
   sqy = 0.0
   for (at, E, F) in train
      sqy += (E/length(at))^2 + norm(mat(F)[:])^2
      Efit = energy(Vfit, at)
      Ffit = forces(Vfit, at)
      sqe += ((E - Efit)/length(at))^2 +
             norm( mat(F)[:] - mat(Ffit)[:] )^2
   end
   return sqrt(sqe / sqy)
end

function trainset(Vref, Ntrain; kwargs...)
   train = []
   for n = 1:Ntrain
      at = rand_config(Vref)
      push!(train, (at = at, E = energy(pot, at), F = forces(pot, at)))
   end
   V0 = get_V0(train)
   @show V0
   train = [ (at = at, E = E - energy(V0, at), F = F)
               for (at, E,  F) in train ]
   return train
end

#set potential
pot = LennardJones(rnn(:W), 1.0) * C2Shift(6.0)

#create train database
train = trainset(pot, 1000)

#set up 2b basis, and perform the least squares problem
basis = get_2b_basis(:W; maxdeg = 14, rcut = 6.0)
c, relrmse = lsq(train, basis, wE=30.0, wF=1.0)

#combine with coefficients to get IP model
IP = JuLIP.MLIPs.combine(basis, c)

E_pot = []
E_ip = []
R = []

for r in range(1,10,length=1000)
   X = [[0.0,0.0,0.0], [0.0, 0.0, r]]
   C = [[100.0,0.0,0.0], [0.0, 100.0, 0.0],[0.0, 0.0, 100.0] ]
   at = Atoms(X, [[0.0,0.0,0.0], [0.0, 0.0, 0.0]], [0.0, 0.0], AtomicNumber.([74, 74]), C, false)
   push!(R,r)
   push!(E_ip, energy(IP, at))
   push!(E_pot, energy(pot,at))
end

plot(R,E_pot, label="LJ")
plot!(R,E_ip, label="FIT")
ylims!((-1,1))
xlims!((2,8))

#FORCE SCATTER PLOT ON TEST SET
F_ip = []
F_2b = []

for i in 1:10
   at = rand_config(pot)
   push!(F_ip, vcat(forces(IP, at)...))
   push!(F_2b, vcat(forces(pot, at)...))
end

scatter(vcat(F_ip...), vcat(F_2b...))


text:
using Weave
# using Fontconfig, Cairo
# using PhyloPlots

set_chunk_defaults(Dict{Symbol, Any}(:results => "hidden", :eval => false))

files_to_weave = ["installation.jmd",
                  "inputdata.jmd",
                  "trait_tree.jmd",
                  "snaq_plot.jmd",
                  "dist_reroot.jmd",
                  "fixednetworkoptim.jmd",
                  "expectedCFs.jmd",
                  "bootstrap.jmd"]

for file in files_to_weave
    println(file)
    weave(Pkg.dir("PhyloNetworks","docs","src", "man", "src", file),
        informat = "markdown",
        out_path = Pkg.dir("PhyloNetworks","docs","src", "man"),
        fig_path = "../assets/figures",
        fig_ext  = ".png",
        doctype  = "github")
end


text:
using Plots
using Reproduce
using JLD2
using FileIO
using Statistics

include("plot_utils.jl")

num_steps = 60000
log_interval = 100

function get_lines(data)
    means = []
    step = 10
    for start in 1:step:size(data)[2]
        print()
        chunk = data[:, start:min(start + step - 1, size(data)[2]),: ]
        # println(size(chunk))
        chunk_mean = mean(chunk, dims=2)
        push!(means, chunk_mean)
    end

    smoothed = cat(means..., dims=2)[:,:, :]
    sum_lines = sum(smoothed, dims=1)
    num_runs = size(sum_lines)[3]
    mean_line = dropdims(mean(sum_lines, dims=3), dims=(1, 3))
    std_err_line = dropdims(std(sum_lines, dims=3), dims=(1, 3)) / sqrt(num_runs)

    xs = 1:log_interval*step:num_steps + 1
    println(size(xs))
    println(size(mean_line))
    println(size(std_err_line))
    return xs, mean_line, std_err_line
end

sweep_params = ["eta", "alpha_init"]

p = plot()

ic = ItemCollection("M:/globus/OneDTMaze_Emphatic_Control_Real")
# ic = ItemCollection("experiment_data/EmphaticTest")
# ic = search(ic, Dict("demon_learner" => "SR"))

demon_learner = "SR"
ic = search(ic, Dict("demon_learner" => demon_learner))
# EmphaticTest
print(diff(ic))

# sdfsdf

metric = :oned_tmaze_start_error
function plot_line(ic, demon_update_algo)
    println(diff(ic))
    filtered_ic = search(ic, Dict("demon_update" => demon_update_algo))
    println(diff(filtered_ic))

    best_ic = get_best(filtered_ic, sweep_params, metric)
    print_params(best_ic, sweep_params, [])

    data = load_results(best_ic, metric)
    println(size(data))

    xs, mean_line, std_err_line = get_lines(data)
    plot!(p, xs, mean_line, ribbons=std_err_line, label=demon_update_algo, ylabel="RMSE", xlabel="time step", ylim=(0,15), title=demon_learner)
end

plot_line(ic, "ESARSA")
plot_line(ic, "EmphESARSA")
plot_line(ic, "TB")
plot_line(ic, "ETB")
plot_line(ic, "PriorESARSA")
plot_line(ic, "PriorTB")

savefig("plotting/chunlok/generated_plots/test_plot_reproduce.svg")

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
# unary # 2^2 - 2 = 2 definitions
@inline fmap(f::F, x::Tuple{X}) where {F,X} = (f(first(x)),)
@inline fmap(f::F, x::NTuple) where {F} = (f(first(x)), fmap(f, Base.tail(x))...)

# binary # 2^3 - 2 = 6 definitions
@inline fmap(f::F, x::Tuple{X}, y::Tuple{Y}) where {F,X,Y} = (f(first(x), first(y)),)
@inline fmap(f::F, x::Tuple{X}, y) where {F,X} = (f(first(x), y),)
@inline fmap(f::F, x, y::Tuple{Y}) where {F,Y} = (f(x, first(y)),)
@inline fmap(f::F, x::Tuple{Vararg{Any,N}}, y::Tuple{Vararg{Any,N}}) where {F,N} = (f(first(x), first(y)), fmap(f, Base.tail(x), Base.tail(y))...)
@inline fmap(f::F, x::Tuple, y) where {F} = (f(first(x), y), fmap(f, Base.tail(x), y)...)
@inline fmap(f::F, x, y::Tuple) where {F} = (f(x, first(y)), fmap(f, x, Base.tail(y))...)

fmap(f::F, x::Tuple{X}, y::Tuple) where {F,X} = throw("Dimension mismatch.")
fmap(f::F, x::Tuple, y::Tuple{Y}) where {F,Y} = throw("Dimension mismatch.")
fmap(f::F, x::Tuple, y::Tuple) where {F} = throw("Dimension mismatch.")

@generated function fmap(f::F, x::Vararg{Any,N}) where {F,N}
    q = Expr(:block, Expr(:meta, :inline))
    t = Expr(:tuple)
    U = 1
    call = Expr(:call, :f)
    syms = Vector{Symbol}(undef, N)
    istup = Vector{Bool}(undef, N)
    for n ∈ 1:N
        syms[n] = xₙ = Symbol(:x_, n)
        push!(q.args, Expr(:(=), xₙ, Expr(:ref, :x, n)))
        istup[n] = ist = (x[n] <: Tuple)
        if ist
            U = length(x[n].parameters)
            push!(call.args, Expr(:ref, xₙ, 1))
        else
            push!(call.args, xₙ)
        end
    end
    push!(t.args, call)
    for u ∈ 2:U
        call = Expr(:call, :f)
        for n ∈ 1:N
            xₙ = syms[n]
            if istup[n]
                push!(call.args, Expr(:ref, xₙ, u))
            else
                push!(call.args, xₙ)
            end
        end
        push!(t.args, call)
    end
    push!(q.args, t); q
end

for op ∈ [:vsub, :vabs, :vfloor, :vceil, :vtrunc, :vround, :vsqrt, :vnot, :vleading_zeros, :vtrailing_zeros, :vsub_fast]
    @eval @inline $op(v1::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap($op, v1.data))
end
# only for `Float32` and `Float64`
@inline inv_approx(v::VecUnroll{N,W,T}) where {N,W,T<:Union{Float32,Float64}} = VecUnroll(fmap(inv_approx, v.data))

@inline vreinterpret(::Type{T}, v::VecUnroll) where {T<:Number} = VecUnroll(fmap(vreinterpret, T, v.data))
for op ∈ [:vadd,:vsub,:vmul,:vand,:vor,:vxor,:vlt,:vle,:vgt,:vge,:veq,:vne,:vadd_fast,:vsub_fast,:vmul_fast]
    @eval begin
        @inline $op(v1::VecUnroll, v2::VecUnroll) = VecUnroll(fmap($op, v1.data, v2.data))
        @inline $op(v1::VecUnroll{N,W,T,V}, ::StaticInt{M}) where {N,W,T,V,M} = VecUnroll(fmap($op, v1.data, vbroadcast(Val{W}(), T(M))))
        @inline $op(::StaticInt{M}, v1::VecUnroll{N,W,T,V}) where {N,W,T,V,M} = VecUnroll(fmap($op, vbroadcast(Val{W}(), T(M)), v1.data))

        @inline $op(v1::VecUnroll{N,W,T,V}, ::StaticInt{1}) where {N,W,T,V} = VecUnroll(fmap($op, v1.data, One()))
        @inline $op(::StaticInt{1}, v1::VecUnroll{N,W,T,V}) where {N,W,T,V} = VecUnroll(fmap($op, One(), v1.data))

        @inline $op(v1::VecUnroll{N,W,T,V}, ::StaticInt{0}) where {N,W,T,V} = VecUnroll(fmap($op, v1.data, Zero()))
        @inline $op(::StaticInt{0}, v1::VecUnroll{N,W,T,V}) where {N,W,T,V} = VecUnroll(fmap($op, Zero(), v1.data))
    end
end
for op ∈ [:vmax,:vmax_fast,:vmin,:vmin_fast,:vcopysign]
    @eval begin
        @inline $op(v1::VecUnroll, v2::VecUnroll) = VecUnroll(fmap($op, v1.data, v2.data))
    end
end
for op ∈ [:vgt,:vge,:vlt,:vle,:veq,:vne,:vmax,:vmax_fast,:vmin,:vmin_fast]
    @eval begin
        @inline $op(v::VecUnroll{N,W,T,V}, s::Union{NativeTypes,AbstractSIMDVector}) where {N,W,T,V} = VecUnroll(fmap($op, v.data, vbroadcast(Val{W}(), s)))
        @inline $op(s::Union{NativeTypes,AbstractSIMDVector}, v::VecUnroll{N,W,T,V}) where {N,W,T,V} = VecUnroll(fmap($op, vbroadcast(Val{W}(), s), v.data))
    end
end

for op ∈ [:vrem, :vshl, :vashr, :vlshr, :vdiv, :vfdiv, :vrem_fast, :vfdiv_fast]
    @eval begin
        @inline $op(v1::VecUnroll, v2::VecUnroll) = VecUnroll(fmap($op, v1.data, v2.data))
    end
end
for op ∈ [:vrem, :vand, :vor, :vxor, :vshl, :vashr, :vlshr, :vlt, :vle,:vgt,:vge,:veq,:vne]
    @eval begin
        @inline $op(vu::VecUnroll, i::MM) = $op(vu, Vec(i))
        @inline $op(i::MM, vu::VecUnroll) = $op(Vec(i), vu)
    end
end
for op ∈ [:vshl, :vashr, :vlshr]
    @eval @inline $op(m::Mask, vu::VecUnroll) = $op(Vec(m), vu)
end
for op ∈ [:rotate_left,:rotate_right,:funnel_shift_left,:funnel_shift_right]
    @eval begin
        @inline $op(v1::VecUnroll{N,W,T,V}, v2::R) where {N,W,T,V,R<:IntegerTypes} = VecUnroll(fmap($op, v1.data, promote_type(V, R)(v2)))
        @inline $op(v1::R, v2::VecUnroll{N,W,T,V}) where {N,W,T,V,R<:IntegerTypes} = VecUnroll(fmap($op, promote_type(V, R)(v1), v2.data))
        @inline $op(v1::VecUnroll, v2::VecUnroll) = VecUnroll(fmap($op, v1.data, v2.data))
    end
end

for op ∈ [:vfma, :vmuladd, :vfma_fast, :vmuladd_fast, :vfnmadd, :vfmsub, :vfnmsub, :vfnmadd_fast, :vfmsub_fast, :vfnmsub_fast,
          :vfmadd231, :vfnmadd231, :vfmsub231, :vfnmsub231, :ifmahi, :ifmalo]
    @eval begin
        @inline $op(v1::VecUnroll{N,W}, v2::VecUnroll{N,W}, v3::VecUnroll{N,W}) where {N,W} = VecUnroll(fmap($op, v1.data, v2.data, v3.data))
    end
end
@inline vifelse(v1::VecUnroll{N,W,<:Boolean}, v2::T, v3::T) where {N,W,T<:NativeTypes} = VecUnroll(fmap(vifelse, v1.data, Vec{W,T}(v2), Vec{W,T}(v3)))
@inline vifelse(v1::VecUnroll{N,W,<:Boolean}, v2::T, v3::T) where {N,W,T<:Real} = VecUnroll(fmap(vifelse, v1.data, v2, v3))
@inline vifelse(v1::Vec{W,Bool}, v2::VecUnroll{N,W,T}, v3::Union{NativeTypes,AbstractSIMDVector,StaticInt}) where {N,W,T} = VecUnroll(fmap(vifelse, Vec{W,T}(v1), v2.data, Vec{W,T}(v3)))
@inline vifelse(v1::Vec{W,Bool}, v2::Union{NativeTypes,AbstractSIMDVector,StaticInt}, v3::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap(vifelse, Vec{W,T}(v1), Vec{W,T}(v2), v3.data))
@inline vifelse(v1::VecUnroll{N,W,<:Boolean}, v2::VecUnroll{N,W,T}, v3::Union{NativeTypes,AbstractSIMDVector,StaticInt}) where {N,W,T} = VecUnroll(fmap(vifelse, v1.data, v2.data, Vec{W,T}(v3)))
@inline vifelse(v1::VecUnroll{N,W,<:Boolean}, v2::Union{NativeTypes,AbstractSIMDVector,StaticInt}, v3::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap(vifelse, v1.data, Vec{W,T}(v2), v3.data))
@inline vifelse(v1::Vec{W,Bool}, v2::VecUnroll{N,W,T}, v3::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap(vifelse, Vec{W,T}(v1), v2.data, v3.data))
@inline vifelse(v1::VecUnroll{N,W,<:Boolean}, v2::VecUnroll{N,W,T}, v3::VecUnroll{N,W,T}) where {N,W,T} = VecUnroll(fmap(vifelse, v1.data, v2.data, v3.data))
@inline function vifelse(v1::VecUnroll{N,W,<:Boolean}, v2::VecUnroll{N,W}, v3::VecUnroll{N,W}) where {N,W}
    v4, v5 = promote(v2, v3)
    VecUnroll(fmap(vifelse, v1.data, v4.data, v5.data))
end

@inline veq(v::VecUnroll{N,W,T}, x::AbstractIrrational) where {N,W,T} = v == vbroadcast(Val{W}(), T(x))
@inline veq(x::AbstractIrrational, v::VecUnroll{N,W,T}) where {N,W,T} = vbroadcast(Val{W}(), T(x)) == v

@inline vunsafe_trunc(::Type{T}, v::VecUnroll) where {T<:Real} = VecUnroll(fmap(vunsafe_trunc, T, v.data))
@inline vrem(v::VecUnroll, ::Type{T}) where {T<:Real} = VecUnroll(fmap(vrem, v.data, T))
@inline vrem(v::VecUnroll{N,W}, ::Type{VecUnroll{N,W,T,V}}) where {N,W,T,V} = VecUnroll(fmap(vrem, v.data, V))

@inline (::Type{VecUnroll{N,W,T,V}})(vu::VecUnroll{N,W,T,V}) where {N,W,T,V<:AbstractSIMDVector{W,T}} = vu
@inline (::Type{VecUnroll{N,W,T,VT}})(vu::VecUnroll{N,W,S,VS})  where {N,W,T,VT<:AbstractSIMDVector{W,T},S,VS<:AbstractSIMDVector{W,S}} = VecUnroll(fmap(convert, Vec{W,T}, vu.data))

function collapse_expr(N, op)
    N += 1
    t = Expr(:tuple); s = Vector{Symbol}(undef, N)
    for n ∈ 1:N
        s_n = s[n] = Symbol(:v_, n)
        push!(t.args, s_n)
    end
    q = quote
        $(Expr(:meta,:inline))
        $t = data(vu)
    end
    while N > 1
        for n ∈ 1:N >>> 1
            push!(q.args, Expr(:(=), s[n], Expr(:call, op, s[n], s[n + (N >>> 1)])))
        end
        isodd(N) && push!(q.args, Expr(:(=), s[1], Expr(:call, op, s[1], s[N])))
        N >>>= 1
    end
    q
end
@generated collapse_add(vu::VecUnroll{N}) where {N} = collapse_expr(N, :vadd)
@generated collapse_mul(vu::VecUnroll{N}) where {N} = collapse_expr(N, :vmul)
@generated collapse_max(vu::VecUnroll{N}) where {N} = collapse_expr(N, :max)
@generated collapse_min(vu::VecUnroll{N}) where {N} = collapse_expr(N, :min)
@generated collapse_and(vu::VecUnroll{N}) where {N} = collapse_expr(N, :&)
@generated collapse_or(vu::VecUnroll{N}) where {N} = collapse_expr(N, :|)
@inline vsum(vu::VecUnroll) = vsum(collapse_add(vu))
@inline vsum(s, vu::VecUnroll) = vsum(s, collapse_add(vu))
@inline vprod(vu::VecUnroll) = vprod(collapse_mul(vu))
@inline vprod(s, vu::VecUnroll) = vprod(s, collapse_mul(vu))
@inline vmaximum(vu::VecUnroll) = vmaximum(collapse_max(vu))
@inline vminimum(vu::VecUnroll) = vminimum(collapse_min(vu))
@inline vall(vu::VecUnroll) = vall(collapse_and(vu))
@inline vany(vu::VecUnroll) = vany(collapse_or(vu))


text:
# Autogenerated wrapper script for SLICOT_jll for x86_64-linux-gnu-libgfortran4
export libslicot

using libblastrampoline_jll
using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("SLICOT")
JLLWrappers.@declare_library_product(libslicot, "libslicot.so")
function __init__()
    JLLWrappers.@generate_init_header(libblastrampoline_jll, CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libslicot,
        "lib/libslicot.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using SparseArrays
using SparseUtils
using Serialization

function create_serialize_sparse_array(nr=30, nc=70, density=0.01)
    sparse_array = SparseArrays.sprandn(nr, nc, density)
    open("sparse_array.dat", "w") do io
        serialize(io, sparse_array)
    end
    return nothing
end

function deserialize_sparse_array()
    return open("sparse_array.dat", "r") do io
        res = deserialize(io)
    end
end