text:
module LeNet5

using Base.Test

# lecun-98 p.41
random_weight(inputs) = (rand() * 2 - 1) * 2.4 / inputs
random_weight(inputs, size...) = (rand(size...) .* 2 - 1) .* 2.4 ./ inputs
valid_weight(w::Number, inputs) = (max = 2.4/inputs; w <= max && w >= -max)
valid_weight(w::Array, inputs) = (max = 2.4/inputs; all((x) -> valid_weight(x, inputs), w) && abs(std(w) - max / 2) <= 0.2 )

# lecun-98 p.8 and 41
const A = 1.7159
const S = 2 / 3
function squash(a)
    A * tanh(S * a)
end
@test_approx_eq_eps 1.0 squash(1) 1e-5
@test_approx_eq_eps -1.0 squash(-1) 1e-5
@test_approx_eq_eps 0.140843782 squash(0.1234) 1e-5

abstract Layer

type C1 <: Layer
    weights
    biases

    output
end
C1() = C1(random_weight(5*5, 6, 5,5), random_weight(5*5, 6), zeros(6,28,28))

function run(input, layer::C1)
    @assert size(input) == (32,32)
    for feature_map in 1:6
        for i in 1:28
            for j in 1:28
                layer.output[feature_map, i, j] = squash(sum(input[i:i+5-1, j:j+5-1] .* layer.weights[feature_map]) + layer.biases[feature_map])
            end
        end
    end
    layer.output
end

parameters(c1::C1) = [c1.weights..., c1.biases...]

function test_c1()
    srand(123)
    input = rand(32,32)
    c1 = C1()
    @test valid_weight(c1.weights, 5*5)
    @test valid_weight(c1.biases, 5*5)
    @test size(parameters(c1)) == (156,)

    output = run(input, c1)
    @test size(output) == (6,28,28)
    @test output[1, 1, 1] == squash(sum(input[1:5, 1:5] .* c1.weights[1]) + c1.biases[1])
    @test output[3, 5, 2] == squash(sum(input[5:9, 2:6] .* c1.weights[3]) + c1.biases[3])
end

type S2 <: Layer
    coefficients
    biases

    output
end
S2() = S2(random_weight(2*2,6), random_weight(2*2,6), zeros(6, 14, 14))

function run(input, layer::S2)
    for f in 1:6, i in 1:14, j in 1:14
        i2 = i*2-1
        j2 = j*2-1
        layer.output[f, i, j] = squash((input[f, i2, j2] + input[f, i2+1, j2] + input[f, i2, j2+1] + input[f, i2+1, j2+1]) * layer.coefficients[f] + layer.biases[f])
    end
    layer.output
end
parameters(s2::S2) = [s2.coefficients..., s2.biases...]

function test_s2()
    srand(123)
    input = rand(6, 28,28)
    s2 = S2()
    @test valid_weight(s2.coefficients, 2*2)
    @test valid_weight(s2.biases, 2*2)
    @test size(parameters(s2)) == (12,)

    output = run(input, s2)
    @test size(output) == (6, 14, 14)
    @test output[1, 1, 1] == squash((input[1,1,1] + input[1,1,2] + input[1,2,1] + input[1,2,2]) * s2.coefficients[1] + s2.biases[1])
    @test output[3, 5, 7] == squash((input[3,9,13] + input[3,9,14] + input[3,10,13] + input[3,10,14]) * s2.coefficients[3] + s2.biases[3])
end

type C3FeatureMap
    input_map_indexes
    weights
    bias
end

type C3 <: Layer
    feature_maps

    output
end
function C3()
    maps = Array(C3FeatureMap, (16))
    index_maps = {
        [0, 1, 2],
        [1, 2, 3],
        [2, 3, 4],
        [3, 4, 5],
        [0, 4, 5],
        [0, 1, 5],
        [0, 1, 2, 3],
        [1, 2, 3, 4],
        [2, 3, 4, 5],
        [0, 3, 4, 5],
        [0, 1, 4, 5],
        [0, 1, 2, 5],
        [0, 1, 3, 4],
        [1, 2, 4, 5],
        [0, 2, 3, 5],
        [0, 1, 2, 3, 4, 5],
    }
    for i in 1:size(index_maps)[1]
        index_map = index_maps[i] + 1
        map = C3FeatureMap(index_map, random_weight(length(index_map)*5*5, length(index_map), 5, 5), random_weight(length(index_map)*5*5))
        maps[i] = map
    end
    C3(maps, zeros(16, 10, 10))
end
function parameters(layer::C3)
    result = {}
    for map in layer.feature_maps
        append!(result, [map.weights..., map.bias])
    end
    result
end
function run(input, layer::C3)
    @assert size(input) == (6, 14, 14)
    fill!(layer.output, 0)
    for feature_map_index in 1:16
        feature_map = layer.feature_maps[feature_map_index]
        for i in 1:10
            for j in 1:10
                for input_index in 1:length(feature_map.input_map_indexes)
                    input_map_index = feature_map.input_map_indexes[input_index]
                    layer.output[feature_map_index, i, j] += sum(input[input_map_index, i:i+5-1, j:j+5-1] .* feature_map.weights[input_index, 1:5, 1:5])
                end
                layer.output[feature_map_index, i, j] = squash(layer.output[feature_map_index, i, j] + feature_map.bias)
            end
        end
    end
    layer.output
end

function test_c3()
    c3 = C3()
    @test size(parameters(c3)) == (1516,)
    @test valid_weight(c3.feature_maps[1].weights, 3*5*5)
    @test valid_weight(c3.feature_maps[8].weights, 4*5*5)

    input = rand(6, 14, 14)
    output = run(input, c3)

    @test size(output) == (16, 10, 10)
    @test output[1, 1, 1] == squash(
        sum(input[1, 1:5, 1:5] .* c3.feature_maps[1].weights[1, 1:5, 1:5]) +
        sum(input[2, 1:5, 1:5] .* c3.feature_maps[1].weights[2, 1:5, 1:5]) +
        sum(input[3, 1:5, 1:5] .* c3.feature_maps[1].weights[3, 1:5, 1:5]) +
        c3.feature_maps[1].bias)
    @test output[8, 3, 5] == squash(
        sum(input[2, 3:7, 5:9] .* c3.feature_maps[8].weights[1, 1:5, 1:5]) +
        sum(input[3, 3:7, 5:9] .* c3.feature_maps[8].weights[2, 1:5, 1:5]) +
        sum(input[4, 3:7, 5:9] .* c3.feature_maps[8].weights[3, 1:5, 1:5]) +
        sum(input[5, 3:7, 5:9] .* c3.feature_maps[8].weights[4, 1:5, 1:5]) +
        c3.feature_maps[8].bias)
end

type S4 <: Layer
    coefficients
    biases

    output
end
S4() = S4(random_weight(2*2,16), random_weight(2*2,16), zeros(16, 5, 5))

function run(input, layer::S4)
    for f in 1:16, i in 1:5, j in 1:5
        i2 = i*2-1
        j2 = j*2-1
        layer.output[f, i, j] = squash((input[f, i2, j2] + input[f, i2+1, j2] + input[f, i2, j2+1] + input[f, i2+1, j2+1]) * layer.coefficients[f] + layer.biases[f])
    end
    layer.output
end
parameters(s4::S4) = [s4.coefficients..., s4.biases...]

function test_s4()
    srand(123)
    input = rand(16,10,10)
    layer = S4()
    @test valid_weight(layer.coefficients, 2*2)
    @test valid_weight(layer.biases, 2*2)
    @test size(parameters(layer)) == (32,)

    output = run(input, layer)
    @test size(output) == (16, 5, 5)
    @test output[1, 1, 1] == squash((input[1,1,1] + input[1,1,2] + input[1,2,1] + input[1,2,2]) * layer.coefficients[1] + layer.biases[1])
    @test output[3, 5, 2] == squash((input[3,9,3] + input[3,9,4] + input[3,10,3] + input[3,10,4]) * layer.coefficients[3] + layer.biases[3])
end

type C5FeatureMap
    weights
    bias
end

type C5 <: Layer
    feature_maps

    output
end
function C5()
    maps = Array(C5FeatureMap, (120,))
    for i in 1:120
        map = C5FeatureMap(random_weight(16*5*5, 16, 5, 5), random_weight(16*5*5))
        maps[i] = map
    end
    C5(maps, zeros(120, 1, 1))
end
function parameters(layer::C5)
    result = {}
    for map in layer.feature_maps
        append!(result, [map.weights..., map.bias])
    end
    result
end

function run(input, layer::C5)
    @assert size(input) == (16, 5, 5)
    fill!(layer.output, 0)
    for feature_map_index in 1:120
        feature_map = layer.feature_maps[feature_map_index]
        for i in 1:1
            for j in 1:1
                for input_index in 1:16
                    layer.output[feature_map_index, i, j] += sum(input[input_index, i:i+5-1, j:j+5-1] .* feature_map.weights[input_index, 1:5, 1:5])
                end
                layer.output[feature_map_index, i, j] = squash(layer.output[feature_map_index, i, j] + feature_map.bias)
            end
        end
    end
    layer.output
end

function test_c5()
    srand(1223)
    input = rand(16, 5, 5)
    layer = C5()
    @test length(parameters(layer)) == 48120
    @test valid_weight(layer.feature_maps[1].weights, 16*5*5)
    @test valid_weight(layer.feature_maps[75].weights, 16*5*5)

    output = run(input, layer)
    @test size(output) == (120,1,1)
    @test_approx_eq_eps output[1, 1, 1] squash(sum(input .* layer.feature_maps[1].weights) + layer.feature_maps[1].bias) 1e-10
    @test_approx_eq_eps output[75, 1, 1] squash(sum(input .* layer.feature_maps[75].weights) + layer.feature_maps[75].bias) 1e-10
end

type F6 <: Layer
    weights
    biases

    output
end
F6() = F6(random_weight(120, 84, 120), random_weight(120, 84), zeros(84))

function parameters(layer::F6)
    [layer.weights..., layer.biases...]
end

function run(input, layer::F6)
    for i in 1:84
        weighted_sum = sum(reshape(input,120) .* reshape(layer.weights[i,:],120)) + layer.biases[i]
        layer.output[i] = squash(weighted_sum)
    end
    layer.output
end

function test_f6()
    srand(123)
    input = rand(120,1,1)
    layer = F6()
    @test valid_weight(layer.weights, 120)
    @test valid_weight(layer.biases, 120)
    @test length(parameters(layer)) == 10164

    output = run(input, layer)
    @test size(output) == (84,)
    expected = 0
    for i in 1:120
        expected += input[i,1,1] * layer.weights[1,i]
    end
    expected += layer.biases[1]
    expected = squash(expected)
    @test output[1] == expected
    @test output[79] == squash(sum(reshape(input,120) .* reshape(layer.weights[79,:],120)) + layer.biases[79])
end

type Output <: Layer
    weights

    output
end

using Images

function Output()
    weights = zeros(96, 12, 7)

    for i in 1:96
        filename = "output_rbfs/$(lpad(i,3,'0')).ppm"
        image = imread(filename)
        array = convert(Array, image)[:,:,1]'
        for y in 1:size(array)[1]
            for x in 1:size(array)[2]
                weights[i,y,x] = array[y, x]  == 0 ? 1 : -1
            end
        end
    end

    Output(weights, zeros(96))
end

function run(input, layer::Output)
    fill!(layer.output, 0)
    for i in 1:96
        for j in 1:84
            layer.output[i] += ((input[j] - layer.weights[i,j]) .^ 2)
        end
    end
    layer.output
end

function test_output()
    srand(123)
    input = rand(84)
    layer = Output()
    @test reshape(layer.weights[ 2, 1, :], 7) == [-1, -1, -1,  1,  1, -1, -1]
    @test reshape(layer.weights[23, 9, :], 7) == [ 1,  1,  1, -1, -1,  1,  1]

    output = run(input, layer)
    @test size(output) == (96,)
    @test_approx_eq_eps output[ 1] sum((input - reshape(layer.weights[ 1,:,:], 84)).^2) 1e-10
    @test_approx_eq_eps output[21] sum((input - reshape(layer.weights[21,:,:], 84)).^2) 1e-10
end

function test_precise_output()
    srand(123)
    input = ones(84)
    layer = Output()

    output = run(input, layer)
    @test size(output) == (96,)
    @test_approx_eq_eps output[ 2] ((7*12-18)*4 + (18*0)) 1e-10
end

type NeuralNetwork <: Layer
    c1
    s2
    c3
    s4
    c5
    f6
    output
end
NeuralNetwork() = NeuralNetwork(C1(), S2(), C3(), S4(), C5(), F6(), Output())

function run(input, network::NeuralNetwork)
    run(input, network.c1)
    run(network.c1.output, network.s2)
    run(network.s2.output, network.c3)
    run(network.c3.output, network.s4)
    run(network.s4.output, network.c5)
    run(network.c5.output, network.f6)
    run(network.f6.output, network.output)
    network.output.output
end

import Base.|>
(|>)(input, layer::Layer) = run(input, layer)

function test_lenet5()
    srand(123)
    input = rand(32, 32)
    lenet5 = NeuralNetwork()

    output = run(input, lenet5)

    expected_output = run(input, lenet5.c1)
    expected_output = run(expected_output, lenet5.s2)
    expected_output = run(expected_output, lenet5.c3)
    expected_output = run(expected_output, lenet5.s4)
    expected_output = run(expected_output, lenet5.c5)
    expected_output = run(expected_output, lenet5.f6)
    expected_output = run(expected_output, lenet5.output)
    @test output == expected_output
end

exp_minus_J = exp(-5)

function maximum_a_posteriori(output)
    exp_output = map((x) -> exp(-x), output)
    log(exp_minus_J + sum(exp_output))
end

function single_loss(output, desired_class)
  error = output[desired_class] + maximum_a_posteriori(output)
  error
end

function loss(outputs, desired_classes)
    training_samples = length(desired_classes)
    @assert size(outputs)[1] == training_samples
    errors = zeros(training_samples)
    for i in 1:training_samples
        errors[i] += single_loss(reshape(outputs[i, :], size(outputs)[2]), desired_classes[i])
    end
    mean(errors)
end

function test_loss()
    training_samples = 2
    network_outputs = ones(training_samples, 96)
    desired_output_classes = [6, 87]
    network_outputs[1, 6] = 0.1
    network_outputs[1, 7] = 0.5
    network_outputs[2, 87] = 0.5
    network_outputs[2, 12] = 0.2

    expected_error = 0.1 + log(exp_minus_J + exp(-0.1) + exp(-0.5) + 94 * exp(-1))
    expected_error += 0.5 + log(exp_minus_J + exp(-0.5) + exp(-0.2) + 94 * exp(-1))
    expected_error /= 2

    @test_approx_eq_eps loss(network_outputs, desired_output_classes) expected_error 1e-10
end

function test_loss_generic()
    srand(123)
    training_samples = 6
    network_outputs = rand(training_samples, 96)
    desired_output_classes = rand(1:96, training_samples)

    expected_error = 0
    for p in 1:training_samples
        expected_error += network_outputs[p, desired_output_classes[p]]
        expected_error += log(exp_minus_J + sum([exp(-output) for output in network_outputs[p,:]]))
    end
    expected_error /= training_samples
    @test loss(network_outputs, desired_output_classes) == expected_error
end

function test_all()
    test_c1()
    test_s2()
    test_c3()
    test_s4()
    test_c5()
    test_f6()
    test_output()
    test_lenet5()
    test_loss()
    test_loss_generic()
end

end


text:
include("data.jl")
include("technical_indicators.jl")

assets = []

push!(assets, "USDJPY")
USDJPY = market("USDJPY")
push!(assets, "USDCAD")
USDCAD = market("USDCAD")

x = []
i = 1
while i < 1000
   global i
   value = rand()
   push!(x, value)
   i += 1
end

create_historic(USDJPY, x, 500, "Daily")

for i in 501:length(x)
   create_event(USDJPY, i, "Daily")
end

ext: jl

max_stars_repo_name: SRelieved/JuliaBacktesting

max_issues_repo_name: MartinGenesse/JuliaBacktesting

max_forks_repo_name: MartinGenesse/JuliaBacktesting


text:
# squared L2 norm (times a constant, or weighted)

export SqrNormL2

"""
**Squared Euclidean norm (weighted)**

    SqrNormL2(λ=1)

With a nonnegative scalar `λ`, returns the function
```math
f(x) = \\tfrac{λ}{2}\\|x\\|^2.
```
With a nonnegative array `λ`, returns the function
```math
f(x) = \\tfrac{1}{2}∑_i λ_i x_i^2.
```
"""
struct SqrNormL2{T <: Union{Real, AbstractArray}} <: ProximableFunction
    lambda::T
    function SqrNormL2{T}(lambda::T) where {T <: Union{Real,AbstractArray}}
        if any(lambda .< 0)
            error("coefficients in λ must be nonnegative")
        else
            new(lambda)
        end
    end
end

is_convex(f::SqrNormL2) = true
is_smooth(f::SqrNormL2) = true
is_separable(f::SqrNormL2) = true
is_quadratic(f::SqrNormL2) = true
is_strongly_convex(f::SqrNormL2) = all(f.lambda .> 0)

SqrNormL2(lambda::T=1) where {T <: Real} = SqrNormL2{T}(lambda)

SqrNormL2(lambda::T) where {T <: AbstractArray} = SqrNormL2{T}(lambda)

function (f::SqrNormL2{S})(x::AbstractArray{T}) where {S <: Real, R <: Real, T <: RealOrComplex{R}}
    return f.lambda / R(2) * norm(x)^2
end

function (f::SqrNormL2{S})(x::AbstractArray{T}) where {S <: AbstractArray, R <: Real, T <: RealOrComplex{R}}
    sqnorm = R(0)
    for k in eachindex(x)
        sqnorm += f.lambda[k] * abs2(x[k])
    end
    return sqnorm / R(2)
end

function gradient!(y::AbstractArray{T}, f::SqrNormL2{S}, x::AbstractArray{T}) where {S <: Real, R <: Real, T <: RealOrComplex{R}}
    sqnx = R(0)
    for k in eachindex(x)
        y[k] = f.lambda * x[k]
        sqnx += abs2(x[k])
    end
    return f.lambda / R(2) * sqnx
end

function gradient!(y::AbstractArray{T}, f::SqrNormL2{S}, x::AbstractArray{T}) where {S <: AbstractArray, R <: Real, T <: RealOrComplex{R}}
    sqnx = R(0)
    for k in eachindex(x)
        y[k] = f.lambda[k] * x[k]
        sqnx += f.lambda[k] * abs2(x[k])
    end
    return sqnx / R(2)
end

function prox!(y::AbstractArray{T}, f::SqrNormL2{S}, x::AbstractArray{T}, gamma::R=R(1)) where {S <: Real, R <: Real, T <: RealOrComplex{R}}
    gl = gamma * f.lambda
    sqny = R(0)
    for k in eachindex(x)
        y[k] = x[k] / (1 + gl)
        sqny += abs2(y[k])
    end
    return f.lambda / R(2) * sqny
end

function prox!(y::AbstractArray{T}, f::SqrNormL2{S}, x::AbstractArray{T}, gamma::R=R(1)) where {S <: AbstractArray, R <: Real, T <: RealOrComplex{R}}
    wsqny = R(0)
    for k in eachindex(x)
        y[k] = x[k] / (1 + gamma * f.lambda[k])
        wsqny += f.lambda[k] * abs2(y[k])
    end
    return wsqny / R(2)
end

function prox!(y::AbstractArray{T}, f::SqrNormL2{S}, x::AbstractArray{T}, gamma::AbstractArray{R}) where {S <: Real, R <: Real, T <: RealOrComplex{R}}
    sqny = R(0)
    for k in eachindex(x)
        y[k] = x[k] / (1 + gamma[k] * f.lambda)
        sqny += abs2(y[k])
    end
    return f.lambda / R(2) * sqny
end

function prox!(y::AbstractArray{T}, f::SqrNormL2{S}, x::AbstractArray{T}, gamma::AbstractArray{R}) where {S <: AbstractArray, R <: Real, T <: RealOrComplex{R}}
    wsqny = R(0)
    for k in eachindex(x)
        y[k] = x[k] / (1 + gamma[k] * f.lambda[k])
        wsqny += f.lambda[k] * abs2(y[k])
    end
    return wsqny / R(2)
end

fun_name(f::SqrNormL2) = "weighted squared Euclidean norm"
fun_dom(f::SqrNormL2) = "AbstractArray{Real}, AbstractArray{Complex}"
fun_expr(f::SqrNormL2{T}) where {T <: Real} = "x ↦ (λ/2)||x||^2"
fun_expr(f::SqrNormL2{T}) where {T <: AbstractArray} = "x ↦ (1/2)sum( λ_i (x_i)^2 )"
fun_params(f::SqrNormL2{T}) where {T <: Real} = "λ = $(f.lambda)"
fun_params(f::SqrNormL2{T}) where {T <: AbstractArray} = string("λ = ", typeof(f.lambda), " of size ", size(f.lambda))

function prox_naive(f::SqrNormL2, x::AbstractArray{T}, gamma=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    y = x./(R(1) .+ f.lambda .* gamma)
    return y, real(dot(f.lambda .* y, y)) / R(2)
end


text:
# test for reso.jl
using Test
#include("unify.jl")
#include("reso.jl")

#==
idea:
unify([:x,:y], "P(x,y)", "P(x,f(x))")
unify([:x,:y], :(P(x,y)), :(P(x,f(x))))

:(P(x,y)) == Meta.Meta.parse("P(x,y)")
youcan readline(from console) and Meta.parse it to get the Expr
==#

@testset "genvars" begin
 @test genvars([]) == []
 gv = genvars([:x,:y])
 @test length(gv) == 2
 @test gv != [:x,:y]

end

@testset "equalclause" begin
 @test equalclause(([:x], [Meta.parse("+Q(x,b)")]), ([:x],[Meta.parse("+Q(x,b)")])) == true
 @test equalclause(([:x], [Meta.parse("+Q(x,c)")]), ([:x],[Meta.parse("+Q(x,b)")])) == false
 @test equalclause(([:y], [Meta.parse("+Q(y,b)")]), ([:x],[Meta.parse("+Q(x,b)")])) == true
 @test equalclause(([:y], [Meta.parse("+Q(y,b)"),Meta.parse("+P(c)")]), ([:x],[Meta.parse("+Q(x,b)"),Meta.parse("-P(a)")])) == false
 @test equalclause(([:x,:y], [Meta.parse("+Q(x,y)"),Meta.parse("+P(a)")]), ([:w,:y],[Meta.parse("+Q(x,y)"),Meta.parse("-P(b)")])) == false
 @test equalclause(([:x,:y], [Meta.parse("+Q(x,y)"),Meta.parse("+P(x)")]), ([:w,:y],[Meta.parse("+Q(x,y)"),Meta.parse("-P(w)")])) == false
 @test equalclause(([:y], [Meta.parse("+Q(y,b)")]), ([:y],[Meta.parse("+Q(x,b)")])) == false
 @test equalclause(([:x], [Meta.parse("+Q(y,x,b)")]), ([:y],[Meta.parse("+Q(y,x,b)")])) == false
end

@testset "pure resolution" begin
 @test resolution(([],Meta.parse("[+P(a)]").args),1, ([],Meta.parse("[+P(a)]").args),1) == :NAP
 @test resolution(([],Meta.parse("[+P(a)]").args),1, ([],Meta.parse("[-P(a)]").args),1) == (EmptyCForm,[])
 @test resolution(([],Meta.parse("[-P(a)]").args),1, ([],Meta.parse("[+P(a)]").args),1) == (EmptyCForm,[])
 @test resolution(([],Meta.parse("[+P(a)]").args),1, ([:x],Meta.parse("[-P(x)]").args),1) == (EmptyCForm,[])
 @test resolution(([:x],Meta.parse("[+P(x)]").args),1, ([:x],Meta.parse("[-P(x)]").args),1) == (EmptyCForm,[])

 @test resolution(([],Meta.parse("[+P(a)]").args),1, ([],Meta.parse("[-P(b)]").args),1) == :NAP
 @test resolution(([],Meta.parse("[+P(a)]").args),1, ([:x,:y],Meta.parse("[-P(b),+Q(x,y),+R(y)]").args),1) == :NAP
 r1 = resolution(([:x],Meta.parse("[+(P(x))]").args),1, ([:x,:y],Meta.parse("[-P(x),+Q(x, y),+R(y)]").args),1)
 @test equalclause(r1[1],([:x,:y],Meta.parse("[+Q(x,y),+R(y)]").args)) == true

 r2 = resolution(([],Meta.parse("[+(P(a))]").args), 1, ([:x,:y],Meta.parse("[-P(x),+Q(x, y),+R(y)]").args),1)
 @test equalclause(r2[1], ([:y],Meta.parse("[+Q(a,y),+R(y)]").args))

 r3 = resolution(([],Meta.parse("[-Q(a,b),+P(a)]").args),1, ([:x,:y],Meta.parse("[-P(x),+Q(x, y),+R(y)]").args),2)
 @test equalclause(r3[1], ([],Meta.parse("[+P(a),-P(a),+R(b)]").args))

 r4 = resolution(([],Meta.parse("[-Q(a,b),+P(a)]").args),2, ([:x,:y],Meta.parse("[+Q(x, y),+R(y),-P(x)]").args),3)
 @test equalclause(r4[1], ([:y],Meta.parse("[-Q(a,b),+Q(a,y),+R(y)]").args))
end

@testset "real resolution" begin
 v,r = resolution([:y], Meta.parse("[-Q(a,y),+P(a)]").args, 2, [:x,:y], Meta.parse("[+Q(x, y),+R(y),-P(x)]").args,3)
 @test rename(v[1],v[2],[:x,:y,:z]) == Meta.parse("[-Q(a,x),+Q(a,z),+R(z)]").args
 v,r = resolution([], Meta.parse("[-Q(a,b),+P(a)]").args, 2, [:x,:y], Meta.parse("[+Q(x, y),+R(y),-P(x)]").args,3)
 @test rename(v[1],v[2],[:x,:y]) ==  Meta.parse("[-Q(a,b),+Q(a,y),+R(y)]").args
end

@testset "reduction" begin
 r = reduction(([:y], Meta.parse("[-Q(a,y),+P(a)]").args), 1)
 @test equalclause(r[1], ([:y],Meta.parse("[-Q(a,y),+P(a)]").args)) == true

 r = reduction(([:y], Meta.parse("[-Q(a,y),+P(a),+P(a)]").args), 2)
 @test equalclause(r[1], ([:y],Meta.parse("[-Q(a,y),+P(a)]").args)) == true

 r = reduction(([:x,:y],Meta.parse("[-Q(a,y),+P(x),+P(a)]").args), 2)
 @test equalclause(r[1], ([:y], Meta.parse("[-Q(a,y),+P(a)]").args)) == true

 r = reduction(([:y], Meta.parse("[-Q(a,y),+P(b),+P(a)]").args), 2)
 @test equalclause(r[1], ([:y],Meta.parse("[-Q(a,y),+P(b),+P(a)]").args)) == true

 r = reduction(([:x], Meta.parse("[-Q(a,x),+P(x),+P(a)]").args), 2)
 @test equalclause(r[1], ([],Meta.parse("[-Q(a,a),+P(a)]").args)) == true

 r = reduction(([:x], Meta.parse("[-Q(a,x),+P(x),+P(a),+R(x)]").args), 2)
 @test equalclause(r[1],([],Meta.parse("[-Q(a,a),+P(a),+R(a)]").args))== true

 r = reduction(([:x], Meta.parse("[-Q(a,x),+P(x),+P(a),+R(x)]").args), 2)
 @test equalclause(r[1],([], Meta.parse("[-Q(a,a),+P(a),+R(a)]").args)) == true

 r = reduction(([:x], Meta.parse("[-Q(a,x),+P(x),+P(f(x)),+R(x)]").args), 2)
 @test equalclause(r[1], ([:x],Meta.parse("[-Q(a,x),+P(x),+P(f(x)),+R(x)]").args)) == true
end

@testset "satisfiable" begin
 @test satisfiable([:x], Meta.parse("[-Q(a,x),+P(x),-P(y),+R(x)]").args) == true
 @test satisfiable([], Meta.parse("[-Q(a,x),+P(x),-P(y),+R(x)]").args) == false

 @test satisfiable([:x,:y], Meta.parse("[-Q(a,x),+P(x),-P(y),+R(x)]").args) == true
 @test satisfiable([:x,:y], Meta.parse("[-Q(a,x),+P(x,a),-P(y,a),+R(x)]").args) == true
 @test satisfiable([:x,:y], Meta.parse("[-Q(a,x),-P(x,a),-P(y,a),+R(x)]").args) == false
end

@testset "readablevars" begin
 @test length(readablevars(1))==1
 @test length(readablevars(10))==10
 newv = readablevars(30)
 @test length(intersect(newv,newv)) == 30
end

@testset "reducesub" begin
 @test reducesub([:x,:y], [:x], [:a,:b]) == [:a]
 @test reducesub([:x,:y,:z,:w], [:x,:w], [:a,:b,:c,[:f :x]]) == [:a, [:f :x]]
end

@testset "fitvars" begin
 @test fitvars([:x,:y], Meta.parse("[+Q(x,b)]").args, [:a,:b]) == (([:x],Meta.parse("[+Q(x,b)]").args), [:a])
 @test fitvars([:x,:y], Meta.parse("[+Q(a,b)]").args, [:a,:b]) == (([],Meta.parse("[+Q(a,b)]").args), [])
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
function test_closure_instantiation(T, closurename)
    closure = getproperty(TurbulenceClosures, closurename)(T)
    return eltype(closure) == T
end

function test_function_interpolation(T=Float64)
    grid = RegularCartesianGrid(T, (3, 3, 3), (3, 3, 3))
    ϕ = rand(T, 3, 3, 3)
    ϕ² = ϕ.^2

    ▶x_ϕ_f = (ϕ²[2, 2, 2] + ϕ²[1, 2, 2]) / 2
    ▶x_ϕ_c = (ϕ²[3, 2, 2] + ϕ²[2, 2, 2]) / 2

    ▶y_ϕ_f = (ϕ²[2, 2, 2] + ϕ²[2, 1, 2]) / 2
    ▶y_ϕ_c = (ϕ²[2, 3, 2] + ϕ²[2, 2, 2]) / 2

    ▶z_ϕ_f = (ϕ²[2, 2, 2] + ϕ²[2, 2, 1]) / 2
    ▶z_ϕ_c = (ϕ²[2, 2, 3] + ϕ²[2, 2, 2]) / 2

    f(i, j, k, grid, ϕ) = ϕ[i, j, k]^2

    return (
        ▶x_caa(2, 2, 2, grid, f, ϕ) == ▶x_ϕ_c &&
        ▶x_faa(2, 2, 2, grid, f, ϕ) == ▶x_ϕ_f &&

        ▶y_aca(2, 2, 2, grid, f, ϕ) == ▶y_ϕ_c &&
        ▶y_afa(2, 2, 2, grid, f, ϕ) == ▶y_ϕ_f &&

        ▶z_aac(2, 2, 2, grid, f, ϕ) == ▶z_ϕ_c &&
        ▶z_aaf(2, 2, 2, grid, f, ϕ) == ▶z_ϕ_f
        )
end

function test_function_differentiation(T=Float64)
    grid = RegularCartesianGrid(T, (3, 3, 3), (3, 3, 3))
    ϕ = rand(T, 3, 3, 3)
    ϕ² = ϕ.^2

    ∂x_ϕ_f = ϕ²[2, 2, 2] - ϕ²[1, 2, 2]
    ∂x_ϕ_c = ϕ²[3, 2, 2] - ϕ²[2, 2, 2]

    ∂y_ϕ_f = ϕ²[2, 2, 2] - ϕ²[2, 1, 2]
    ∂y_ϕ_c = ϕ²[2, 3, 2] - ϕ²[2, 2, 2]

    # Note reverse indexing here!
    ∂z_ϕ_f = ϕ²[2, 2, 1] - ϕ²[2, 2, 2]
    ∂z_ϕ_c = ϕ²[2, 2, 2] - ϕ²[2, 2, 3]

    f(i, j, k, grid, ϕ) = ϕ[i, j, k]^2

    return (
        ∂x_caa(2, 2, 2, grid, f, ϕ) == ∂x_ϕ_c &&
        ∂x_faa(2, 2, 2, grid, f, ϕ) == ∂x_ϕ_f &&

        ∂y_aca(2, 2, 2, grid, f, ϕ) == ∂y_ϕ_c &&
        ∂y_afa(2, 2, 2, grid, f, ϕ) == ∂y_ϕ_f &&

        ∂z_aac(2, 2, 2, grid, f, ϕ) == ∂z_ϕ_c &&
        ∂z_aaf(2, 2, 2, grid, f, ϕ) == ∂z_ϕ_f
        )
end

function test_constant_isotropic_diffusivity_basic(T=Float64; ν=T(0.3), κ=T(0.7))
    closure = ConstantIsotropicDiffusivity(T, κ=κ, ν=ν)
    return (                    closure.ν == ν &&
                                closure.κ == κ &&
        κ_ccc(nothing, nothing, nothing, nothing, closure) == κ
    )
end

function test_constant_isotropic_diffusivity_fluxdiv(TF=Float64; ν=TF(0.3), κ=TF(0.7))
    closure = ConstantIsotropicDiffusivity(TF, κ=κ, ν=ν)
    grid = RegularCartesianGrid(TF, (3, 1, 1), (3, 1, 1))
    eos = LinearEquationOfState()
    g = 1.0

    u = zeros(TF, 3, 1, 1); v = zeros(TF, 3, 1, 1); w = zeros(TF, 3, 1, 1)
    T = zeros(TF, 3, 1, 1); S = zeros(TF, 3, 1, 1)

    u[:, 1, 1] .= [0, -1, 0]
    v[:, 1, 1] .= [0, -2, 0]
    w[:, 1, 1] .= [0, -3, 0]
    T[:, 1, 1] .= [0, -1, 0]

    return (∇_κ_∇ϕ(2, 1, 1, grid, T, closure, eos, g, u, v, w, T, S) == 2κ &&
            ∂ⱼ_2ν_Σ₁ⱼ(2, 1, 1, grid, closure, eos, g, u, v, w, T, S) == 2ν &&
            ∂ⱼ_2ν_Σ₂ⱼ(2, 1, 1, grid, closure, eos, g, u, v, w, T, S) == 4ν &&
            ∂ⱼ_2ν_Σ₃ⱼ(2, 1, 1, grid, closure, eos, g, u, v, w, T, S) == 6ν
            )
end

function test_anisotropic_diffusivity_fluxdiv(TF=Float64; νh=TF(0.3), κh=TF(0.7), νv=TF(0.1), κv=TF(0.5))
    closure = ConstantAnisotropicDiffusivity(TF, κh=κh, νh=νh, κv=κv, νv=νv)
    grid = RegularCartesianGrid(TF, (3, 1, 3), (3, 1, 3))
    eos = LinearEquationOfState()
    g = 1.0

    u = zeros(TF, 3, 1, 3); v = zeros(TF, 3, 1, 3); w = zeros(TF, 3, 1, 3)
    T = zeros(TF, 3, 1, 3); S = zeros(TF, 3, 1, 3)

    u[:, 1, 1] .= [0,  1, 0]
    u[:, 1, 2] .= [0, -1, 0]
    u[:, 1, 3] .= [0,  1, 0]

    v[:, 1, 1] .= [0,  1, 0]
    v[:, 1, 2] .= [0, -2, 0]
    v[:, 1, 3] .= [0,  1, 0]
    w[:, 1, 1] .= [0,  1, 0]
    w[:, 1, 2] .= [0, -3, 0]
    w[:, 1, 3] .= [0,  1, 0]

    T[:, 1, 1] .= [0,  1, 0]
    T[:, 1, 2] .= [0, -4, 0]
    T[:, 1, 3] .= [0,  1, 0]

    return (∇_κ_∇ϕ(2, 1, 2, grid, T, closure, eos, g, u, v, w, T, S) == 8κh + 10κv &&
            ∂ⱼ_2ν_Σ₁ⱼ(2, 1, 2, grid, closure, eos, g, u, v, w, T, S) == 2νh + 4νv &&
            ∂ⱼ_2ν_Σ₂ⱼ(2, 1, 2, grid, closure, eos, g, u, v, w, T, S) == 4νh + 6νv &&
            ∂ⱼ_2ν_Σ₃ⱼ(2, 1, 2, grid, closure, eos, g, u, v, w, T, S) == 6νh + 8νv
            )
end

function test_smag_divflux_finiteness(TF=Float64)
    closure = ConstantSmagorinsky(TF)
    grid = RegularCartesianGrid(TF, (3, 3, 3), (3, 3, 3))
    eos = LinearEquationOfState()
    g = 1.0
    u, v, w = rand(TF, size(grid)...), rand(TF, size(grid)...), rand(TF, size(grid)...)
    T, S = rand(TF, size(grid)...), rand(TF, size(grid)...)

    return (
        isfinite(∇_κ_∇ϕ(2, 1, 2, grid, T, closure, eos, g, u, v, w, T, S)) &&
        isfinite(∂ⱼ_2ν_Σ₁ⱼ(2, 1, 2, grid, closure, eos, g, u, v, w, T, S)) &&
        isfinite(∂ⱼ_2ν_Σ₂ⱼ(2, 1, 2, grid, closure, eos, g, u, v, w, T, S)) &&
        isfinite(∂ⱼ_2ν_Σ₃ⱼ(2, 1, 2, grid, closure, eos, g, u, v, w, T, S))
        )
end


text:
point_count = 50
flow = gen_tiled_flow((200, 200), 30, 40)
inds = gen_rand_points(flow, point_count, "Semigridded")

### scatter

points = reshape([ind[i] for i=1:length(inds[1]) for ind=inds], length(inds), length(inds[1]))'
samples = [flow[I] for I in inds]

meshgrid(x, y) = [repeat(x, outer=length(y)) repeat(y, inner=length(x))]
meshgrid(x::Real, y::Real) = [repeat(1:x, outer=y) repeat(1:y, inner=x)]

function test_scatter(flow_size, points, samples; method)
    gridPoints = meshgrid(flow_size...)'
    itp = interpolate(method, points, samples);
    interpolated = ScatteredInterpolation.evaluate(itp, gridPoints)
    gridded = reshape(interpolated, flow_size)
    return gridded
end

bench = @benchmark test_scatter(size(flow), points, samples; method = Polyharmonic(2))

### Dier

x = [ind[1] for ind in inds]
y = [ind[2] for ind in inds]
z1 = [real(flow[ind]) for ind in inds]
z2 = [imag(flow[ind]) for ind in inds]

function test_dier(flow_size, x, y, z1, z2)
    storage = 1000
    spl1 = Spline2D(x, y, z1; w=ones(length(x)), kx=1, ky=1, s=storage)
    real1 = evalgrid(spl1, 1:flow_size[1], 1:flow_size[2])

    spl2 = Spline2D(x, y, z2; w=ones(length(x)), kx=1, ky=1, s=storage)
    imag1 = evalgrid(spl2, 1:flow_size[1], 1:flow_size[2])

    gridded = real1 .+ imag1 .* im
    return gridded
end

@benchmark test_dier(size(flow), x, y, z1, z2)

### generate indices

"""
    gen_rand_points(flow, point_count, mode; max_skew_const=4)

Generate `point_count` random indices over `flow` depending on selection `mode`.
"""
function gen_rand_points(flow, point_count, mode; max_skew_const=4)
    if mode == "Random_points"
        ret_inds = rand(CartesianIndices(flow), point_count)
    elseif mode == "Gridded"
        points_in_dim = round(Int64, sqrt(point_count))
        point_spacing = floor.(size(flow) ./ points_in_dim)

        ret_inds = Array{CartesianIndex{2}, 1}(undef, 0)
        for ind in CartesianIndices(flow)
            if all((Tuple(ind)) .% point_spacing .== 0)
                add_ind = CartesianIndex(Tuple(ind) .- floor.(Int64, point_spacing ./ 2))
                push!(ret_inds, add_ind)
            end
        end
    elseif mode == "Semigridded"
        points_in_dim = round(Int64, sqrt(point_count))
        point_spacing = floor.(size(flow) ./ points_in_dim)

        ret_inds = Array{CartesianIndex{2}, 1}(undef, 0)
        for ind in CartesianIndices(flow)
            if all((Tuple(ind)) .% point_spacing .== 0)
                skew = (rand((-floor(Int64, point_spacing[1]/max_skew_const):floor(Int64, point_spacing[1]/max_skew_const))),
                        rand((-floor(Int64, point_spacing[1]/max_skew_const):floor(Int64, point_spacing[2]/max_skew_const))))
                add_ind = CartesianIndex(Tuple(ind) .- floor.(Int64, point_spacing ./ 2) .+ skew)
                @assert add_ind in CartesianIndices(flow)

                push!(ret_inds, add_ind)
            end
        end
    end
    return ret_inds
end

## add to df
index = 1

push!(df, Dict(:Index => index
               ,:Method => method
               ,:Mode => point_mode
               ,:Point_count => point_count
               ,:Median_speed => median(bench).time
               ,:Magnitude_of_MSE => vec_len(mse(method_flow, flow))
               ,:Angle_RMS => angle_rmse(method_flow, flow)
               ,:Angle_mean => angle_mae(method_flow, flow)
               ,:Benchmark => bench
               ,:Truth_flow => showflow(flow, figtitle="Truth flow", ret="PyObject")
               ,:Method_flow => showflow(method_flow, figtitle=(method * " flow"), ret="PyObject")
               ,:Scatter_points => scatter(x, y, marker=:x)
               ))

bench = @benchmark test_scatter($flow_size, $points, $samples; method = $Polyharmonic(2))

typeof(bench)
## run

# parameters:
test_count = 3
flow_size = (200, 200)

Methods = ["Dier",
           "Scatter-Thin_plate_spline",
           "Scatter-Multiquadratic"]

Modes = ["Random_points",
         "Semigridded",
         "Gridded"]

Point_counts = [9,
                16,
                25,
                49,
                81,
                100,
                144]
import LAP_julia: helpers.mse, helpers.vec_len, helpers.angle_rmse, helpers.angle_mae

let
    index = 0

    for test in test_count
        flow = gen_tiled_flow(flow_size, 30, 40)
        for point_count in Point_counts
            for point_mode in Modes
                inds = gen_rand_points(flow, point_count, point_mode)
                for method in Methods
                        x = [ind[1] for ind in inds]
                        y = [ind[2] for ind in inds]
                    if method == "Dier"
                        # set inputs
                        z1 = [real(flow[ind]) for ind in inds]
                        z2 = [imag(flow[ind]) for ind in inds]
                        # do a benchmark
                        bench = @benchmark test_dier($x, $y, $z1, $z2)
                        # eval with current input
                        method_flow = test_dier(x, y, z1, z2)
                    elseif method == "Scatter-Thin_plate_spline"
                        # set inputs
                        points = reshape([ind[i] for i=1:length(inds[1]) for ind=inds], length(inds), length(inds[1]))'
                        samples = [flow[I] for I in inds]
                        # do a benchmark
                        bench = @benchmark test_scatter($flow_size, $points, $samples; method = $Polyharmonic(2))
                        # eval with current input
                        method_flow = test_scatter(flow_size, points, samples; method = Polyharmonic(2))
                    elseif method == "Scatter-Multiquadratic"
                        # set inputs
                        points = reshape([ind[i] for i=1:length(inds[1]) for ind=inds], length(inds), length(inds[1]))'
                        samples = [flow[I] for I in inds]
                        # do a benchmark
                        bench = @benchmark test_scatter($flow_size, $points, $samples; method = $Multiquadratic(2))
                        # eval with current input
                        method_flow = test_scatter(flow_size, points, samples; method = Multiquadratic(2))
                    end
                    # ADD TO DATA FRAME
                    index = index + 1
                    push!(df, Dict(:Index => index,
                                   :Method => method,
                                   :Mode => point_mode,
                                   :Point_count => point_count,
                                   :Median_speed => median(bench).time,
                                   :Magnitude_of_MSE => vec_len(mse(method_flow, flow)),
                                   :Angle_RMS => angle_rmse(method_flow, flow),
                                   :Angle_mean => angle_mae(method_flow, flow),
                                   :Benchmark => bench,
                                   :Truth_flow => showflow(flow, figtitle="Truth flow", ret="PyObject"),
                                   :Method_flow => showflow(method_flow, figtitle=(method * " flow"), ret="PyObject"),
                                   :Scatter_points => scatter(x, y, marker=:x)))
                    println("at index", index)
                end
            end
        end
    end
end # let


text:
# ---
# title: 1381. Design a Stack With Increment Operation
# id: problem1381
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Stack, Design
# link: <https://leetcode.com/problems/design-a-stack-with-increment-operation/description/>
# hidden: true
# ---
#
# Design a stack which supports the following operations.
#
# Implement the `CustomStack` class:
#
#   * `CustomStack(int maxSize)` Initializes the object with `maxSize` which is the maximum number of elements in the stack or do nothing if the stack reached the `maxSize`.
#   * `void push(int x)` Adds `x` to the top of the stack if the stack hasn't reached the `maxSize`.
#   * `int pop()` Pops and returns the top of stack or **-1** if the stack is empty.
#   * `void inc(int k, int val)` Increments the bottom `k` elements of the stack by `val`. If there are less than `k` elements in the stack, just increment all the elements in the stack.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["CustomStack","push","push","pop","push","push","push","increment","increment","pop","pop","pop","pop"]
#     [[3],[1],[2],[],[2],[3],[4],[5,100],[2,100],[],[],[],[]]
#     **Output**
#     [null,null,null,2,null,null,null,null,null,103,202,201,-1]
#     **Explanation**
#     CustomStack customStack = new CustomStack(3); // Stack is Empty []
#     customStack.push(1);                          // stack becomes [1]
#     customStack.push(2);                          // stack becomes [1, 2]
#     customStack.pop();                            // return 2 --> Return top of the stack 2, stack becomes [1]
#     customStack.push(2);                          // stack becomes [1, 2]
#     customStack.push(3);                          // stack becomes [1, 2, 3]
#     customStack.push(4);                          // stack still [1, 2, 3], Don't add another elements as size is 4
#     customStack.increment(5, 100);                // stack becomes [101, 102, 103]
#     customStack.increment(2, 100);                // stack becomes [201, 202, 103]
#     customStack.pop();                            // return 103 --> Return top of the stack 103, stack becomes [201, 202]
#     customStack.pop();                            // return 202 --> Return top of the stack 102, stack becomes [201]
#     customStack.pop();                            // return 201 --> Return top of the stack 101, stack becomes []
#     customStack.pop();                            // return -1 --> Stack is empty return -1.
#
#
#
#
# **Constraints:**
#
#   * `1 <= maxSize <= 1000`
#   * `1 <= x <= 1000`
#   * `1 <= k <= 1000`
#   * `0 <= val <= 100`
#   * At most `1000` calls will be made to each method of `increment`, `push` and `pop` each separately.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using KitBase, Flux, Plots
using KitBase.JLD2
using Flux: @epochs
import PyPlot as plt

###
# model
###

cd(@__DIR__)
#@load "../nn.jld2" nn
#@load "../sampler/nn_rif.jld2" nn
@load "nn_pro.jld2" nn

@load "kn3_ma4_aux.jld2" X Y
X1, Y1 = deepcopy(X), deepcopy(Y)
@load "kn3_ma5_aux.jld2" X Y
X2, Y2 = deepcopy(X), deepcopy(Y)
@load "kn2_ma4_aux.jld2" X Y
X3, Y3 = deepcopy(X), deepcopy(Y)
@load "kn2_ma5_aux.jld2" X Y
X4, Y4 = deepcopy(X), deepcopy(Y)
X = hcat(X1, X2, X3, X4)
Y = hcat(Y1, Y2, Y3, Y4)

device = cpu
data = Flux.Data.DataLoader((X, Y), shuffle = true) |> device
ps = params(nn)
#sqnorm(x) = sum(abs2, x)
#loss(x, y) = sum(abs2, nn(x) - y) / size(x, 2) #+ 1e-6 * sum(sqnorm, ps)
loss(x, y) = Flux.binarycrossentropy(nn(x), y)
cb = () -> println("loss: $(loss(X, Y))")
opt = ADAM()

@epochs 10 Flux.train!(loss, ps, data, opt, cb = Flux.throttle(cb, 1))

include("../common.jl")
accuracy(nn, X, Y)

@save "nn_pro.jld2" nn

###
# analysis
###

@load "kn3ref.jld2" ctr
#@load "kn2ref.jld2" ctr
include("tools.jl")
begin
    set = Setup(
        case = "cylinder",
        space = "2d2f2v",
        boundary = ["maxwell", "extra", "mirror", "mirror"],
        limiter = "minmod",
        cfl = 0.5,
        maxTime = 15.0, # time
    )
    ps = CSpace2D(1.0, 6.0, 60, 0.0, π, 50, 1, 1)
    vs = VSpace2D(-10.0, 10.0, 48, -10.0, 10.0, 48)
    gas = Gas(Kn = 1e-3, Ma = 5.0, K = 1.0)

    prim0 = [1.0, 0.0, 0.0, 1.0]
    prim1 = [1.0, gas.Ma * sound_speed(1.0, gas.γ), 0.0, 1.0]
    fw = (args...) -> prim_conserve(prim1, gas.γ)
    ff = function(args...)
        prim = conserve_prim(fw(args...), gas.γ)
        h = maxwellian(vs.u, vs.v, prim)
        b = h .* gas.K / 2 / prim[end]
        return h, b
    end
    bc = function(x, y)
        if abs(x^2 + y^2 - 1) < 1e-3
            return prim0
        else
            return prim1
        end
    end
    ib = IB2F(fw, ff, bc)

    ks = SolverSet(set, ps, vs, gas, ib)
end

rmap = zeros(ks.ps.nr, ks.ps.nθ)
@inbounds Threads.@threads for j = 1:ks.ps.nθ
    rmap[1, j] = 1
    for i = 2:ks.ps.nr
        swx1 = (ctr[i+1, j].w - ctr[i-1, j].w) / (1e-6 + ks.ps.x[i+1, j] - ks.ps.x[i-1, j])
        swy1 = (ctr[i+1, j].w - ctr[i-1, j].w) / (1e-6 + ks.ps.y[i+1, j] - ks.ps.y[i-1, j])
        swx2 = (ctr[i, j+1].w - ctr[i, j-1].w) / (1e-6 + ks.ps.x[i, j+1] - ks.ps.x[i, j-1])
        swy2 = (ctr[i, j+1].w - ctr[i, j-1].w) / (1e-6 + ks.ps.y[i, j+1] - ks.ps.y[i, j-1])
        swx = (swx1 + swx2) ./ 2
        swy = (swy1 + swy2) ./ 2
        #sw = sqrt.(sw1.^2 + sw2.^2)

        rmap[i, j] = judge_regime(ks, ctr[i, j].h, ctr[i, j].prim, swx, swy)
    end
end

@inbounds Threads.@threads for j = 1:ks.ps.nθ
    rmap[1, j] = 1
    for i = 2:ks.ps.nr
        swx1 = (ctr[i+1, j].w - ctr[i-1, j].w) / (1e-6 + ks.ps.x[i+1, j] - ks.ps.x[i-1, j])
        swy1 = (ctr[i+1, j].w - ctr[i-1, j].w) / (1e-6 + ks.ps.y[i+1, j] - ks.ps.y[i-1, j])
        swx2 = (ctr[i, j+1].w - ctr[i, j-1].w) / (1e-6 + ks.ps.x[i, j+1] - ks.ps.x[i, j-1])
        swy2 = (ctr[i, j+1].w - ctr[i, j-1].w) / (1e-6 + ks.ps.y[i, j+1] - ks.ps.y[i, j-1])
        swx = (swx1 + swx2) ./ 2
        swy = (swy1 + swy2) ./ 2
        sw = sqrt.(swx.^2 + swy.^2)
        tau = vhs_collision_time(ctr[i, j].prim, ks.gas.μᵣ, ks.gas.ω)

        rmap[i, j] = nn([ctr[i, j].w; sw; tau])[1] |> round |> Int
    end
end

begin
    contourf(
        ps.x[1:ks.ps.nr, 1:ks.ps.nθ],
        ps.y[1:ks.ps.nr, 1:ks.ps.nθ],
        rmap[:, :],
        ratio = 1,
    )
end

savefig("nn.pdf")

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
"""
load all content from shaders located in shaders folder in root dir
"""
function loadShaders()

CONTENT_SH = fileGetContents("shaders/global.glsl")
CONTENT_INST_VSH = fileGetContents("shaders/inst_vsh.glsl")
CONTENT_INST_VSH_GSH = fileGetContents("shaders/inst_vsh_gsh.glsl")
CONTENT_INST_GSH = fileGetContents("shaders/inst_gsh.glsl")
CONTENT_INST_FSH = fileGetContents("shaders/inst_fsh.glsl")
CONTENT_VSH_TEXTURE = fileGetContents("shaders/vsh_texture.glsl")
CONTENT_VSH = fileGetContents("shaders/vsh.glsl")
CONTENT_FSH = fileGetContents("shaders/fsh.glsl")
CONTENT_GSH = fileGetContents("shaders/gsh.glsl")

GLOBAL_CONTENT_SH = """
$(get_glsl_version_string())
$CONTENT_SH
"""

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

global INST_VSH = (:INST_VSH,"""
$GLOBAL_CONTENT_SH
$CONTENT_INST_VSH
""")

global INST_VSH_GSH = (:INST_VSH_GSH,"""
$GLOBAL_CONTENT_SH
$CONTENT_INST_VSH_GSH
""")

global INST_GSH = (:INST_GSH,"""
$GLOBAL_CONTENT_SH
$CONTENT_INST_GSH
""")

global INST_FSH = (:INST_FSH,"""
$GLOBAL_CONTENT_SH
$CONTENT_INST_FSH
""")

global VSH_TEXTURE = (:VSH,"""
$GLOBAL_CONTENT_SH
$CONTENT_VSH_TEXTURE
""")

global VSH = (:VSH,"""
$GLOBAL_CONTENT_SH
$CONTENT_VSH
""")

global FSH = (:FSH,"""
$GLOBAL_CONTENT_SH
$CONTENT_FSH
""")

global GSH = (:GSH,"""
$GLOBAL_CONTENT_SH
$CONTENT_GSH
""")

(INST_VSH, INST_VSH_GSH, INST_GSH, INST_FSH, VSH_TEXTURE, VSH, FSH, GSH)
end

max_stars_repo_name: Gilga/JuliaOpenGL

max_issues_repo_name: Gilga/JuliaOpenGL

max_forks_repo_name: Gilga/JuliaOpenGL


text:
module ReagentsTests

include("utils.jl")
include("test_immutablelists.jl")
include("test_internal_bags.jl")
include("test_try.jl")
include("test_treiberstack.jl")
include("test_msqueue.jl")
include("test_hmlist.jl")
include("test_dualcontainers.jl")
include("test_channels.jl")
include("test_nack.jl")
include("test_cancellablecontainers.jl")
include("test_blocking.jl")
include("test_anchors.jl")
include("test_locks.jl")
include("test_promises.jl")
include("test_dissolve.jl")
include("test_doctest.jl")
include("test_aqua.jl")

end  # module ReagentsTests


text:
# # Modal analysis of a free-floating steel circle

# Source code: [`circle_modal_conv_tut.jl`](circle_modal_conv_tut.jl)

# ## Description

# Vibration analysis of a free-floating steel ring. This is a
# benchmark from the NAFEMS Selected Benchmarks for Natural Frequency Analysis,
# publication: Test VM09: Circular Ring --  In-plane and Out-of-plane
# Vibration.

# The results can be compared with analytical expressions, but the main purpose
# is to compute data for extrapolation to the limit to predict the true natural
# frequencies.

# ## Reference frequencies

# There will be 6 rigid body modes (zero natural frequencies).

# The numerical results are due to the publication:
# NAFEMS Finite Element Methods & Standards, Abbassian, F., Dawswell, D. J., and
# Knowles, N. C. Selected Benchmarks for Natural Frequency Analysis, Test No.
# 6. Glasgow: NAFEMS, Nov., 1987.

# The reference values were analytically determined (Blevins, FORMULAS FOR
# DYNAMICS, ACOUSTICS AND VIBRATION, Table 4.16). Note that shear flexibility
# was neglected when computing the reference values.

# | Mode       |         Reference Value (Hz)  |  NAFEMS Target Value (Hz) |
# | -------   |     -------  |  ------- |
# | 7, 8 | (out of plane)   |        51.85          |         52.29  |
# | 9, 10 |  (in plane)       |       53.38         |          53.97  |
# | 11, 12 |  (out of plane)   |     148.8          |         149.7  |
# | 13, 14 |  (in plane)       |     151.0          |         152.4  |
# | 15, 16 |  (out of plane)   |     287.0          |         288.3  |
# |  17, 18 |  (in plane)      |      289.5         |          288.3  |

# ## Goals

# - Show convergence relative to reference values.
# - Compute data for extrapolation to the limit to predict the true natural
#   frequencies.
#

##
# ## Definition of the basic inputs

# Include the needed packages and modules.
using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexBeams.MeshFrameMemberModule: frame_member
using FinEtoolsFlexBeams.CrossSectionModule: CrossSectionCircle
using FinEtoolsFlexBeams.FEMMCorotBeamModule
CB = FEMMCorotBeamModule
using FinEtoolsFlexBeams.FEMMCorotBeamModule: FEMMCorotBeam
using FinEtoolsFlexBeams.FESetCorotBeamModule: MASS_TYPE_CONSISTENT_NO_ROTATION_INERTIA,
    MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA,
    MASS_TYPE_LUMPED_DIAGONAL_NO_ROTATION_INERTIA,
    MASS_TYPE_LUMPED_DIAGONAL_WITH_ROTATION_INERTIA

# The material parameters may be defined with the specification of the units.
# The elastic properties are:
E = 200.0 * phun("GPa")
nu = 0.3;

# The mass density is
rho = 8000 * phun("kg/m^3")
# Here are the cross-sectional dimensions and the length of the beam between supports.
radius = 1.0 * phun("m"); diameter = 0.1 * phun("m");

# We shall calculate these eigenvalues, but we are mostly interested in the
# first three  natural frequencies.
neigvs = 18;

# The mass shift needs to be applied since the structure is free-floating.
oshift = (2*pi*15)^2

# Here we get to choose the model: Bernoulli or Timoshenko
shear_correction_factor = 6/7 # Timoshenko
# shear_correction_factor = Inf # Bernoulli
cs = CrossSectionCircle(s -> diameter/2, s -> [1.0, 0.0, 0.0], shear_correction_factor)
@show cs.parameters(0.0)

# Here we can choose the mass- matrix type:

mtype = MASS_TYPE_CONSISTENT_NO_ROTATION_INERTIA
mtype = MASS_TYPE_CONSISTENT_WITH_ROTATION_INERTIA
mtype = MASS_TYPE_LUMPED_DIAGONAL_NO_ROTATION_INERTIA
mtype = MASS_TYPE_LUMPED_DIAGONAL_WITH_ROTATION_INERTIA

# Here are the formulas for the first two natural frequencies, obtained
# analytically with the shear flexibility  neglected. The parameters of the
# structure:
R = radius
Im = cs.parameters(0.0)[4]
m = rho * cs.parameters(0.0)[1]

# For instance the the first out of plane mode is listed in this table as
J = cs.parameters(0.0)[2]
G = E/2/(1+nu)
i = 2 # the first non-rigid body mode
@show i*(i^2-1)/(2*pi*R^2)*sqrt(E*Im/m/(i^2+E*Im/G/J))

# The first "ovaling" (in-plane) mode is:
i=2 # the first ovaling mode
@show i*(i^2-1)/(2*pi*R^2*(i^2+1)^(1/2))*sqrt(E*Im/m)

material = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

# We will generate this many elements  along the length of the circular ring.
results = let
    results = []
    for i in 1:3
        n = 10*2^i
        # beam elements along the member.
        tolerance = radius/n/1000;
        # Generate the mesh of a straight member.
        fens, fes = frame_member([0 0 0; 2*pi 0 0], n, cs)
        # Twist the straight member into a ring.
        for i in 1:count(fens)
            a = fens.xyz[i, 1]
            fens.xyz[i, :] .= (radius+radius*cos(a), radius*sin(a), 0)
        end
        # Merge the nodes of the bases, which involves renumbering the connectivity.
        fens, fes = mergenodes(fens, fes, tolerance, [1, n+1])
        # Generate the discrete model.
        geom0 = NodalField(fens.xyz)
        u0 = NodalField(zeros(size(fens.xyz, 1), 3))
        using FinEtoolsFlexBeams.RotUtilModule: initial_Rfield
        Rfield0 = initial_Rfield(fens)
        dchi = NodalField(zeros(size(fens.xyz, 1), 6))
        applyebc!(dchi)
        numberdofs!(dchi);
        femm = FEMMCorotBeam(IntegDomain(fes, GaussRule(1, 2)), material);
        K = CB.stiffness(femm, geom0, u0, Rfield0, dchi);
        M = CB.mass(femm, geom0, u0, Rfield0, dchi; mass_type = mtype);
        # Solve the free vibration problem.
        evals, evecs, nconv = eigs(K + oshift * M, M; nev=neigvs, which=:SM, ncv = 3*neigvs, maxiter = 2000);
        # Correct for the mass shift.
        evals = evals .- oshift;
        sigdig(n) = round(n * 10000) / 10000
        fs = real(sqrt.(complex(evals)))/(2*pi)
        println("Eigenvalues: $(sigdig.(fs)) [Hz]")
        push!(results, (fs[6+1], fs[6+3], fs[6+5]))
    end
    results # return these values from the block
end

@show results

##
# ## Richardson extrapolation

# Here we will use Richardson extrapolation from the three sets of data. This
# will allow us to predict the convergence rate and the true solution for each
# of the three frequencies (or rather the pairs of frequencies, 7 and 8, 9 and
# 10, and 11 and 12).

# We will immediately set up the convergence plots. We will extrapolate and then
# compute from that the normalized error to be plotted with respect to the
# refinement factors (which in this case are 4, 2, and 1). We use the
# refinement factor as a convenience: we will calculate the element size by
# dividing the circumference of the ring with a number of elements generated
# circumferentially.

using FinEtools.AlgoBaseModule: richextrapol

using Gnuplot
@gp  "set terminal windows 1 "  :-

# Modes 7 and 8
sols = [r[1] for r in results]
resextrap = richextrapol(sols, [4.0, 2.0, 1.0])
print("Predicted frequency 7 and 8: $(resextrap[1])\n")
errs = abs.(sols .- resextrap[1])./resextrap[1]
@gp  :- 2*pi*radius./[80, 160, 320] errs " lw 2 lc rgb 'red' with lp title 'Mode 7, 8' "  :-

# Modes 9 and 10
sols = [r[2] for r in results]
resextrap = richextrapol(sols, [4.0, 2.0, 1.0])
print("Predicted frequency 9 and 10: $(resextrap[1])\n")
errs = abs.(sols .- resextrap[1])./resextrap[1]
@gp  :- 2*pi*radius./[80, 160, 320] errs " lw 2 lc rgb 'green' with lp title 'Mode 9, 10' "  :-

# Modes 11 and 12
sols = [r[3] for r in results]
resextrap = richextrapol(sols, [4.0, 2.0, 1.0])
print("Predicted frequency 11 and 12: $(resextrap[1])\n")
errs = abs.(sols .- resextrap[1])./resextrap[1]
@gp  :- 2*pi*radius./[80, 160, 320] errs " lw 2 lc rgb 'blue' with lp title 'Mode 11, 12' "  :-

@gp  :- "set xrange [0.01:0.1]" "set logscale x" :-
@gp  :- "set logscale y" :-
@gp  :- "set xlabel 'Element size'" :-
@gp  :- "set ylabel 'Normalized error [ND]'" :-
@gp  :- "set title 'Beam: Convergence of modes 7, ..., 12'"


text:
"""
    zi = interpola(z, ns, nt)

Interpolate field `z` to size `(ns,ns,nt)` using nearest neighbors.
"""
function interpola(z, ns, nt)
    nss = size(z)
    if (length(nss)>=3)
        nas = nss[1]; nat = nss[3]
    else
        nas = nss[1]; nat = 1
    end

    sdim = div(ns,nas)
    tdim = div(nt,nat)
    ir = 1:sdim; jr = ir; kr = 1:tdim
    rs = 0:sdim:(ns-sdim)
    rt = 0:tdim:(nt-tdim)
    zi = zeros(ns, ns, nt)
    for i = 1:sdim
        for j = 1:sdim
            for k = 1:tdim
                zi[i.+rs,j.+rs,k.+rt] = z
            end
        end
    end

    return zi
end


text:
# [test/interpd_tests.jl]
using Test

@testset "Interpd Tests" begin

    using JLD2
    import ODEHybrid: interp1, interpd

    function generate_random_times(matrix_size = 1)
        """
            Convenience function used to generate ind and dep variables,
                as well as interpolation points
        """
        start = round(10.0 * rand())
        step  = round(2.0 * rand(), digits = 1)
        if abs(step) < 0.1
            step = 0.1
        end
        fin   = start + round(15.0 * rand())

        x = Array(start:step:fin) # Convert to an array for MATLAB
        y = zeros(length(x), matrix_size)
        for i = 1:matrix_size
            dx = 5.0 * rand()
            y_int = 10.0 * rand()
            y[:, i] .= dx .* x .+ y_int
        end

        if matrix_size == 1
            y = y[:, 1]
        end

        xq = Array(start:step:fin)
        xq = round.(Array(xq .+ randn(length(xq))), digits = 2)
        xq[ xq .< x[1] ] .= x[1]
        xq[ xq .> x[end]] .= x[end]

        return x, y, xq
    end;

    # Run a bunch of random interpolations to see if they match
    @testset "Interp1" begin # verbose = true
        @testset "Vector dependent variable" begin
            # data = []
            # for i = 1:50
            #     push!(data, generate_random_times())
            # end

            # r_mats = []
            # for i = 1:50
            #     x, y, xq = data[i]
            #     if length(x) < 3
            #         push!(r_mats, [])
            #         continue
            #     end
            #     r = mat"interp1($x, $y, $xq)";
            #     push!(r_mats, r)
            # end

            @load "solutions/interp1_vec_dep.jld2" data  r_mats

            for i = 1:50
                x, y, xq = data[i]
                if length(x) < 3
                    continue
                end

                # r_og = mat"interp1($x, $y, $xq)"
                r_jl = interp1((x,), y, xq)

                @test r_jl == r_mats[i]
            end
        end

        # Test for a non-column vector dependent variable
        @testset "Matrix dependent variable" begin
            # data = []
            # for i = 1:50
            #     N = Int(round(5.0 * rand())) + 1;
            #     if N ≤ 1
            #         N = 2
            #     end
            #     push!(data, generate_random_times(N))
            # end

            # r_mats = []
            # for i = 1:50
            #     x, y, xq = data[i]
            #     if length(x) < 3
            #         push!(r_mats, [])
            #         continue
            #     end
            #     r = mat"interp1($x, $y, $xq)";
            #     push!(r_mats, r)
            # end

            @load "solutions/interp1_mat_dep.jld2" data r_mats

            for i = 1:50
                x, y, xq = data[i]
                if length(x) < 3
                    continue
                end

                r_jl = interp1((x,), y, xq)

                @test r_jl == r_mats[i]
            end
        end
    end;

    # Run the full Interpd in a variety of cases to verify the julia and matlab implementations match
    @testset "Interpd" begin

        # Test when there are doubled values
        @testset "Doubled values" begin
            x = [1.0, 2, 3, 3, 3, 4, 5, 5, 6, 6, 7, 8, 10]
            y = 5.0 .* x .- 16
            xi = [1.5, 2.75, 5.5, 8.9]
            # og_ans = mat"interpd($x, $y, $xi)"
            @load "solutions/interpd_doubled_values_1.jld2" og_ans;
            jl_ans = interpd(x, y, xi)
            @test og_ans == jl_ans

            x = [5.0, 5.0, 6.0, 10.0, 10.0, 15.0, 15.0, 15.0, 15.0]
            y = 2.0 .* x .+ 1.0
            xi = [5.5, 6.5, 14.14, 14.14, 14.5 ]
            # og_ans = mat"interpd($x, $y, $xi)"
            @load "solutions/interpd_doubled_values_2.jld2" og_ans
            jl_ans = interpd(x, y, xi)
            @test og_ans == jl_ans
        end

        # Test with an unequal step size
        @testset "Unequal step size" begin
            # data = []
            # for i = 1:25
            #     x = sort(10.0 * round.(rand(10), digits = 2));
            #     dist = x[end] - x[1];
            #     y = 5.0 * rand() .* x .+ 5.0 * rand();
            #     xq = Array(range(x[1], x[end], step = dist/10));

            #     push!(data, [x, y, xq])
            # end

            # r_mats = []
            # for i = 1:25
            #     x, y, xq = data[i]
            #     rp = mat"""interpd($x, $y, $xq, "+")""";
            #     rn = mat"""interpd($x, $y, $xq, "-")""";
            #     push!(r_mats, (rp, rn))
            # end

            @load "solutions/interpd_unequal_step.jld2" data r_mats

            for i = 1:25
                x, y, xq = data[i]
                rp_mat, rn_mat = r_mats[i]

                mode = "-";
                jl_ans = interpd(x, y, xq, mode);
                @test jl_ans == rn_mat

                mode = "+"
                jl_ans = interpd(x, y, xq, mode)
                @test rp_mat == jl_ans
            end
        end

        # Test with an [n × m] dependent variable (doubled, mode, step_size)
        @testset "Doubled with matrix dependent" begin
            x = [5.0, 5.0, 6.0, 10.0, 10.0, 15.0, 15.0, 15.0, 15.0]
            y1 = 2.0 .* x .+ 1.0
            y2 = 5.5 .* x .- 7.2
            y3 = -14.2 .* x .+ 17.8
            y = [y1 y2 y3]
            xq = [5.5, 6.5, 14.14, 14.14, 14.5 ]

            mode = "+";
            # og_ans = mat"interpd($x, $y, $xq, $mode)";
            @load "solutions/interpd_doubled_with_matrix_1.jld2" og_ans
            jl_ans = interpd(x, y, xq, mode);
            @test og_ans == jl_ans

            mode = "-";
            # og_ans = mat"interpd($x, $y, $xq, $mode)";
            @load "solutions/interpd_doubled_with_matrix_2.jld2" og_ans
            jl_ans = interpd(x, y, xq, mode);
            @test og_ans == jl_ans
        end

    end;

    @testset "Interpolating Instantaneous Jump" begin

        # mat"""
        #     t = [2.76, 2.91, 3,   3,   3.12].';
        #     x = [0.2,  0.3,  0.4, 1.1, 1.2].';

        #     ti = (2.8:0.1:3.1).';
        #     $xi_p_mat = interpd(t, x, ti, '+');
        #     $xi_m_mat = interpd(t, x, ti, '-');
        # """

        t = [2.76; 2.91; 3; 3; 3.12];
        x = [0.2;  0.3;  0.4; 1.1; 1.2];
        ti = (2.8:0.1:3.1);
        xi_p_jul = interpd(t, x, ti, "+");
        xi_m_jul = interpd(t, x, ti, "-");

        @load "solutions/interpd_instant_jump.jld2" xi_p_mat xi_m_mat

        @test xi_p_mat ≈ xi_p_jul
        @test xi_m_mat ≈ xi_m_jul
    end

end;


text:
using DiscreteElasticRods
using Test
include("unit_tests.jl")
include("gradient_tests.jl")


text:
"""Compute the KL Divergence between the GP Prior and the variational distribution for the variational full batch model"""
function GaussianKL(model::VGP)
    return sum(broadcast(GaussianKL,model.μ,model.μ₀,model.Σ,model.invKnn))
end

"""Compute the KL Divergence between the Sparse GP Prior and the variational distribution for the sparse variational model"""
function GaussianKL(model::SVGP)
    return sum(broadcast(GaussianKL,model.μ,model.μ₀,model.Σ,model.invKmm))
end

function GaussianKL(μ::AbstractVector{T},μ₀::MeanPrior,Σ::Symmetric{T,Matrix{T}},invK::Symmetric{T,Matrix{T}}) where {T<:Real}
    0.5*(-logdet(Σ)-logdet(invK)+opt_trace(invK,Σ)+dot(μ-μ₀,invK*(μ-μ₀))-length(μ))
end

""" Compute the equivalent of KL divergence between an improper prior p(λ) (``1_{[0,\\infty]}``) and a variational Gamma distribution"""
function GammaEntropy(model::AbstractGP)
    return model.inference.ρ*(-sum(model.likelihood.α)+sum(log,model.likelihood.β[1])-sum(lgamma,model.likelihood.α)-dot(1.0.-model.likelihood.α,digamma.(model.likelihood.α)))
end

"""KL(q(ω)||p(ω)), where q(ω) = IG(α,β) and p(ω) = IG(α_p,β_p)"""
function InverseGammaKL(α,β,α_p,β_p)
    sum((α_p-α).*digamma(α_p) .- log.(gamma.(α_p)).+log.(gamma.(α)) .+  α.*(log.(β_p).-log.(β)).+α_p.*(β.-β_p)./β_p)
end

"""KL(q(ω)||p(ω)), where q(ω) = Po(ω|γ) and p(ω) = Po(ω|λ) with ν = E[λ] and ψ = E[log(λ)]"""
function PoissonKL(γ::AbstractVector{<:Real},λ::Real)
    λ*length(γ)-(1.0+log(λ))*sum(γ)+sum(xlogx,γ)
end

"""KL(q(ω)||p(ω)), where q(ω) = Po(ω|γ) and p(ω) = Po(ω|λ) with ν = E[λ] and ψ = E[log(λ)]"""
function PoissonKL(γ::AbstractVector{<:Real},ν::AbstractVector{<:Real},ψ::AbstractVector{<:Real})
    sum(ν)-sum(γ)+sum(xlogx,γ)-dot(γ,ψ)
end

"""KL(q(ω)||p(ω)), where q(ω) = PG(b,c) and p(ω) = PG(b,0). θ = 𝑬[ω]"""
function PolyaGammaKL(b,c,θ)
    dot(b,logcosh.(0.5*c))-0.5*dot(abs2.(c),θ)
end

"""Compute Entropy for Generalized inverse Gaussian latent variables (BayesianSVM)"""
function GIGEntropy(model::AbstractGP{<:BayesianSVM})
    return model.inference.ρ*sum(broadcast(b->0.5*sum(log.(b))+sum(log.(2.0*besselk.(0.5,sqrt.(b))))-0.5*sum(sqrt.(b)),model.likelihood.ω))
end

"""Entropy of GIG variables with parameters a,b and p and omitting the derivative d/dpK_p cf <https://en.wikipedia.org/wiki/Generalized_inverse_Gaussian_distribution#Entropy>"""
function GIGEntropy(a,b,p)
    sqrtab = sqrt.(a.*b)
    return sum(0.5*log.(a./b))+sum(log.(2*besselk.(p,sqrtab)))+ sum(0.5*sqrtab./besselk.(p,sqrtab).*(besselk.(p+1,sqrtab)+besselk.(p-1,sqrtab)))
end


text:
"""
Equivalent to `aws s3 sync <from> <to> <tags>` in terminal.

# Example

```julia
import Util
Util.s3sync(from="resultsdir", to="s3://my-bucket", tags=`--exclude '*.nfs'`)
```
"""
s3sync(; from, to, tags=``) = run(`aws s3 sync $(from) $(to) $(tags)`)


text:
using EzPugiXML
using Test

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


text:
using ASTInterpreter

# Simple evaling of function argument
function evalfoo1(x,y)
    x+y
end
interp = ASTInterpreter.enter_call_expr(nothing, :($(evalfoo1)(1,2)))
state = dummy_state(interp)
ok, res = ASTInterpreter.eval_code(state, "x")
@assert res == 1

ok, res = ASTInterpreter.eval_code(state, "y")
@assert res == 2

# Evaling with sparams
function evalsparams{T}(x::T)
    x
end
interp = ASTInterpreter.enter_call_expr(nothing, :($(evalsparams)(1)))
state = dummy_state(interp)
ok, res = ASTInterpreter.eval_code(state, "x")
@assert res == 1

ok, res = ASTInterpreter.eval_code(state, "T")
@assert res == Int


text:
abstract type AbstractTrajectory end

"""
Bezier curve
θ: a vector whose length is (N+1); each element can also be vector,
e.g., θ = [[1, 2, 3], [2, 3, 4]]
"""
struct Bezier <: AbstractTrajectory
    θ::AbstractVector
    t0::Real
    tf::Real
end

function (bezier::Bezier)(t::Real)
    @unpack t0, tf, θ = bezier
    N = length(θ) -1
    @assert t <= tf && t >= t0
    s = (t-t0)/(tf-t0)
    basis = 0:N |> Map(i -> binomial(N, i) * s^i * (1-s)^(N-i)) |> collect
    basis .* θ |> sum  # sum of N_C_i * s^(N-i) * (1-s)^i * θ_i
end


text:
using StaticArrays, LinearAlgebra
using JuLIP: JVec
import JuLIP
using JuLIP.MLIPs: IPBasis
using JuLIP.Potentials: SZList, ZList
import JuLIP: alloc_temp, alloc_temp_d

using SHIPs.SphericalHarmonics: SHBasis, sizeY, cart2spher, index_y
using SHIPs.Rotations: ClebschGordan
using SparseArrays: SparseMatrixCSC, sparse
using LinearAlgebra: svd

import Base: Dict, convert, ==

export SHIPBasis, bodyorder

"""
`struct SHIPBasis` : one of the two main types provided by `SHIPs.jl`;
represents a SHIP basis.

The standard constructor is
```
SHIPBasis(spec::BasisSpec, trans::DistanceTransform, fcut::PolyCutoff; kwargs...)
```
Keyword arguments are
* filter = true
* pure = false
* T = Float64
"""
struct SHIPBasis{T, NZ, TJ} <: IPBasis
   J::TJ               # specifies the radial basis
   SH::SHBasis{T}      # specifies the angular basis
   # ------------------------------------------------------------------------
   bgrps::NTuple{NZ, Vector{Tuple}}  # specification of basis functions
                                     # each group def. by a tuple (izz, kk, ll)
   zlist::SZList{NZ}                 # list of species (S=static)
   # ------------------------------------------------------------------------
   alists::NTuple{NZ, AList}         # datastructure to assemble A
   aalists::NTuple{NZ, AAList}       # datastructure to assemble AA
   A2B::NTuple{NZ, SparseMatrixCSC{Complex{T}, IntS}}   # convert AA -> B
   firstb::NTuple{NZ, Vector{IntS}}  # for each bgrp, store the zeroth index of
                                     # in the B vector; firstb[iz0][end] is the
                                     # last index
end

# ---------------(de-)dictionisation---------------------------------
Dict(shipB::SHIPBasis) = Dict(
      "__id__" => "SHIPs_SHIPBasis_v3",
      "J" => Dict(shipB.J),
      "zlist" => Dict(shipB.zlist),
      "bgrps" => bgrp2vecvec.(shipB.bgrps)
   )
convert(::Val{:SHIPs_SHIPBasis_v3}, D::Dict) = SHIPBasis(D)
SHIPBasis(D::Dict) = SHIPBasis(TransformedJacobi(D["J"]),
                               SZList(D["zlist"]),
                               vecvec2bgrps(D["bgrps"]))
bgrp2vecvec(bgrp) = [ Vector{Int}(vcat(b...)) for b in bgrp ]
vecvec2bgrps(Vbs) = convert.(Vector{Tuple}, tuple(vecvec2bgrp.(Vbs)...))
vecvec2bgrp(Vb) = _vec2b.(Vb)
_vec2b(v::Vector) = (
   N = length(v) ÷ 3; ( SVector{N, Int16}(v[1:N]...),
                        SVector{N,  IntS}(v[N+1:2*N]...),
                        SVector{N,  IntS}(v[2*N+1:3*N]...) ) )

==(B1::SHIPBasis, B2::SHIPBasis) =
      (B1.J == B2.J) && (B1.bgrps == B2.bgrps) && (B1.zlist == B2.zlist)
# ------------------------------------------------------------------------

function SHIPBasis(spec::BasisSpec, trans::DistanceTransform, fcut::PolyCutoff;
                   kwargs...)
   J = TransformedJacobi(maxK(spec), trans, fcut)
   return SHIPBasis(spec, J; kwargs...)
end

function SHIPBasis(spec::BasisSpec, J; T=Float64, kwargs...)
   # precompute the rotation-coefficients
   Bcoefs = Rotations.CoeffArray(T)
   # instantiate the basis specification
   allKL, NuZ = generate_ZKL_tuples(spec, Bcoefs)
   # # get the Ylm basis coefficients
   # rotcoefs = precompute_rotcoefs(allKL, NuZ, Bcoefs)
   # convert to new format ...
   bgrps = convert_basis_groups(NuZ, allKL) # zkl tuples
   # z-list to get i2z and z2i maps
   Zs = spec.Zs
   @assert issorted(Zs)
   zlist = ZList([Zs...]; static=true)

   return SHIPBasis(J, zlist, bgrps; T=T, Bcoefs=Bcoefs, kwargs...)
end

function SHIPBasis(J, zlist::SZList, bgrps::NTuple{NZ, Vector{Tuple}};
                   filter = true, T = Float64, pure = false,
                   Bcoefs = Rotations.CoeffArray(T)
                   ) where {NZ}
   @assert pure == false
   SH = SHBasis(get_maxL(bgrps))
   alists, aalists = alists_from_bgrps(bgrps)        # zklm tuples, A, AA
   A2B, firstb = A2B_matrices(bgrps, alists, aalists, Bcoefs, T)
   preB = SHIPBasis( J, SH,
                     bgrps, zlist,
                     alists, aalists, A2B, firstb )
   if filter
      return alg_filter_rpi_basis(preB)
   end
   return preB
end

function A2B_matrices(bgrps, alists, aalists, Bcoefs, T=Float64)
   NZ = length(alists)
   A2B = Vector{Any}(undef, NZ)
   firstb = Vector{Vector{IntS}}(undef, NZ)
   idx0 = 0
   for iz0  = 1:NZ
      A2B_iz0, firstb_iz0 = A2B_matrix(bgrps[iz0], alists[iz0], aalists[iz0],
                                       Bcoefs, T)
      A2B[iz0] = A2B_iz0
      # firstb_iz0 gives indexing only within the iz0 groups; adding
      # idx0 gives the global basis indices
      firstb[iz0] = idx0 .+ firstb_iz0
      idx0 = firstb[iz0][end]
   end
   return ntuple(i->A2B[i], NZ), ntuple(i->firstb[i], NZ)
end

function A2B_matrix(bgrps, alist, aalist, Bcoefs, T=Float64)
   # allocate triplet format
   Irow, Jcol, vals = IntS[], IntS[], Complex{T}[]
   firstb = IntS[]
   idxB = 0
   # loop through all (zz, kk, ll) tuples; each specifies 1 to several B
   for (izz, kk, ll) in bgrps
      # store the zeroth index for this basis group
      push!(firstb, idxB)
      # get the rotation-coefficients for this basis group
      Ull = SHIPs.Rotations.basis(Bcoefs, ll)
      # loop over the columns of Ull -> each specifies a basis function
      for ibasis = 1:size(Ull, 2)
         idxB += 1
         # next we loop over the list of admissible mm to write the
         # CG-coefficients into the A2B matrix
         for (im, mm) in enumerate(_mrange(ll))
            # the (izz, kk, ll, mm) tuple corresponds to an entry in the
            # AA vector (for the species iz0) at index idxA:
            idxA = aalist[(izz, kk, ll, mm)]
            push!(Irow, idxB)
            push!(Jcol, idxA)
            push!(vals, Ull[im, ibasis])
         end
      end
   end
   push!(firstb, idxB)
   # create CSC: [   triplet    ]  nrows   ncols
   return sparse(Irow, Jcol, vals, idxB, length(aalist)), firstb
end

_zkl(νz, ZKL) = (νz.izz, _kl(νz.ν, νz.izz, ZKL)...)

function convert_basis_groups(NuZ, ZKL)
   BO = size(NuZ, 1)
   NZ = size(NuZ, 2)
   @assert NZ == length(ZKL)
   bgrps = ntuple(iz0 -> Tuple[], NZ)
   for iz0 = 1:NZ, νz in vcat(NuZ[:, iz0]...)
      izz, kk, ll = _zkl(νz, ZKL)
      push!(bgrps[iz0], (izz, kk, ll))
   end
   return bgrps
end

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

nspecies(B::SHIPBasis{T, NZ}) where {T, NZ} = NZ

bodyorder(ship::SHIPBasis) = maximum(bodyorder, ship.aalists)

# get_maxL(allKL) = maximum( maximum( kl.l for kl in allKL_ )
#                            for allKL_ in allKL )

get_maxL(bgrps::NTuple{NZ, Vector{Tuple}}) where {NZ} = maximum(get_maxL, bgrps)
get_maxL(bgrp::Vector{Tuple}) = maximum(maximum(zkl[3]) for zkl in bgrp)

# the length of the basis depends on how many RI-coefficient sets there are
# so we have to be very careful how we define this.
Base.length(ship::SHIPBasis) = sum(size(A2B, 1) for A2B in ship.A2B)

# ----------------------------------------------
#      Computation of the B-basis
# ----------------------------------------------

alloc_B(ship::SHIPBasis, args...) = zeros(Float64, length(ship))
alloc_dB(ship::SHIPBasis, N::Integer) = zeros(JVec{Float64}, N, length(ship))
alloc_dB(ship::SHIPBasis, Rs::AbstractVector, args...) = alloc_dB(ship, length(Rs))

alloc_temp(ship::SHIPBasis{T, NZ}, args...) where {T, NZ} = (
      A = [ alloc_A(ship.alists[iz0])  for iz0 = 1:NZ ],
      AA = [ alloc_AA(ship.aalists[iz0])  for iz0 = 1:NZ ],
      Bc = [ zeros(Complex{T}, size(ship.A2B[iz0], 1)) for iz0=1:NZ],
      J = alloc_B(ship.J),
      Y = alloc_B(ship.SH),
      tmpJ = alloc_temp(ship.J),
      tmpY = alloc_temp(ship.SH)
   )

alloc_temp_d(shipB::SHIPBasis, Rs::AbstractVector{<:JVec}, args...) =
      alloc_temp_d(shipB, length(Rs))

function alloc_temp_d(ship::SHIPBasis{T, NZ}, N::Integer) where {T, NZ}
   J1 = alloc_B(ship.J)
   dJ1 = alloc_dB(ship.J)
   Y1 = alloc_B(ship.SH)
   dY1 = alloc_dB(ship.SH)
   return (
         A = [ alloc_A(ship.alists[iz0])  for iz0 = 1:NZ ],
        dA = [ zeros(JVec{Complex{T}}, N, length(ship.alists[iz0])) for iz0 = 1:NZ ],
        AA = [ alloc_AA(ship.aalists[iz0])  for iz0 = 1:NZ ],
       dBc = [zeros(JVec{Complex{T}}, size(ship.A2B[iz0], 1)) for iz0=1:NZ],
      dAAj = [ zeros(JVec{Complex{T}}, length(ship.aalists[iz0])) for iz0 = 1:NZ ],
         JJ = zeros(eltype(J1), N, length(J1)),
        dJJ = zeros(eltype(dJ1), N, length(dJ1)),
         YY = zeros(eltype(Y1), N, length(Y1)),
        dYY = zeros(eltype(dY1), N, length(dY1)),
        J = J1,
       dJ = dJ1,
        Y = Y1,
       dY = dY1,
      tmpJ = alloc_temp_d(ship.J, N),
      tmpY = alloc_temp_d(ship.SH, N)
      )
end

"""
return all basis indices corresponding to species index iz0
"""
function _get_I_iz0(ship::SHIPBasis, iz0)
   idx0 = 0
   for iz = 1:iz0-1
      idx0 += size(ship.A2B[iz], 1)
   end
   idxend = idx0 + size(ship.A2B[iz0], 1)
   I_iz0 = (idx0+1):idxend
   @assert I_iz0 == (ship.firstb[iz0][1]+1):ship.firstb[iz0][end]
   return I_iz0
end

# _get_I_iz0(ship::SHIPBasis, iz0) = (ship.firstb[iz0][1]+1):ship.firstb[iz0][end]

function eval_basis!(B, tmp, ship::SHIPBasis{T},
                     Rs::AbstractVector{<: JVec},
                     Zs::AbstractVector{<: Integer},
                     z0::Integer ) where {T}
   fill!(B, 0)
   iz0 = z2i(ship, z0)
   precompute_A!(tmp, ship, Rs, Zs, iz0)
   precompute_AA!(tmp, ship, iz0)
   _my_mul!(tmp.Bc[iz0], ship.A2B[iz0], tmp.AA[iz0])
   Iz0 = _get_I_iz0(ship, iz0)
   B[Iz0] .= real.(tmp.Bc[iz0])
   return B
end

function eval_basis_d!(dB, tmp, ship::SHIPBasis{T},
                       Rs::AbstractVector{<: JVec},
                       Zs::AbstractVector{<: Integer},
                       z0::Integer ) where {T}
   fill!(dB, zero(JVec{T}))
   iz0 = z2i(ship, z0)
   len_AA = length(ship.aalists[iz0])
   precompute_dA!(tmp, ship, Rs, Zs, iz0)
   for j = 1:length(Rs)
      dAAj = grad_AA_Rj!(tmp, ship, j, Rs, Zs, iz0)  # writes into tmp.dAAj[iz0]
      _my_mul!(tmp.dBc[iz0], ship.A2B[iz0], dAAj)
      Iz0 = _get_I_iz0(ship, iz0)
      @inbounds for i = 1:length(tmp.dBc[iz0])
         dB[j, Iz0[i]] = real(tmp.dBc[iz0][i])
      end
   end
   return dB
end

# -------------------------------------------------------
# linking to the functions implemented in `Alist.jl`

precompute_A!(tmp, ship::SHIPBasis, Rs, Zs, iz0) =
   precompute_A!(tmp.A[iz0], tmp, ship.alists[iz0], Rs, Zs, ship)

precompute_dA!(tmp, ship::SHIPBasis,
                    Rs::AbstractVector{<:JVec},
                    Zs::AbstractVector{<:Integer}, iz0 ) =
   precompute_dA!(tmp.A[iz0], tmp.dA[iz0], tmp, ship.alists[iz0],
                  Rs, Zs, ship)

precompute_AA!(tmp, ship::SHIPBasis, iz0) =
   precompute_AA!(tmp.AA[iz0], tmp.A[iz0], ship.aalists[iz0])

# -------------------------------------------------------------
#       Filtering the extra basis functions
# -------------------------------------------------------------

len_bgrp(shpB::SHIPBasis, igrp, iz0) =
      shpB.firstb[iz0][igrp+1] - shpB.firstb[iz0][igrp]

alllen_bgrp(shpB::SHIPBasis, iz0) =
      [len_bgrp(shpB, i, iz0) for i = 1:length(shpB.bgrps[iz0])]

I_bgrp(shpB::SHIPBasis, igrp, iz0) =
      (shpB.firstb[iz0][igrp]+1):shpB.firstb[iz0][igrp+1]

function _algebraic_gramian(ship, zkl, Igr, U, iz0)
   izz, kk, ll = zkl
   N = length(kk)
   n = length(Igr)
   @assert size(U, 1) == n
   G = zeros(n, n)
   for σ in permutations(1:N)
      if (izz[σ] != izz) || (kk[σ] != kk) || (ll[σ] != ll); continue; end
      for mm1 in SHIPs._mrange(ll), mm2 in SHIPs._mrange(ll)
         if mm1[σ] == mm2
            iU1 = ship.aalists[iz0][(izz, kk, ll, mm1)]
            iU2 = ship.aalists[iz0][(izz, kk, ll, mm2)]
            for i1 = 1:n, i2 = 1:n
               G[i1, i2] += conj(U[i1, iU1]) * U[i2, iU2]
            end
         end
      end
   end
   return G
end

function _alg_filter_group(ship::SHIPBasis{T}, zkl, Igr, U, iz0) where {T}
   G = _algebraic_gramian(ship, zkl, Igr, U, iz0)
   S = svd(G)
   rk = rank(G)
   UT = convert(SparseMatrixCSC{Complex{T},IntS}, S.U[:, 1:rk]')
   return Diagonal(sqrt.(S.S[1:rk])) * UT * U
end

function alg_filter_rpi_basis(preB::SHIPBasis{T, NZ}) where {T, NZ}

   new_A2B = [ sparse(IntS[], IntS[], Complex{T}[], 0, size(preB.A2B[1], 2))
               for iz0 = 1:NZ ]
   new_firstb = [ IntS[] for iz0 = 1:NZ ]
   bidx0 = 0

   for iz0 = 1:NZ
      for igrp = 1:length(preB.bgrps[iz0])
         zkl = preB.bgrps[iz0][igrp]
         Ib_grp =  I_bgrp(preB, igrp, iz0) .- preB.firstb[iz0][1]
         U = preB.A2B[iz0][Ib_grp, :]
         Ufiltered = _alg_filter_group(preB, zkl, Ib_grp, U, iz0)
         new_A2B[iz0] = vcat(new_A2B[iz0], Ufiltered)
         push!(new_firstb[iz0], bidx0)
         bidx0 += size(Ufiltered, 1)
      end
      # and then add one more to get the total length of the basis
      push!(new_firstb[iz0], bidx0)
      # double-check that we have exactly the right number of firstb entries
      @assert length(new_firstb[iz0]) == length(preB.firstb[1])
   end

   return SHIPBasis(preB.J, preB.SH, preB.bgrps, preB.zlist,
                    preB.alists, preB.aalists,
                    ntuple(i->new_A2B[i], NZ),
                    ntuple(i->new_firstb[i], NZ) )
end

# -------------------------------------------------------------
#       JuLIP Calculators: energies and forces
# -------------------------------------------------------------

# TODO: move all of this into JuLIP.MLIPs

using NeighbourLists: max_neigs, neigs
using JuLIP: Atoms, sites, neighbourlist
using JuLIP.Potentials: neigsz!
import JuLIP: energy, forces, virial, cutoff, site_energy, site_energy_d

cutoff(shipB::SHIPBasis) = cutoff(shipB.J)

function energy(shipB::SHIPBasis, at::Atoms{T}) where {T}
   E = zeros(length(shipB))
   B = alloc_B(shipB)
   nlist = neighbourlist(at, cutoff(shipB))
   maxnR = maxneigs(nlist)
   tmp = alloc_temp(shipB, maxnR)
   tmpRZ = (R = zeros(JVec{T}, maxnR), Z = zeros(Int16, maxnR))
   for i = 1:length(at)
      j, R, Z = neigsz!(tmpRZ, nlist, at, i)
      eval_basis!(B, tmp, shipB, R, Z, at.Z[i])
      E[:] .+= B[:]
   end
   return E
end

function forces(shipB::SHIPBasis, at::Atoms{T}) where {T}
   # precompute the neighbourlist to count the number of neighbours
   nlist = neighbourlist(at, cutoff(shipB))
   maxR = max_neigs(nlist)
   # allocate space accordingly
   F = zeros(JVec{T}, length(at), length(shipB))
   dB = alloc_dB(shipB, maxR)
   tmp = alloc_temp_d(shipB, maxR)
   tmpRZ = (R = zeros(JVec{T}, maxR), Z = zeros(Int16, maxR))
   # assemble site gradients and write into F
   for i = 1:length(at)
      j, R, Z = neigsz!(tmpRZ, nlist, at, i)
      eval_basis_d!(dB, tmp, shipB, R, Z, at.Z[i])
      for a = 1:length(R)
         F[j[a], :] .-= dB[a, :]
         F[i, :] .+= dB[a, :]
      end
   end
   return [ F[:, iB] for iB = 1:length(shipB) ]
end

function virial(shipB::SHIPBasis, at::Atoms{T}) where {T}
   # precompute the neighbourlist to count the number of neighbours
   nlist = neighbourlist(at, cutoff(shipB))
   maxR = max_neigs(nlist)
   # allocate space accordingly
   V = zeros(JMat{T}, length(shipB))
   dB = alloc_dB(shipB, maxR)
   tmp = alloc_temp_d(shipB, maxR)
   tmpRZ = (R = zeros(JVec{T}, maxR), Z = zeros(Int16, maxR))
   # assemble site gradients and write into F
   for i = 1:length(at)
      j, R, Z = neigsz!(tmpRZ, nlist, at, i)
      eval_basis_d!(dB, tmp, shipB, R, Z, at.Z[i])
      for iB = 1:length(shipB)
         V[iB] += JuLIP.Potentials.site_virial(dB[:, iB], R)
      end
   end
   return V
end

function _get_neigs(at::Atoms{T}, i0::Integer, rcut) where {T}
   nlist = neighbourlist(at, rcut)
   maxR = maxneigs(nlist)
   tmpRZ = (R = zeros(JVec{T}, maxR), Z = zeros(Int16, maxR))
   j, R, Z = neigsz!(tmpRZ, nlist, at, i0)
   return j, R, Z
end

function site_energy(basis::SHIPBasis, at::Atoms, i0::Integer)
   j, Rs, Zs = _get_neigs(at, i0, cutoff(basis))
   return eval_basis(basis, Rs, Zs, at.Z[i0])
end

function site_energy_d(basis::SHIPBasis, at::Atoms{T}, i0::Integer) where {T}
   Ineigs, Rs, Zs = _get_neigs(at, i0, cutoff(basis))
   dEs = [ zeros(JVec{T}, length(at)) for _ = 1:length(basis) ]
   dB = alloc_dB(basis, length(Rs))
   tmp = alloc_temp_d(basis, length(Rs))
   eval_basis_d!(dB, tmp, basis, Rs, Zs, at.Z[i0])
   @assert dB isa Matrix{JVec{T}}
   @assert size(dB) == (length(Rs), length(basis))
   for iB = 1:length(basis), n = 1:length(Ineigs)
      dEs[iB][Ineigs[n]] += dB[n, iB]
      dEs[iB][i0] -= dB[n, iB]
   end
   return dEs
end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
import RCall.sexp

struct JuliaObject
    id :: Ptr{RCall.ExtPtrSxp}
    typ :: String
    JuliaObject(id :: Ptr{RCall.ExtPtrSxp}, typ = "Regular") = new(id, typ)
end

function getID(x :: JuliaObject)
    x.id
end

function getType(x :: JuliaObject)
    x.typ
end

function new_obj(obj, typ = "Regular")
    # wrap in a `Ref`
    refj = Ref(obj)
    jptr = pointer_from_objref(refj)
    s = RCall.makeExternalPtr(jptr)
    RCall.jtypExtPtrs[s] = refj
    RCall.registerCFinalizerEx(s)

    JuliaObject(s, typ)
end

## debugging for freeing
# function rm_obj(id)
#     RCall.decref_extptr(RCall.sexp(x))
# end

JuliaObject(x :: JuliaObject, typ = "Regular") = x
JuliaObject(x :: RObject, typ = "Regular") = new_obj(rcopy(x), typ)
JuliaObject(x :: RCall.Sxp, typ = "Regular") = new_obj(rcopy(x), typ)
JuliaObject(x, typ = "Regular") = new_obj(x, typ)

## Conversion related to JuliaObject

const makeJuliaObjectInR = reval("JuliaCall:::juliaobjectnew")

function sexp(x :: JuliaObject)
    RCall.rcall_p(makeJuliaObjectInR, getID(x), getType(x))
end

import RCall.rcopy

function rcopy(::Type{JuliaObject}, x::Ptr{EnvSxp})
    try
        # get(julia_object_stack, rcopy(Int32, RObject(x)[:getID]()))
        # get(julia_object_stack, rcopy(Int32, RObject(x)[:id]))
        RCall.jtypExtPtrs[RCall.sexp(x[:id])].x
    catch e
        nothing
    end
end

import RCall: RClass, rcopytype

rcopytype(::Type{RClass{:JuliaObject}}, x::Ptr{EnvSxp}) = JuliaObject

## Fallback conversions

@static if julia07
    import RCall.sexpclass
    sexpclass(x) = RClass{:JuliaObject}
    sexp(::Type{RClass{:JuliaObject}}, x) = sexp(JuliaObject(x))
else
    sexp(x) = sexp(JuliaObject(x))
end

## Regarding to issue #12, #13 and #16,
## we should use JuliaObject for general AbstractArray
@static if julia07
    @suppress_err begin
        JuliaCall.sexpclass(x :: AbstractArray{T}) where {T} = RClass{:JuliaObject}
    end

    ## AbstractArray{Any} should be converted to R List
    sexpclass(x :: AbstractArray{Any}) = RClass{:list}
else
    @suppress_err begin
        JuliaCall.sexp(x :: AbstractArray{T}) where {T} = sexp(JuliaObject(x))
    end

    ## AbstractArray{Any} should be converted to R List
    sexp(x :: AbstractArray{Any}) = sexp(VecSxp, x)
end

## Preserve BigFloat precision,
## as the design decision in issue #16
# sexp(x::AbstractArray{BigFloat}) = sexp(JuliaObject(x))
# sexp(x::BigFloat) = sexp(JuliaObject(x))

function setfield1!(object, name, value1)
    value = fieldtype(typeof(object), name)(value1)
    setfield!(object, name, value)
end

max_stars_repo_name: KevCaz/JuliaCall

max_issues_repo_name: KevCaz/JuliaCall

max_forks_repo_name: KevCaz/JuliaCall


text:
"""
    newton(x, xd, yd)
adalah fungsi yang digunakan untuk mencari nilai interpolasi pada titik/vektor `x`, jika
diketahui suatu himpunan pasangan terurut `(xd,yd)`.

# Example
```jl
julia> xd = [1,2,3,5];

julia> yd = [1.06 1.12 1.34 1.78];

julia> y,D = newton(4,xd,yd);

julia> y
1.6

julia> D
4×4 Array{Float64,2}:
 1.06  0.0   0.0    0.0
 1.12  0.06  0.0    0.0
 1.34  0.22  0.08   0.0
 1.78  0.22  0.0   -0.02
```
return solusi hampiran interpolasi `y` dan matriks beda-terbagi `D`.
"""
function newton(x, xd, yd)
  m=length(x);
  y=zeros(m)
  if m > 1
    for i=1:m
      # proses satu titik demi satu titik
      y[i],D =newton(x[i], xd, yd);
    end
    return y,D
  end
  #% periksa jumlah titik dan tentukan derajat polinom
  ntitik = length(xd);
  #% hitung tabel beda-terbagi (divided-difference)
  D = zeros(ntitik,ntitik)
  D[:,1]=yd;          #% kolom pertama
  for j=2:ntitik      #% kolom ke-2 dan seterusnya
    for k=j:ntitik
        D[k,j] = (D[k,j-1]-D[k-1,j-1])/(xd[k]-xd[k-j+1]);
    end
  end
  #% hitung interpolasi Newton
  y = D[1,1]; s = 1;
  for i=2:ntitik
    s = s * (x-xd[i-1]);
    y = y + D[i,i]*s;
  end
  return y, D
end