text:
using ReservoirComputing

const res_size = 20
const radius = 1.0
const sparsity = 0.1
const weight = 0.2
const jump_size = 3

#testing RandSparseReservoir implicit and esplicit constructors
reservoir_constructor = RandSparseReservoir(res_size, radius, sparsity)
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)

reservoir_constructor = RandSparseReservoir(res_size, radius=radius, sparsity=sparsity)
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)

#testing PseudoSVDReservoir implicit and esplicit constructors
reservoir_constructor = PseudoSVDReservoir(res_size, radius, sparsity)
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)
@test maximum(reservoir_matrix) <= radius

reservoir_constructor = PseudoSVDReservoir(res_size, max_value=radius, sparsity=sparsity)
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)
@test maximum(reservoir_matrix) <= radius

#testing DelayLineReservoir implicit and esplicit constructors
reservoir_constructor = DelayLineReservoir(res_size, weight)
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)
@test maximum(reservoir_matrix) == weight

reservoir_constructor = DelayLineReservoir(res_size, weight=weight)
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)
@test maximum(reservoir_matrix) == weight

#testing DelayLineReservoir implicit and esplicit constructors
reservoir_constructor = DelayLineBackwardReservoir(res_size, weight, weight)
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)
@test maximum(reservoir_matrix) == weight

reservoir_constructor = DelayLineBackwardReservoir(res_size, weight=weight, fb_weight=weight)
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)
@test maximum(reservoir_matrix) == weight

#testing SimpleCycleReservoir implicit and esplicit constructors
reservoir_constructor = SimpleCycleReservoir(res_size, weight)
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)
@test maximum(reservoir_matrix) == weight

reservoir_constructor = SimpleCycleReservoir(res_size, weight=weight)
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)
@test maximum(reservoir_matrix) == weight

#testing CycleJumpsReservoir implicit and esplicit constructors
reservoir_constructor = CycleJumpsReservoir(res_size, weight, weight, jump_size)
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)
@test maximum(reservoir_matrix) == weight

reservoir_constructor = CycleJumpsReservoir(res_size, cycle_weight=weight, jump_weight=weight, jump_size=jump_size)
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)
@test maximum(reservoir_matrix) == weight

#testing NullReservoir constructors
reservoir_constructor = NullReservoir()
reservoir_matrix = create_reservoir(reservoir_constructor, res_size)
@test size(reservoir_matrix) == (res_size, res_size)

max_stars_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_issues_repo_name: MartinuzziFrancesco/EchoStateNetwork

max_forks_repo_name: MartinuzziFrancesco/EchoStateNetwork


text:
using ENN.Neurons, ENN.TimePetriNets, GLMakie, GraphMakie

n1 = Neuron((((((),[:A1,:A2],:A),),[:B1,:B2],:B),),[:C1,:C2],:C)
n2 = Neuron((((),[:B1,:B2],:B),),[:trig],:A)
#tpn=TPN(neuron_sequential)

net = NeuralNetwork(
    Dict(
        :n1 => n1,
        :n2 => n2
    ),
    Dict(
        :n1 => [(:n1, :A, :A1), (:n2, :A, :trig)],
        :n2 => [(:n1, :B, :inh), (:n1, :B, :B2)],
    )
)
net_tpn = TPN(net)

f,ax,p = graphplot(net_tpn); f

ext: jl


text:
@testset "219.contains-duplicate-ii.jl" begin
    @test contains_nearby_duplicate([1, 2, 3, 1], 3) == true
    @test contains_nearby_duplicate([1, 0, 1, 1], 1) == true
    @test contains_nearby_duplicate([1, 2, 3, 1, 2, 3], 2) == false
    @test contains_nearby_duplicate([1, 1, 3, 2], 0) == false
end


text:
using EconLH.ProductionFunctionsLH
using Random, Test
# Tests that output is continuous around elasticity 1
function elast_close_to_one_test()
  n = 3;
  T = 5;
  xM = LinRange(0.5, 5.0, T) * LinRange(2.0, 0.1, n)';

  elastV = LinRange(0.95, 1.05, 50);
  y1V = zeros(length(elastV));
  for (j, substElast) in enumerate(elastV)
    fS = pf.make_test_ces(T, n, substElast);
    yV = output(fS, xM);
    y1V[j] = yV[1]
  end
  @test all(abs.(diff(y1V)) .< 0.01)
end

@testset "CES" begin
  elast_close_to_one_test();
end

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


text:
import SolveIPlusOrMinusAkronB: IPlusAtKronBWs, solvi_real_eliminate!,
    QuasiUpperTriangular, solve1!, transformation1, transform2,
    generalized_sylvester_solver!, solvi_complex_eliminate!, solviip_complex_eliminate!,
    solviip, solviip2, solviip_real_eliminate!

using Random
using Test

Random.seed!(123)

if false
n = 3
x = randn(n,n*n)
a = randn(n,n)
b = randn(n,n)
c = randn(n,n)
d = a*x + b*x*kron(c,c)
dold = d
ws = EyePlusAtKronBWs(n,n,n,2)

index = 1
depth = 1
nd = n^depth
t = QuasiUpperTriangular(lu(a\b)[2])
s = QuasiUpperTriangular(lu(c)[2])
d = zeros(n^2)
r = 1.0
d[1:n] = randn(n)
y = copy(d[1:n])
solvi_real_eliminate!(index,n,nd,1:nd,depth-1,r,t,s,d,ws)
d_target = -kron(s[[1],:]',t)*y
@test d[n+1:n^2] ≈ d_target[n+1:n^2]

index = 1
depth = 3
ws = EyePlusAtKronBWs(n,n,n,depth)
nd = n^depth
t = QuasiUpperTriangular(lu(a\b)[2])
s = QuasiUpperTriangular(lu(c)[2])
d = zeros(n^(depth+1))
r = 1.0
d[1:nd] = randn(nd)
y = copy(d[1:nd])
solvi_real_eliminate!(index,n,nd,1:nd,depth-1,r,t,s,d,ws)
d_target = -kron(s[1,:],kron(kron(s',s'),t))*y
@test d[nd+1:n^(depth+1)] ≈ d_target[nd+1:n^(depth+1)]

index = 2
d = zeros(n^(depth+1))
r = 1.0
d[nd+(1:nd)] = randn(nd)
y = copy(d[nd+(1:nd)])
solvi_real_eliminate!(index,n,nd,nd+(1:nd),depth-1,r,t,s,d,ws)
d_target = -kron(s[2,:],kron(s',s'),t)*y
@test d[2*nd+1:n^(depth+1)] ≈ d_target[2*nd+1:n^(depth+1)]

aa = a\b
t = QuasiUpperTriangular(lu(aa)[2])
s = QuasiUpperTriangular(lu(c)[2])

r = 1.0
d = randn(n)
d_orig = copy(d)
td = similar(d)
order = 0
depth = order
t2 = t*t
s2 = s*s
solve1!(r,depth,t,t2,s,s2,d,ws)
d_target = (I(n) + t)\d_orig
@test d_target ≈ d

r = 1.0
d = randn(n*n*n)
d_orig = copy(d)
td = similar(d)
order = 2
depth = order
solve1!(r,depth,t,t2,s,s2,d,ws)

d_target = (I(n^3) + kron(kron(s',s'),t))\d_orig

@test d ≈ d_target

a = 0.5
b1 = 0.1
b2 = 1.3
depth = 1
tt = QuasiUpperTriangular(t)
ss = QuasiUpperTriangular(s)
d = randn(2*n^(depth+1))
d_orig = copy(d)
transformation1(a,b1,b2,depth,tt,ss,d,ws)
d_target = d_orig + kron([a -b1; -b2 a],kron(s',t))*d_orig

@test d ≈ d_target

a = 0.5
b1 = 0.1
b2 = 1.3
depth = 2
tt = QuasiUpperTriangular(t)
ss = QuasiUpperTriangular(s)
d = randn(2*n^(depth+1))
d_orig = copy(d)
transformation1(a,b1,b2,depth,tt,ss,d,ws)
d_target = d_orig + kron([a -b1; -b2 a],kron(kron(s',s'),t))*d_orig

@test d ≈ d_target

index = 1
depth = 3
nd = n^depth
t = QuasiUpperTriangular(lu(aa)[2])
s = QuasiUpperTriangular(lu(c)[2])
s[1,2] = 0.5
s2 = QuasiUpperTriangular(s*s)
t2 = QuasiUpperTriangular(t*t)
d = zeros(n^(depth+1))
d[1:nd] = randn(nd)
d_orig = copy(d)
r1 = 1.0
r2 = 0.8
td = []
solviip_real_eliminate!(index,n,nd,1:nd,depth-1,r1,r2,t,t2,s,s2,d,ws)
d_target = -2*r1*kron(s[index,:],kron(kron(s',s'),t))*d_orig[1:nd] - (r1*r1+r2*r2)*kron(s2[index,:],kron(kron(s2',s2'),t2))*d_orig[1:nd]
@test d[nd+1:n^(depth+1)] ≈ d_target[nd+1:n^(depth+1)]

depth = 2
nd = n^depth
gamma = 0.3
delta1 = -0.4
delta2 = 0.6
d = randn(2*nd)
d_orig = copy(d)
transform2(r1, r2, gamma, delta1, delta2, nd, depth, t, t2, s, s2, d, ws)
G = [gamma delta1; delta2 gamma]
d_target = (I(2*nd) + 2*r1*kron(G,kron(s',t)) + (r1*r1 + r2*r2)*kron(G*G,kron(s2',t2)))*d_orig
@test d ≈ d_target

index = 1
depth = 3
nd = n^depth
t = QuasiUpperTriangular(schur(a\b)[1])
s = QuasiUpperTriangular(schur(c)[1])
s2 = QuasiUpperTriangular(s*s)
t2 = QuasiUpperTriangular(t*t)
d = zeros(n^(depth+1))
d[1:2*nd] = randn(2*nd)
d_orig = copy(d)
r1 = 1.0
r2 = 0.8
td = []
solviip_complex_eliminate!(index,n,nd,1:nd,depth-1,r1,r2,t,t2,s,s2,d,ws)
d_target = (-2*r1*kron(kron(s[1,:],kron(s',s')),t)*d_orig[1:nd] - (r1*r1+r2*r2)*kron(kron(s2[1,:],kron(s2',s2')),t2)*d_orig[1:nd]
            -2*r1*kron(kron(s[2,:],kron(s',s')),t)*d_orig[nd+1:2*nd] - (r1*r1+r2*r2)*kron(kron(s2[2,:],kron(s2',s2')),t2)*d_orig[nd+1:2*nd])
@test d[2*nd+1:n^3] ≈ d_target[2*nd+1:n^3]

depth = 1
n = 2
nd = n^depth
t = QuasiUpperTriangular(schur([1 2; -5 4])[1])
s = QuasiUpperTriangular(schur([1 -3; 2 3])[1])
s2 = QuasiUpperTriangular(s*s)
t2 = QuasiUpperTriangular(t*t)
d = randn(n^(depth+1))
d_orig = copy(d)
alpha = s[1,1]
# s is transposed
beta1 = s[2,1]
beta2 = s[1,2]
beta = sqrt(-beta1*beta2)
td = similar(d)
solviip(alpha,beta,depth,t,t2,s,s2,d,ws)
d_target = (I(n^2) + 2*alpha*kron(s',t) + (alpha*alpha + beta*beta)*kron(s2',t2))\d_orig
@test d ≈ d_target

depth = 1
n = 3
srand(1) #first diag block 2x2
#srand(2) # upper triangular
#srand(3) # second diag block 2x2
a = randn(n,n)
b = randn(n,n)
c = randn(n,n)
ws = EyePlusAtKronBWs(n,n,n,2)
nd = n^depth
t = QuasiUpperTriangular(schur(a\b)[1])
s = QuasiUpperTriangular(schur(c)[1])
s2 = QuasiUpperTriangular(s*s)
t2 = QuasiUpperTriangular(t*t)
d = randn(n^(depth+1))
d_orig = copy(d)
alpha = s[1,1]
# s is transposed
beta1 = s[2,1]
beta2 = s[1,2]
r1 = alpha
r2 = sqrt(-beta1*beta2)
td = similar(d)
solviip(r1,r2,depth,t,t2,s,s2,d,ws)
d_target = (I(n^2) + 2*alpha*kron(s',t) + (alpha*alpha + r2*r2)*kron(s2',t2))\d_orig
@test d ≈ d_target

depth = 2
nd = n^depth
t = QuasiUpperTriangular(schur(a\b)[1])
s = QuasiUpperTriangular(schur(c)[1])
s2 = QuasiUpperTriangular(s*s)
t2 = QuasiUpperTriangular(t*t)
d = randn(n^3)
d_orig = copy(d)
alpha = s[1,1]
# s is transposed
beta1 = s[2,1]
beta2 = s[1,2]
r1 = alpha
r2 = sqrt(-beta1*beta2)
td = similar(d)
solviip(r1,r2,depth,t,t2,s,s2,d,ws)
d_target = (I(n^3) + 2*alpha*kron(s',kron(s',t)) + (alpha*alpha + r2*r2)*kron(kron(s2',s2'),t2))\d_orig
@test d ≈ d_target

n = 4
a = randn(n,n)
b = randn(n,n)
c = randn(n,n)
t = QuasiUpperTriangular(schur(a\b)[1])
s = QuasiUpperTriangular(schur(c)[1])
s2 = QuasiUpperTriangular(s*s)
t2 = QuasiUpperTriangular(t*t)
order = 1
d = randn(n^(order+1))
d_orig = copy(d)
ws = EyePlusAtKronBWs(n,n,n,order)
println("begin test")
solver!(t,s,d,order,ws)
d_target = (I(n^(order+1)) + kron(s',t))\d_orig
#d_target = (I(n^(order+1)) + kron(kron(s',s'),t))\d_orig
@test d ≈ d_target

n = 4
a = randn(n,n)
b = randn(n,n)
c = randn(n,n)
t = QuasiUpperTriangular(schur(a\b)[1])
s = QuasiUpperTriangular(schur(c)[1])
s2 = QuasiUpperTriangular(s*s)
t2 = QuasiUpperTriangular(t*t)
order = 3
d = randn(n^(order+1))
d_orig = copy(d)
ws = EyePlusAtKronBWs(n,n,n,order)
println("begin test")
solver!(t,s,d,order,ws)
d_target = (I(n^(order+1)) + kron(s',kron(kron(s',s'),t)))\d_orig
@test d ≈ d_target

n = 4
a = randn(n,n)
b = randn(n,n)
c = randn(n,n)
t = QuasiUpperTriangular(schur(a\b)[1])
s = QuasiUpperTriangular(schur(c)[1])
s2 = QuasiUpperTriangular(s*s)
t2 = QuasiUpperTriangular(t*t)
order = 3
ws = EyePlusAtKronBWs(n,n,n,order)
d = randn(n^(order+1))
d_orig = copy(d)
solver!(t,s,d,order,ws)
d_target = (I(n^(order+1)) + kron(s',kron(kron(s',s'),t)))\d_orig
@test d ≈ d_target

n1 = 4
n2 = 3
a = randn(n1,n1)
b = randn(n1,n1)
c = randn(n2,n2)
order = 2
ws = EyePlusAtKronBWs(n1,n1,n2,order)
d = randn(n1,n2^order)
a_orig = copy(a)
b_orig = copy(b)
c_orig = copy(c)
d_orig = copy(d)

generalized_sylvester_solver!(a,b,c,vec(d),2,ws)
@test a_orig*d + b_orig*d*kron(c_orig,c_orig) ≈ d_orig
@test d ≈ reshape((kron(I(n2^order),a_orig) + kron(kron(c_orig',c_orig'),b_orig))\vec(d_orig),n1,n2^order)
end

function f(t,s,d,order,ws)
    for i = 1:100
        solver!(t,s,d,order,ws)
    end
end

#@profile  f(t,s,d,order,ws)
#Profile.print(combine=true,sortedby=:count)

function kron_power(x,order)
    if order == 0
        y = 1
    else
        y = x
        for i = 1:(order-1)
            y = kron(y, x)
        end
    end
    y
end

function test1(m, n, order)

    a = randn(m,m)
    b = randn(m,m)
    c = randn(n,n)
    d = randn(m*n^order)
    alpha = randn()
    beta = randn()
    t = lu(a\b).U
    s = lu(c).U
    t2 = QuasiUpperTriangular(t*t)
    s2 = QuasiUpperTriangular(s*s)
    t = QuasiUpperTriangular(t)
    s = QuasiUpperTriangular(s)
    r = 1.0
    nd = m*n^(order-1)

    ws = IPlusAtKronBWs(m, m, n, order)

    @test kron_power(a,1) == a
    @test kron_power(a,3) ≈ kron(kron(a,a),a)

    println("SOLVI_REAL_ELIMINATE")
    drange = 1:nd
    for index = 1:(n-1)
        d_orig = copy(d)
        @time solvi_real_eliminate!(index, n, nd, drange, order-1, r, t, s, d, ws)
        d_target = d_orig[(index*nd + 1): m*n^order]  - r*kron(s[index,(index+1):n],kron(kron_power(s',order-1),t))*d_orig[(index-1)*nd .+ (1:m*n^(order-1))]
        @test d_target ≈ d[(index*nd + 1): m*n^order]
        drange = drange .+ nd
    end

    println("SOLVI_COMPLEX_ELIMINATE")
    drange = 1:nd
    for index = 1:2:(n-2)
        println(index)
        d_orig = copy(d)
        @time solvi_complex_eliminate!(index, n, nd, drange, order-1, r, t, s, d, ws)
        d_target = (d_orig[(index+1)*nd .+ 1: m*n^order]  - r*kron(s[index,(index+2):n],kron(kron_power(s',order-1),t))*d_orig[(index-1)*nd .+ (1:m*n^(order-1))]
                    - r*kron(s[index + 1,(index+2):n],kron(kron_power(s',order-1),t))*d_orig[(index-1)*nd + m*n^(order-1) .+ (1:m*n^(order-1))])
        @test d_target ≈ d[(index+1)*nd .+ 1: m*n^order]
        drange = drange .+ 2*nd
    end

    println("SOLVIIP_REAL_ELIMINATE")
    drange = 1:nd
    for index = 1:(n-1)
        d_orig = copy(d)
        @time solviip_real_eliminate!(index, n, nd, drange, order-1, alpha, beta, t, t2, s, s2, d, ws)
        d_target = (d_orig[(index*nd + 1): m*n^order]
                    - 2*alpha*kron(s[index,(index+1):n],kron(kron_power(s',order-1),t))*d_orig[(index-1)*nd .+ (1:m*n^(order-1))]
                    - (alpha^2 + beta^2)*kron(s2[index,(index+1):n],kron(kron_power(s2',order-1),t2))*d_orig[(index-1)*nd .+ (1:m*n^(order-1))])
        @test d_target ≈ d[(index*nd + 1): m*n^order]
        drange = drange .+ nd
    end

    println("SOLVIIP_COMPLEX_ELIMINATE")
    drange = 1:nd
    for index = 1:2:(n-2)
        d_orig = copy(d)
        @time solviip_complex_eliminate!(index, n, nd, drange, order-1, alpha, beta, t, t2, s, s2, d, ws)
        d_target = (d_orig[(index+1)*nd + 1: m*n^order]
                    - 2*alpha*kron(s[index,(index+2):n],kron(kron_power(s',order-1),t))*d_orig[(index-1)*nd .+ (1:m*n^(order-1))]
                    - (alpha^2 + beta^2)*kron(s2[index,(index+2):n],kron(kron_power(s2',order-1),t2))*d_orig[(index-1)*nd .+ (1:m*n^(order-1))]
                    - 2*alpha*kron(s[index + 1,(index+2):n],kron(kron_power(s',order-1),t))*d_orig[(index-1)*nd + m*n^(order-1) .+ (1:m*n^(order-1))]
                    - (alpha^2 + beta^2)*kron(s2[index + 1,(index+2):n],kron(kron_power(s2',order-1),t2))*d_orig[(index-1)*nd  + m*n^(order-1) .+ (1:m*n^(order-1))])
        @test d_target ≈ d[(index+1)*nd .+ 1: m*n^order]
        drange = drange .+ 2*nd
    end

    println("TRANSFORMATION1")
    a = randn()
    b1 = -rand()
    b2 = rand()
    d_orig = copy(d[1:2*m*n^(order-1)])
    @time transformation1(a,b1,b2,order-1,t,s,d,ws)
    d_target = d_orig + kron([a -b1; -b2 a],kron(kron_power(s',order-1),t))*d_orig
    @test d[1:2*m*n^(order-1)] ≈ d_target

    if order > 1
        println("TRANSFORMATION2")
        # check handling of order and nd
        a = randn()
        b1 = -rand()
        b2 = rand()
        r1 = randn()
        r2 = randn()
        d1 = randn(2*m*n^(order-1))
        nd1 = m*n^(order-2)
        d_orig = copy(d1[1:2*m*n^(order-2)])
        @time transform2(r1, r2, a, b1, b2, nd1, order-1, t, t2, s, s2, d1, ws)
        d_target = (I(2*m*n^(order-2)) + 2*r1*kron([a b1; b2 a],kron(kron_power(s',(order-2)),t))
                    + (r1*r1 + r2*r2)*kron([a b1; b2 a]*[a b1; b2 a],kron(kron_power(s2',(order-2)),t2)))*d_orig
        @test d1[1:2*m*n^(order-2)] ≈ d_target
    end

    println("SOLVE1 real eigenvalues case")
    d_orig = copy(d)
    sreal = QuasiUpperTriangular(triu(s))
    ts = triu(s)
    sreal2 = QuasiUpperTriangular(ts*ts)
    @time solve1!(r,order,t,t2,sreal,sreal2,d,ws)
    d_target = (I(m*n^order) + kron(kron_power(sreal',order),t))\d_orig
    @test d ≈ d_target

    if n == 3
        alpha = randn()
        beta1 = rand()
        beta2 = rand()
        d_orig = copy(d)

        println("SOLVEIIP beta == 0 s, with real eigenvalues")
        beta1 = 0.0
        @time solviip(alpha, beta1, order, t, t2, sreal, sreal2, d, ws)
        d_target = (I(m*n^order) + 2*alpha*kron(kron_power(sreal',order),t) + (alpha*alpha + beta1*beta1)*kron(kron_power(sreal2',order),t2))\d_orig
        @test d ≈ d_target

        scplx = copy(s)
        scplx[1:2,1:2] = [alpha -beta1; beta2 alpha]
        scplx[2,3] = 0
        scplx2 = QuasiUpperTriangular(scplx*scplx)

        println("SOLVE1 complex eigenvalues case")
        d = randn(m*n^order)
        d_orig = copy(d)
        @time solve1!(r,order,t,t2,scplx,scplx2,d,ws)
        d_target = (I(m*n^order) + kron(kron_power(scplx',order),t))\d_orig
        @test d ≈ d_target
    end

    println("SOLVEIIP with real eigenvalue")
    d_orig = copy(d)
    @time solviip(alpha, beta, order,
                                       t, t2, sreal, sreal2,
                                       d, ws)
d_target = (I(m*n^order)
            + 2*alpha*kron(kron_power(sreal',order),t)
            + (alpha*alpha + beta*beta)*kron(kron_power(sreal2',order),t2))\d_orig
@test d ≈ d_target

println("SOLVEIIP2")
d = randn(2*m*n^(order-1))
d_orig = copy(d)
a = randn()
b1 = -rand()
b2 = rand()
G = [a b1; b2 a]
@time solviip2(alpha,beta,a,b2,b1,order,t,t2,s,s2,d,ws)
d_target = (I(2*m*n^(order-1)) + 2*alpha*kron(kron(G',kron_power(s',order-1)),t)
            + (alpha*alpha + beta*beta)*kron(kron(G'*G',kron_power(s2',order-1)),t2))\d_orig
@test d ≈ d_target

println("SOLVE1 general case")
d = randn(m*n^order)
d_orig = copy(d)
@time solve1!(r,order,t,t2,s,s2,d,ws)
d_target = (I(m*n^order) + kron(kron_power(s',order),t))\d_orig
@test d ≈ d_target

end

n=9
test1(n,n,2)
test1(n,n,2)


text:
"""
T is for dispatch.
"""

struct Token{T}
    lineno :: Int64
    colno  :: Int32
    offset :: Int64
    str    :: String
    span   :: Int32
end

"""
    Token{T}(str::String; lineno::Int=0, colno::Int=0, offset::Int=0, span::Int=0)

Create a Token of type `T` with content of given `str::String`.
"""
Token{T}(str; lineno::Int=0, colno::Int=0, offset::Int=0, span::Int=0) where T =
    Token{T}(lineno, colno, offset, String(str), span)

@as_record Token


text:
struct ZCube
  ar
end


text:
# Seasonality and unobserved components
# models: an overview
# Andrew Harvey
# December 4, 2006

# Are Seasonal Patterns Constant Over Time?
# A Test for Seasonal Stability
# Fabio CANOVA
# Bruce E. HANSEN
# Journal of Business & Economic Statistics, July 1995, Vol. 13, No. 3

λ(j::Int, s::Int) = (2*pi*j)/s

const CRIT_VALS_05 = [0.470, 0.749, 1.01, 1.24, 1.47, 1.68, 1.9, 2.11, 2.32, 2.54, 2.75, 2.96]

function crit_val_05_generalized_von_mises_distribution(s::Int)
    return s > 12 ? 0.269 * s^0.928 : CRIT_VALS_05[s]
end

function π_statistic(ε::Vector{Fl}, σ::Fl, T::Int) where Fl
    return (T^(-2)*σ^(-2))*
            sum(
                sum(ε[i]*(-1)^(i) for i = 1:t) for t = 1:T
            )
end

function std_statistic(ε::Vector{Fl}, σ::Fl, T::Int, s::Int, j::Int) where Fl
    return (2*T^(-2)*σ^(-2))*
            sum(
                sum(ε[i]*cos(λ(j, s)*i) for i = 1:t)^2
                +
                sum(ε[i]*sin(λ(j, s)*i) for i = 1:t)^2 for t = 1:T
            )
end

function test_statistic(ε::Vector{Fl}, σ::Fl, T::Int, s::Int) where Fl
    # joint statistical test for all frequencies in
    # trigonometric pattern
    ω = 0
    # for each frequency in the trigonometric pattern
    for j = 1:floor(Int, s/2)
        ω += j == s/2 ? π_statistic(ε, σ, T) : std_statistic(ε, σ, T, s, j)
    end
    return ω
end

function build_std(T::Int, s::Int, j::Int)
    # add cossine and sine with frequency (2*pi*j)/s
    # to exogenous matrix X
    return hcat(cos.(λ(j, s)*(1:T)), sin.(λ(j, s)*(1:T)))
end

function build_π(T::Int, s::Int, j::Int)
    # when j = s/2
    # add cossine with frequency pi and intercept
    # to exogenous matrix X
    return hcat(cos.(λ(j, s)*(1:T)), ones(T))
end

function build_X(s::Int, T::Int, Fl::Type)
    # create the exogenous matrix
    X = Array{Fl, 2}(undef, T, s)
    # for each frequency in the trigonometric pattern
    for j = 1:floor(Int, s/2)
        X[:,[2*j-1, 2*j]] = j == s/2 ? build_π(T, s, j) : build_std(T, s, j)
    end
    return X
end

function seasonal_regression(y::Vector{Fl}, s::Int, T::Int) where Fl
    X = build_X(s, T, Fl)
    return y - X * (X \ y)
end

function seasonal_stationarity_test(y::Vector{Fl}, s::Int) where Fl
    T  = length(y)
    ε  = seasonal_regression(y, s, T)
    σ  = std(ε)
    ω  = test_statistic(ε, σ, T, s)
    crit_val = crit_val_05_generalized_von_mises_distribution(s)
    res = ω - crit_val
    res > 0 ? println("Rejected seasonal stationarity at 5% significance level"*
            " Test Statistic: $ω - Critical Value: $crit_val") :
              println("Didn't reject seasonal stationarity at 5% significance level"*
            " Test Statistic: $ω - Critical Value: $crit_val")
    return res > 0 ? false : true
end


text:
function demoPlotDAG()
#    try
        # variables:
        D=Dict{Integer,VariablePlot}()
        D[1]=VariablePlot("1")
        D[2]=VariablePlot("2")
        D[3]=VariablePlot("3")
        D[4]=VariablePlot("4")
#        D[5]=VariablePlot("5")
#        D[6]=VariablePlot("6")

        # connections:
        L=Dict{Integer,Array{Integer}}()

        if false
        L[1]=[2,4]
        L[2]=[3,5]
        L[3]=[6]
            end

        L[1]=[2,3]
        L[4]=[3]
        L[3]=[4]

        if false
        L[1]=[2,5,4] # connection from variable 1 to 2
        L[2]=[1,5,3]
        L[3]=[2,5,4]
        L[4]=[3,5,1]
        L[5]=[1,2,3,4]
        end

        if false
        L[1]=[2,3,4] # connection from variable 1 to 2
        L[2]=[4]
        L[3]=[4]
        L[4]=[5]
        L[5]=[2,1]
            end

        PlotGraph(PlaceVertices(D,L,scale=2),L)
        #PlotGraph(PlaceVertices(D,L,scale=2,UseEnergyApproach=true),L,arrowsize=0.02) # rerun this if the vertex placement is not good
#    end
end


text:
#!/usr/bin/env julia
# convert date to day of year
using Dates

doy = Dates.dayofyear(DateTime(ARGS[1], "yyyy-mm-dd"))

println(doy)

max_stars_repo_name: scienceopen/sciencedates

max_issues_repo_name: geospace-code/sciencedates

max_forks_repo_name: geospace-code/sciencedates


text:
module DiscreteExteriorCalculus

using StaticArrays: SVector, SMatrix
using LinearAlgebra: Symmetric
using UniqueVectors: UniqueVector

export Point
"""
    Point{N}(coords::SVector{N, Float64}) where N
    Point(coords::AbstractVector{<:Real})
    Point(coords::Vararg{<:Real})
    Point(b::Barycentric)
    Point(b::SimpleBarycentric)

A point in `N` dimensional space.
"""
struct Point{N}
    coords::SVector{N, Float64}
end

export Simplex
"""
    Simplex{N, K}(points::SVector{K, Point{N}}) where {N, K}
    Simplex(points::AbstractVector{Point{N}}) where N
    Simplex(points::Vararg{Point{N}}) where N
    Simplex(c::Cell)
    Simplex(s::SimpleSimplex)

A simplex of dimension `K-1` embedded in `N` dimensional space.
"""
struct Simplex{N, K}
    points::SVector{K, Point{N}}
    function Simplex{N, K}(points::SVector{K, Point{N}}) where {N, K}
        @assert 1 <= K <= N+1
        return new{N, K}(points)
    end
end

export SimpleSimplex
"""
    SimpleSimplex{N}(Vector{Point{N}}) where N
    SimpleSimplex(points::AbstractVector{Point{N}}) where N
    SimpleSimplex(s::Simplex)
    SimpleSimplex(s::Vector{SimpleBarycentric{N}}) where N

Like a `simplex` but not parametrized by dimension.
"""
struct SimpleSimplex{N}
    points::Vector{Point{N}}
end

export Metric
"""
    Metric{N, T<:SMatrix{N, N, Float64}}(mat::Symmetric{Float64, T})
    Metric(mat::T) where {N, T<:SMatrix{N, N, Float64}}
    Metric(mat::AbstractMatrix{<:Real})

A symmetric but not necessarily positive semi-definite metric.
"""
struct Metric{N, T<:SMatrix{N, N, Float64}}
    mat::Symmetric{Float64, T}
end

export Wedge
"""
    Wedge{N, K}(vectors::SVector{K, SVector{N, Float64}})
    Wedge(vectors::Vector{SVector{N, Float64}}) where N

The wedge product of `K` vectors in `N` dimensional space.
"""
struct Wedge{N, K}
    vectors::SVector{K, SVector{N, Float64}}
end

export Barycentric
"""
    Barycentric{N, K}(s::Simplex{N, K}, coords::SVector{K, Float64}) where {N, K}
    Barycentric(s::Simplex{N, K}, coords::AbstractVector{<:Real}) where {N, K}
    Barycentric(s::Simplex, coords::Vararg{<:Real})
    Barycentric(b::SimpleBarycentric)

A representation of a point with respect to a simplex using barycentric coordinates. If the
simplex `s` has vertices `Aᵢ` and `x=coords` then the point represented by
`b = Barycentric(s, coords)` is `Point(b) = ΣAᵢxᵢ`.
"""
struct Barycentric{N, K}
    simplex::Simplex{N, K}
    coords::SVector{K, Float64}
    function Barycentric{N, K}(s::Simplex{N, K}, coords::SVector{K, Float64}) where {N, K}
        @assert sum(coords) ≈ 1
        if sum(coords) != 1
            coords = SVector{K, Float64}([coords[1:end-1]..., 1 - sum(coords[1:end-1])])
        end
        return new{N, K}(s, coords)
    end
end

export SimpleBarycentric
"""
    SimpleBarycentric{N}(s::SimpleSimplex{N}, coords::Vector{Float64}) where N
    SimpleBarycentric(s::SimpleSimplex{N}, coords::AbstractVector{<:Real}) where N
    SimpleBarycentric(b::Barycentric)

Like a `Barycentric` but not parametrized by dimension.
"""
struct SimpleBarycentric{N}
    simplex::SimpleSimplex{N}
    coords::Vector{Float64}
    function SimpleBarycentric{N}(simplex::SimpleSimplex{N}, coords::Vector{Float64}) where N
        @assert length(simplex.points) == length(coords)
        @assert sum(coords) ≈ 1
        if sum(coords) != 1
            coords = [coords[1:end-1]..., 1 - sum(coords[1:end-1])]
        end
        return new{N}(simplex, coords)
    end
end

export Cell
"""
    Cell{N}(children::Vector{Cell{N}}, parents::Dict{Cell{N}, Bool},
        points::Vector{Point{N}}, K::Int) where N
    Cell(points::AbstractVector{Point{N}}, K::Int) where N
    Cell(s::Simplex{N, K}) where {N, K}

A polytope of dimension `K-1` embedded in `N` dimensional space.
"""
struct Cell{N}
    children::Vector{Cell{N}}
    parents::Dict{Cell{N}, Bool}
    points::Vector{Point{N}}
    K::Int
    function Cell{N}(children::Vector{Cell{N}}, parents::Dict{Cell{N}, Bool},
                     points::Vector{Point{N}}, K::Int) where N
        @assert 1 <= K <= N+1
        for c in children
            @assert c.K == K-1
        end
        for c in keys(parents)
            @assert c.K == K+1
        end
        return new{N}(children, parents, points, K)
    end
end

export CellComplex
"""
    CellComplex{N,K}(cells::SVector{K, UniqueVector{Cell{N}}}) where {N, K}
    CellComplex{N,K}() where {N, K}
    CellComplex{N,K}(cells::AbstractVector{Cell{N}}) where {N, K}
    CellComplex(simplices::AbstractVector{Simplex{N, K}})

A cell complex with cells of dimension `0` through `K-1` embedded in `N` dimensional space.
"""
struct CellComplex{N, K}
    cells::SVector{K, UniqueVector{Cell{N}}}
    function CellComplex{N,K}(cells::SVector{K, UniqueVector{Cell{N}}}) where {N, K}
        @assert 1 <= K <= N+1
        for k in 1:K
            for c in cells[k]
                @assert c.K == k
            end
        end
        return new{N,K}(cells)
    end
end

SignedSimpleSimplex{N} = Tuple{SimpleSimplex{N}, Bool}

export TriangulatedComplex
"""
    TriangulatedComplex{N, K}(complex::CellComplex{N, K},
        simplices::Dict{Cell{N}, Vector{Tuple{SimpleSimplex{N}, Bool}})
    TriangulatedComplex{N,K}() where {N, K}
    TriangulatedComplex(simplices::AbstractVector{Simplex{N, K}}) where {N, K}

A cell complex along with a decomposition of each cell into simplices. The decomposition
comes with signs which are used to compute a signed volume for each cell.
"""
struct TriangulatedComplex{N, K}
    complex::CellComplex{N, K}
    simplices::Dict{Cell{N}, Vector{SignedSimpleSimplex{N}}}
end

export Mesh
"""
    Mesh{N,K}(primal::TriangulatedComplex{N,K}, dual::TriangulatedComplex{N,K})
    Mesh(tcomp::TriangulatedComplex{N}, center::Function) where N

A pair of TriangulatedComplexes representing a primal and dual pair. The second constructor
takes a primal TriangulatedComplex and a function `center` that takes a Simplex{N, K} to a
Barycentric{N, K} (e.g. the circumcenter or centroid) and returns the corresponding dual
TriangulatedComplex.
"""
struct Mesh{N,K}
    primal::TriangulatedComplex{N,K}
    dual::TriangulatedComplex{N,K}
end

include("utils.jl")
include("metrics.jl")
include("simplices.jl")
include("cell_complex.jl")
include("components_and_boundary.jl")
include("orientation.jl")
include("mesh.jl")
include("operators.jl")
include("constraints.jl")
end


text:
"""
SamplerState for preconditioned stochastic gradient Langevin dynamics with Adan following

Li, C., Chen, C., Carlson, D., & Carin., L. (2016). Preconditioned Stochastic Gradient Langevin Dynamics for Deep Neural Networks. AAAI. Retrieved from http://arxiv.org/abs/1512.07666

Fields:

    - x: the current state
    - preconditioner: current preconditioner
    - beta: constant for moving average

    - t: iteration counter
    - stepsize: s.stepsize(s.t) gives the stepsize to be used in iteration s.t
    - lambda: minimum mass
    - betat: current debiasing factor

Constructors:

    - construct using `SGLDAdamState(x,stepsize::Function)` or `SGLDAdamState(x,stepsize::Float64)` for a fixed stepsize
"""
type SGLDAdamState <: SamplerState
  x::Array{Float64,1}
  preconditioner::Array{Float64,1}
  beta::Float64

  t::Int64
  stepsize::Function
  lambda::Float64
  betat::Float64

  function SGLDAdamState(x::Array{Float64,1},
    stepsize::Float64;
    beta = 0.9,
    preconditioner=zeros(size(x)),
    lambda=1e-8)
    @assert size(x) == size(preconditioner)
    s(niters) = stepsize
    new(x,preconditioner,beta,0,s,lambda)
  end
  function SGLDAdamState(x::Array{Float64,1},
    stepsize::Function;
    beta =0.9,
    preconditioner=zeros(size(x)),
    lambda=1e-8)
    @assert size(x) == size(preconditioner)
    new(x,preconditioner,beta, 0,stepsize,lambda)
  end
end

function sample!(s::SGLDAdamState, grad::Function)
    # increment iteration
    s.t += 1
    s.betat *= s.beta

    # improve readability
    t = s.t
    x = s.x
    beta = s.beta
    betat = s.betat
    lambda = s.lambda
    precond = s.preconditioner
    stepsize = s.stepsize(t)

    # update preconditioner
    g = grad(x)
    if t == 1
        precond[:] = g.*g
    else
        precond *= beta
        add_inplace!((1-beta),g.*g,precond)
    end
    # Adam debiasing
    m = lambda + sqrt(precond/(1-betat))
    # gradient step (inplace addition)
    add_inplace!(stepsize,g./m,x)

    # injected noise (inplace addition)
    add_inplace!(sqrt(2.*stepsize),randn(length(x))./sqrt(m),x)
    s
end


text:
# Reactive message-passing in linear Gaussian dynamical system
# Experiment to pass messages without a schedule
#
# Wouter Kouw, BIASlab
# 12-08-2019

using Random
using Revise
using CPUTime
using ForneyLab
using Plots
pyplot()

# Functions for generating data
include(joinpath(@__DIR__, "../gen_data.jl"))

"""
Experiment parameters
"""

# Signal time horizon
T = 30

# Reaction-time clock
TT = 10

# Known transition and observation matrices
transition_coefficient = 0.6
likelihood_coefficient = 1.0

# Noise parameters
process_noise = 0.5
measurement_noise = 0.3

# Clamped parameters
x0 = [0.1, 1.0]
a0 = [1.0, 1.0]

# Generate data
observed, hidden = gendata_LGDS(transition_coefficient,
                                likelihood_coefficient,
                                process_noise,
                                measurement_noise,
                                x0[1],
                                x0[2],
                                time_horizon=T)

"""
Model specification
"""

# Define composite transition node
@composite TransitionGaussian (x_t, a, x_tmin) begin
    @RV x_tmin_ = x_tmin*a
    @RV x_t ~ GaussianMeanVariance(x_tmin_, process_noise)
end

@NaiveVariationalRule(:node_type     => TransitionGaussian,
                      :outbound_type => Message{GaussianMeanPrecision},
                      :inbound_types => (Nothing, Message{Gaussian}, Message{Gaussian}))

# Start timer
CPUtic()

g = FactorGraph()

@RV a ~ GaussianMeanVariance(a0[1], a0[2])
@RV x_0 ~ GaussianMeanVariance(x0[1], x0[2])

global x = Vector{Variable}(undef, T)
global y = Vector{Variable}(undef, T)
global x_tmin = x_0
for t = 1:T

    # @RV x_tmin_ = a*x_tmin
    # @RV x[t] ~ GaussianMeanVariance(x_tmin_, process_noise)

    @RV x[t] ~ TransitionGaussian(x_tmin, a)
    @RV y[t] ~ GaussianMeanVariance(x[t], measurement_noise)

    global x_tmin = x[t]

    placeholder(y[t], :y, index=t)
end

"""
Mean-field algorithm generation
"""

# Define recognition factorization
RecognitionFactorization()

q_a = RecognitionFactor(a)
q_x_0 = RecognitionFactor(x_0)

q_x = Vector{RecognitionFactor}(undef, T)
for t=1:T
    q_x[t] = RecognitionFactor(x[t])
end

# Compile algorithm
algo_a_mf = variationalAlgorithm(q_a, name="AMF")
algo_x_mf = variationalAlgorithm([q_x_0; q_x], name="SMF")
algo_F_mf = freeEnergyAlgorithm(name="MF");

"""
Run inference algorithm
"""

# Load algorithm
eval(Meta.parse(algo_a_mf))
eval(Meta.parse(algo_x_mf))
eval(Meta.parse(algo_F_mf))

# Initialize data
global data = Dict(:y => observed)

# Initial recognition distributions
global marginals_mf = Dict{Symbol, ProbabilityDistribution}()
global marginals_mf[:a] = vague(GaussianMeanPrecision)
for t = 0:T
    global marginals_mf[:s_*t] = vague(GaussianMeanPrecision)
end

# Run algorithm
global estimated_states = zeros(T,2,TT)
global F_mf = Vector{Float64}(undef, TT)
for tt = 1:TT
    stepAMF!(data, marginals_mf)
    stepSMF!(data, marginals_mf)

    for t = 1:T
        estimated_states[t,1,tt] = mean(marginals_mf[:x_*t])
        estimated_states[t,2,tt] = sqrt(var(marginals_mf[:x_*t]))
    end

    global F_mf[tt] = freeEnergyMF(data, marginals_mf)
end

CPUtoc()

"""
Visualize experimental results
"""

# Visualize final estimates over time
plot(hidden[2:end], color="red", label="states")
plot!(estimated_states[:,1,end], color="blue", label="estimates")
plot!(estimated_states[:,1,end],
      ribbon=[estimated_states[:,2,end], estimated_states[:,2,end]],
      linewidth=2,
      color="blue",
      fillalpha=0.2,
      fillcolor="blue",
      label="")
scatter!(observed, color="black", label="observations")
savefig(joinpath(@__DIR__, "viz/state_estimates_baseline.png"))

max_stars_repo_name: wmkouw/schedulefreeMP

max_issues_repo_name: wmkouw/schedulefreeMP

max_forks_repo_name: wmkouw/schedulefreeMP


text:
using JuMP, LinearAlgebra, Suppressor, StatsBase, CSV, DataFrames, Compat, Test, Random#, CPLEX
include("alt_min.jl")
include("convexpenalties.jl")

results_template = DataFrame(n=Int[], r=Int[], p=Real[], err_GD_msk=Real[], err_GD_msk_am=Real[], t_gd_msk=Real[], t_gd_msk_am=Real[])
results_exact=similar(results_template, 0)

for ARG in ARGS
    array_num = parse(Int, ARG)
    for n in [50*array_num] #1-12
        p=0.25
        r=5
        gamma=20/p
        M_L=randn(n, r)
        M_R=randn(n, r)
        M=M_L*M_R'

        numEntriesSampled=Int(round(p*n^2))
        coords=[repeat(1:n, inner=size(1:n,1)) repeat(1:n, outer=size(1:n,1))];
        entries=sample(1:n^2, Int(numEntriesSampled), replace=false)

        theIndices=zeros(n,n)
        for t=1:Int(numEntriesSampled)
            theIndices[coords[entries[t],1], coords[entries[t],2]]=1.0
        end

        k=r
        t_GD_msk=@elapsed X_GD,LB, UB, Y_relax=getMatrixFrobeniusNorm_greedy(coords, entries, M, n, numEntriesSampled, k, gamma)

        @show "SDP lower bound is: " LB
        u,sigma,v,=svd(X_GD)
        @show rnk_GD=sum(sigma.>1e-2);
        @show err_GD_msk = norm(X_GD.-M)/norm(M);

        t_GD_msk_am=@elapsed X_GD=alternatingminimization_matrixcompletion_solvesdprelax(coords, entries, M, n, numEntriesSampled, k, gamma)
        @show t_GD_msk_am
        u,sigma,v,=svd(X_GD)
        @show rnk_GD=sum(sigma.>1e-2);
        @show err_GD_am = norm(X_GD.-M)/norm(M);

        push!(results_exact, [n, r, p, err_GD_msk, err_GD_am,t_GD_msk, t_GD_msk_am])

    end
end

CSV.write("results_evaluatescalability.csv", results_exact, append=true)

max_stars_repo_path: Matrix Completion/createTable3.jl

max_stars_repo_name: ryancorywright/MixedProjectionSoftware

max_issues_repo_path: Matrix Completion/createTable3.jl

max_issues_repo_name: ryancorywright/MixedProjectionSoftware

max_forks_repo_path: Matrix Completion/createTable3.jl

max_forks_repo_name: ryancorywright/MixedProjectionSoftware


text:
module Minerva

using Plots, Juno

include("environment.jl")
include("environments/cartpole.jl")
include("environments/mountaincar.jl")

export RandomWalk,
       CartPole,
       MountainCar,
       reset!,
       step!,
       render

include("agent.jl")
include("agents/deep_q.jl")
include("agents/vanilla_pg.jl")

export DeepQAgent,
       VanillaPG,
       action,
       train!,
       evaluate

include("schedulers.jl")

export StepDecay,
       ExpDecay

include("recorders.jl")

export ExperimentRecorder,
       stamp!

end # module


text:
type NEOSServer
	useragent::ASCIIString
	host::ASCIIString
	contenttype::ASCIIString
	email::ASCIIString
	NEOSServer(;email="") = new("JuliaXMLRPC", "http://neos-server.org:3332", "text/xml", email)
end

function addemail!(s::NEOSServer, email::ASCIIString)
	s.email = email
end

type NEOSJob
	number::Int64
	password::ASCIIString
end

function _method(s::NEOSServer, name::ASCIIString, args...)
	xml = XMLDocument()
	mthd = 	create_root(xml, "methodCall")
	mname = new_child(mthd, "methodName")
	add_text(mname, name)
	if length(args) > 0
		params = new_child(mthd, "params")
		for a in args
			param = new_child(params, "param")
			value = new_child(param, "value")
			if isa(a, ASCIIString)
				add_text(value, a)
			else
				myint = new_child(value, "int")
				add_text(myint, string(a))
			end
		end
	end
	return _send(s, string(xml))
end

function _send(s::NEOSServer, xml::ASCIIString)
	hdrs = Dict{AbstractString, AbstractString}("user-agent" => s.useragent, "host" => s.host,
	"content-type" => s.contenttype, "content-length" => string(length(xml)))
	res = post(s.host; headers=hdrs, data=xml)
	if res.status == 200
		return extractResponse(res.data)
	else
		error("XML-RPC failed")
	end
end

function getValues!(values, c)
	types = ["int", "i4", "string", "double", "base64", "dateTime.iso8601"]
	for d in child_nodes(c)
		if name(d) in types
			push!(values, content(d))
		else
			getValues!(values, d)
		end
	end
end

function extractResponse(s)
	parameters = Array(Any, 0)
	xml = parse_string(convert(ASCIIString, s))
	xroot = root(xml)
	getValues!(parameters, xroot)
	return parameters
end

# ===========================================
#
# 	NEOS API Methods
#
function neosHelp(s::NEOSServer)
	return _method(s, "help")[1]
end

messages = [:emailHelp, :welcome, :version, :ping, :printQueue]
for m in messages
	m_str = string(m)
	@eval 	function ($m)(s::NEOSServer)
				return _method(s, $(m_str))[1]
			end
end

lists = [:listAllSolvers, :listCategories]
for m in lists
	m_str = string(m)
	@eval 	function ($m)(s::NEOSServer)
				return _method(s, $(m_str))
			end
end

function getSolverTemplate(s::NEOSServer, category::Symbol, solvername::Symbol, inputMethod::Symbol)
	_method(s, "getSolverTemplate", string(category), string(solvername), string(inputMethod))[1]
end

function listSolversInCategory(s::NEOSServer, category::Symbol)
	_method(s, "listSolversInCategory", string(category))
end

function submitJob(s::NEOSServer, xmlstring::ASCIIString)
	res = _method(s, "submitJob", xmlstring)
	println("===================")
	println("NEOS Job submitted")
	println("number:\t$(res[1])")
	println("pwd:\t$(res[2])")
	println("==================")
	return NEOSJob(parse(Int, res[1]), res[2])
end

job_methods = [:getJobStatus, :killJob]
for m in job_methods
	m_str = string(m)
	@eval 	function ($m)(s::NEOSServer, j::NEOSJob)
				return _method(s, $(m_str), j.number, j.password)[1]
			end
end

function getJobInfo(s::NEOSServer, j::NEOSJob)
	_method(s, "getJobInfo", j.number, j.password)
end

for s in ["", "NonBlocking"]
	_final_str = "getFinalResults$s"
	_final = symbol(_final_str)
	_intermediate_str = "getIntermediateResults$s"
	_intermediate = symbol(_intermediate_str)
	@eval function ($_final)(s::NEOSServer, j::NEOSJob)
		decode_to_string(_method(s, $(_final_str), j.number, j.password)[1])
	end

	@eval function ($(_intermediate))(s::NEOSServer, j::NEOSJob; offset=0)
		decode_to_string(_method(s, $(_intermediate_str), j.number, j.password, offset)[1])
	end
end

function decode_to_string(s)
	bytestring(decode(Base64, replace(s, "\n", "")))
end


text:
const year_seconds = 31556926

using OrdinaryDiffEq
using DiffEqCallbacks
using Printf
using Plots

function odefun(dψV, ψδ, p, t)
  reject_step = p.reject_step
  Vp = p.Vp
  lop = p.lop
  F = p.F
  u = p.u
  Δτ = p.Δτ
  ge = p.ge
  μshear = p.μshear
  RSa = p.RSa
  RSb = p.RSb
  σn = p.σn
  η = p.η
  RSV0 = p.RSV0
  τz0 = p.τz0
  RSDc = p.RSDc
  RSf0 = p.RSf0
  LFtoB = p.LFtoB

  if reject_step[1]
    return
  end
  δNp = div(length(ψδ), 2)
  ψ  = @view ψδ[        (1:δNp) ]
  δ  = ψδ[ δNp .+ (1:δNp) ]

  bc_Dirichlet = (lf, x, y, e) -> (2-lf)*(δ ./ 2) + (lf-1)*fill(t * Vp/2, size(x))
  bc_Neumann   = (lf, x, y, nx, ny, e) -> zeros(size(x))

  # solve for displacements everywhere in domain
  e = 1
  locbcarray_mod!(ge, lop[e], LFtoB, bc_Dirichlet, bc_Neumann,(e))
  u[:] = F \ ge

  # set up rates of change for  state and slip
  dψ = @view dψV[       (1:δNp) ]
  V  = @view dψV[δNp .+ (1:δNp) ]

  dψ .= 0
  V  .= 0

  # Update the fault data
  Δτ .= 0
  lf1 = 1

  Δτ .= -μshear .* computetraction_mod(lop[1], lf1, u, δ)

  #δex = Vp*t/2
  #τ_ex = μshear * (Vp * t/2 -δex/2)/20
  Vex = zeros(δNp)
  #yf = lop[1].facecoord[2][1]
  #Vtest = Vp/2
  # update velocity
#@show t./31556926

  for n = 1:δNp
    ψn = ψ[n]
    an = RSa[n]

    τn = Δτ[n] + τz0
    if isnan(τn)
      println("τ reject")
      reject_step[1] = true
      return
    end

    VR = abs(τn / η)
    VL = -VR
    Vn = V[n]
    obj_rs(V) = rateandstate(V, ψn, σn, τn, η, an, RSV0)
    (Vn, _, iter) = newtbndv(obj_rs, VL, VR, Vn; ftol = 1e-9,
                                 atolx = 1e-9, rtolx = 1e-9)

    Vex[n] = 2*RSV0*sinh(τn/(σn * an)) * exp(-ψn / an)

    if isnan(Vn) || iter < 0
      println("V reject")
      reject_step[1] = true
      return
          #error()
    end
    V[n] = Vn

    dψ[n] = (RSb * RSV0 / RSDc) * (exp((RSf0 - ψn) / RSb) - abs(Vn) / RSV0)
    if !isfinite(dψ[n])
      println("ψ reject")
      dψ[n] = 0
      reject_step[1] = true
      return
    end
  end
  #@show δ[end]
  #sleep(2)
  #@show [maximum(V) minimum(ψ) maximum(τ)]
  #pyplot()

  #display(plot(yf, δ))
  #sleep(1)
 #@show [minimum(V) maximum(V) t/31556926]
 #sleep(.1)
#  @show t/31556926
#sleep(1)
  nothing
end

max_stars_repo_name: brittany-erickson/HybridSBP

max_issues_repo_name: brittany-erickson/HybridSBP

max_forks_repo_name: brittany-erickson/HybridSBP


text:
using LLVM
using ObjectFile
using Libdl

module FFI
    using LLVM
    module BLASSupport
        # TODO: LAPACK handling
        using LinearAlgebra
        using ObjectFile
        using Libdl
        if VERSION >= v"1.7"
            function __init__()
                if VERSION > v"1.8"
                  global blas_handle = Libdl.dlopen(BLAS.libblastrampoline)
                else
                  global blas_handle = Libdl.dlopen(BLAS.libblas)
                end
            end
            function get_blas_symbols()
                symbols = BLAS.get_config().exported_symbols
                if BLAS.USE_BLAS64
                    return map(n->n*"64_", symbols)
                end
                return symbols
            end

            function lookup_blas_symbol(name)
                Libdl.dlsym(blas_handle::Ptr{Cvoid}, name; throw_error=false)
            end
        else
            function __init__()
                global blas_handle = Libdl.dlopen(BLAS.libblas)
            end
            function get_blas_symbols()
                symbols = Set{String}()
                path = Libdl.dlpath(BLAS.libblas)
                ignoreSymbols = Set(String["", "edata", "_edata", "end", "_end", "_bss_start", "__bss_start", ".text", ".data"])
                for s in Symbols(readmeta(open(path, "r")))
                    name = symbol_name(s)
                    if !Sys.iswindows() && BLAS.vendor() == :openblas64
                        endswith(name, "64_") || continue
                    else
                        endswith(name, "_") || continue
                    end
                    if !in(name, ignoreSymbols)
                        push!(symbols, name)
                    end
                end
                symbols = collect(symbols)
                if Sys.iswindows() &&  BLAS.vendor() == :openblas64
                    return map(n->n*"64_", symbols)
                end
                return symbols
            end

            function lookup_blas_symbol(name)
                Libdl.dlsym(blas_handle::Ptr{Cvoid}, name; throw_error=false)
            end
        end
    end

    const ptr_map = Dict{Ptr{Cvoid},String}()

    function __init__()
        known_names = (
            "jl_alloc_array_1d", "jl_alloc_array_2d", "jl_alloc_array_3d",
            "ijl_alloc_array_1d", "ijl_alloc_array_2d", "ijl_alloc_array_3d",
            "jl_new_array", "jl_array_copy", "jl_alloc_string",
            "jl_in_threaded_region", "jl_enter_threaded_region", "jl_exit_threaded_region", "jl_set_task_tid", "jl_new_task",
            "malloc", "memmove", "memcpy", "jl_array_grow_beg", "jl_array_grow_end", "jl_array_grow_at", "jl_array_del_beg",
            "jl_array_del_end", "jl_array_del_at", "jl_array_ptr", "jl_value_ptr", "jl_get_ptls_states", "jl_gc_add_finalizer_th",
            "jl_symbol_n", "jl_", "jl_object_id"
        )
        for name in known_names
            sym = LLVM.find_symbol(name)
            if sym == C_NULL
                continue
            end
            if haskey(ptr_map, sym)
                # On MacOS memcpy and memmove seem to collide?
                if name == "memcpy"
                    continue
                end
            end
            @assert !haskey(ptr_map, sym)
            ptr_map[sym] = name
        end
        for sym in BLASSupport.get_blas_symbols()
            ptr = BLASSupport.lookup_blas_symbol(sym)
            if ptr !== nothing
                if haskey(ptr_map, ptr)
                    if ptr_map[ptr] != sym
                        @warn "Duplicated symbol in ptr_map" ptr, sym, ptr_map[ptr]
                    end
                    continue
                end
                ptr_map[ptr] = sym
            end
        end
    end

    function memoize!(ptr, fn)
        fn = get(ptr_map, ptr, fn)
        if !haskey(ptr_map, ptr)
            ptr_map[ptr] = fn
        else
            @assert ptr_map[ptr] == fn
        end
        return fn
    end
end

import GPUCompiler: IRError, InvalidIRError

function restore_lookups(mod::LLVM.Module)
    i64 = LLVM.IntType(64; ctx=context(mod))
    for (v, k) in FFI.ptr_map
        if haskey(functions(mod), k)
            f = functions(mod)[k]
            replace_uses!(f, LLVM.Value(LLVM.API.LLVMConstIntToPtr(ConstantInt(i64, convert(Int, v)), llvmtype(f))))
            unsafe_delete!(mod, f)
        end
    end
end

function check_ir(job, mod::LLVM.Module)
    errors = check_ir!(job, IRError[], mod)
    unique!(errors)
    if !isempty(errors)
        throw(InvalidIRError(job, errors))
    end
end

function check_ir!(job, errors, mod::LLVM.Module)
    imported = Set(String[])
    if haskey(functions(mod), "malloc")
        f = functions(mod)["malloc"]
        name!(f, "")
        ctx = context(mod)
        ptr8 = LLVM.PointerType(LLVM.IntType(8; ctx))

        prev_ft = eltype(llvmtype(f)::LLVM.PointerType)::LLVM.FunctionType

        mfn = LLVM.API.LLVMAddFunction(mod, "malloc", LLVM.FunctionType(ptr8, parameters(prev_ft)))
        replace_uses!(f, LLVM.Value(LLVM.API.LLVMConstPointerCast(mfn, llvmtype(f))))
        unsafe_delete!(mod, f)
    end
    for f in collect(functions(mod))
        check_ir!(job, errors, imported, f)
    end

    return errors
end

function check_ir!(job, errors, imported, f::LLVM.Function)
    calls = []
    isInline = API.EnzymeGetCLBool(cglobal((:EnzymeInline, API.libEnzyme))) != 0
    for bb in blocks(f), inst in instructions(bb)
        if isa(inst, LLVM.CallInst)
            push!(calls, inst)
        # remove illegal invariant.load and jtbaa_const invariants
        elseif isInline && isa(inst, LLVM.LoadInst)
            md = metadata(inst)
            if haskey(md, LLVM.MD_tbaa)
                modified = LLVM.Metadata(ccall((:EnzymeMakeNonConstTBAA, API.libEnzyme), LLVM.API.LLVMMetadataRef, (LLVM.API.LLVMMetadataRef,), md[LLVM.MD_tbaa]))
                setindex!(md, modified, LLVM.MD_tbaa)
            end
            if haskey(md, LLVM.MD_invariant_load)
                delete!(md, LLVM.MD_invariant_load)
            end
        end
    end

    while length(calls) > 0
        inst = pop!(calls)
        check_ir!(job, errors, imported, inst, calls)
    end
    return errors
end

const libjulia = Ref{Ptr{Cvoid}}(C_NULL)

import GPUCompiler: DYNAMIC_CALL, DELAYED_BINDING, RUNTIME_FUNCTION, UNKNOWN_FUNCTION, POINTER_FUNCTION
import GPUCompiler: backtrace, isintrinsic
function check_ir!(job, errors, imported, inst::LLVM.CallInst, calls)
    bt = backtrace(inst)
    dest = called_value(inst)
    if isa(dest, LLVM.Function)
        fn = LLVM.name(dest)

        # some special handling for runtime functions that we don't implement
        if fn == "jl_get_binding_or_error"
        elseif fn == "jl_invoke"
        elseif fn == "jl_apply_generic"
        elseif fn == "gpu_malloc"
            ofn = LLVM.parent(LLVM.parent(inst))
            mod = LLVM.parent(ofn)
            ctx = context(mod)

            b = Builder(ctx)
            position!(b, inst)

            mfn = LLVM.API.LLVMGetNamedFunction(mod, "malloc")
            if mfn == C_NULL
                ptr8 = LLVM.PointerType(LLVM.IntType(8; ctx))
                mfn = LLVM.API.LLVMAddFunction(mod, "malloc", LLVM.FunctionType(ptr8, [llvmtype(LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(inst, 0)))]))
            end
            mfn2 = LLVM.Function(mfn)
            nval = ptrtoint!(b, call!(b, mfn2, [LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(inst, 0))]), llvmtype(inst))
            replace_uses!(inst, nval)
            LLVM.API.LLVMInstructionEraseFromParent(inst)
        elseif fn == "jl_load_and_lookup"
            ofn = LLVM.parent(LLVM.parent(inst))
            mod = LLVM.parent(ofn)
            ctx = context(mod)

            flib = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(inst, 0))
            if isa(flib, LLVM.ConstantExpr)
                flib = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(flib, 0))
            end
            if isa(flib, LLVM.GlobalVariable)
                flib = LLVM.initializer(flib)
            end
            if isa(flib, LLVM.ConstantArray) && eltype(llvmtype(flib)) == LLVM.IntType(8; ctx)
                flib = String(map((x)->convert(UInt8, x), collect(flib)[1:(end-1)]))
            end

            fname = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(inst, 1))
            if isa(fname, LLVM.ConstantExpr)
                fname = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(fname, 0))
            end
            if isa(fname, LLVM.GlobalVariable)
                fname = LLVM.initializer(fname)
            end
            if isa(fname, LLVM.ConstantArray) && eltype(llvmtype(fname)) == LLVM.IntType(8; ctx)
                fname = String(map((x)->convert(UInt8, x), collect(fname)[1:(end-1)]))
            end
            hnd = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(inst, 2))
            if isa(hnd, LLVM.GlobalVariable)
                hnd = LLVM.name(hnd)
            end

            if !isa(hnd, String) || !isa(fname, String) || !isa(flib, String)
                push!(errors, ("jl_load_and_lookup", bt, nothing))
                return
            end
            # res = ccall(:jl_load_and_lookup, Ptr{Cvoid}, (Cstring, Cstring, Ptr{Cvoid}), flib, fname, cglobal(Symbol(hnd)))
            push!(errors, ("jl_load_and_lookup", bt, nothing))
        elseif fn == "jl_lazy_load_and_lookup" || fn == "ijl_lazy_load_and_lookup"
            ofn = LLVM.parent(LLVM.parent(inst))
            mod = LLVM.parent(ofn)
            ctx = context(mod)

            flib = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(inst, 0))
            if isa(flib, LLVM.LoadInst)
                op = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(flib, 0))
                if isa(op, LLVM.ConstantExpr)
                    op1 = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(op, 0))
                    if isa(op1, LLVM.ConstantExpr)
                        op2 = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(op1, 0))
                        if isa(op2, ConstantInt)
                            rep = reinterpret(Ptr{Cvoid}, convert(Csize_t, op2)+8)
                            ld = unsafe_load(convert(Ptr{Ptr{Cvoid}}, rep))
                            flib = Base.unsafe_pointer_to_objref(ld)
                        end
                    end
                end
            end

            fname = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(inst, 1))
            if isa(fname, LLVM.ConstantExpr)
                fname = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(fname, 0))
            end
            if isa(fname, LLVM.GlobalVariable)
                fname = LLVM.initializer(fname)
            end
            if isa(fname, LLVM.ConstantArray) && eltype(llvmtype(fname)) == LLVM.IntType(8; ctx)
                fname = String(map((x)->convert(UInt8, x), collect(fname)[1:(end-1)]))
            end

            if !isa(fname, String) || !isa(flib, String)
                push!(errors, ("jl_lazy_load_and_lookup", bt, nothing))
                return
            end

            data = open(flib, "r") do io
                lib = readmeta(io)
                sections = Sections(lib)
                if !(".llvmbc" in sections)
                    return nothing
                end
                llvmbc = read(findfirst(sections, ".llvmbc"))
                return llvmbc
            end

            found = false
            if data !== nothing
                inmod = parse(LLVM.Module, data; ctx)
                found = haskey(functions(inmod), fname)
            end

            if found
                if !(fn in imported)
                    internalize = String[]
                    for fn in functions(inmod)
                        if !isempty(LLVM.blocks(fn))
                            push!(internalize, name(fn))
                        end
                    end
                    for g in globals(inmod)
                        linkage!(g, LLVM.API.LLVMExternalLinkage)
                    end
                    # override libdevice's triple and datalayout to avoid warnings
                    triple!(inmod, triple(mod))
                    datalayout!(inmod, datalayout(mod))
                    GPUCompiler.link_library!(mod, inmod)
                    for n in internalize
                        linkage!(functions(mod)[n], LLVM.API.LLVMInternalLinkage)
                    end
                    push!(imported, fn)
                end
                replaceWith = functions(mod)[fname]

                for u in LLVM.uses(inst)
                    st = LLVM.user(u)
                    if isa(st, LLVM.StoreInst) && LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(st, 0)) == inst
                        ptr = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(st, 1))
                        for u in LLVM.uses(ptr)
                            ld = LLVM.user(u)
                            if isa(ld, LLVM.LoadInst)
                                b = Builder(ctx)
                                position!(b, ld)
                                replace_uses!(ld, LLVM.pointercast!(b, replaceWith, llvmtype(inst)))
                            end
                        end
                    end
                end

                b = Builder(ctx)

                position!(b, inst)
                replace_uses!(inst, LLVM.pointercast!(b, replaceWith, llvmtype(inst)))
                LLVM.API.LLVMInstructionEraseFromParent(inst)

            else
                if fn == "jl_lazy_load_and_lookup"
                    res = ccall(:jl_lazy_load_and_lookup, Ptr{Cvoid}, (Any, Cstring), flib, fname)
                else
                    res = ccall(:ijl_lazy_load_and_lookup, Ptr{Cvoid}, (Any, Cstring), flib, fname)
                end
                replaceWith = LLVM.ConstantInt(LLVM.IntType(64; ctx), reinterpret(UInt64, res))
                for u in LLVM.uses(inst)
                    st = LLVM.user(u)
                    if isa(st, LLVM.StoreInst) && LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(st, 0)) == inst
                        ptr = LLVM.Value(LLVM.LLVM.API.LLVMGetOperand(st, 1))
                        for u in LLVM.uses(ptr)
                            ld = LLVM.user(u)
                            if isa(ld, LLVM.LoadInst)
                                b = Builder(ctx)
                                position!(b, ld)
                                for u in LLVM.uses(ld)
                                    u = LLVM.user(u)
                                    if isa(u, LLVM.CallInst)
                                        push!(calls, u)
                                    end
                                end
                                replace_uses!(ld, LLVM.inttoptr!(b, replaceWith, llvmtype(inst)))
                            end
                        end
                    end
                end

                b = Builder(ctx)
                position!(b, inst)
                replacement = LLVM.inttoptr!(b, replaceWith, llvmtype(inst))
                            for u in LLVM.uses(inst)
                                u = LLVM.user(u)
                                if isa(u, LLVM.CallInst)
                                    push!(calls, u)
                                end
                                if isa(u, LLVM.PHIInst)
                                    if all(x->first(x) == inst || first(x) == replacement, LLVM.incoming(u))

                                        for u in LLVM.uses(u)
                                            u = LLVM.user(u)
                                            if isa(u, LLVM.CallInst)
                                                push!(calls, u)
                                            end
                                            if isa(u, LLVM.BitCastInst)
                                                for u1 in LLVM.uses(u)
                                                    u1 = LLVM.user(u1)
                                                    if isa(u1, LLVM.CallInst)
                                                        push!(calls, u1)
                                                    end
                                                end
                                                replace_uses!(u, LLVM.inttoptr!(b, replaceWith, llvmtype(u)))
                                            end
                                        end
                                    end
                                end
                            end
                replace_uses!(inst, replacement)
                LLVM.API.LLVMInstructionEraseFromParent(inst)
            end
        end

    elseif isa(dest, InlineAsm)
        # let's assume it's valid ASM

    elseif isa(dest, ConstantExpr)
        # Enzyme should be able to handle these
        # detect calls to literal pointers and replace with function name, if possible
        if occursin("inttoptr", string(dest))
            # extract the literal pointer
            ptr_arg = first(operands(dest))
            GPUCompiler.@compiler_assert isa(ptr_arg, ConstantInt) job
            ptr_val = convert(Int, ptr_arg)
            ptr = Ptr{Cvoid}(ptr_val)

            # look it up in the Julia JIT cache
            frames = ccall(:jl_lookup_code_address, Any, (Ptr{Cvoid}, Cint,), ptr, 0)

            if length(frames) >= 1
                if VERSION >= v"1.4.0-DEV.123"
                    fn, file, line, linfo, fromC, inlined = last(frames)
                else
                    fn, file, line, linfo, fromC, inlined, ip = last(frames)
                end

                # Remember pointer in our global map
                fn = FFI.memoize!(ptr, string(fn))

                if length(fn) > 1 && fromC
                    mod = LLVM.parent(LLVM.parent(LLVM.parent(inst)))
                    lfn = LLVM.API.LLVMGetNamedFunction(mod, fn)
                    if lfn == C_NULL
                        lfn = LLVM.API.LLVMAddFunction(mod, fn, LLVM.API.LLVMGetCalledFunctionType(inst))
                    else
                        lfn = LLVM.API.LLVMConstBitCast(lfn, LLVM.PointerType(LLVM.FunctionType(LLVM.API.LLVMGetCalledFunctionType(inst))))
                    end
                    LLVM.API.LLVMSetOperand(inst, LLVM.API.LLVMGetNumOperands(inst)-1, lfn)
                end
            end
        end
    end

    return errors
end


text:
#TODO: define GF type that knows about which dimension stores which variable
using Base.Iterators

function calc_bubble(νGrid::Vector{AbstractArray}, Gνω::SharedArray{Complex{Float64},2}, kGrid::T,
        mP::ModelParameters, sP::SimulationParameters) where T <: Union{ReducedKGrid,Nothing}
    res = SharedArray{Complex{Float64},3}(2*sP.n_iω+1, length(kGrid.kMult), 2*sP.n_iν)
    @sync @distributed for ωi in 1:2*sP.n_iω+1
        for (j,νₙ) in enumerate(νGrid[ωi])
            v1 = view(Gνω, νₙ+sP.n_iω, :)
            v2 = view(Gνω, νₙ+ωi-1, :)
            res[ωi,:,j] .= -mP.β .* conv_fft(kGrid, v1, v2)[:]
        end
    end
    return res
end

"""
Solve χ = χ₀ - 1/β² χ₀ Γ χ
    ⇔ (1 + 1/β² χ₀ Γ) χ = χ₀
    ⇔      (χ⁻¹ - χ₀⁻¹) = 1/β² Γ
    with indices: χ[ω, q] = χ₀[]
"""
function calc_χ_trilex(Γr::SharedArray{Complex{Float64},3}, bubble::SharedArray{Complex{Float64},3},
                       kGrid::T2, νGrid, sumHelper::T1, U::Float64,
                       mP::ModelParameters, sP::SimulationParameters) where {T1 <: SumHelper, T2 <: ReducedKGrid}
    χ = SharedArray{eltype(bubble), 2}((size(bubble)[1:2]...))
    γ = SharedArray{eltype(bubble), 3}((size(bubble)...))
    χ_ω = Array{Float64, 1}(undef, size(bubble,1))  # ωₙ (summed over νₙ and ωₙ)
    ωZero = sP.n_iω

    indh = ceil(Int64, size(bubble,1)/2)
    fixed_ω = typeof(sP.ωsum_type) == Tuple{Int,Int}
    ωindices = if sP.fullChi
            (1:size(bubble,1))
        elseif fixed_ω
            mid_index = Int(ceil(size(bubble,1)/2))
            default_sum_range(mid_index, sP.ωsum_type)
        else
            [(i == 0) ? indh : ((i % 2 == 0) ? indh+floor(Int64,i/2) : indh-floor(Int64,i/2)) for i in 1:size(bubble,1)]
        end

    νIndices = 1:size(bubble,3)
    lower_flag = false
    upper_cut = false
    @sync @distributed for ωi in ωindices

        Γview = view(Γr,ωi,νIndices,νIndices)
        UnitM = Matrix{eltype(Γr)}(I, length(νIndices),length(νIndices))
        for qi in 1:size(bubble, 2)
            bubble_i = view(bubble,ωi, qi, νIndices)
            bubbleD = Diagonal(bubble_i)
            χ_full = (bubbleD * Γview + UnitM)\bubbleD
            @inbounds χ[ωi, qi] = sum_freq_full(χ_full, sumHelper, mP.β)
            #TODO: absor this loop into sum_freq
            tmp = (bubble_i * (1.0 + U * χ[ωi, qi]))
            for νp in νIndices
                @inbounds γ[ωi, qi, νp] = sum_freq_full((@view χ_full[νp,:]), sumHelper, 1.0) / tmp[νp]
            end
            if sP.tc_type_f != :nothing
                extend_γ!(view(γ,ωi, qi, :), 2*π/mP.β)
            end
        end
        if (!sP.fullChi && !fixed_ω)
            @warn "Deactivating the fullChi option can lead to issues and is not recomended for this version."
            #usable = find_usable_interval(real(χ_ω), sum_type=sP.ωsum_type, reduce_range_prct=sP.usable_prct_reduction)
            #(first(usable) > ωi) && (lower_flag = true)
            #(last(usable) < ωi) && (upper_flag = true)
            #(lower_flag && upper_flag) && break
        end
    end

    if sP.ω_smoothing != :full
        for ωi in ωindices
            χ_ω[ωi] = real.(kintegrate(kGrid, χ[ωi,:])[1])
        end
    end
    if sP.ω_smoothing == :full
        for qi in 1:size(bubble, 2)
            filter_MA!(χ[1:ωZero,qi],3,χ[1:ωZero,qi])
            filter_MA!(χ[ωZero:end,qi],3,χ[ωZero:end,qi])
        end
        for ωi in ωindices
            χ_ω[ωi] = real.(kintegrate(kGrid, χ[ωi,:])[1])
        end
    elseif sP.ω_smoothing == :range
        filter_MA!(χ_ω[1:ωZero],3,χ_ω[1:ωZero])
        filter_MA!(χ_ω[ωZero:end],3,χ_ω[ωZero:end])
    end

    usable = !fixed_ω ? find_usable_interval(χ_ω, sum_type=sP.ωsum_type, reduce_range_prct=sP.usable_prct_reduction) : ωindices
    return NonLocalQuantities(χ, γ, usable, 0.0)
end

function Σ_internal!(tmp, ωindices, bubble::BubbleT, FUpDo, sumHelper::T) where T <: SumHelper
    @sync @distributed for ωi in 1:length(ωindices)
        ωₙ = ωindices[ωi]
        for qi in 1:size(bubble,2)
            for νi in 1:size(tmp,3)
                val = bubble[ωₙ,qi,:] .* FUpDo[ωₙ,νi,:]
                @inbounds tmp[ωi, qi, νi] = sum_freq_full(val, sumHelper, 1.0)
            end
        end
    end
end

function calc_Σ_ω!(Σ::SharedArray{Complex{Float64},3}, ωindices::AbstractArray{Int,1}, ωZero::Int, νZero::Int, shift::Bool, χsp, χch, γsp, γch, Gνω,
                     tmp::SharedArray{Complex{Float64},3}, U::Float64, kGrid::ReducedKGrid, sP)
    @sync @distributed for ωi in 1:length(ωindices)
        ωₙ = ωindices[ωi]
        fsp = 1.5 .* (1 .+ U*χsp[ωₙ, :])
        fch = 0.5 .* (1 .- U*χch[ωₙ, :])
        νZerop = νZero + shift*(trunc(Int64,(ωₙ - ωZero - 1)/2) + 1)
        for νi in 1:size(Σ,3)
            Kνωq = γsp[ωₙ, :, νZerop+νi] .* fsp .- γch[ωₙ, :, νZerop+νi] .* fch .- 1.5 .+ 0.5 .+ tmp[ωi,:,νZerop+νi]
            @inbounds Σ[ωi,:, νi] = conv_fft1(kGrid, Kνωq, view(Gνω,νZero + νi + ωₙ - 1,:))
        end
    end
end

function calc_Σ(Q_sp::NonLocalQuantities, Q_ch::NonLocalQuantities, bubble::BubbleT,
                Gνω::GνqT, FUpDo::SharedArray{Complex{Float64},3}, kGrid::T1,
                sumHelper_f::T2, mP::ModelParameters, sP::SimulationParameters) where {T1 <:  ReducedKGrid, T2 <: SumHelper}
    #TODO: move transform stuff to Dispersions.jl
    ωindices = (sP.dbg_full_eom_omega) ? (1:size(bubble,1)) : intersect(Q_sp.usable_ω, Q_ch.usable_ω)
    sh_b = Naive() #get_sum_helper(ωindices, sP, :b)

    tmp = SharedArray{Complex{Float64},3}(length(ωindices), size(bubble,2), size(bubble,3))
    Σ_ladder_ω = SharedArray{Complex{Float64},3}(length(ωindices), size(bubble,2), sP.n_iν-sP.shift*(trunc(Int,sP.n_iω/2) + 1))

    Σ_internal!(tmp, ωindices, bubble, FUpDo, sumHelper_f)
    calc_Σ_ω!(Σ_ladder_ω, ωindices, sP.n_iω, sP.n_iν, sP.shift, Q_sp.χ, Q_ch.χ, Q_sp.γ, Q_ch.γ,Gνω, tmp, mP.U, kGrid, sP)
    res = permutedims( mP.U .* sum_freq(Σ_ladder_ω, [1], Naive(), mP.β)[1,:,:], [2,1])
    return  res
end


text:
"""
	FixedPlane{N,T}

Fixed `N`-dimensional hyperplane.

- `p`: point on the plane
- `n`: normal vector
"""
struct FixedPlane{N,T}
    p::SVector{N,T}
    n::SVector{N,T}
end

normal(p::FixedPlane, x) = p.n

"""
"Rectangular" box in N dimensions.

`lower` and `upper` contain the lower and upper bounds in each dimension.
"""
struct RectangularBox{N,T}
    lower::SVector{N,T}
    upper::SVector{N,T}
    walls::Vector{FixedPlane{N,T}}

    function RectangularBox{N,T}(lower::SVector{N,T}, upper::SVector{N,T}) where {N,T}

        z = zero(SVector{N,T})

        walls = FixedPlane{N,T}[]

        for i in 1:N
            n = setindex(z, -1, i)  # unit vector with one non-zero component

            push!(walls, FixedPlane(lower, n))
            push!(walls, FixedPlane(upper, -n))
        end

        new{N,T}(lower, upper, walls)
    end
end

RectangularBox(lower::SVector{N,T}, upper::SVector{N,T}) where {N,T} = RectangularBox{N,T}(lower, upper)

function unit_hypercube(N, T)
	return RectangularBox(-0.5 .* ones(SVector{N,T}),
						  +0.5 .* ones(SVector{N,T})
						  )
end

max_stars_repo_name: dpsanders/HardSphereDynamics

max_forks_repo_name: dpsanders/HardSphereDynamics


text:
include("combination_models.jl")
include("turbines.jl")

abstract type AbstractWakeModel end

struct Jensen{TF} <: AbstractWakeModel
    alpha::TF
end

struct Multizone{TF} <: AbstractWakeModel
    me::TF
    ke::TF
    MU::TF
    aU::TF
    bU::TF
end

struct Gauss <: AbstractWakeModel
    version
    # for version==2014
    k_star
    # for version==2016
    turbulence_intensity
    horizontal_spread_rate
    vertical_spread_rate
    alpha_star
    beta_star
end

function wake_model(loc, deflection, model::Jensen, turbine::Turbine)
    """the original jensen wake model, from the paper: "A Note on Wind
    Generator Interaction" by N.O. Jensen (1983)"""
    # pull out the deflection distances in y (cross stream) and z (up and down)
    deflection_y = deflection[1]
    deflection_z = deflection[2]

    # find delta x, y, and z. dx is the downstream distance from the turbine to
    # the point of interest. dy and dz are the distances from the point of interest
    # and the wake center (in y and z)
    dx = loc[1]-turbine.coord.x
    dy = loc[2]-(turbine.coord.y+deflection_y)
    dz = loc[3]-(turbine.coord.z+turbine.hub_height+deflection_z)

    r0 = turbine.rotor_diameter/2.0 #turbine rotor radius
    del = sqrt(dy^2+dz^2) #distance from wake center to the point of interest
    r = model.alpha*dx + r0 #figure (1) from the paper

    if dx < 0.
        loss = 0.0
    else
        if del > r #if you're outside the wake
            loss = 0.0
        else # if you're inside the wake
            loss = 2.0*turbine.aI*(r0/(r0+model.alpha*dx))^2 #equation (2) from the paper
        end
    end
    return loss
end

function wake_model(loc, deflection, model::Multizone, turbine::Turbine)
    """The original multizone "FLORIS" wake model, from the paper:
    "Wind plant power optimization through yaw control using a parametric model
    for wake effects—a CFD simulation study" by Gebraad et al. (2014)"""

    Dt = turbine.rotor_diameter
    # extract model parameters
    ke = model.ke
    me = model.me
    MU = model.MU
    aU = model.aU
    bU = model.bU
    # pull out the deflection distances in y (cross stream) and z (up and down)
    deflection_y = deflection[1]
    deflection_z = deflection[2]

    # find delta x, y, and z. dx is the downstream distance from the turbine to
    # the point of interest. dy and dz are the distances from the point of interest
    # and the wake center (in y and z)
    dx = loc[1]-turbine.coord.x
    dy = loc[2]-(turbine.coord.y+deflection_y)
    dz = loc[3]-(turbine.coord.z+turbine.hub_height+deflection_z)

    if dx < 0.
        c = 0.0
    else
        del = sqrt(dy^2+dz^2) #distance from wake center to the point of interest

        # calculate the diameter of the wake in each of the three zones (at the specified dx)
        Dw = zeros(3)
        for i = 1:3
            Dw[i] = max(Dt+2*ke*me[i]*dx,0) # equation (13) from the paper
        end
        Rw = Dw./2 # radius of the wake zones
        # calculate the coefficient c of the appropriate wake zone:
        # equations (15, 16, and 17) from the paper
        if del < Rw[1]
            mU = MU[1]/(cos(aU+bU*turbine.gamma))
            c = (Dt/(Dt+2.0*ke*mU*dx))^2
        elseif del < Rw[2]
            mU = MU[2]/(cos(aU+bU*turbine.gamma))
            c = (Dt/(Dt+2.0*ke*mU*dx))^2
        elseif del < Rw[3]
            mU = MU[3]/(cos(aU+bU*turbine.gamma))
            c = (Dt/(Dt+2.0*ke*mU*dx))^2
        else
            c = 0.0
        end
    end
    loss = 2.0*turbine.aI*c # calculate the wake loss. Equation (14) fro mthe paper
    return loss
end

function wake_model(loc, deflection, model::Gauss, turbine::Turbine)

    deflection_y = deflection[1]
    deflection_z = deflection[2]

    dx = loc[1]-turbine.coord.x
    dy = loc[2]-(turbine.coord.y+deflection_y)
    dz = loc[3]-(turbine.coord.z+turbine.hub_height+deflection_z)

    # extract turbine properties
    Dt = turbine.rotor_diameter
    yaw = turbine.yaw
    ct = turbine.ct

    as = model.alpha_star
    bs = model.beta_star
    TI = model.turbulence_intensity

    # extract model parameters
    ks = model.k_star       # wake spread rate (k* in 2014 paper)
    TI = model.turbulence_intensity
    ky = model.horizontal_spread_rate
    kz = model.vertical_spread_rate
    as = model.alpha_star
    bs = model.beta_star

    if model.version == 2014
      """The Gaussian wake model presented by Bastankhah and Porte-Agel in
      the paper: "A new analytical model for wind-turbine wakes" (2014)"""

      # calculate beta (paper eq: 6)
      beta = 0.5*(1.0+sqrt(1.0-ct))/sqrt(1.0-ct)

      # calculate the length of the potential core (2016 paper eq: 7.3)
      x0 = ((1.0+sqrt(1.0+ct)))/(sqrt(2.0)*as*TI+bs*(1.0-sqrt(1.0-ct)))

      # calculate loss (paper eq: 23)
      enum =((dz/Dt)^2+(dy/Dt)^2)
      if dx > x0
        denom = (ks*dx/Dt+0.2*sqrt(beta))^2
      else
        denom = (ks*x0/Dt+0.2*sqrt(beta))^2
      end
      loss = (1.0 - sqrt(1.0-ct/(8.0*denom)))*exp(-enum/(2.0*denom))

    elseif model.version == 2016
      """The Gaussian wake model presented by Bastankhah and Porte-Agel in
      the paper: "Experimental and theoretical study of wind turbine wakes
      in yawed conditions" (2016)"""

      # calculate the length of the potential core (paper eq: 7.3)
      x0 = (cos(yaw)*(1+sqrt(1+ct)))/(sqrt(2)*as*TI+bs*(1-sqrt(1-ct)))

      # calculate wake spread
      if dx > x0

        # calculate horizontal wake spread (paper eq: 7.2)
        sigma_y = Dt*(ky*dx/Dt+cos(yaw)/sqrt(8))

        # calculate vertical wake spread (paper eq: 7.2)
        sigma_z = Dt*(kz*dx/Dt+1/sqrt(8))

      else # linear interpolation in the near wakes

        # calculate horizontal wake spread
        sigma_y = Dt*(ky*x0/Dt+cos(yaw)/sqrt(8))

        # calculate vertical wake spread
        sigma_z = Dt*(kz*x0/Dt+1/sqrt(8))

      end

      # calculate velocity deficit
      ey = exp(-0.5*(dy/sigma_y)^2)
      ez = exp(-0.5*(dz/sigma_z)^2)

      loss = (1-sqrt(1-ct*cos(yaw)/(8*(sigma_y*sigma_z/Dt^2))))*ey*ez

    end
end


text:
import Registrator.RegEdit: write_registry
using Pkg.TOML

using Test

@testset "RegEdit" begin

@testset "write_registry" begin
    registry = """
        name = "General"
        uuid = "23338594-aafe-5451-b93e-139f81909106"
        repo = "https://github.com/JuliaRegistries/General.git"

        description = \"\"\"
        Official general Julia package registry where people can
        register any package they want without too much debate about
        naming and without enforced standards on documentation or
        testing. We nevertheless encourage documentation, testing and
        some amount of consideration when choosing package names.
        \"\"\"

        [packages]
        00701ae9-d1dc-5365-b64a-a3a3ebf5695e = { name = "BioAlignments", path = "B/BioAlignments" }
        00718b61-6157-5045-8849-3d4c4093d022 = { name = "Convertible", path = "C/Convertible" }
        0087ddc6-3964-5e57-817f-9937aefb0357 = { name = "MathOptInterfaceMosek", path = "M/MathOptInterfaceMosek" }
        """

    registry_data = TOML.parse(registry)
    written_registry = sprint(write_registry, registry_data)
    written_registry_data = TOML.parse(written_registry)

    @test written_registry_data == registry_data
    @test written_registry == registry
end

end


text:
@using LinearAlgebraicRepresentation as LAR
@using LinearAlgebra as LA
@using ViewerGL as GL
include("../../src/Lar.jl")
LarA = Lar.Arrangement

expl = 2.0

store = [];
scaling = 1.20;
V,(VV,EV,FV,CV) = LAR.cuboid([0.25,0.25,0.25],true,[-0.25,-0.25,-0.25]);
mybox = (V,CV,FV,EV);
for k=1:20
	size = rand()*scaling
	scale = LAR.s(size,size,size)
	transl = LAR.t(rand(3)...)
	alpha = 2*pi*rand()
	rx = LAR.r(alpha,0,0); ry = LAR.r(0,alpha,0); rz = LAR.r(0,0,alpha)
	rot = rx * ry * rz
	str = LAR.Struct([ transl, scale, rot, mybox ])
	obj = LAR.struct2lar(str)
	vs = obj[1]
	diag = LinearAlgebra.norm(vs[:,8]-vs[:,1])
	if diag > 1/5
		push!(store, obj)
	end
end

str = LAR.Struct(store);
V,CV,FV,EV = LAR.struct2lar(str);

GL.VIEW([ GL.GLPol(V,CV, GL.COLORS[1]) ]);

cop_EV = LAR.coboundary_0(EV::Lar.Cells);
cop_FE = LAR.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);

model = Lar.Model(V);
Lar.addModelCells!(model, 1, cop_EV)
Lar.addModelCells!(model, 2, cop_FE)

sp_idx = LarA.spaceIndex(model, 2)

de_models = [
	LarA.face_decomposition(model, face_idx, sp_idx[face_idx])
	for face_idx = 1 : size(model, 2, 1)
]

m = deepcopy(de_models[1])
for i = 2 : length(de_models)
	Lar.uniteModels!(m, de_models[i])
end

triangulated_faces = LAR.triangulate2D(convert(Lar.Points, m.G'), m.T[1:2])
FVs = convert(Array{Lar.Cells}, triangulated_faces);
GL.VIEW( GL.GLExplode(m.G,FVs,expl,expl,expl,99,1) );

#-------------------------------------------------------------------------------
#
# V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement(
# 	 W::Lar.Points, cop_EV::Lar.ChainOp, cop_FE::Lar.ChainOp);
#
#
# triangulated_faces = Lar.triangulate2D(V, [copEV, copFE]);
# FVs = convert(Array{Lar.Cells}, triangulated_faces);
# V = convert(Lar.Points, V');
# GL.VIEW( GL.GLExplode(V,FVs,1.5,1.5,1.5,99,1) );
#
#
# EVs = Lar.FV2EVs(copEV, copFE); # polygonal face fragments
# GL.VIEW( GL.GLExplode(V,EVs,1.5,1.5,1.5,99,1) );


text:
"""
    Multi-loop iterators
"""
struct Range
    cols::Array{Int}
    rows::Array{Int}
end
const Ranges{N} = Tuple{Range,N} where N<:Int
const RangeProduct{N} = Tuple{Array{Tuple{Int, Int}}, N} where N<:Int end

Iterators.product(r::Range) = Iterators.product(r.cols,r.rows)
Base.filter(r::Ranges{N}) where N = filter()
Base.collect(r::Range) = collect(Iterators.product(r))

Base.filter(f::Function,r::Tuple{Range,N})::RangeProduct where N::Int = filter(f,Iterators.product(r...))

function rangefilter(ranges::Ranges{N2}, predicate::function) :: RangeProd{N2} where N2::Int
    collect(filter(predicate, Iterators.product(ranges...)))
end

tworangefilter(r :: Range) = s -> (s.rbegin - 1 <= 0) && (s.rend - r.rend) == 0
tworange(r :: Range, pred = _ -> true) = tworangefilter(Ranges(r), tworangefilter(Ranges(r)) ○ pred))
lowrange(r :: Range, Λ) = tworangefilter(Ranges(r), _ -> 0 <= r.rend - 1 <= Λ))
highrange(r :: Range, Λ) = tworangefilter(Ranges(r), s -> Λ + 1 <= s.rbegin - 1 <= r.rbegin - 1, r))

function cprange(r :: Range) :: RangeProd{2}
    m1lim(r1 :: Range, r2 :: Range) = r1.rend - 1 >= 1
    m2lim(r1 :: Range, r2 :: Range) = 0 <= r2.rned-1 <= min(r1.rend - 1, r.rbegin - r1.rend)
    n2lim(r1 :: Range, r2 :: Range) = 1 - r2.rbegin <= r1.rbegin - r.rend <= min(r.rend - 1, 1 - r2.r.rbegin)
    rangefilter(Ranges(r, r), m1lim ∘ m2lim ∘ n2lim)
end

function cmrange(r :: Range) :: RangeProd{2}
    m1lim(_ :: Range, r2 :: Range) = 2 <= r2.rend <= rend
    m2lim(r1 :: Range, r2:: Range) = 1 <= r1.rbegin <= r2.rend
    function n2lim(r1 :: Range, r2 :: Range)
        r2max = r1.rend == r2.rend ? r2.rend- r.rend - 1 : r.end -1
        r2.rbegin - 2r.rend + 1 <= r1.rbegin - r.rend <= min(r2max, r2.rend - 1)
    end
    rangefilter(Ranges(r, r), m1lim ∘ m2lim ∘ n2lim)
end

function iterate(rprod::RangeProd, state = eachindex(rprod.itr))
    y = iterate(state...)
    y === nothing && return y
    idx, itrs = y
    (rprod.itr[idx], (state[1], itrs))
end


text:
# Add
+(x::Posit8,y::Posit8) = ccall((:p8_add, SoftPositPath), Posit8, (Posit8,Posit8),x,y)
+(x::Posit16,y::Posit16) = ccall((:p16_add, SoftPositPath), Posit16, (Posit16,Posit16),x,y)
+(x::Posit32,y::Posit32) = ccall((:p32_add, SoftPositPath), Posit32, (Posit32,Posit32),x,y)

+(x::Posit8_1,y::Posit8_1) = ccall((:pX1_add, SoftPositPath), Posit8_1, (Posit8_1,Posit8_1,Int64),x,y,8)
+(x::Posit16_1,y::Posit16_1) = ccall((:pX1_add, SoftPositPath), Posit16_1, (Posit16_1,Posit16_1,Int64),x,y,16)
+(x::Posit24_1,y::Posit24_1) = ccall((:pX1_add, SoftPositPath), Posit24_1, (Posit24_1,Posit24_1,Int64),x,y,24)

+(x::Posit8_2,y::Posit8_2) = ccall((:pX2_add, SoftPositPath), Posit8_2, (Posit8_2,Posit8_2,Int64),x,y,8)
+(x::Posit16_2,y::Posit16_2) = ccall((:pX2_add, SoftPositPath), Posit16_2, (Posit16_2,Posit16_2,Int64),x,y,16)
+(x::Posit24_2,y::Posit24_2) = ccall((:pX2_add, SoftPositPath), Posit24_2, (Posit24_2,Posit24_2,Int64),x,y,24)

# Subtract
-(x::Posit8,y::Posit8) = ccall((:p8_sub, SoftPositPath), Posit8, (Posit8,Posit8),x,y)
-(x::Posit16,y::Posit16) = ccall((:p16_sub, SoftPositPath), Posit16, (Posit16,Posit16),x,y)
-(x::Posit32,y::Posit32) = ccall((:p32_sub, SoftPositPath), Posit32, (Posit32,Posit32),x,y)

-(x::Posit8_1,y::Posit8_1) = ccall((:pX1_sub, SoftPositPath), Posit8_1, (Posit8_1,Posit8_1,Int64),x,y,8)
-(x::Posit16_1,y::Posit16_1) = ccall((:pX1_sub, SoftPositPath), Posit16_1, (Posit16_1,Posit16_1,Int64),x,y,16)
-(x::Posit24_1,y::Posit24_1) = ccall((:pX1_sub, SoftPositPath), Posit24_1, (Posit24_1,Posit24_1,Int64),x,y,24)

-(x::Posit8_2,y::Posit8_2) = ccall((:pX2_sub, SoftPositPath), Posit8_2, (Posit8_2,Posit8_2,Int64),x,y,8)
-(x::Posit16_2,y::Posit16_2) = ccall((:pX2_sub, SoftPositPath), Posit16_2, (Posit16_2,Posit16_2,Int64),x,y,16)
-(x::Posit24_2,y::Posit24_2) = ccall((:pX2_sub, SoftPositPath), Posit24_2, (Posit24_2,Posit24_2,Int64),x,y,24)

# Multiply
*(x::Posit8,y::Posit8) = ccall((:p8_mul, SoftPositPath), Posit8, (Posit8,Posit8),x,y)
*(x::Posit16,y::Posit16) = ccall((:p16_mul, SoftPositPath), Posit16, (Posit16,Posit16),x,y)
*(x::Posit32,y::Posit32) = ccall((:p32_mul, SoftPositPath), Posit32, (Posit32,Posit32),x,y)

*(x::Posit8_1,y::Posit8_1) = ccall((:pX1_mul, SoftPositPath), Posit8_1, (Posit8_1,Posit8_1,Int64),x,y,8)
*(x::Posit16_1,y::Posit16_1) = ccall((:pX1_mul, SoftPositPath), Posit16_1, (Posit16_1,Posit16_1,Int64),x,y,16)
*(x::Posit24_1,y::Posit24_1) = ccall((:pX1_mul, SoftPositPath), Posit24_1, (Posit24_1,Posit24_1,Int64),x,y,24)

*(x::Posit8_2,y::Posit8_2) = ccall((:pX2_mul, SoftPositPath), Posit8_2, (Posit8_2,Posit8_2,Int64),x,y,8)
*(x::Posit16_2,y::Posit16_2) = ccall((:pX2_mul, SoftPositPath), Posit16_2, (Posit16_2,Posit16_2,Int64),x,y,16)
*(x::Posit24_2,y::Posit24_2) = ccall((:pX2_mul, SoftPositPath), Posit24_2, (Posit24_2,Posit24_2,Int64),x,y,24)

# Divide
/(x::Posit8,y::Posit8) = ccall((:p8_div, SoftPositPath), Posit8, (Posit8,Posit8),x,y)
/(x::Posit16,y::Posit16) = ccall((:p16_div, SoftPositPath), Posit16, (Posit16,Posit16),x,y)
/(x::Posit32,y::Posit32) = ccall((:p32_div, SoftPositPath), Posit32, (Posit32,Posit32),x,y)

/(x::Posit8_1,y::Posit8_1) = ccall((:pX1_div, SoftPositPath), Posit8_1, (Posit8_1,Posit8_1,Int64),x,y,8)
/(x::Posit16_1,y::Posit16_1) = ccall((:pX1_div, SoftPositPath), Posit16_1, (Posit16_1,Posit16_1,Int64),x,y,16)
/(x::Posit24_1,y::Posit24_1) = ccall((:pX1_div, SoftPositPath), Posit24_1, (Posit24_1,Posit24_1,Int64),x,y,24)

/(x::Posit8_2,y::Posit8_2) = ccall((:pX2_div, SoftPositPath), Posit8_2, (Posit8_2,Posit8_2,Int64),x,y,8)
/(x::Posit16_2,y::Posit16_2) = ccall((:pX2_div, SoftPositPath), Posit16_2, (Posit16_2,Posit16_2,Int64),x,y,16)
/(x::Posit24_2,y::Posit24_2) = ccall((:pX2_div, SoftPositPath), Posit24_2, (Posit24_2,Posit24_2,Int64),x,y,24)

# Square root
sqrt(x::Posit8) = ccall((:p8_sqrt, SoftPositPath), Posit8, (Posit8,),x)
sqrt(x::Posit16) = ccall((:p16_sqrt, SoftPositPath), Posit16, (Posit16,),x)
sqrt(x::Posit32) = ccall((:p32_sqrt, SoftPositPath), Posit32, (Posit32,),x)

sqrt(x::Posit8_1) = ccall((:pX1_sqrt, SoftPositPath), Posit8_1, (Posit8_1,Int64),x,8)
sqrt(x::Posit16_1) = ccall((:pX1_sqrt, SoftPositPath), Posit16_1, (Posit16_1,Int64),x,16)
sqrt(x::Posit24_1) = ccall((:pX1_sqrt, SoftPositPath), Posit24_1, (Posit24_1,Int64),x,24)

sqrt(x::Posit8_2) = ccall((:pX2_sqrt, SoftPositPath), Posit8_2, (Posit8_2,Int64),x,8)
sqrt(x::Posit16_2) = ccall((:pX2_sqrt, SoftPositPath), Posit16_2, (Posit16_2,Int64),x,16)
sqrt(x::Posit24_2) = ccall((:pX2_sqrt, SoftPositPath), Posit24_2, (Posit24_2,Int64),x,24)


text:
#!/usr/bin/env julia
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.

if isempty(ARGS)
  if Pkg.installed("CUDArt") == nothing
    exit(0)
  end
  const file,io = mktemp()
else
  const file = ARGS[1]
  const io = open(file, "w")
end

using Compat, CUDArt, JLD, SHA

include(Pkg.dir("Seep", "src", "cudahacks_common.jl"))

macro cuda(x) end
macro llvm(x) end

macro cuda_text(text)
   quote
       println(io, $text)
   end
end

macro cuda_include(filename)
   quote
       println(io, "#include \"$($filename)\"")
   end
end

function ctype(ex::Symbol)
    if ex == :Void return "void" end
    if string(ex)[1] == 'C'
        return Symbol(string(ex)[2:end])
    end
    string(ex)
end

function ctype(ex::Expr)
    if (ex.head == :&)
        return string(ctype(ex.args[1]), "* __restrict__")
    end
    error("??")
end

function cuda_kernel(name, signature, templates, text)
    if (!isempty(templates))
        template_args = templates |> x->map(x->"typename $(first(x))", x) |> x->join(x, ", ")
        println(io, "template <$template_args>")
    end

    hashname = "$(name)_$(hexstring(sha1(text)))"
    fullsig = signature |> x->map(x->"$(ctype(x[1])) $(x[2])", x) |> x->join(x, ", ")
    argnames = signature |> x->map(x->string(x[2]), x) |> x->join(x, ", ")

    println(io, "__device__ void $hashname($fullsig) {")
    println(io, text)
    println(io, "}")
    println(io)

    if (!isempty(templates))
        println(io, "extern \"C\" {")
        for type_map in TemplateIterator(templates)
            typed_sig = replace_type_parameters(signature, type_map)
            fullsig = typed_sig |> x->map(x->"$(ctype(x[1])) $(x[2])", x) |> x->join(x, ", ")

            print(io, "__global__ void $(hashname)_$(type_names(templates, type_map))($fullsig) { ")
            print(io, "$hashname($argnames);")
            println(io, " }")
        end
        println(io, "}")
        println(io)
    end
end

include(Pkg.dir("Seep", "src", "sources.jl"))

close(io)
compiler = haskey(ENV, "CXX") ? [string("-ccbin=", ENV["CXX"])] : []
run(`nvcc -ptx -O3 --use_fast_math -o $(Pkg.dir("Seep", "kernels.ptx")) $compiler -x cu -I$(Pkg.dir("Seep", "src")) $file`)

if isempty(ARGS)
    rm(file)
end


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ f4f1a9d4-fcde-11ea-24d9-efff04ac07bc
using Pkg, DrWatson

# ╔═╡ f4f1e034-fcde-11ea-08da-b7f09891f0a5
begin
	#@quickactivate "StatisticalRethinkingStan"
	using StanSample
	using StatisticalRethinking
	using GLM
end

# ╔═╡ 1bfdf3ee-fcde-11ea-0161-539e0e2b0932
md"## Clip-05-13-14s.jl"

# ╔═╡ 4fa2fe54-81dd-11eb-00c8-89b0092adf2c
begin
	df = CSV.read(sr_datadir("WaffleDivorce.csv"), DataFrame)
	df = DataFrame(
		:A => df[:, :MedianAgeMarriage],
		:M => df[:, :Marriage],
		:D => df[:, :Divorce]
	)
	scale!(df, [:M, :A, :D])
end;

# ╔═╡ 6ba505f2-81dd-11eb-311d-5f1944f30ca6
# Define the Stan language model

stan5_4_AM = "
data {
  int N;
  vector[N] A;
  vector[N] M;
}
parameters {
  real a;
  real bAM;
  real<lower=0> sigma;
}
model {
  vector[N] mu = a + bAM * M;
  a ~ normal( 0 , 0.2 );
  bAM ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  A ~ normal( mu , sigma );
}
";

# ╔═╡ 6bc53098-81dd-11eb-17d8-23c55ba9bfee
begin
	m5_4_AMs = SampleModel("m5.4.AM", stan5_4_AM)
	m5_4_data = Dict(
		"N" => size(df, 1),
		"M" => df[:, :M_s],
		"A" => df[:, :A_s]
	)
	rc5_4_AMs = stan_sample(m5_4_AMs, data=m5_4_data)
	if success(rc5_4_AMs)
		part5_4_AMs = read_samples(m5_4_AMs, :particles)
		part5_4_AMs
	end
end

# ╔═╡ a3421202-81dd-11eb-2583-d9cc66c08963
stan5_4_MA = "
data {
  int N;
  vector[N] A;
  vector[N] M;
}
parameters {
  real a;
  real bMA;
  real<lower=0> sigma;
}
model {
  vector[N] mu = a + bMA * A;
  a ~ normal( 0 , 0.2 );
  bMA ~ normal( 0 , 0.5 );
  sigma ~ exponential( 1 );
  M ~ normal( mu , sigma );
}
";

# ╔═╡ a1eb539c-81de-11eb-2f32-7dfa617502a4
begin
	m5_4_MAs = SampleModel("m5.4", stan5_4_MA);
	rc5_4_MAs = stan_sample(m5_4_MAs, data=m5_4_data);
	if success(rc5_4_MAs)

	  # Rethinking results

	  rethinking_results = "
			   mean   sd  5.5% 94.5%
		a      0.00 0.09 -0.14  0.14
		bMA   -0.69 0.10 -0.85 -0.54
		sigma  0.68 0.07  0.57  0.79
	  ";

	  part5_4_MAs = read_samples(m5_4_MAs, :particles)
	  part5_4_MAs |> display
	end
end

# ╔═╡ f501e93e-fcde-11ea-2bc3-256fb8778233
if success(rc5_4_AMs)
	begin
		post5_4_MAs_df = read_samples(m5_4_MAs, :dataframe)
		post5_4_AMs_df = read_samples(m5_4_AMs, :dataframe)

		pMA = plotbounds(df, :M, :A, post5_4_MAs_df,
			[:a, :bMA, :sigma]; ylims=(10, 30))
		pAM = plotbounds(df, :A, :M, post5_4_AMs_df,
			[:a, :bAM, :sigma]; ylims=(10, 30))
		plot(pAM, pMA, layout=(1, 2))
	end
end

# ╔═╡ f502a090-fcde-11ea-37d0-233bf56ce068
md"##### Compute standardized residuals."

# ╔═╡ f50ffec0-fcde-11ea-2670-534a1a8a9725
if success(rc5_4_MAs)
	begin
		figs = Vector{Plots.Plot{Plots.GRBackend}}(undef, 4)
		a = -2.5:0.1:3.0
		mu_MA = mean(part5_4_MAs.a) .+ mean(part5_4_MAs.bMA)*a
		figs[1] = plot(xlab="Age at marriage (std)", ylab="Marriage rate (std)", leg=false)
		plot!(a, mu_MA)
		scatter!(df[:, :A_s], df[:, :M_s])
		annotate!([(df[9, :A_s]-0.1, df[9, :M_s], Plots.text("DC", 6, :red, :right))])
	end
end

# ╔═╡ f510ac4e-fcde-11ea-05a8-fdda93bbfc2a
if success(rc5_4_AMs)
	begin
		m = -2.0:0.1:3.0
		mu_AM = mean(part5_4_AMs.a) .+ mean(part5_4_AMs.bAM)*m
		figs[2] = plot(ylab="Age at marriage (std)", xlab="Marriage rate (std)", leg=false)
		plot!(m, mu_AM)
		scatter!(df[:, :M_s], df[:, :A_s])
		annotate!([(df[9, :M_s]+0.2, df[9, :A_s], Plots.text("DC", 6, :red, :left))])
	end
end

# ╔═╡ f51db3b0-fcde-11ea-1be0-eb23b24429e1
if success(rc5_4_MAs)
	begin
		mu_MA_obs = mean(part5_4_MAs.a) .+ mean(part5_4_MAs.bMA)*df[:, :A_s]
		res_MA = df[:, :M_s] - mu_MA_obs

		df2 = DataFrame(
			:d => df[:, :D_s],
			:r => res_MA
		)

		m1 = lm(@formula(d ~ r), df2)
		#coef(m1) |> display

		figs[3] = plot(xlab="Marriage rate residuals", ylab="Divorce rate (std)", leg=false)
		plot!(m, coef(m1)[1] .+ coef(m1)[2]*m)
		scatter!(res_MA, df[:, :D_s])
		vline!([0.0], line=:dash, color=:black)
		annotate!([(res_MA[9], df[9, :D_s]+0.1, Plots.text("DC", 6, :red, :bottom))])
	end
end

# ╔═╡ f52450f0-fcde-11ea-066b-153ac07ad60d
if success(rc5_4_AMs)
	begin
		mu_AM_obs = mean(part5_4_AMs.a) .+ mean(part5_4_AMs.bAM)*df[:, :M_s]
		res_AM = df[:, :A_s] - mu_AM_obs
		df3 = DataFrame(
			:d => df[:, :D_s],
			:r => res_AM
		)

		m2 = lm(@formula(d ~ r), df3)
		#coef(m2) |> display

		figs[4] = plot(xlab="Age at marriage residuals", ylab="Divorce rate (std)", leg=false)
		plot!(a, coef(m2)[1] .+ coef(m2)[2]*a)
		scatter!(res_AM, df[:, :D_s])
		vline!([0.0], line=:dash, color=:black)
		annotate!([(res_AM[9]-0.1, df[9, :D_s], Plots.text("DC", 6, :red, :right))])
	end
end

# ╔═╡ f533f69a-fcde-11ea-0181-e5e75979c929
plot(figs..., layout=(2,2))

# ╔═╡ f53b5e30-fcde-11ea-05c9-41a05f46fa1d
md"## End of clip-05-13-14s.jl"

# ╔═╡ Cell order:
# ╟─1bfdf3ee-fcde-11ea-0161-539e0e2b0932
# ╠═f4f1a9d4-fcde-11ea-24d9-efff04ac07bc
# ╠═f4f1e034-fcde-11ea-08da-b7f09891f0a5
# ╠═4fa2fe54-81dd-11eb-00c8-89b0092adf2c
# ╠═6ba505f2-81dd-11eb-311d-5f1944f30ca6
# ╠═6bc53098-81dd-11eb-17d8-23c55ba9bfee
# ╠═a3421202-81dd-11eb-2583-d9cc66c08963
# ╠═a1eb539c-81de-11eb-2f32-7dfa617502a4
# ╠═f501e93e-fcde-11ea-2bc3-256fb8778233
# ╟─f502a090-fcde-11ea-37d0-233bf56ce068
# ╠═f50ffec0-fcde-11ea-2670-534a1a8a9725
# ╠═f510ac4e-fcde-11ea-05a8-fdda93bbfc2a
# ╠═f51db3b0-fcde-11ea-1be0-eb23b24429e1
# ╠═f52450f0-fcde-11ea-066b-153ac07ad60d
# ╠═f533f69a-fcde-11ea-0181-e5e75979c929
# ╟─f53b5e30-fcde-11ea-05c9-41a05f46fa1d