text:
using LinearAlgebra

useVerbose  = [true,true,false]

macro verbose(level::Integer, string::String)
    @assert(0 < level < 4)
    if level > 0; padding = ""  ; end
    if level > 1; padding = "--"; end
    if level > 2; padding = "--"; end
    string = padding * string
    if useVerbose[level]
        return :( println($string) )
    end
end

function get_matrix_and_copy(in)
    return in, deepcopy(in), matrix_double_t(in)
end

function skew(u)
    S = zeros(3,3)
    S[1,2] = -u[3];
    S[1,3] = +u[2];
    S[2,1] = +u[3];
    S[2,3] = -u[1];
    S[3,1] = -u[2];
    S[3,2] = +u[1];
    return S;
end

function SO3_hat(u)
    return skew(u);
end

function SO3_vee(gfrak)
    out = zeros(3,1);
    out[1] = gfrak[3,2];
    out[2] = gfrak[1,3];
    out[3] = gfrak[2,1];
    return out;
end

function SU2_hat(u)
    L1 = [ 0.0   im;  im  0.0];
    L2 = [ 0.0 -1.0; 1.0  0.0];
    L3 = [  im  0.0; 0.0  -im];
    gfrak = L1*u[1] + L2*u[2] + L3*u[3];
    return gfrak
end

function SU2_vee(gfrak)
    out = [ imag(gfrak[1,2] + gfrak[2,1])/2.0;
            real(gfrak[2,1] - gfrak[1,2])/2.0;
            imag(gfrak[1,1] - gfrak[2,2])/2.0];
    return out
end

function quat_2_SU2(qX::Array{Float64, 2})
    X = zeros(2,2) + zeros(2,2)*im;
    X[1,1] = +qX[1]+qX[4]*im
    X[1,2] = -qX[3]+qX[2]*im
    X[2,1] = +qX[3]+qX[2]*im
    X[2,2] = +qX[1]-qX[4]*im
    return X
end

function SU2_2_quat(X::Array{Complex{Float64}, 2})
    q = zeros(4,1);
    q[1,1] = real(X[1,1])
    q[2,1] = imag(X[2,1])
    q[3,1] = real(X[2,1])
    q[4,1] = imag(X[1,1])
    return q
end

function quat_2_SO3(q::Array{Float64, 2})
    a, b, c, d = q[1,1], q[2,1], q[3,1], q[4,1];
    R = zeros(3,3);
    R[1,1] = a^2+b^2-c^2-d^2;
    R[1,2] = 2*(b*c-a*d);
    R[1,3] = 2*(b*d+a*c);
    R[2,1] = 2*(b*c+a*d);
    R[2,2] = a^2-b^2+c^2-d^2;
    R[2,3] = 2*(c*d-a*b);
    R[3,1] = 2*(b*d-a*c);
    R[3,2] = 2*(c*d+a*b);
    R[3,3] = a^2-b^2-c^2+d^2;
    return R
end

function quaternion_left_product(q)
    qim = q[2:4]
    qre = q[1]
    QL = zeros(4,4)
    QL[1,2:4]   = -qim;
    QL[2:4,1]   = +qim;
    QL[2:4,2:4] = skew(qim);
    QL = QL + I*qre;
    return QL
end

function quaternion_right_product(q)
    qim = q[2:4]
    qre = q[1]
    QR = zeros(4,4);
    QR[1,2:4]   = -qim;
    QR[2:4,1]   = +qim;
    QR[2:4,2:4] = -skew(qim);
    QR = QR + I*qre;
    return QR
end

function attitude_dynamics(q, w, J, tau; delta = zeros(3,))
    we      = zeros(4,1);
    we[2:4] = w;
    dq      = quaternion_left_product(q) * we / 2.0;
    dw      = (J) \ (-skew(w)*(J*w) + tau + delta);
    return dq, dw
end

function ntuple_2_mat(tup, Nrow, Ncol)
    mat = zeros(Nrow, Ncol)
    for i = 1:Nrow
        for j = 1:Ncol
            mat[i,j] = tup[i][j];
        end
    end
    return mat;
end

function rand_PSD(N)
    J  = rand(N,N); J  = J + J'; J  = J + I*maximum(abs.(eigvals(J)));
    return J;
end

function maximal_feasible_kc_SU2(kX, kw, J)
    ev   = eigvals(J);
    maxJ = maximum(ev);
    minJ = minimum(ev);
    bounds = [4.0 * kw, 2.0 * sqrt(kw * minJ), 4.0 * kw * kX * minJ * minJ / (maxJ * kw * kw + minJ * minJ * kX)]
    return minimum(bounds);
end

function maximal_feasible_kc_SO3(kR, kw, J)
    ev   = eigvals(J);
    maxJ = maximum(ev);
    minJ = minimum(ev);
    bounds = [kw, sqrt(kR * minJ), 4.0 * kw * kR * minJ * minJ / (maxJ * kw * kw + 4.0 * minJ * minJ * kR)]
    return minimum(bounds);
end

function allFieldsEqual(
    C1::con_state_qw_fsf_t,
    C2::con_state_qw_fsf_t,
    checkStatic::Bool
)
    @test isapprox(C1.gain_kR, C2.gain_kR);
    @test isapprox(C1.gain_kc, C2.gain_kc);
    @test isapprox(C1.gain_kw, C2.gain_kw);
    @test isapprox(C1.param_a, C2.param_a);
    @test isapprox(C1.param_b, C2.param_b);
    @test isapprox(C1.param_c, C2.param_c);
    @test isapprox(C1.param_d, C2.param_d);
    @test isapprox(C1.thrust, C2.thrust);
    for i = 1:3
        for j = 1:3
            @test isapprox(C1.inertia[i][j], C2.inertia[i][j]);
        end
    end
    if checkStatic
        for i = 1:3
            @test isapprox(C1.torques[i],  C2.torques[i] );
        end
        @test isapprox(C1.dist_Psi,       C2.dist_Psi );
        @test isapprox(C1.dist_Gamma,     C2.dist_Gamma );
        @test isapprox(C1.dist_lyapunov,  C2.dist_lyapunov );
    end
end

function allFieldsEqual(
    R1::ref_state_qw_t,
    R2::ref_state_qw_t,
)
    for i = 1:4
        @test  isapprox(R1.quaternion[i], R2.quaternion[i]);
    end
    for i = 1:3
        @test  isapprox(R1.attrates[i], R2.attrates[i]);
        @test  isapprox(R1.attaccel[i], R2.attaccel[i] );
    end
end

function allFieldsEqual(
    S1::dyn_state_qw_t,
    S2::dyn_state_qw_t,
)
    for i = 1:4
        @test  isapprox(S1.quaternion[i], S2.quaternion[i]);
    end
    for i = 1:3
        @test  isapprox(S1.attrates[i], S2.attrates[i]);
    end
end

function get_M1M2W_SO3(C, phi)
    # Extract parameters
    kR, kc, kw  = C.gain_kR, C.gain_kc, C.gain_kw;
    J           = ntuple_2_mat(C.inertia, 3, 3);
    minJ, maxJ  = minimum(eigvals(J)), maximum(eigvals(J));

    # Evaluate matrices
    M1 = [[+kR, -kc] [-kc, +minJ]] / 2.0;
    M2 = [[2*kR/(2 - phi), kc] [kc, maxJ]] / 2.0;
    W  = [[+kR*kc/maxJ, -kc*kw/(2 * minJ)] [-kc*kw/(2 * minJ), kw - kc]];
    return M1, M2, W
end

function get_M1M2W_SU2(C, phi)
    # Extract parameters
    kX, kc, kw  = C.gain_kR, C.gain_kc, C.gain_kw;
    J           = ntuple_2_mat(C.inertia, 3, 3);
    minJ, maxJ  = minimum(eigvals(J)), maximum(eigvals(J));

    # Evaluate matrices
    M1 = [[+4.0*kX, -kc] [-kc, +minJ]] / 2.0;
    M2 = [[+8.0*kX/(2.0 - phi), kc] [kc, maxJ]]/2.0;
    W  = [[+kX*kc/maxJ, -kc*kw/(2 * minJ)] [-kc*kw/(2 * minJ), kw - kc/4.0]];
    return M1, M2, W
end

function numerical_differentiation(in, time)
    NT  = size(in, 2);
    out = zeros(size(in));
    for t = 2:NT;
        dt = time[t] - time[t - 1];
        out[:,t] = (in[:,t] - in[:,t-1])/dt;
    end
    return out
end

function maximal_feasible_eps_SO3(
    C::con_state_qw_fsf_t,
    R0::Array{T,2},
    Rr0::Array{T,2},
    w0::Array{T,2},
    wr0::Array{T,2}
) where {T<:Number}

    # Extract parameters
    kR, kc, kw  = C.gain_kR, C.gain_kc, C.gain_kw;
    J           = ntuple_2_mat(C.inertia, 3, 3);

    # Compute maximum permissible phi
    Psi         = tr(I-Rr0'*R0)/2.0;
    eR0         = SO3_vee(Rr0'*R0 - R0'*Rr0) / 2.0;
    ew0         = w0 -  R0'*Rr0*wr0;
    attpart     = kR * Psi;
    crosspart   = kc * ew0'*eR0;
    ratepart    = ew0'*J*ew0/2.0;
    Lyap0       = attpart + crosspart[1,1] + ratepart[1,1];
    phi         = Lyap0/kR;

    M1, M2, W   = get_M1M2W_SO3(C, phi)

    M1max, M1min = maximum(eigvals(M1)), minimum(eigvals(M1));
    M2max, M2min = maximum(eigvals(M2)), minimum(eigvals(M2));
    Wmax,  Wmin  = maximum(eigvals(W)),  minimum(eigvals(W));

    mindecayrate = Wmin/M2max;
    maxeps       = ((M1min*Wmin) / M2max)*phi*(2.0 - phi);
    uniboundgain = (M2max / (M1min*Wmin));

    return maxeps, mindecayrate, uniboundgain, phi, M1, M2, W
end


text:
using LinearAlgebra
using Random

Random.seed!(843)

include("feast.jl")
using ..feastLinear
using IterativeSolvers

n=100

#generate B-matrix eigenvectors and eigenvalues
#vectors:
R=rand(n,n)
(Xb,R)=qr(R)
#values:
Lb=rand(n)*0.5 .+ 0.5 #random positive numbers

#generate A-matrix eigenvectors:
R=rand(n,n)
(X1,R)=qr(R)

#generate A-matrix eigenvalues
L1=zeros(n)
for i in 1:n
    L1[i]=i
end

#generate eigenvalue problem matrices
A=Xb*diagm(0 => sqrt.(Lb))*X1*diagm(0 => L1)*X1'*diagm(0 => sqrt.(Lb))*Xb'
B=Xb*diagm(0 => Lb)*Xb'

#FEAST parameters
emid=1.0+0.0*im #contour center
ra=0.5 #contour radius 1
rb=0.5 #contour radius 2
nc=6 #number of contour points
m0=2 #subspace dimension
eps=1e-5 #residual convergence tolerance
maxit=100 #maximum FEAST iterations
x0=rand(n,m0) #eigenvector initial guess

println("Standard FEAST")
(lf,xf)=feast_linear(A,B,x0,nc,emid,ra,rb,eps,maxit)
p=sortperm(real(lf))
println(lf[p])

println("Iterative FEAST BiCGSTAB")
alpha=0.01 #relative accuracy of linear system solves
isMaxit=220  #maximum number of linear system iterations per FEAST iteration
(lf,xf,data)=ifeast_linear(A,B,x0,alpha,isMaxit,nc,emid,ra,rb,eps,maxit;log=true)
p=sortperm(real(lf))
println(lf[p])

ext: jl

max_stars_repo_name: rleegates/feastjl

max_issues_repo_name: rleegates/feastjl

max_forks_repo_name: rleegates/feastjl


text:
using Pkg
packages = ["ArgParse"]
Pkg.add(packages)
using ArgParse


text:
using DynamicGrids, Test, Dates, Unitful

# life glider sims

init =  [0 0 0 0 0 0;
         0 0 0 0 0 0;
         0 0 0 0 0 0;
         0 0 0 1 1 1;
         0 0 0 0 0 1;
         0 0 0 0 1 0]

test2 = [0 0 0 0 0 0;
         0 0 0 0 0 0;
         0 0 0 0 1 0;
         0 0 0 0 1 1;
         0 0 0 1 0 1;
         0 0 0 0 0 0]

test3 = [0 0 0 0 0 0;
         0 0 0 0 0 0;
         0 0 0 0 1 1;
         0 0 0 1 0 1;
         0 0 0 0 0 1;
         0 0 0 0 0 0]

test5 = [0 0 0 0 0 0;
         0 0 0 0 0 0;
         1 0 0 0 1 1;
         1 0 0 0 0 0;
         0 0 0 0 0 1;
         0 0 0 0 0 0]

ruleset = Ruleset(Life(); init=init, overflow=WrapOverflow(), timestep=Day(2))
output = ArrayOutput(init, 5)
sim!(output, ruleset; tspan=(Date(2001, 1, 1), Date(2001, 1, 10)))

@testset "stored results match glider behaviour" begin
    @test output[2] == test2
    @test output[3] == test3
    @test output[5] == test5
end

@testset "converted results match glider behaviour" begin
    output2 = ArrayOutput(output)
    @test output2[2] == test2
    @test output2[3] == test3
    @test output2[5] == test5
end

@testset "REPLOutput block works, in Unitful.jl seconds" begin
    ruleset = Ruleset(Life(); init=init, overflow=WrapOverflow(), timestep=5u"s")
    output = REPLOutput(init; style=Block(), fps=100, store=true)
    sim!(output, ruleset; tspan=(0u"s", 10u"s"))
    resume!(output, ruleset; tstop=25u"s")
    @test output[2] == test2
    @test output[3] == test3
    @test output[5] == test5
    # replay(output, ruleset)
end

@testset "REPLOutput braile works, in Months" begin
    ruleset = Ruleset(Life(); init=init, overflow=WrapOverflow(), timestep=Month(1))
    output = REPLOutput(init; style=Braile(), fps=100, store=true)
    sim!(output, ruleset; tspan=(Date(2010, 4), Date(2010, 7)))
    resume!(output, ruleset; tstop=Date(2010, 9))
    @test output[2] == test2
    @test output[3] == test3
    @test output[5] == test5
    # replay(output, ruleset)
end


text:
# ---
# title: CropRatio
# cover: cropratio.gif
# ---

# Crop centered window to fit given aspect ratio

using Augmentor
using ImageShow, ImageCore

img_in = testpattern(RGB, ratio=0.5)
img_out = augment(img_in, CropRatio()) # crop out a square window

mosaicview(img_in, img_out; nrow=1)

# `RCropRatio` is a random version that randomly choose a crop center -- not necessarily the center
# of the input image.

augment(img_in, RCropRatio())
#md nothing #hide

# ![](cropratio.gif)

# ## Reference

#md # ```@docs
#md # CropRatio
#md # RCropRatio
#md # ```

## save covers #src
using ImageMagick #src
using FileIO #src
include(joinpath("..", "assets", "utilities.jl")) #src
cover = make_gif(img_in, RCropRatio(), 5) #src
ImageMagick.save("cropratio.gif", cover; fps=1) #src


text:
"""
Finds the earliest bus that occur after `time`, then return
its ID given by (wait time) * (bus number).
"""
function find_earliest(time::Int, buses::Array{String})
  lowest_wait = nothing
  earliest_bus = nothing

  # Remove x's, parse into Int
  buses = filter(x -> x != "x", buses)
  buses = map(x -> parse(Int, x), buses)

  for bus in buses
    div = time ÷ bus
    wait = ((div + 1) * bus) - time
    bid = wait * bus

    if lowest_wait == nothing || wait < lowest_wait
      lowest_wait = wait
      earliest_bus = bus
    end
  end

  return lowest_wait * earliest_bus
end

"""
Finds the time where buses will begin arriving according to their index in the input.
"""
function find_increasing_arrival_time(buses::Array{String})
  mods = Dict{Int, Int}()

  for i in 1:length(buses)
    bus = tryparse(Int, buses[i])
    if (!isnothing(bus))
      push!(mods, bus => mod(-(i - 1), bus))
    end
  end

  bus_product = nothing
  value = nothing

  for (bus, mod) in mods
    if isnothing(bus_product) || isnothing(value)
      # Seed initial value
      bus_product = bus
      value = mod
      continue
    end

    while value % bus != mod
      value += bus_product
    end

    bus_product *= bus
  end

  return value
end

# Entry point, starts reading the input file and solves the puzzles.
lines = readlines("var/aoc13_input.txt")
time = parse(Int, lines[1])
buses = split(lines[2], ",")

# Convert Array{SubString{String}} to Array{String}
buses = map(String, buses)

bid = find_earliest(time, buses)
println("ID for earliest bus: $bid")

time = find_increasing_arrival_time(buses)
println("Time where schedules align: $time")


text:
## Panel helpers
#============================================#

# Source singularities

constant_source_stream(σ, panel :: AbstractPanel2D, xi, yi) = panel_scalar(constant_source_stream, σ, panel, xi, yi)

linear_source_stream_1(σ1, panel :: AbstractPanel2D, xi, yi) =  panel_scalar(linear_source_stream_minus, σ1, panel, xi, yi) - panel_scalar(linear_source_stream_plus, σ1, panel, xi, yi)

linear_source_stream_2(σ2, panel :: AbstractPanel2D, xi, yi) = panel_scalar(linear_source_stream_plus, σ2, panel, xi, yi)

linear_source_stream(σ1, σ2, panel :: AbstractPanel2D, xi, yi) =
linear_source_stream_1(σ1, panel, xi, yi), linear_source_stream_2(σ2, panel, xi, yi)

source_stream(σ1, σ2, panel :: AbstractPanel2D, xi, yi) = linear_source_stream_1(σ1, panel, xi, yi) + linear_source_stream_2(σ2, panel, xi, yi)

source_stream(σ, panel :: AbstractPanel2D, xi, yi) = panel_scalar(constant_source_stream, σ, panel, xi, yi)

linear_source_velocity(σ1, σ2, panel :: AbstractPanel2D, xi, yi) = panel_velocity(linear_source_velocity_a, σ1, panel, xi, yi), panel_velocity(linear_source_velocity_b, σ2, panel, xi, yi)

source_velocity(σ, panel :: AbstractPanel2D, xi, yi) = panel_velocity(constant_source_velocity, σ, panel, xi, yi)

source_velocity(σ1, σ2, panel :: AbstractPanel2D, xi, yi) = panel_velocity(linear_source_velocity_a, linear_source_velocity_b, σ1, σ2, panel, xi, yi)

source_velocity(σ1, σ2, panel :: AbstractPanel2D, panel_i :: AbstractPanel2D) = panel_velocity(linear_source_velocity_a, linear_source_velocity_b, σ1, σ2, panel, panel_i)

source_velocity(σ, panel :: AbstractPanel2D, panel_i :: AbstractPanel2D) = panel_velocity(constant_source_velocity, σ, panel, panel_i)

# Vortex singularities
linear_vortex_stream_1(γ1, panel :: AbstractPanel2D, xi, yi) =  panel_scalar(linear_vortex_stream_minus, γ1, panel, xi, yi) - panel_scalar(linear_vortex_stream_plus, γ1, panel, xi, yi)

linear_vortex_stream_2(γ2, panel :: AbstractPanel2D, xi, yi) = panel_scalar(linear_vortex_stream_plus, γ2, panel, xi, yi)

linear_vortex_stream(σ1, σ2, panel :: AbstractPanel2D, xi, yi) =
linear_vortex_stream_1(σ1, panel, xi, yi), linear_vortex_stream_2(σ2, panel, xi, yi)

vortex_stream(γ1, γ2, panel :: AbstractPanel2D, xi, yi) = linear_vortex_stream_1(γ1, panel, xi, yi) + linear_vortex_stream_2(γ2, panel, xi, yi)

vortex_stream(γ, panel :: AbstractPanel2D, xi, yi) = vortex_stream(γ, γ, panel :: AbstractPanel2D, xi, yi)

vortex_velocity(γ1, γ2, panel :: AbstractPanel2D, x, y) = panel_velocity(linear_vortex_velocity_a, linear_vortex_velocity_b, γ1, γ2, panel, x, y)

vortex_velocity(γ, panel :: AbstractPanel2D, x, y) = vortex_velocity(γ, γ, panel :: AbstractPanel2D, x, y)

vortex_velocity(γ1, γ2, panel :: AbstractPanel2D, panel_i :: AbstractPanel2D) = panel_velocity(linear_vortex_velocity_a, linear_vortex_velocity_b, γ1, γ2, panel, panel_i)

linear_vortex_velocity(γ1, γ2, panel :: AbstractPanel2D, xi, yi) = panel_velocity(linear_vortex_velocity_a, γ1, panel, xi, yi), panel_velocity(linear_vortex_velocity_b, γ2, panel, xi, yi)

# Total velocity

total_velocity(σ_1j :: Real, σ_2j :: Real, γ_1j :: Real, γ_2j :: Real, panel_j :: AbstractPanel2D, x, y) = source_velocity(σ_1j, σ_2j, panel_j, x, y) + vortex_velocity(γ_1j, γ_2j, panel_j, x, y)

total_velocity(σ_1j, σ_2j, γ_1j, γ_2j, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D) = source_velocity(σ_1j, σ_2j, panel_j, panel_i) .+ vortex_velocity(γ_1j, γ_2j, panel_j, panel_i)

# Influence coefficient computations
neumann_influence_coefficient(velocity_func :: F, angle_func, panel_j :: AbstractPanel2D, panel_i :: AbstractPanel2D) where F <: Function = dot(panel_velocity(velocity_func, 1., panel_j, panel_i), angle_func(panel_i))

# dirichlet_influence_coefficient()

## Matrix construction
#============================================#

function two_point_neumann_matrix(vel_a :: F, vel_b :: G, panels_1, panels_2) where F where G
    N     = min(length(panels_1), length(panels_2))
    inf_a = neumann_influence_matrix(vel_a, normal_vector, panels_1, panels_2)
    inf_b = neumann_influence_matrix(vel_b, normal_vector, panels_1, panels_2)
    [ inf_a zeros(N) ] + [ zeros(N) inf_b ]
end

linear_source_neumann_matrix(panels_1, panels_2) = two_point_neumann_matrix(linear_source_velocity_a, linear_source_velocity_b, panels_1, panels_2)
linear_vortex_neumann_matrix(panels_1, panels_2) = two_point_neumann_matrix(linear_vortex_velocity_a, linear_vortex_velocity_b, panels_1, panels_2)

kutta_condition(panels_1, panels_2) = [ 1 zeros(length(panels_1) - 1)' 1 zeros(length(panels_2) - length(panels_1))' ]
source_influence_matrix(panels_1, panels_2) = linear_source_matrix(panels_1, panels_2)
vortex_influence_matrix(panels_1, panels_2) = linear_vortex_matrix(panels_1, panels_2)

# Diagonal cases
kutta_condition(panels) = [ 1 zeros(length(panels) - 1)' 1 ]
source_influence_matrix(panels) = [ linear_source_matrix(panels, panels); kutta_condition(panels) ]
vortex_influence_matrix(panels) = [ linear_vortex_matrix(panels, panels); kutta_condition(panels) ]

constant_source_influence_coefficient(panel_j, panel_i) = ifelse(panel_i === panel_j, 0.5, influence_coefficient(constant_source_velocity, normal_vector, panel_j, panel_i))

constant_source_matrix(panels) = [ constant_source_influence_coefficient(panel_j, panel_i) for (panel_j, panel_i) in product(panels, panels)]

neumann_influence_matrix(func, angle_func, panel_is, panel_js) = [ influence_coefficient(func, angle_func, panel_j, panel_i) for (panel_i, panel_j) in product(panel_is, panel_js) ]

function two_point_dirichlet_matrix(src_a :: F, src_b :: G, panels_1, panels_2) where {F,G}
    N     = min(length(panels_1), length(panels_2))
    inf_a = dirichlet_influence_matrix(src_a, normal_vector, panels_1, panels_2)
    inf_b = dirichlet_influence_matrix(src_b, normal_vector, panels_1, panels_2)
    [ inf_a zeros(N) ] + [ zeros(N) inf_b ]
end

# Boundary conditions
#============================================#

constant_source_boundary_condition(panels, u) = -map(pan -> dot(u, normal_vector(pan)), panels)

neumann_boundary_condition(panels, u) = [ constant_source_boundary_condition(panels, u); 0 ]

## Block matrix setup
#===============================================#

# linear_vortex_matrix()

function influence_matrix(panels, pts)
    N   = length(pts[:,1])

    # Linear vorticity-streamfunction matrices
    A1 = [ linear_vortex_stream_1(1., pj, xi[1], xi[2]) for xi in pts, pj in panels ]
    A2 = [ linear_vortex_stream_2(1., pj, xi[1], xi[2]) for xi in pts, pj in panels ]

    AIC_pts = [ A1 zeros(N) ] + [ zeros(N) A2 ]

    ## Trailing edge
    te_panel, h_TE, tcp, tdp, _ = trailing_edge_info(panels)

    # Constant-strength vortex and source terms
    A_TE = [
            linear_vortex_stream_1(1., te_panel, xi[1], xi[2]) + linear_vortex_stream_2(1., te_panel, xi[1], xi[2])
            for xi in pts
            ]
    B_TE = [
            constant_source_stream(1., te_panel, xi[1], xi[2]) for xi in pts ]

    AB_TE = (A_TE * tdp - B_TE * tcp) / 2

    AIC_TE = [ AB_TE zeros(N, N - 2) -AB_TE ]

    AIC_surf = AIC_pts + AIC_TE

    # Check trailing edge thickness and augment last equation
    if abs(h_TE) < 1e-10
        AIC[end-1,:] .= 0
        AIC[end-1,[1,2,3,end-3,end-2,end-1]] = [1,-2,1,-1,2,-1]
    end

    # Kutta condition
    kutta = [ 1. zeros(1, N - 2) 1. 0. ]

    # Assembly
    AIC  = [ AIC_surf  fill(-1., N);
               kutta            ]

    AIC
end

boundary_condition(pts) = @views [ -getindex.(pts, 2) getindex.(pts, 1);
                                    zeros(1, 2) ]

function solve_linear(panels)
    # Panel points
    pts  = panel_points(panels)
    AIC  = influence_matrix(panels, pts)
    boco = boundary_condition(pts)
    ψs   = AIC \ boco

    ψs, AIC, boco
end

## Preallocated setup
#===============================================#

function solve_linear!(A, b, panels :: AbstractVector{<: Panel2D{T}}) where T <: Real
    assemble_system!(A, b, panels)
    ψ = A \ b
end

function assemble_system!(AIC, boco, panels)
    # Panel points
    pts = panel_points(panels)

    # Trailing edge info
    te_panel, h_TE, tcp, tdp, _ = trailing_edge_info(panels)

    # Construct everything in one loop
    for i in eachindex(pts)
        x, y = pts[i]
        for j in eachindex(panels)
            AIC[i,j]   += linear_vortex_stream_1(1., panels[j], x, y)
            AIC[i,j+1] += linear_vortex_stream_2(1., panels[j], x, y)
        end

        # Trailing edge influence coefficients

        # Constant-strength vortex panel
        a_TE = linear_vortex_stream_1(1., te_panel, x, y) + linear_vortex_stream_2(1., te_panel, x, y)
        AIC[i,1]     += a_TE * ( tdp / 2)
        AIC[i,end-1] += a_TE * (-tdp / 2)

        # Constant-strength source panel
        b_TE = constant_source_stream(1., te_panel, x, y)
        AIC[i,1]     += b_TE * (-tcp / 2)
        AIC[i,end-1] += b_TE * ( tcp / 2)

        # Constant streamfunction value on the surface
        AIC[i,end] = -1

        # Check trailing edge thickness and augment last equation accordingly
        if abs(h_TE) < 1e-10
            AIC[end-1,:] .= 0
            AIC[end-1,[1,2,3,end-3,end-2,end-1]] = [1,-2,1,-1,2,-1]
        end

        # RHS
        boco[i,1] = -y
        boco[i,2] =  x
    end

    # Kutta condition
    AIC[end,1]     = 1
    AIC[end,end-1] = 1

    nothing
end


text:
"""
    gradientflops(layer, input)::Float64

Calculates the number of non-embedding Floating Point Operations `neFLOPs` for most layers.
`layer` should be any of the [Flux](https://github.com/FluxML/Flux.jl) model layers except
`GlobalMaxPool` and `GlobalMeanPool`.

# Example
```julia
layer = Conv((2, 2), 1=>16, relu)
input = (4,4)
layerflops(Conv((2, 2), 1=>16, relu),(4,4))
```
"""
function gradientflops(layer::Dense, input::Tuple)
    x = 0; y = 0; Fm = 1
    if length(input) == 3
        x,y,Fn = input
        xo, yo = outdims(layer,(x,y))
    elseif length(input) == 2
        x,y,_ = input
        xo, yo = outdims(layer,(x,y))
    else
        error("Not a valid Input size, expected (::Int,::Int) or (::Int,::Int,::Int)")
    end
    bi = length(layer.b)
    noofoper = gradientoperations(layer)
    return convert(Float64, ((xo*yo+bi)*noofoper)*Fm), (xo,yo)
end

function gradientflops(layer::Maxout,input::Tuple)
    x = 0; y = 0; Fm = 1
    if length(input) == 3
        x,y,Fm = input
        xo, yo = outdims(layer,(x,y))
    elseif length(input) == 2
        x,y,_ = input
        xo, yo = outdims(layer,(x,y))
    else
        error("Not a valid Input size, expected (::Int,::Int) or (::Int,::Int,::Int)")
    end
    a,b = layer.over
    return convert(Float64,xo*yo*Fm), (x,1)
end

function gradientflops(layer::Conv, input::Tuple)
    _,_,_,outputfeaturemaps = size(layer.weight)
    x = 0; y = 0; Fm = 1
    if length(input) == 3
        x,y,Fm = input
        xo, yo = outdims(layer,(x,y))
    elseif length(input) == 2
        x,y,_ = input
        xo, yo = outdims(layer,(x,y))
    else
        error("Not a valid Input size, expected (::Int,::Int) or (::Int,::Int,::Int)")
    end
    noofoper = gradientoperations(layer)
    if size(layer.bias) == ()
        return convert(Float64,xo*yo*noofoper*Fm), (xo,yo,outputfeaturemaps)
    else
        return convert(Float64,((xo*yo+size(layer.bias))*noofoper)*Fm), (xo,yo,outputfeaturemaps)
    end
end

function gradientflops(layer::MaxPool, input::Tuple)
    kernelsizex,kernelsizey = layer.k
    x = 0; y = 0; Fm = 1
    if length(input) == 3
        x,y,Fm = input
        xo, yo = outdims(layer,(x,y))
    elseif length(input) == 2
        x,y,_ = input
        xo, yo = outdims(layer,(x,y))
    else
        error("Not a valid Input size, expected (::Int,::Int) or (::Int,::Int,::Int)")
    end
    outx, outy = outdims(layer,(x,y))
    return convert(Float64,outx*outy*Fm), (outx,outy)
end

gradientflops(layer::GlobalMaxPool, input::Tuple) = error("Must be implemented in a future release")

function gradientflops(layer::MeanPool, input::Tuple)
    kernelsizex,kernelsizey = layer.k
    x = 0; y = 0; Fm = 1
    if length(input) == 3
        x,y,Fm = input
        xo, yo = outdims(layer,(x,y))
    elseif length(input) == 2
        x,y,_ = input
        xo, yo = outdims(layer,(x,y))
    else
        error("Not a valid Input size, expected (::Int,::Int) or (::Int,::Int,::Int)")
    end
    outx, outy = outdims(layer,(x,y))
    return convert(Float64,outx*outy*Fm), (outx,outy)
end

gradientflops(layer::GlobalMeanPool, input::Tuple) = error("Must be implemented in a future release")

function gradientflops(layer::DepthwiseConv, input::Tuple)
    kernelsizex,kernelsizey,inputfeaturemaps,outputfeaturemaps = size(layer.weight)
    x = 0; y = 0; Fm = 1
    if length(input) == 3
        x,y,Fm = input
        xo, yo = outdims(layer,(x,y))
    elseif length(input) == 2
        x,y,_ = input
        xo, yo = outdims(layer,(x,y))
    else
        error("Not a valid Input size, expected (::Int,::Int) or (::Int,::Int,::Int)")
    end
    outx, outy = outdims(layer,(x,y))
    noofoper = gradientoperations(layer)
    if size(layer.bias) == ()
        return convert(Float64,kernelsizex*kernelsizey*noofoper*Fm), (outx,outy)
    else
        return convert(Float64,kernelsizex*kernelsizey*noofoper*Fm + (length(layer.bias)-1)), (outx,outy)
    end
end

function gradientflops(layer::ConvTranspose, input::Tuple)
    kernelsizex,kernelsizey,inputfeaturemaps,outputfeaturemaps = size(layer.weight)
    x = 0; y = 0; Fm = 1
    if length(input) == 3
        x,y,Fm = input
        xo, yo = outdims(layer,(x,y))
    elseif length(input) == 2
        x,y,_ = input
        xo, yo = outdims(layer,(x,y))
    else
        error("Not a valid Input size, expected (::Int,::Int) or (::Int,::Int,::Int)")
    end
    outx, outy = outdims(layer,(x,y))
    noofoper = gradientoperations(layer)
    if size(layer.bias) == ()
        return convert(Float64,kernelsizex*kernelsizey*noofoper*Fm), (outx,outy)
    else
        return convert(Float64,kernelsizex*kernelsizey*noofoper*Fm + (length(layer.bias)-1)), (outx,outy)
    end
end

function gradientflops(layer::CrossCor, input::Tuple)
    kernelsizex,kernelsizey,inputfeaturemaps,outputfeaturemaps = size(layer.weight)
    x = 0; y = 0; Fm = 1
    if length(input) == 3
        x,y,Fm = input
        xo, yo = outdims(layer,(x,y))
    elseif length(input) == 2
        x,y,_ = input
        xo, yo = outdims(layer,(x,y))
    else
        error("Not a valid Input size, expected (::Int,::Int) or (::Int,::Int,::Int)")
    end
    outx, outy = outdims(layer,(x,y))
    noofoper = gradientoperations(layer)
    if size(layer.bias) == ()
        return convert(Float64,kernelsizex*kernelsizey*noofoper*Fm), (outx,outy)
    else
        return convert(Float64,kernelsizex*kernelsizey*noofoper*Fm + (length(layer.bias)-1)), (outx,outy)
    end
end

function gradientflops(layer::Recur{T}, input::Tuple) where {T <: RNNCell}
    inM,inN = input
    WhM,WhN = size(layer.cell.Wh)
    WiM,WiN = size(layer.cell.Wi)
    hM = length(layer.cell.h)
    noofoper = activationoperations(layer)
    if size(layer.cell.b) == ()
        return convert(Float64,((2*WhN*WhM-hM + 2*WiM*WiN-inM)*noofoper)+((2*WhN*WhM-hM)*noofoper)), (hM,1)
    else
        bM = length(layer.cell.b)
        return convert(Float64,((2*WhN*WhM-hM + 2*WiM*WiN-inM + 2*bM)*noofoper)+((2*WhN*WhM-hM + bM)*noofoper)), (bM,1)
    end
end

function gradientflops(layer::Recur{T}, input::Tuple) where {T <: LSTMCell}
    inM,inN = input
    WhM,WhN = size(layer.cell.Wh)
    WiM,WiN = size(layer.cell.Wi)
    hM = length(layer.cell.h)
    noofoper = 3
    if size(layer.cell.b) == ()
        return convert(Float64,((2*WhN*WhM-hM + 2*WiM*WiN-inM)*noofoper)+((2*WhN*WhM-hM + 2*WiM*WiN-inM)*noofoper)+((2*WhN*WhM-hM + 2*WiM*WiN-inM)*noofoper)+((2*WhN*WhM-hM + 2*WiM*WiN-inM)*noofoper)+(3*hM)+((hM*noofoper)+hM)), (hM,1)
    else
        bM = length(layer.cell.b)
        return convert(Float64,((2*WhN*WhM-hM + 2*WiM*WiN-inM + 2*bM)*noofoper)+((2*WhN*WhM-hM + 2*WiM*WiN-inM + 2*bM)*noofoper)+((2*WhN*WhM-hM + 2*WiM*WiN-inM + 2*bM)*noofoper)+((2*WhN*WhM-hM + 2*WiM*WiN-inM + 2*bM)*noofoper)+(3*bM)+((bM*noofoper)+bM)), (bM,1)
    end
end

function gradientflops(layer::Recur{T}, input::Tuple) where {T <: GRUCell}
    inM,inN = input
    WhM,WhN = size(layer.cell.Wh)
    WiM,WiN = size(layer.cell.Wi)
    hM = length(layer.cell.h)
    noofoper = 3
    if size(layer.cell.b) == ()
        return convert(Float64,((2*WhN*WhM-hM + 2*WiM*WiN-inM)*noofoper)+((2*WhN*WhM-hM + 2*WiM*WiN-inM)*noofoper)+((2*WhN*WhM-hM + 2*WiM*WiN-inM)*noofoper)+((2*WhN*WhM-hM + 2*WiM*WiN-inM)*noofoper)+(3*hM)+((hM*noofoper)+hM)), (hM,1)
    else
        bM = length(layer.cell.b)
        return convert(Float64,((2*WhN*WhM-hM + 2*WiM*WiN-inM + 2*bM)*noofoper)+((2*WhN*WhM-hM + 2*WiM*WiN-inM + 2*bM)*noofoper)+((2*WhN*WhM-hM + 2*WiM*(WiN+hM)-inM + 2*bM)*noofoper)+(4bM)), (bM,1)
    end
end


text:
load("src/init.jl")

# Given inputs, ouputs, lambda and current weights,
# Assess the ridge loss function.
# Assume first entry of w is an intercept and do not penalize.

function ridge_error(x, y, l, w)
  (1.0 / 2.0) * sum(map(z -> z^2, y - x * w)) + (l / 2.0) * sum(map(z -> z^2, w[2:length(w)]))
end

ridge_error([1.0 2; 3 4], [1, 2], 1, [1, 1])

function ridge_gradient(x, y, l, w)
	unshift(l * w[2:length(w)], 0) - x' * (y - x * w)
end

ridge_gradient([1.0 1 2; 1 3 4], [1, 2], 1, [1, 1, 1])

function ridge_regression(x, y, l)
	function f(w)
    ridge_error(x, y, l, w)
  end

  function g(w)
    ridge_gradient(x, y, l, w)
  end

	w0 = zeros(size(x, 2))

	solution = gradient_descent(f, g, w0)

	solution.minimum
end

x = [1 1 2; 1 3 3; 1 5 6]
y = [1, 2, 2]

w = ridge_regression(x, y, 0.0)
x * w
norm(w)

w = ridge_regression(x, y, 1.0)
x * w
norm(w)

w = ridge_regression(x, y, 10.0)
x * w
norm(w)

w = ridge_regression(x, y, 100.0)
x * w
norm(w)

p = 10

x = hcat([1, 1], vcat([1:p]', [(p + 1):(2p)]'))
y = [1, 2]

w = ridge_regression(x, y, 0)
x * w
norm(w)

w = ridge_regression(x, y, 1)
x * w
norm(w)

w = ridge_regression(x, y, 10)
x * w
norm(w)


text:
function Element(message::Message)
    element_handle = blpapi_Message_elements(message.handle)
    Element(element_handle, message)
end

function MessageIterator(event::Event)
    return MessageIterator(blpapi_MessageIterator_create(event.handle))
end

function next_message(iter::MessageIterator) :: Union{Nothing, Message}
    result_msg_handle_ref = Ref{Ptr{Cvoid}}(C_NULL)
    result = blpapi_MessageIterator_next(iter.handle, result_msg_handle_ref)

    if result == 0
        return Message(result_msg_handle_ref[])
    else
        return nothing
    end
end

each_message(event::Event) = MessageIterator(event)

function Base.iterate(iter::MessageIterator, state=nothing)
    iter_result = next_message(iter)
    if iter_result == nothing
        return nothing
    else
        return (iter_result, nothing)
    end
end

@inline function has_correlation_id(message::Message, corr_id::CorrelationId) :: Bool
    return corr_id ∈ message.correlation_ids
end

@inline function check_has_correlation_id(message::Message, corr_id::CorrelationId)
    @assert has_correlation_id(message, corr_id) "Got message with unexpected correlation id: $(message.correlation_ids) (expected $corr_id)."
end

max_stars_repo_licenses: OpenSSL MIT

max_issues_repo_licenses: OpenSSL MIT

max_forks_repo_licenses: OpenSSL MIT


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMCoupling.jl/blob/master/LICENSE

# using JuliaFEM
using FEMCoupling
using FEMBase.Test

# This test is decoupled from JuliaFEM by assembling stiffness matrix and
# saving it as a test resource. Stiffness matrix can be read back using
# readdlm.

# Creating 2d geometry for square shaped element and one independent node.
d = 3.0
X = Dict(1 => [0.0, 0.0],
         2 => [d, 0.0],
         3 => [d, d],
         4 => [0.0, d],
         5 => [2*d, d/2])

# Creating one plane_stress element
element1 = Element(Quad4, [1, 2, 3, 4])
update!(element1, "geometry", X)
update!(element1, "youngs modulus", 210.0e8)
update!(element1, "poissons ratio", 1/3)

# Creating Elasticity problem named "test problem" with 2 dimensions.
# problem1 = Problem(Elasticity, "test problem", 2)
# problem1.properties.formulation = :plane_stress
# add_elements!(problem1, [element1])
# assemble!(problem1, 0.0)
# K = full(problem1.assembly.K)
# fn = Pkg.dir("FEMCoupling", "test", "test_dcoupling_2d", "stiffness_matrix")
# writedlm(fn, K)

# Creating Dirichlet problem for boundary conditions with Seg2
# element which fixes nodes 4 and 1.
element2 = Element(Seg2, [4, 1])
update!(element2, "geometry", X)
update!(element2, "displacement 1", 0.0)
update!(element2, "displacement 2", 0.0)
# bc = Problem(Dirichlet, "fixed", 2, "displacement")
# add_elements!(bc, [element2])

###############################################################################
# Creating coupling problem.
coupling1 = Problem(Coupling, "test", 2, "displacement")

# Creating Poi1 elements for coupling nodes and numbering
# them with original node numbers.
coupling_node2 = Element(Poi1, [2])
coupling_node3 = Element(Poi1, [3])
update!([coupling_node2, coupling_node3], "geometry",X)

# Creating Poi1 element for the reference node and
# applying loads for it.
reference_node = Element(Poi1, [5])
update!(reference_node, "geometry", X)
update!(reference_node, "point force 1", 10.0e5)
update!(reference_node, "point force 2", 35.0e5)
update!(reference_node, "point moment 3", 80.0e5)

# Pointing out that coupling nodes and reference nodes
# are for the coupling problem.
add_coupling_nodes!(coupling1, [coupling_node2, coupling_node3])
add_reference_node!(coupling1, reference_node)

# Making step for nonlinear analysis and adding all three problems
# to the step. Running the analysis.
# step = Analysis(Nonlinear)
# add_problems!(step, [problem1, bc, coupling1])
# run!(step)

# Testing results by comparing them against abaqus results.
# Global force vector and plane stress element's displacements
# are compared.

assemble!(coupling1, 0.0)
f = full(sparse(coupling1.assembly.f, 8, 1))
f_expected = [0.0, 0.0, 6.66667e6, 1.75e6, -5.66667e6, 1.75e6, 0.0, 0.0]
# println("f = $f")
# println("f_expected = $f_expected")
@test isapprox(f, f_expected, rtol=1.0e-3)

K = readdlm(Pkg.dir("FEMCoupling", "test", "test_dcoupling_2d", "stiffness_matrix"))
free_dofs = [3, 4, 5, 6]
u = zeros(8)
u[free_dofs] = K[free_dofs, free_dofs] \ f[free_dofs]
u2 = u[3:4]
u3 = u[5:6]
u2_expected = [0.00155379, 0.00196296]
u3_expected = [-0.00146208, 0.0019418]
@test isapprox(u2, u2_expected, rtol=1.0e-5)
@test isapprox(u3, u3_expected, rtol=1.0e-5)

# If reference node is not defined, error is thrown:

coupling2 = Problem(Coupling, "test coupling", 2, "displacement")
add_coupling_nodes!(coupling2, [coupling_node2, coupling_node3])
@test_throws Exception assemble!(coupling2, 0.0)


text:
#!/usr/bin/env julia

using Luxor, Random

using Test

using Random
Random.seed!(42)

function test_hypotrochoid(fname)
    Drawing(1400, 1400, fname)
    origin()
    background("grey20")
    sethue("antiquewhite")
    setopacity(0.8)
    tiles = Tiler(1200, 1200, 6, 6)
    for (pos, n) in tiles
        gsave()
        setline(1.5)
        translate(pos)

        R = 45; r = n; d = 20
        sethue("antiquewhite")
        poly(hypotrochoid(R, n, d, :stroke, stepby=pi/300, vertices=true), :stroke, close=true)
        text(string("R = $(R), r = $(r), d = $(d)"), halign=:center, O.x, O.y+tiles.tileheight/2)

        # test perimeterlength
        ptslist = hypotrochoid(R, n, d, :stroke, stepby=pi/300, vertices=true)
        polyperimeter(ptslist)

        R = 2n; r = 12; d = 5
        sethue("orange")
        (R == r) && (r += 1)
        p = hypotrochoid(R, r, d, :none, stepby=pi/305, vertices=true)
        polydistances(p) # just testing
        poly(p, :stroke, close=true)
        text(string("R = $(R), r = $(r), d = $(d)"), halign=:center, O.x, O.y+tiles.tileheight/2 + 10)

        # testing portions of polys
        portion = rand()
        setline(.25)
        sethue("magenta")
        section = polyportion(p, portion)
        poly(section, :stroke, close=false)
        sethue("green")
        circle(section[end], 2, :stroke)
        poly(polyremainder(p, portion), :stroke, close=false)

        text(string("portion ", round(portion, digits=2)), halign=:center, O.x, O.y+tiles.tileheight/2 + 30)

        sethue("cyan")
        setline(0.3)
        poly(offsetpoly(p, 15), :stroke, close=true)
        text("offset 15", halign=:center, O.x, O.y+tiles.tileheight/2 + 20)

        grestore()
    end
    @test finish() == true
end

using Random

Random.seed!(42)

fname = "hypotrochoid-test.pdf"
test_hypotrochoid(fname)
println("...finished test: output in $(fname)")


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     cell_metadata_filter: -all
#     formats: jl:light,notebooks//ipynb
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.5'
#       jupytext_version: 1.6.0
#   kernelspec:
#     display_name: Julia 1.5.1
#     language: julia
#     name: julia-1.5
# ---

# # Example Julia
#
# ## Documentation de fonction

"""
    bspline(p, j, x)

Return the value at x in [0,1[ of the B-spline with integer nodes of degree p with support starting at j.
Implemented recursively using the [De Boor's Algorithm](https://en.wikipedia.org/wiki/De_Boor%27s_algorithm)

```math
B_{i,0}(x) := \\left\\{
\\begin{matrix}
1 & \\mathrm{if}  \\quad t_i ≤ x < t_{i+1} \\\\
0 & \\mathrm{otherwise}
\\end{matrix}
\\right.
```

```math
B_{i,p}(x) := \\frac{x - t_i}{t_{i+p} - t_i} B_{i,p-1}(x)
+ \\frac{t_{i+p+1} - x}{t_{i+p+1} - t_{i+1}} B_{i+1,p-1}(x).
```
"""
function bspline(p::Int, j::Int, x::Float64)
   if p == 0
       if j == 0
           return 1.0
       else
           return 0.0
       end
   else
       w = (x - j) / p
       w1 = (x - j - 1) / p
   end
   return (w * bspline(p - 1, j, x) + (1 - w1) * bspline(p - 1, j + 1, x))
end

# Switch the next cell from `Markdown` to `Code` format

# ?bspline

# +
using Random, Plots

function generate_data( n = 2000, seed = 1234 )
    seuil = 0.25
    rng = MersenneTwister(seed)
    X1 = rand( rng, n)
    X2 = rand( rng, n)
    U  = rand( rng, n)
    Y  = zeros(Int,n)
    Y[(X1 .<= 0.25) .& (U  .<= seuil)] .= 1
    Y[(X1 .>  0.25) .& (X2 .>= 0.75) .& (U .<= seuil)] .= 1
    Y[(X1 .>  0.25) .& (X2 .<  0.75) .& (U .>  seuil)] .= 1
    return X1, X2, Y
end

X1, X2, Y = generate_data()
scatter(X1,X2, marker_z = Y)
# -
using DataFrames, StatsPlots

data = DataFrame( X1=X1, X2=X2, Y=Y)
head(data)

@df data scatter(:X1,:X2, zcolor= :Y, xaxis = "X1", yaxis="X2", lab="Y")

"""
[x1,x2,x1^2,x1x2,x2^2……x2^6]
"""
function map_features(X1,X2)
    degree = 6
    out = ones(size(X1[:,1]))
    for i=1:6
        for j=0:i
            out = hcat(out,(X1.^(i-j)).*(X2.^j))
        end
    end
    return out
end

X = map_features(X1,X2)


text:
struct NetketJsonDataType <: LazyDataType end

function _lazyload(m::MetricsData, ::NetketJsonDataType)
    m.loaded && return m

    path = m.path
    data_json = JSON.parse(read(path, String))["Output"]

    for iter_data in data_json
        ks = keys(iter_data)
        step = iter_data["Iteration"]

        for k in ks
            k == "Iteration" && continue
            data = iter_data[k]

            if data isa Dict
                expval = data["Mean"] ± data["Sigma"]
                variance = data["Variance"]
                rhat = data["R_hat"]
                taucorr = data["TauCorr"]

                push!(m.data, Symbol(k), step, expval)
                push!(m.data, Symbol(k*"_var"), step, variance)
                push!(m.data, Symbol(k*"_rhat"), step, rhat)
                push!(m.data, Symbol(k*"_tau"), step, taucorr)
            else
                if data isa Integer
                    data = convert(Float64, data)
                end

                push!(m.data, Symbol(k), step, data)
            end
        end
    end

    m.loaded = true
    return m
end


text:
# Propriedades dos materiais
AE = 0.72e6
EI = 24e6

# Propriedades dos elementos
l1 = 600
l2 = 600
A = 1

# Elemento 1
k1 = [
    (4*EI)/l1 (2*EI)/l1 (6*EI)/l1^2 -(6*EI)/l1^2 0 0
    (2*EI)/l1 (4*EI)/l1 (6*EI)/l1^2 -(6*EI)/l1^2 0 0
    (6*EI)/l1^2 (6*EI)/l1^2 (12*EI)/l1^3 -(12*EI)/l1^3 0 0
    -(6*EI)/l1^2 -(6*EI)/l1^2 -(12*EI)/l1^3 (12*EI)/l1^3 0 0
    0 0 0 0 (AE)/l1 -(AE)/l1
    0 0 0 0 -(AE)/l1 (AE)/l1
]

f = [0, 0, 20, 3000]

@show k1 / 100
@show k1[1:4, 1:4]
@show f \ k1

max_stars_repo_name: rsd-consultorias/scripts-engenharia

max_issues_repo_name: rsd-consultorias/scripts-engenharia

max_forks_repo_name: rsd-consultorias/scripts-engenharia


text:
import Winston

function Winston.plot(data::Vector{Float64} , changepoints::Vector{Int})
    p = Winston.plot(data)
    # then put a vertical line everywhere where there is a changept
    for i in changepoints
        Winston.add(p,Winston.LineX(i, color = "red"))
    end
    return(p)
end

function Winston.plot( crops_object::Dict{ASCIIString, Array} )
	# on x axis plot no of chpts
	# on y axis constrained likelihood
	x = crops_object["number"]
	y = crops_object["constrained"]
	p = Winston.FramedPlot(title="Elbow plot for range of penalties",
                               ylabel="Cost",
                               xlabel="Number of segments")
	a = Winston.Points(x, y, kind="circle")
	Winston.add(p, a, Winston.Curve(x, y))
	return(p)
end


text:
# avoid to unzip manually and save the storage
# refer to https://stackoverflow.com/questions/44876852/julia-extract-zip-files-within-a-zip-file
using ZipFile
using DelimitedFiles

using GradientBoost
using GradientBoost.Util
using GradientBoost.ML
using Statistics

r = ZipFile.Reader("houses.zip")
# the first 27 lines are description
data = readdlm(f, skipstart=27)
# extract the 15-th line manually
line15 = "The file cadata.txt contains all the the variables. Specifically, it contains median house value, median income, housing median age, total rooms, total bedrooms, population, households, latitude, and longitude in that order. "
colnames = split(line15, ", ")[2:end]
colnames[1] = split(colnames[1], "contains ")[2]
colnames[end] = split(colnames[end])[2]
colnames = Array{String, 1}(colnames)

# notice the unit is $100,000
Y = data[:, 1] / 100_000
X = data[:, 2:end]

function Housing_GB(X::Matrix{Float64}, Y::Vector{Float64}; maxdepth = 6, nsubfeatures = 0, sampling_rate = 0.8, ν = 0.1)
    # divide into training and test set
    train_ind, test_ind = holdout(size(X, 1), 0.2)
    train_X = X[train_ind, :]
    train_Y = Y[train_ind]
    test_X = X[test_ind, :]
    test_Y = Y[test_ind]

    err_train = Float64[]
    err_test = Float64[]
    # train under number of iterations
    for niter in 100:100:800
        println("current niter = $niter")
        # train
        gbdt = GBDT(; loss_function=LeastSquares(),
                      sampling_rate=sampling_rate,
                      learning_rate=ν,
                      num_iterations=niter,
                      tree_options = Dict(
                        :maxdepth => maxdepth,
                        :nsubfeatures => nsubfeatures
                      ))
        gbl = GBLearner(gbdt, :regression)
        fit!(gbl, train_X, train_Y)

        train_Y_hat = predict!(gbl, train_X)
        test_Y_hat = predict!(gbl, test_X)
        push!(err_train, mean(abs.(train_Y - train_Y_hat)))
        push!(err_test, mean(abs.(test_Y - test_Y_hat)))
    end
    return err_train, err_test
end

err_train, err_test = Housing_GB(X, Y)

using Plots
niters = 100:100:800
plot(niters, err_train, label="Train Error", xlab = "Iterations M", ylab = "Absolute Error", title = "Training and Test Absolute Error (with square loss)")
plot!(niters, err_test, label="Test Error")
savefig("fig10.13_square_loss.png")

max_stars_repo_name: IanFla/Statistical-Learning

max_issues_repo_name: IanFla/Statistical-Learning

max_forks_repo_name: IanFla/Statistical-Learning


text:
using Test, Base.CoreLogging
import Base.CoreLogging: Info

using CUDAnative, CUDAdrv
using GPUCompiler
import Adapt, LLVM

include("util.jl")

# TODO: also run the tests with CuArrays.jl
const CuArray = CUDAnative.CuHostArray

@testset "CUDAnative" begin

@test CUDAnative.functional(true)

CUDAnative.version()
CUDAnative.release()

# ensure CI is using the requested version
if haskey(ENV, "CI") && haskey(ENV, "JULIA_CUDA_VERSION")
  @test CUDAnative.release() == VersionNumber(ENV["JULIA_CUDA_VERSION"])
end

length(devices()) > 0 || error("The CUDAnative.jl test suite requires a CUDA device")
include("init.jl")
include("pointer.jl")
include("codegen.jl")

capability(device()) >= v"2.0" || error("The CUDAnative.jl test suite requires a CUDA device with compute capability 2.0 or higher")
include("device/codegen.jl")
include("device/execution.jl")
include("device/pointer.jl")
include("device/array.jl")
include("device/cuda.jl")
include("device/wmma.jl")

include("nvtx.jl")

include("examples.jl")

end


text:
let
    # Include inside this quote any rules that should have FastMath versions
    # IMPORTANT:
    # Do not add any rules here for functions that do not have varients in Base.FastMath
    # e.g. do not add `foo` unless `Base.FastMath.foo_fast` exists.
    fastable_ast = quote
        #  Trig-Basics
        ## use `sincos` to compute `sin` and `cos` at the same time
        ## for the rules for `sin` and `cos`
        ## See issue: https://github.com/JuliaDiff/ChainRules.jl/issues/291
        ## sin
        function rrule(::typeof(sin), x::Number)
            sinx, cosx = sincos(x)
            sin_pullback(Δy) = (NoTangent(), cosx' * Δy)
            return (sinx, sin_pullback)
        end

        function frule((_, Δx), ::typeof(sin), x::Number)
            sinx, cosx = sincos(x)
            return (sinx, cosx * Δx)
        end

        ## cos
        function rrule(::typeof(cos), x::Number)
            sinx, cosx = sincos(x)
            cos_pullback(Δy) = (NoTangent(), -sinx' * Δy)
            return (cosx, cos_pullback)
        end

        function frule((_, Δx), ::typeof(cos), x::Number)
            sinx, cosx = sincos(x)
            return (cosx, -sinx * Δx)
        end

        @scalar_rule tan(x) 1 + Ω ^ 2

        # Trig-Hyperbolic
        @scalar_rule cosh(x) sinh(x)
        @scalar_rule sinh(x) cosh(x)
        @scalar_rule tanh(x) 1 - Ω ^ 2

        # Trig- Inverses
        @scalar_rule acos(x) -(inv(sqrt(1 - x ^ 2)))
        @scalar_rule asin(x) inv(sqrt(1 - x ^ 2))
        @scalar_rule atan(x) inv(1 + x ^ 2)

        # Trig-Multivariate
        @scalar_rule atan(y, x) @setup(u = x ^ 2 + y ^ 2) (x / u, -y / u)
        @scalar_rule sincos(x) @setup((sinx, cosx) = Ω) cosx -sinx

        # exponents
        @scalar_rule cbrt(x) inv(3 * Ω ^ 2)
        @scalar_rule inv(x) -(Ω ^ 2)
        @scalar_rule sqrt(x) inv(2Ω)  # gradient +Inf at x==0
        @scalar_rule exp(x) Ω
        @scalar_rule exp10(x) Ω * log(oftype(x, 10))
        @scalar_rule exp2(x) Ω * log(oftype(x, 2))
        @scalar_rule expm1(x) exp(x)
        @scalar_rule log(x) inv(x)
        @scalar_rule log10(x) inv(x) / log(oftype(x, 10))
        @scalar_rule log1p(x) inv(x + 1)
        @scalar_rule log2(x) inv(x) / log(oftype(x, 2))

        # Unary complex functions
        ## abs
        function frule((_, Δx), ::typeof(abs), x::Union{Real, Complex})
            Ω = abs(x)
            # `ifelse` is applied only to denominator to ensure type-stability.
            signx = x isa Real ? sign(x) : x / ifelse(iszero(x), one(Ω), Ω)
            return Ω, realdot(signx, Δx)
        end

        function rrule(::typeof(abs), x::Union{Real, Complex})
            Ω = abs(x)
            function abs_pullback(ΔΩ)
                signx = x isa Real ? sign(x) : x / ifelse(iszero(x), one(Ω), Ω)
                return (NoTangent(), signx * real(ΔΩ))
            end
            return Ω, abs_pullback
        end

        ## abs2
        function frule((_, Δz), ::typeof(abs2), z::Union{Real, Complex})
            return abs2(z), 2 * realdot(z, Δz)
        end

        function rrule(::typeof(abs2), z::Union{Real, Complex})
            function abs2_pullback(ΔΩ)
                Δu = real(ΔΩ)
                return (NoTangent(), 2Δu*z)
            end
            return abs2(z), abs2_pullback
        end

        ## conj
        function frule((_, Δz), ::typeof(conj), z::Union{Real, Complex})
            return conj(z), conj(Δz)
        end
        function rrule(::typeof(conj), z::Union{Real, Complex})
            function conj_pullback(ΔΩ)
                return (NoTangent(), conj(ΔΩ))
            end
            return conj(z), conj_pullback
        end

        ## angle
        function frule((_, Δx), ::typeof(angle), x)
            Ω = angle(x)
            # `ifelse` is applied only to denominator to ensure type-stability.
            n = ifelse(iszero(x), one(real(x)), abs2(x))
            ∂Ω = _imagconjtimes(x, Δx) / n
            return Ω, ∂Ω
        end

        function rrule(::typeof(angle), x::Real)
            function angle_pullback(ΔΩ::Real)
                return (NoTangent(), ZeroTangent())
            end
            function angle_pullback(ΔΩ)
                Δu, Δv = reim(ΔΩ)
                return (NoTangent(), im*Δu/ifelse(iszero(x), one(x), x))
                # `ifelse` is applied only to denominator to ensure type-stability.
            end
            return angle(x), angle_pullback
        end
        function rrule(::typeof(angle), z::Complex)
            function angle_pullback(ΔΩ)
                x,  y  = reim(z)
                Δu, Δv = reim(ΔΩ)
                # `ifelse` is applied only to denominator to ensure type-stability.
                n = ifelse(iszero(z), one(real(z)), abs2(z))
                return (NoTangent(), (-y + im*x)*Δu/n)
            end
            return angle(z), angle_pullback
        end

        # Binary functions

        ## `hypot`
        function frule(
            (_, Δx, Δy),
            ::typeof(hypot),
            x::T,
            y::T,
        ) where {T<:Union{Real,Complex}}
            Ω = hypot(x, y)
            n = ifelse(iszero(Ω), one(Ω), Ω)
            ∂Ω = (realdot(x, Δx) + realdot(y, Δy)) / n
            return Ω, ∂Ω
        end

        function rrule(::typeof(hypot), x::T, y::T) where {T<:Union{Real,Complex}}
            Ω = hypot(x, y)
            function hypot_pullback(ΔΩ)
                c = real(ΔΩ) / ifelse(iszero(Ω), one(Ω), Ω)
                return (NoTangent(), c * x, c * y)
            end
            return (Ω, hypot_pullback)
        end

        @scalar_rule x + y (true, true)
        @scalar_rule x - y (true, -1)
        @scalar_rule x / y (one(x) / y, -(Ω / y))

        ## power
        # literal_pow is in base.jl
        function frule((_, Δx, Δp), ::typeof(^), x::Number, p::Number)
            y = x ^ p
            _dx = _pow_grad_x(x, p, float(y))
            if iszero(Δp)
                # Treat this as a strong zero, to avoid NaN, and save the cost of log
                return y, _dx * Δx
            else
                # This may do real(log(complex(...))) which matches ProjectTo in rrule
                _dp = _pow_grad_p(x, p, float(y))
                return y, muladd(_dp, Δp, _dx * Δx)
            end
        end

        function rrule(::typeof(^), x::Number, p::Number)
            y = x^p
            project_x = ProjectTo(x)
            project_p = ProjectTo(p)
            function power_pullback(dy)
                _dx = _pow_grad_x(x, p, float(y))
                return (
                    NoTangent(),
                    project_x(conj(_dx) * dy),
                    # _pow_grad_p contains log, perhaps worth thunking:
                    @thunk project_p(conj(_pow_grad_p(x, p, float(y))) * dy)
                )
            end
            return y, power_pullback
        end

        ## `rem`
        @scalar_rule(
            rem(x, y),
            @setup((u, nan) = promote(x / y, NaN16), isint = isinteger(x / y)),
            (ifelse(isint, nan, one(u)), ifelse(isint, nan, -trunc(u))),
        )
        ## `min`, `max`
        @scalar_rule max(x, y) @setup(gt = x > y) (gt, !gt)
        @scalar_rule min(x, y) @setup(gt = x > y) (!gt, gt)

        # Unary functions
        @scalar_rule +x true
        @scalar_rule -x -1

        ## `sign`
        function frule((_, Δx), ::typeof(sign), x)
            n = ifelse(iszero(x), one(real(x)), abs(x))
            Ω = x isa Real ? sign(x) : x / n
            ∂Ω = Ω * (_imagconjtimes(Ω, Δx) / n) * im
            return Ω, ∂Ω
        end

        function rrule(::typeof(sign), x)
            n = ifelse(iszero(x), one(real(x)), abs(x))
            Ω = x isa Real ? sign(x) : x / n
            function sign_pullback(ΔΩ)
                ∂x = Ω * (_imagconjtimes(Ω, ΔΩ) / n) * im
                return (NoTangent(), ∂x)
            end
            return Ω, sign_pullback
        end

        # product rule requires special care for arguments where `mul` is non-commutative
        function frule((_, Δx, Δy), ::typeof(*), x::Number, y::Number)
            # Optimized version of `Δx .* y .+ x .* Δy`. Also, it is potentially more
            # accurate on machines with FMA instructions, since there are only two
            # rounding operations, one in `muladd/fma` and the other in `*`.
            ∂xy = muladd(Δx, y, x * Δy)
            return x * y, ∂xy
        end
        frule((_, Δx), ::typeof(*), x::Number) = x, Δx

        function rrule(::typeof(*), x::Number, y::Number)
            function times_pullback2(Ω̇)
                ΔΩ = unthunk(Ω̇)
                return (NoTangent(), ProjectTo(x)(ΔΩ * y'), ProjectTo(y)(x' * ΔΩ))
            end
            return x * y, times_pullback2
        end
        # While 3-arg * calls 2-arg *, this is currently slow in Zygote:
        # https://github.com/JuliaDiff/ChainRules.jl/issues/544
        function rrule(::typeof(*), x::Number, y::Number, z::Number)
            function times_pullback3(Ω̇)
                ΔΩ = unthunk(Ω̇)
                return (
                    NoTangent(),
                    ProjectTo(x)(ΔΩ * y' * z'),
                    ProjectTo(y)(x' * ΔΩ * z'),
                    ProjectTo(z)(x' * y' * ΔΩ),
                )
            end
            return x * y * z, times_pullback3
        end
        # Instead of this recursive rule for N args, you could write the generic case
        # directly, by nesting ntuples, but this didn't infer well:
        # https://github.com/JuliaDiff/ChainRules.jl/pull/547/commits/3558951c9f1b3c70e7135fd61d29cc3b96a68dea
        function rrule(::typeof(*), x::Number, y::Number, z::Number, more::Number...)
            Ω3, back3 = rrule(*, x, y, z)
            Ω4, back4 = rrule(*, Ω3, more...)
            function times_pullback4(Ω̇)
                Δ4 = back4(unthunk(Ω̇))  # (0, ΔΩ3, Δmore...)
                Δ3 = back3(Δ4[2])       # (0, Δx, Δy, Δz)
                return (Δ3..., Δ4[3:end]...)
            end
            return Ω4, times_pullback4
        end
        rrule(::typeof(*), x::Number) = rrule(identity, x)
    end  # fastable_ast

    # Rewrite everything to use fast_math functions, including the type-constraints
    fast_ast = Base.FastMath.make_fastmath(fastable_ast)

    # Guard against mistakenly defining something as fast-able when it isn't.
    # NOTE: this check is not infallible, it will be tricked if a function itself is not
    # fastmath_able but it's pullback uses something that is. So manual check should also be
    # done.
    non_transformed_definitions = intersect(fastable_ast.args, fast_ast.args)
    filter!(expr->!(expr isa LineNumberNode), non_transformed_definitions)
    if !isempty(non_transformed_definitions)
        error(
            "Non-FastMath compatible rules defined in fastmath_able.jl. \n Definitions:\n" *
            join(non_transformed_definitions, "\n")
        )
        # This error() may not play well with Revise. But a wanring @error does:
        # @error "Non-FastMath compatible rules defined in fastmath_able.jl." non_transformed_definitions
    end

    eval(fast_ast)
    eval(fastable_ast)  # Get original definitions
    # we do this second so it overwrites anything we included by mistake in the fastable
end

## power
# Thes functions need to be defined outside the eval() block.
# The special cases they aim to hit are in POWERGRADS in tests.
_pow_grad_x(x, p, y) = (p * y / x)
function _pow_grad_x(x::Real, p::Real, y)
    return if !iszero(x) || p < 0
        p * y / x
    elseif isone(p)
        one(y)
    elseif iszero(p) || p > 1
        zero(y)
    else
        oftype(y, Inf)
    end
end

_pow_grad_p(x, p, y) = y * log(complex(x))
function _pow_grad_p(x::Real, p::Real, y)
    return if !iszero(x)
        y * real(log(complex(x)))
    elseif p > 0
        zero(y)
    else
        oftype(y, NaN)
    end
end


text:
module FEniCS

using PyCall
using Requires

# determine if we can include mshr
include_mshr = false
try
    pyimport_conda("mshr", "mshr=2019.1.0", "conda-forge")
    global include_mshr = true
catch ee
    print("mshr has not been included")
end

# global PyObject constants that get initialized at runtime.  We
# initialize them here (rather than via "global foo = ..." in __init__)
# so that their type is known at compile-time.
const fenics = PyCall.PyNULL()
const ufl = PyCall.PyNULL()
const mshr = PyCall.PyNULL()

function __init__()
	@require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" include("jplot.jl")
	@require ProgressMeter="92933f4c-e287-5a05-a399-4b506db050ca" begin
	  using ProgressMeter
    end
    copy!(fenics, pyimport_conda("fenics", "fenics=2019.1.0", "conda-forge"))
    copy!(ufl, pyimport_conda("ufl", "ufl=2019.1.0", "conda-forge"))
    include_mshr && copy!(mshr, pyimport_conda("mshr", "mshr=2019.1.0", "conda-forge"))

    # initialize constants at loadtime
    global dx = Measure(fenics.dx)
    global ds = Measure(fenics.ds)
    global dS = Measure(fenics.dS)
    global dP = Measure(fenics.dP)
    global tetrahedron = fenics.tetrahedron
    global hexahedron = fenics.hexahedron #matplotlib cannot handle hexahedron elements
    global triangle = fenics.triangle
    global quadrilateral = fenics.quadrilateral
    global CellType = fenics.CellType
end
#the below code is an adaptation of aleadev.FEniCS.jl
import Base: size, length, show, *, +, -,/, repr, div, sqrt,split,write
abstract type
  fenicsobject
end #creates placeholder for the fenicsobject type

fenicspycall(object::fenicsobject, func::Union{Symbol,String}, args...) = getproperty(object.pyobject,func)(args...)
export fenicspycall

macro fenicsclass(name::Symbol, base1::Symbol=:fenicsobject)
  impl = Symbol(name, "Impl")
  esc(quote
    abstract type
       $name <: $base1
     end
    struct $impl <: $name
      pyobject::PyObject
    end
    $(name)(pyobject::PyObject) = $impl(pyobject)
  end)
end
export fenicsclass

@enum LOGLEVEL begin
    CRITICAL  = 50
    ERROR     = 40
    WARNING   = 30
    INFO      = 20
    PROGRESS  = 16
    TRACE     = 13
    DBG       = 10
end
set_log_level(lvl::LOGLEVEL) = set_log_level(Int(lvl))
set_log_level(lvl::Int) = fenics.set_log_level(lvl)

str(obj::fenicsobject) = fenicspycall(obj, :__str__)
repr(obj::fenicsobject) = fenicspycall(obj, :__repr__)
show(io::IO, obj::fenicsobject) = show(io, repr(obj))
Docs.getdoc(obj::fenicsobject) = obj.pyobject.__doc__
export str, repr

include("jmesh.jl") #this file contains the mesh functions
include("jfem.jl") #this file contains the fem functions
include("jmisc.jl") #this file contains various miscallaneous functions to assist with solving etc
include("jsolve.jl") #this file contains the solver functions/routines
include("jinterface.jl")
include_mshr && include("fmshr.jl") #this file contains various geometrical objects using the mshr package

end #module


text:
"""
Main module for SimHPC.jl – a discrete event process oriented simulation framework for Julia.
"""
module SimHPC

  using DataStructures
  using Dates
  using ResumableFunctions

  import Base.run, Base.isless, Base.show, Base.yield, Base.get
  import Base.(&), Base.(|)
  import Dates.now

  export AbstractEvent, Environment, value, state, environment
  export Event, succeed, fail, @callback, remove_callback
  export timeout
  export Operator, (&), (|)
  export @resumable, @yield
  export AbstractProcess, Simulation, run, now, active_process, StopSimulation
  export Process, @process, interrupt
  export Container, Resource, Store, put, get, request, release, cancel
  export nowDatetime

  include("base.jl")
  include("events.jl")
  include("operators.jl")
  include("simulations.jl")
  include("processes.jl")
  include("resources/base.jl")
  include("resources/containers.jl")
  include("resources/stores.jl")
  include("utils/time.jl")
end

max_stars_repo_name: zhangliye/SimHPC

max_issues_repo_name: zhangliye/SimHPC

max_forks_repo_name: zhangliye/SimHPC


text:
; Maze generate in JL
; Joe Wingbermuehle

(define width 39)       ; Width of the maze (must be odd)
(define height 23)      ; Height of the maze (must be odd)
(define seed 17)        ; Random number seed.

; Generate a random number.
(define get-rand (lambda (s) (mod (+ (* 19 s) 1) 16383)))

; Set an element in a list.
(define set-element (lambda (lst n v)
   (if n
      (cons (head lst) (set-element (rest lst) (- n 1) v))
      (cons v (rest lst)))))

; Get a particular list element.
(define get-element (lambda (lst n)
   (if n
      (get-element (rest lst) (- n 1))
      (head lst))))

; Initialize the maze matrix.
(define init-maze (lambda (x y)
   (if (= y 1)
      (if (= x 1) (list 0) (cons 0 (init-maze (- x 1) y)))
      (if (= y height)
         (if (= x 1)
            (cons 0 (init-maze width (- y 1)))
            (cons 0 (init-maze (- x 1) y)))
         (if (= x 1)
            (cons 0 (init-maze width (- y 1)))
            (if (= x width)
               (cons 0 (init-maze (- x 1) y))
               (cons 1 (init-maze (- x 1) y))))))))

(define update-x (lambda (x d) (+ x (get-element (list 1 -1 0 0) (mod d 4)))))

(define update-y (lambda (y d) (+ y (get-element (list 0 0 1 -1) (mod d 4)))))

; Carve a maze.
(define carve-maze (lambda (maze rand x y c)
   (define maze (set-element maze (+ (* y width) x) 0))
   (define x1 (update-x x rand))
   (define y1 (update-y y rand))
   (define x2 (update-x x1 rand))
   (define y2 (update-y y1 rand))
   (if (and (> x2 0) (< x2 width) (> y2 0) (< y2 height)
            (get-element maze (+ (* y1 width) x1))
            (get-element maze (+ (* y2 width) x2)))
      (begin
         (define maze (set-element maze (+ (* y1 width) x1) 0))
         (define maze (carve-maze maze (get-rand rand) x2 y2 0))
         (if (< c 4) (carve-maze maze (+ rand 1) x y (+ c 1)) maze))
      (if (< c 4) (carve-maze maze (+ rand 1) x y (+ c 1)) maze))))

; Initialize and carve a maze.
(define generate-maze (lambda ()
   (define init (init-maze width height))
   (define carved (carve-maze init seed 2 2 0))
   (define temp (set-element carved (+ (* 1 width) 2) 0))
   (set-element temp (+ (* (- height 2) width) (- width 3)) 0)))

; Show a maze.
(define show-maze (lambda (maze x y)
   (define element (get-element maze (+ (* y width) x)))
   (if (= element 1) (print "[]") (print "  "))
   (if (= (+ x 1) width)
      (begin
         (print "\n")
         (if (< (+ y 1) height) (show-maze maze 0 (+ y 1))))
      (show-maze maze (+ x 1) y))))

(show-maze (generate-maze) 0 0)

max_stars_repo_name: ambotaku/Pico-JL

max_issues_repo_name: ambotaku/Pico-JL

max_forks_repo_name: ambotaku/Pico-JL


text:
module MLJModels

import MLJModelInterface
import MLJModelInterface: MODEL_TRAITS

using ScientificTypes
using MLJBase
import MLJBase: @load
import MLJBase: Table, Continuous, Count, Finite, OrderedFactor, Multiclass

using Requires, Pkg, Pkg.TOML, OrderedCollections, Parameters
using Tables, CategoricalArrays, StatsBase, Statistics, Dates
import Distributions

# for administrators to update Metadata.toml:
export @update, check_registry

# from loading.jl:
export load, @load, info

# from model_search:
export models, localmodels

# from model/Constant
export ConstantRegressor, ConstantClassifier,
        DeterministicConstantRegressor, DeterministicConstantClassifier

# from model/ThresholdPredictors
export BinaryThresholdPredictor

# from model/Transformers
export FeatureSelector, StaticTransformer, UnivariateDiscretizer,
    UnivariateStandardizer, Standardizer, UnivariateBoxCoxTransformer,
    OneHotEncoder, ContinuousEncoder, FillImputer, UnivariateFillImputer,
    UnivariateTimeTypeToContinuous

const srcdir = dirname(@__FILE__) # the directory containing this file

if VERSION < v"1.3"
    nonmissingtype = ScientificTypes.nonmissing
end
nonmissing = nonmissingtype

include("metadata.jl")
include("model_search.jl")
include("loading.jl")
include("registry/src/Registry.jl")
include("registry/src/check_registry.jl")
import .Registry.@update

# load built-in models:
include("builtins/Constant.jl")
include("builtins/Transformers.jl")
include("builtins/ThresholdPredictors.jl")

const INFO_GIVEN_HANDLE = Dict{Handle,Any}()
const PKGS_GIVEN_NAME   = Dict{String,Vector{String}}()
const AMBIGUOUS_NAMES   = String[]
const NAMES             = String[]

metadata_file = joinpath(srcdir, "registry", "Metadata.toml")

merge!(INFO_GIVEN_HANDLE, info_given_handle(metadata_file))
merge!(PKGS_GIVEN_NAME, pkgs_given_name(INFO_GIVEN_HANDLE))
append!(AMBIGUOUS_NAMES, ambiguous_names(INFO_GIVEN_HANDLE))
append!(NAMES, model_names(INFO_GIVEN_HANDLE))
@info "Model metadata loaded from registry. "

# lazily load in strap-on model interfaces for external packages:
function __init__()
    @require(DecisionTree="7806a523-6efd-50cb-b5f6-3fa6f1930dbb",
             include("DecisionTree.jl"))
    @require(GLM="38e38edf-8417-5370-95a0-9cbb8c7f171a",
             include("GLM.jl"))
    @require(Clustering="aaaa29a8-35af-508c-8bc3-b662a17a0fe5",
             include("Clustering.jl"))
    @require(NaiveBayes="9bbee03b-0db5-5f46-924f-b5c9c21b8c60",
             include("NaiveBayes.jl"))
    @require(XGBoost = "009559a3-9522-5dbb-924b-0b6ed2b22bb9",
             include("XGBoost.jl"))
    @require(LIBSVM="b1bec4e5-fd48-53fe-b0cb-9723c09d164b",
             include("LIBSVM.jl"))
    @require(NearestNeighbors="b8a86587-4115-5ab1-83bc-aa920d37bbce",
             include("NearestNeighbors.jl"))
    @require(MultivariateStats = "6f286f6a-111f-5878-ab1e-185364afe411",
             include("MultivariateStats.jl"))
end

end # module


text:
module SIR

using Test

include("../../example/sir.jl")

g_result = run_gillespie(MersenneTwister(1), (0.0, T))[1]
g_result_rec = run_gillespie_record(MersenneTwister(1), T)[1]
m_result = run_manually(MersenneTwister(1), T)
m_result_rec = run_manually_record(MersenneTwister(1), T)

for sym in (:S, :I, :R)
    @eval @test g_result.$sym == g_result_rec.$sym == m_result.$sym == m_result_rec.$sym
    @eval @test getts(getentries(g_result_rec.$sym)) == getts(getentries(m_result_rec.$sym))
    @eval @test getvs(getentries(g_result_rec.$sym)) == getvs(getentries(m_result_rec.$sym))
end

end # module SIRTest


text:
function makesingleplot(vec1::AbstractVector, vec2::AbstractVector,
                                        E::Int, tau_s::Int, tau_p::Int,
                                        var1name::ASCIIString, var2name::ASCIIString;
                                        lib_start=1, xmin=false, xmax=false, ymin=false,
                                        ymax=false, nboots=0, plot=true,
                                        libsizemin=0, libsizemax=0,
                                        npred=0, pred_start=0)
    """
    vec1: Time series 1
    vec2: Time series 2
    libsizemin: Minimum library size
    ibsizemax: Maximum library size
    E: System dimensionality
    tau_s: Lag for manifold reconstruction
    tau_p: Lag for causal effect
    npred: Number of points to predict
    pred_start: Start for prediction
    var1name: Name of first variable
    var2name: Name of second variable
    lib_start: Start of library
    """
    @pyimport matplotlib.pyplot as plt
    shadowmat_dict_vec1, distmat_dict_vec1 = precalc_manif_dists(E:E, tau_s:tau_s, vec1)
    shadowmat_dict_vec2, distmat_dict_vec2 = precalc_manif_dists(E:E, tau_s:tau_s, vec2)

    librange12, yval_12 = calcCCM(vec1, vec2,
                                                        shadowmat_dict_vec1, distmat_dict_vec1,
                                                        E, tau_s, tau_p; libsizemin=libsizemin, libsizemax=libsizemax,
                                                        npred=npred, pred_start=pred_start,
                                                        lib_start=lib_start, nboots=nboots
                                                    )

    librange21, yval_21 = calcCCM(vec2, vec1,
                                                        shadowmat_dict_vec2, distmat_dict_vec2,
                                                        E, tau_s, tau_p; libsizemin=libsizemin, libsizemax=libsizemax,
                                                        npred=npred, pred_start=pred_start,
                                                        lib_start=lib_start, nboots=nboots
                                                    )

    plt.plot(librange12, yval_12, label = "$var2name influences $var1name?")
    plt.plot(librange21, yval_21, label = "$var1name influences $var2name?")
    plt.xlabel("L")
    plt.ylabel("\$\\rho_{ccm}\$")
    plt.legend(loc="lower right")
    ax = plt.gca()
    ax[:xaxis][:set_ticks_position]("bottom")
    ax[:yaxis][:set_ticks_position]("left")
    if xmin != false
        ax[:set_xlim]((xmin, xmax))
    end
    if ymin != false
        ax[:set_ylim]((ymin, ymax))
    end
    if plot
        plt.show()
    else
        plt.close()
    end
end

function makeoptimizationplots(vec1::AbstractVector, vec2::AbstractVector, E_vals::AbstractVector,
                                                    tau_s_vals::AbstractVector, tau_p_vals::AbstractVector,
                                                    var1name::ASCIIString, var2name::ASCIIString;
                                                    nreps=5, b_offset=1, ncols=28, left_E=false, left_tau_p=false,
                                                    right_E=false, right_tau_p=false, lagunit=1, unit=false,
                                                    imfont="medium", nboots=0, legend=true,
                                                    lib_start=1, plot=true, libsizemin=0, libsizemax=0,
                                                    pred_start=0, npred=0, show_tau_s=true)

    println("Calculating manifolds")
    shadowmat_dict_vec1, distmat_dict_vec1 = precalc_manif_dists(E_vals, tau_s_vals, vec1)
    shadowmat_dict_vec2, distmat_dict_vec2 = precalc_manif_dists(E_vals, tau_s_vals, vec2)

    println("\nCoordDescent opt1")
    res12 = CoordDescentOpt(vec1, vec2,
                                                shadowmat_dict_vec1, distmat_dict_vec1,
                                                E_vals, tau_s_vals, tau_p_vals;
                                                libsizemin=libsizemin, libsizemax=libsizemax,
                                                npred=npred, pred_start=pred_start,
                                                nreps=nreps, nboots=nboots)
    println("\nCoordDescent opt2")
    res21 =CoordDescentOpt(vec2, vec1,
                                                shadowmat_dict_vec2, distmat_dict_vec2,
                                                E_vals, tau_s_vals, tau_p_vals;
                                                libsizemin=libsizemin, libsizemax=libsizemax,
                                                npred=npred, pred_start=pred_start,
                                                nreps=nreps, nboots=nboots)

    libsizemin_12 = max(res12["E"] + b_offset + 1, 10)
    libsizemin_21 = max(res21["E"] + b_offset + 1, 10)

    println("\nstarting calcCCM1")
    librange12, yval_12 = calcCCM(vec1, vec2,
                                                        shadowmat_dict_vec1, distmat_dict_vec1,
                                                        res12["E"], res12["tau_s"], res12["tau_p"];
                                                        libsizemin=libsizemin_12, libsizemax=libsizemax,
                                                        npred=npred, pred_start=pred_start,
                                                        lib_start=1, nboots=nboots)
    println("starting calcCCM2")
    librange21, yval_21 = calcCCM(vec2, vec1,
                                                        shadowmat_dict_vec2, distmat_dict_vec2,
                                                        res12["E"], res12["tau_s"], res12["tau_p"];
                                                        libsizemin=libsizemin_21, libsizemax=libsizemax,
                                                        npred=npred, pred_start=pred_start,
                                                        lib_start=1, nboots=nboots)

    @pyimport matplotlib.pyplot as plt
    fig = plt.figure()
    label1 = "\n$(ucfirst(var2name)) influences $var1name?"
    label2 = "\n$(ucfirst(var1name)) influences $var2name?"

    rho1   = @sprintf("%.2f", res12["rho"])
    rho2   = @sprintf("%.2f", res21["rho"])

    if show_tau_s
        statsuff1 = "$(res12["tau_s"])"
        statsuff2 = "$(res21["tau_s"])"
    else
        statsuff1 = ""
        statsuff2 = ""
    end
    stats1 = "   E = $(res12["E"]), \$\\tau_{p}\$ = $(int(res12["tau_p"] * lagunit))" * statsuff1
    stats2 = "   E = $(res21["E"]), \$\\tau_{p}\$ = $(int(res21["tau_p"] * lagunit))" * statsuff2

    ax1 = plt.subplot2grid((2,ncols), (0,0), rowspan=2, colspan=ifloor(ncols/3))

    ax1[:plot](librange12, yval_12, label="$label1\n$stats1")
    ax1[:plot](librange21, yval_21, label="$label2\n$stats2")

    ax1[:yaxis][:set_ticks_position]("left")
    ax1[:xaxis][:set_ticks_position]("bottom")

    ax1[:set_ylabel]("\$\\rho_{ccm}\$", fontsize="x-large")
    ax1[:set_xlabel]("L", fontsize="large", labelpad=12)

    if legend == true
        ax1[:legend](loc=4, fontsize="x-small")
    end

    mat12, Es12, taus12 = get_E_taupcurves(vec1, vec2,
                                                                        shadowmat_dict_vec1, distmat_dict_vec1,
                                                                        libsizemax, res12["E"], res12["tau_s"], res12["tau_p"],
                                                                        E_vals, tau_p_vals, npred, pred_start;
                                                                        left_E=left_E, left_tau_p=left_tau_p,
                                                                        right_E=right_E, right_tau_p=right_tau_p, nboots=nboots
                                                                    )

    mat21, Es21, taus21 = get_E_taupcurves(vec2, vec1,
                                                                        shadowmat_dict_vec2, distmat_dict_vec2,
                                                                        libsizemax, res21["E"], res21["tau_s"], res21["tau_p"],
                                                                        E_vals, tau_p_vals, npred, pred_start;
                                                                        left_E=left_E, left_tau_p=left_tau_p,
                                                                        right_E=right_E, right_tau_p=right_tau_p, nboots=nboots
                                                                    )

    catmat = vcat(mat12, mat21)
    if minimum(catmat) < 0
        vmin = 0
        llab = "\$\\leq 0\$"
    else
        vmin = minimum(catmat)
        llab = @sprintf("%.2f", vmin)
    end

    vmax = 1  # maximum correlation

    if unit != false
        xlabel = "\$\\tau_p\$\ ($unit)"
    else
        xlabel =  "\$\\tau_p\$"
    end

    ax2 = plt.subplot2grid((2,ncols), (0,ifloor(ncols/3)), colspan=2*ifloor(ncols/3))
    imax12 = plotheatmap(ax2, mat12, vmin, vmax, label1, taus12, Es12, imfont, lagunit, xlabel)

    ax3 = plt.subplot2grid((2,ncols), (1,ifloor(ncols/3)), colspan=2*ifloor(ncols/3))
    imax21 = plotheatmap(ax3, mat21, vmin, vmax, label2, taus21, Es21, imfont, lagunit, xlabel)

    ax4 = plt.subplot2grid((2,ncols), (0,ncols-1), rowspan=2)  # axis for colorbar

    cbar = fig[:colorbar](imax21, cax=ax4, ticks=[vmin, vmax], format = "%.2f")
    cbar[:ax][:set_yticklabels]([llab, "\$1\$"], fontsize="large")
    cbar[:set_label]("\$max\\ \\rho_{ccm}\$", fontsize="large", labelpad=2)
    cbar[:ax][:get_yaxis]()[:labelpad]=0

    if plot
        plt.show()
    else
        plt.close()
    end
end

function get_totest(val::Int, numtotest::Int, minval::Int, maxval::Int)
    if maxval - minval - 1 < numtotest
        println("Maxval and minval for E_taup curves are inconsistent with number to test")
        println("Fixing maxval")
        maxval = minval + numtotest - 1
    end

    left::Int  = val-iceil(numtotest/2)
    right::Int = val+ifloor(numtotest/2)-1

    if left < minval
        right += (minval-left)
        left   = minval
    end

    if right >  maxval
        left -= right%maxval
        right = maxval
    end

    if right-left+1 != numtotest
        error("start and stop not equal to npred")
    elseif left < minval
        println("left: $left, right: $right, Max: $maxval, min: $minval.")
        error("Left value less than minval")
    end

    return left, right
end

function get_E_taupcurves(source_series::AbstractVector, target_series::AbstractVector, shadowmat_dict::Dict,
                                            distmat_dict::Dict, libsizemax::Int, E::Int, tau_s::Int, tau_p::Int,
                                            E_vals::AbstractVector, tau_p_vals::AbstractVector, npred::Int, pred_start::Int;
                                            nlag=10, b_offset=1, num_E=6, Estart=2, num_tau_p=14, tau_p_start=0,
                                            left_E=false, left_tau_p=false, right_E=false, right_tau_p=false, nboots=0)

    num_E = min(maximum(E_vals) - Estart + 1, num_E)
    num_tau_p = min(maximum(tau_p_vals-tau_p_start+1), num_tau_p)

    if (typeof(left_E) == Bool) && (left_E == false)
        left_E, right_E = get_totest(E, num_E, Estart, maximum(E_vals))
    end

    if (typeof(left_tau_p) == Bool) && (left_tau_p == false)
        left_tau_p, right_tau_p = get_totest(tau_p, num_tau_p, tau_p_start, maximum(tau_p_vals))
    end

    res = fill(NaN, num_E, num_tau_p)

    for (count_E, Eval) in enumerate(left_E:right_E)
        println("\nCalculating CCM max for E: $Eval")
        for (count_tau_p, tau_p_val) in enumerate(left_tau_p:right_tau_p)
            # libsizemin will be ignored. libsizemax will be reduced if it is too large
            println("\ttau_p: $tau_p_val")

            librange, rhos = calcCCM(source_series, target_series,
                                                        shadowmat_dict, distmat_dict,
                                                        Eval, tau_s, tau_p_val;
                                                        libsizemin=libsizemax - nlag, libsizemax=libsizemax,
                                                        npred=npred, pred_start=pred_start,
                                                        quick=true, nboots=nboots)
            rhos = rhos[!isnan(rhos)]

            if length(rhos) > 0
                rho = median(rhos)
            else
                rho = NaN
            end

            res[count_E, count_tau_p] = rho
        end
    end
    return res, left_E:right_E, left_tau_p:right_tau_p
end

function plotheatmap(ax, mat, vmin, vmax, label, taus, Es, imfont, lagunit, xlabel)
    imax = ax[:imshow](mat, vmin=vmin, vmax=vmax, origin="lower")
    ax[:set_title](label, fontsize="medium")

    ax[:xaxis][:set_ticks_position]("bottom")
    ax[:yaxis][:set_ticks_position]("right")

    xticks = 0:(length(taus)-1)
    ax[:set_xticks](xticks[1:2:end]) # take every other tick
    ax[:set_yticks](0:(length(Es)-1))

    ax[:set_xticklabels](
                                    format_ticks(
                                                        (taus * lagunit)[1:2:end]
                                                        ), fontsize=imfont)
    ax[:set_yticklabels](Es, fontsize=imfont)

    ax[:set_xlim]((0, length(taus)-1))
    ax[:set_ylim]((0, length(Es)-1))
    ax[:set_xlabel](xlabel, fontsize=imfont)
    ax[:set_ylabel]("E", fontsize=imfont, labelpad = 12)
    return imax
end

function format_ticks(xticks)

    try
        xticks = convert(Vector{Float64}, xticks) # float range to float vector (otherwise subtraction fails)
    catch
        return xticks # This may fail if you have an int range
    end

    xticks_int = fill(0, length(xticks))
    try
        xticks_int = convert(Vector{Int}, xticks)  # If values are too far from ints, this raises an InexactError
    catch
        return xticks
    end

    variance = var(xticks)
    if all(abs(xticks - xticks_int) .< (variance * 1e-12))  # possibly more stringent than InexactError if values are small
        return xticks_int
    else
        return xticks
    end
end


text:
using Test, WAVI

@testset "Fields" begin
    @info "Testing Fields...."

    @testset "Testing HGrid" begin
    @info "Testing HGrid construction..."
    hgrid = WAVI.HGrid(nxh = 10, nyh = 10, mask = trues(10,10), b = ones(10,10), h = ones(10,10), ηav = ones(10,10))
    @test hgrid isa WAVI.HGrid
    end

    @testset "Testing HGrid errors" begin
        @info "Testing HGrid size input errors..."
        @test_throws DimensionMismatch WAVI.HGrid(nxh = 11, nyh = 10, mask = trues(10,10), b = ones(10,10), h = ones(10,10), ηav = ones(10,10))
        @test_throws DimensionMismatch WAVI.HGrid(nxh = 10, nyh = 11, mask = trues(10,10), b = ones(10,10), h = ones(10,10), ηav = ones(10,10))
        @test_throws DimensionMismatch WAVI.HGrid(nxh = 10, nyh = 10, mask = trues(11,10), b = ones(10,10), h = ones(10,10), ηav = ones(10,10))
        @test_throws DimensionMismatch WAVI.HGrid(nxh = 10, nyh = 10, mask = trues(10,10), b = ones(11,10), h = ones(10,10), ηav = ones(10,10))
        @test_throws DimensionMismatch WAVI.HGrid(nxh = 10, nyh = 10, mask = trues(10,10), b = ones(10,10), h = ones(11,10), ηav = ones(10,10))
        @test_throws DimensionMismatch WAVI.HGrid(nxh = 10, nyh = 10, mask = trues(10,10), b = ones(10,10), h = ones(10,10), ηav = ones(11,10))
    end

    @testset "Testing UGrid" begin
        @info "Testing UGrid construction..."
        ugrid = WAVI.UGrid(nxu = 10, nyu = 10, mask = trues(10,10), levels = 5, dx = 10., dy = 10.)
        @test ugrid isa WAVI.UGrid
        end

    @testset "Testing UGrid errors" begin
        @info "Testing UGrid size input errors..."
        @test_throws DimensionMismatch WAVI.UGrid(nxu = 11, nyu = 10, mask = trues(10,10), levels = 5, dx = 10., dy = 10.)
        @test_throws DimensionMismatch WAVI.UGrid(nxu = 10, nyu = 11, mask = trues(10,10), levels = 5, dx = 10., dy = 10.)
        @test_throws DimensionMismatch WAVI.UGrid(nxu = 10, nyu = 10, mask = trues(11,10), levels = 5, dx = 10., dy = 10.)
    end

    @testset "Testing VGrid" begin
        @info "Testing VGrid construction..."
        vgrid = WAVI.VGrid(nxv = 10, nyv = 10, mask = trues(10,10), levels = 5, dx = 10., dy = 10.)
        @test vgrid isa WAVI.VGrid
        end

    @testset "Testing UGrid errors" begin
        @info "Testing UGrid size input errors..."
        @test_throws DimensionMismatch WAVI.VGrid(nxv = 11, nyv = 10, mask = trues(10,10), levels = 5, dx = 10., dy = 10.)
        @test_throws DimensionMismatch WAVI.VGrid(nxv = 10, nyv = 11, mask = trues(10,10), levels = 5, dx = 10., dy = 10.)
        @test_throws DimensionMismatch WAVI.VGrid(nxv = 10, nyv = 10, mask = trues(11,10), levels = 5, dx = 10., dy = 10.)
    end

    @testset "Testing CGrid" begin
        @info "Testing CGrid construction..."
        cgrid = WAVI.CGrid(nxc = 10, nyc = 10, mask = trues(10,10))
        @test cgrid isa WAVI.CGrid
        end

    @testset "Testing CGrid errors" begin
        @info "Testing CGrid size input errors..."
        @test_throws DimensionMismatch WAVI.CGrid(nxc = 11, nyc = 10, mask = trues(10,10))
        @test_throws DimensionMismatch WAVI.CGrid(nxc = 10, nyc = 11, mask = trues(10,10))
        @test_throws DimensionMismatch WAVI.CGrid(nxc = 10, nyc = 10, mask = trues(11,10))
    end

    @testset "Testing SigmaGrid" begin
        @info "Testing SigmaGrid construction..."
        params = Params()
        grid = Grid()
        sigmagrid = WAVI.SigmaGrid(nxs = 10, nys = 10, nσs = 10, σ = collect(range(0., 1., length = 10)), η = ones(10,10,10), θ = ones(10,10,10), Φ = ones(10,10,10), glen_b = fill(WAVI.glen_b(params.default_temperature,params.default_damage,params.glen_a_ref,params.glen_n, params.glen_a_activation_energy, params.glen_temperature_ref, params.gas_const),10,10,10))
        @test sigmagrid isa WAVI.SigmaGrid
        end

end


text:
using Test,
    Documenter,
    ExperimentalDesign,
    Random,
    Distributions,
    StatsModels

# tests = ["kl_exchange.jl", "doctests.jl"]
tests = ["kl_exchange.jl"]

@testset "ExperimentalDesign" begin
    for test in tests
        include(test)
    end
end


text:
using Manifolds, Manopt, Test, ManifoldsBase

using Dates

struct TestProblem <: Problem{AllocatingEvaluation} end
struct TestOptions <: Options end

@testset "generic Options test" begin
    p = TestProblem()
    o = TestOptions()
    a = ArmijoLinesearch(1.0, ExponentialRetraction(), 0.99, 0.1)
    @test get_last_stepsize(p, o, a) == 1.0
    @test get_initial_stepsize(a) == 1.0
end
@testset "Decresaing Stepsize" begin
    ds = DecreasingStepsize(; length=10.0, factor=1.0, subtrahend=0.0, exponent=1.0)
    @test get_initial_stepsize(ds) == 10.0
    @test ds(TestProblem(), TestOptions(), 1) == 10.0
    @test ds(TestProblem(), TestOptions(), 2) == 5.0
end

@testset "Decorator Options test" begin
    o = TestOptions()
    r = RecordOptions(o, RecordIteration())
    d = DebugOptions(o, DebugIteration())
    dr = DebugOptions(r, DebugIteration())

    @test has_record(o) == has_record(d)
    @test !has_record(o)
    @test has_record(r)
    @test has_record(r) == has_record(dr)

    @test is_options_decorator(r)
    @test is_options_decorator(dr)
    @test is_options_decorator(d)
    @test !is_options_decorator(o)

    @test dispatch_options_decorator(r) === Val(true)
    @test dispatch_options_decorator(dr) === Val(true)
    @test dispatch_options_decorator(d) === Val(true)
    @test dispatch_options_decorator(o) === Val(false)

    @test get_options(r) == o
    @test get_options(dr) == o
    @test get_options(d) == o
    @test get_options(o) == o
    @test get_options(o, Val(false)) == o

    @test Manopt._extract_val(Val(true))
    @test !Manopt._extract_val(Val(false))
end


text:
export initAddNoise
export addNoise

"""
    addNoise(speech::Array, noise::Array, dB::Number)
add `noise` in `speech` according to dB
"""
function addNoise(speech::Array{T}, noise::Array{T}, dB::Number) where T
    Ls = veclen(speech)
    Ln = veclen(noise)
    if Ln<Ls
        # noise is shoter, should be repeated
        noise = repeat(noise, div(Ls,Ln)+1)
        Ln = veclen(noise)
    end
    s = floor(Int,rand()*(Ln-Ls) + 1);   # a random start index
    e = s + Ls - 1;                      # a random end index
    Es = sum(speech.^2);                 # energy of speech
    En = sum(noise[s:e].^2) + eps(T);    # energy of noise
    K  = sqrt(Es/En * 10^(-dB/10));
    speech .+= K * noise[s:e]
    return speech
end

"""
    initAddNoise(path::String, period::Int, dBSpan::NTuple{2,Number}) -> addnoise(speech::Array)
init a adding noise function
+ `path` dir only having noise audios
+ `period` how often the addnoise function would change a background noise audio
+ `dBSpan` e.g. (dBMin, dBMax)
"""
function initAddNoise(path::String, period::Int, dBSpan::NTuple{2,Number})
    counter = 1
    noise = nothing
    FILES = readdir(path)
    dBMin, dBMax = dBSpan
    @assert dBMin <= dBMax
    function addnoise(speech::Array)
        if counter == 1
            file = rand(FILES,1)[1]
            @assert endswith(file, "wav") "$path should only keep *.wav files"
            noise, fs = wavread(joinpath(path, file))
        end # every period we read another noise
        (counter == period) ? (counter=1) : (counter+=1)
        return addNoise(speech, noise, rand()*(dBMax - dBMin) + dBMin)
    end
    return addnoise
end

max_stars_repo_name: sonosole/SpeechAugment

max_issues_repo_name: sonosole/SpeechAugment

max_forks_repo_name: sonosole/SpeechAugment


text:
module Test

export init

"""
# Deprecated
    init(n)

Generate two `Vector{Float64}` with some values for use in
array operation testing.

# Arguments
    `n::Int64`: length of the generated arrays

### Returns
    `(x::Vector{Float64}, y::Vector{Float64})`: Two arrays in a tuple
"""
function init(n::Int64)
    x = Array{Float64, 1}(n)
    y = Array{Float64, 1}(n)
    for i in 1:n
        x[i] = n
        y[i] = n
    end

    return x, y
end

"""
    init!(x, y)

Populate `x, y` with random numbers.

# Arguments
    `x::Vector{T}`: vector to be populated
    `y::Vector{T}`: vector to be populated

"""
function init!{T}(x::Vector{T}, y::Vector{T})
    n::Int64 = length(x)

    srand(17)
    rand!(x)
    rand!(y)
    # for i in 1:n
    #     x[i] = n
    #     y[i] = n
    # end

    return nothing
end

"""
    init!(x, y)

Populate `x, y` with random numbers between 1 and 10.

# Arguments
    `x::Vector{T<:Integer}`: vector to be populated
    `y::Vector{T<:Integer}`: vector to be populated

"""
function init!{T <: Integer}(x::Vector{T}, y::Vector{T})
    n::Int64 = length(x)

    srand(17)
    rand!(x, 1:10)
    rand!(y, 1:10)
    # for i in 1:n
    #     x[i] = n
    #     y[i] = n
    # end

    return nothing
end

"""
# Deprecated
    init()
Generate two `Vector{Float64}` with a size of `1,000,000`.

### returns
    `(x, y)`: Two arrays in a tuple
"""
function init()
    n::Int64 = 1000000
    return init(n)
end

end


text:
# TESTIMAGES ==========================================================================================================

# https://github.com/timholy/TestImages.jl

using TestImages

# List of images in package.
#
readdir(joinpath(homedir(), ".julia/v0.4/TestImages/images/"))

# Load Lena.
#
# The archetypical test image (Lena Söderberg, https://en.wikipedia.org/wiki/Lenna)
#
lena = testimage("lena_color_256.tif")

# IMAGEVIEW ===========================================================================================================

# https://github.com/timholy/ImageView.jl

using ImageView

view(lena)

# IMAGES ==============================================================================================================

# https://github.com/timholy/Images.jl

using Images
using Colors

earth = imread(joinpath(homedir(), ".julia/v0.4/TestImages/images/earth_apollo17.jpg"))
#
#
# Functions for writing images are imwrite() and writemime().

colorspace(earth)
height(earth)
width(earth)

earth_data = data(earth)

earth_rgb = separate(earth)
#
# Channels can be recombined as follows:
#
convert(Image{RGB}, earth_rgb);

earth_array = convert(Array, earth)

# Color conversion.
#
lena_gray = convert(Image{Gray}, lena)

# Kernel functions.
#
gaussian2d(1, [3, 5])
imaverage([3, 3])                                           # Like manual kernel above
#
# Numerous other kernels available.

# Filtering.
#
lena_smooth = imfilter(lena, imaverage([3, 3]))             # Kernel filter via FIR
#
# There is an optional argument which determines how the boundaries are handled. Options are "replicate" (default),
# "circular", "reflect" and "symmetric".
#
lena_very_smooth = imfilter_fft(lena, ones(10, 10) / 100)   # Kernel filter via FFT
#
lena_gauss_smooth = imfilter_gaussian(lena, [1, 2])         # Filter with multivariate Gaussian

# Edge detection. Various algorithms available: sobel, prewitt, ando3 etc.
#
(lena_sobel_x, lena_sobel_y) = imgradients(lena, "sobel");
#
lena_sobel_x                                                # Horizontal gradients
lena_sobel_y                                                # Vertical gradients

magnitude(lena_sobel_x, lena_sobel_y);                      # Magnitude of gradient images
phase(lena_sobel_x, lena_sobel_y);                          # Rotation angle of gradient images
orientation(lena_sobel_x, lena_sobel_y);
#
# imedge() is a convenience function which calculates many of these in a single call.

# Morphology.
#
lena_dilate = dilate(lena);
lena_erode = erode(lena);
#
# opening() and closing() are combinations of dilate() and erode().

# Need to clamp the values in the two Sobel arrays because some of them extend beyond the range [0, 1] which is used
# for image content.
#
imwrite(lena_sobel_x, "lena-sobel-x.png", mapi = mapinfo(Clamp, lena_sobel_x))
imwrite(lena_sobel_y, "lena-sobel-y.png", mapi = mapinfo(Clamp, lena_sobel_y))
imwrite(lena_dilate, "lena-dilate.png")
imwrite(lena_erode, "lena-erode.png")

# IMAGEMAGICK =========================================================================================================

# https://github.com/JuliaIO/ImageMagick.jl

using ImageMagick

# PIECEWISE AFFINE TRANSFORMATIONS ====================================================================================

# https://github.com/dfdx/PiecewiseAffineTransforms.jl

using PiecewiseAffineTransforms