text:
using ImageFiltering, ImageCore, Colors, FixedPointNumbers, OffsetArrays, Base.Test

isapprox_const(A::AbstractArray, n::Number) = isapprox(A, fill(n, size(A)))

@testset "Deprecated" begin
    @testset "Array padding" begin
        A = [1 2; 3 4]
        @test padarray(A, (0,0), (0,0), "replicate") == A
        @test padarray(A, (1,2), (2,0), "replicate") == OffsetArray([1 1 1 2; 1 1 1 2; 3 3 3 4; 3 3 3 4; 3 3 3 4], (-1,-2))
        @test padarray(A, [2,1], [0,2], "circular") == OffsetArray([2 1 2 1 2; 4 3 4 3 4; 2 1 2 1 2; 4 3 4 3 4], (-2,-1))
        @test padarray(A, (1,2), (2,0), "symmetric") == OffsetArray([2 1 1 2; 2 1 1 2; 4 3 3 4; 4 3 3 4; 2 1 1 2], (-1,-2))
        @test padarray(A, (1,2), (2,0), "value", -1) == OffsetArray([-1 -1 -1 -1; -1 -1 1 2; -1 -1 3 4; -1 -1 -1 -1; -1 -1 -1 -1], (-1,-2))
        A = [1 2 3; 4 5 6]
        @test padarray(A, (1,2), (2,0), "reflect") == OffsetArray([6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6; 3 2 1 2 3; 6 5 4 5 6], (-1,-2))
        A = [1 2; 3 4]
        @test padarray(A, (1,1)) == OffsetArray([1 1 2 2; 1 1 2 2; 3 3 4 4; 3 3 4 4], (-1,-1))
        @test padarray(A, (1,1), "replicate", "both") == OffsetArray([1 1 2 2; 1 1 2 2; 3 3 4 4; 3 3 4 4], (-1,-1))
        @test padarray(A, (1,1), "circular", "pre") == OffsetArray([4 3 4; 2 1 2; 4 3 4], (-1,-1))
        @test padarray(A, (1,1), "symmetric", "post") == [1 2 2; 3 4 4; 3 4 4]
        A = ["a" "b"; "c" "d"]
        @test padarray(A, (1,1)) == OffsetArray(["a" "a" "b" "b"; "a" "a" "b" "b"; "c" "c" "d" "d"; "c" "c" "d" "d"], (-1,-1))
        @test_throws ArgumentError padarray(A, (1,1), "unknown")
        # issue #292
        A = trues(3,3)
        @test typeof(parent(padarray(A, (1,2), (2,1), "replicate"))) == BitArray{2}
        @test eltype(padarray(grayim(A), (1,2), (2,1), "replicate")) == Gray{Bool}
        # issue #525
        A = falses(10,10,10)
        B = view(A,1:8,1:8,1:8)
        @test isa(parent(padarray(A, ones(Int,3), ones(Int,3), "replicate")), BitArray{3})
        @test isa(parent(padarray(B, ones(Int,3), ones(Int,3), "replicate")), BitArray{3})
    end

    @testset "Filtering" begin
        EPS = 1e-14
        imgcol = colorim(rand(3,5,6))
        imgcolf = convert(Array{RGB{UFixed8}}, imgcol)
        for T in (Float64, Int)
            A = zeros(T,3,3); A[2,2] = 1
            kern = rand(3,3)
            @test maximum(abs(imfilter(A, kern) - rot180(kern))) < EPS
            kern = rand(2,3)
            @test maximum(abs(imfilter(A, kern)[1:2,:] - rot180(kern))) < EPS
            kern = rand(3,2)
            @test maximum(abs(imfilter(A, kern)[:,1:2] - rot180(kern))) < EPS
        end
        kern = zeros(3,3); kern[2,2] = 1
        @test maximum(map(abs, imgcol - imfilter(imgcol, kern))) < EPS
        @test maximum(map(abs, imgcolf - imfilter(imgcolf, kern))) < EPS
        for T in (Float64, Int)
            # Separable kernels
            A = zeros(T,3,3); A[2,2] = 1
            kern = rand(3).*rand(3)'
            @test maximum(abs(imfilter(A, kern) - rot180(kern))) < EPS
            kern = rand(2).*rand(3)'
            @test maximum(abs(imfilter(A, kern)[1:2,:] - rot180(kern))) < EPS
            kern = rand(3).*rand(2)'
            @test maximum(abs(imfilter(A, kern)[:,1:2] - rot180(kern))) < EPS
        end
        A = zeros(3,3); A[2,2] = 1
        kern = rand(3,3)
        @test maximum(abs(imfilter_fft(A, kern) - rot180(kern))) < EPS
        kern = rand(2,3)
        @test maximum(abs(imfilter_fft(A, kern)[1:2,:] - rot180(kern))) < EPS
        kern = rand(3,2)
        @test maximum(abs(imfilter_fft(A, kern)[:,1:2] - rot180(kern))) < EPS
        kern = zeros(3,3); kern[2,2] = 1
        @test maximum(map(abs, imgcol - imfilter_fft(imgcol, kern))) < EPS
        @test maximum(map(abs, imgcolf - imfilter_fft(imgcolf, kern))) < EPS

        @test isapprox_const(imfilter(ones(4,4), ones(3,3)), 9.0)
        @test isapprox_const(imfilter(ones(3,3), ones(3,3)), 9.0)
        @test isapprox_const(imfilter(ones(3,3), [1 1 1;1 0.0 1;1 1 1]), 8.0)
        img = ones(4,4)
        @test isapprox_const(imfilter(img, ones(3,3)), 9.0)
        A0 = zeros(5,5,3); A0[3,3,[1,3]] = 1
        A = colorview(RGB, permuteddimsview(A0, (3,1,2)))
        kern = rand(3,3)
        kernpad = zeros(5,5); kernpad[2:4,2:4] = kern
        Af = permuteddimsview(channelview(imfilter(A, kern)), (2,3,1))

        @test_approx_eq cat(3, rot180(kernpad), zeros(5,5), rot180(kernpad)) Af
        @test isapprox_const(imfilter(ones(4,4),ones(1,3),"replicate"), 3.0)

        A = zeros(5,5); A[3,3] = 1
        kern = rand(3,3)
        Af = imfilter(A, kern, Inner())
        @test Af == OffsetArray(rot180(kern), (1,1))
        Afft = imfilter_fft(A, kern, "inner")
        @test_approx_eq Af Afft
        h = [0.24,0.87]
        hfft = imfilter_fft(eye(3), h, "inner")
        hfft[abs(hfft) .< 3eps()] = 0
        @test_approx_eq imfilter(eye(3), h, Inner()) hfft  # issue #204

        # circular
        A = zeros(3, 3)
        A[3,2] = 1
        kern = rand(3,3)
        @test_approx_eq imfilter_fft(A, kern, "circular") kern[[1,3,2],[3,2,1]]

        A = zeros(5, 5)
        A[5,3] = 1
        kern = rand(3,3)
        @test_approx_eq imfilter_fft(A, kern, "circular")[[1,4,5],2:4] kern[[1,3,2],[3,2,1]]

        A = zeros(5, 5)
        A[5,3] = 1
        kern = rand(3,3)
        @test_approx_eq imfilter(A, kern, "circular")[[1,4,5],2:4] kern[[1,3,2],[3,2,1]]

        @test isapprox_const(imfilter_gaussian(ones(4,4), [5,5]), 1.0)
        A = fill(convert(Float32, NaN), 4, 4)
        A[1:4,1] = 1:4
        @test isequal(imfilter_gaussian(A, [0,0]), A)
        @test_approx_eq imfilter_gaussian(A, [0,3]) A
        B = copy(A)
        B[isfinite(B)] = 2.5
        @test_approx_eq imfilter_gaussian(A, [10^3,0]) B
        @test maximum(map(abs, imfilter_gaussian(imgcol, [10^3,10^3]) - mean(imgcol))) < 1e-4
        @test maximum(map(abs, imfilter_gaussian(imgcolf, [10^3,10^3]) - mean(imgcolf))) < 1e-4
        A = rand(4,5)
        img = reinterpret(Gray{Float64}, A)
        imgf = imfilter_gaussian(img, [2,2])
        @test_approx_eq reinterpret(Float64, data(imgf)) imfilter_gaussian(A, [2,2])
        A = rand(3,4,5)
        img = colorim(A)
        imgf = imfilter_gaussian(img, [2,2])
        @test_approx_eq channelview(data(imgf)) imfilter_gaussian(A, [0,2,2])

        # A = zeros(Int, 9, 9); A[5, 5] = 1
        # @test maximum(abs(imfilter_LoG(A, [1,1]) - imlog(1.0))) < EPS
        # @test maximum(imfilter_LoG([0 0 0 0 1 0 0 0 0], [1,1]) - sum(imlog(1.0),1)) < EPS
        # @test maximum(imfilter_LoG([0 0 0 0 1 0 0 0 0]', [1,1]) - sum(imlog(1.0),2)) < EPS

        # @test imaverage() == fill(1/9, 3, 3)
        # @test imaverage([3,3]) == fill(1/9, 3, 3)
        # @test_throws ErrorException imaverage([5])
    end
end

nothing


text:
using Random

function is_sorted(a::Vector{Float64})
    for i = 1:length(a)-1
        if (a[i] > a[i + 1])
            return false
        end
    end
    return true
end

function bogo_sort!(a::Vector{Float64})
    while(!is_sorted(a))
        shuffle!(a)
    end
end

function main()
    a = [1.0, 3.0, 2.0, 4.0]
    bogo_sort!(a)
    println(a)
end

main()

ext: jl


text:
module escape_dimerB
export escape_exit_num_dimer
using DifferentialEquations
using PolynomialRoots
using Printf
# using MATLAB
# using ProgressMeter

include("Dimer_Eq_of_M.jl")
include("leap_frog_definitions.jl")
H=.25
max_hit_q1 = 11
max_hit_p2 = 1000
barrier = 5
condition_max_hits_back_forward(u, t, integrator) = u[5] > max_hit_q1 || u[6] > max_hit_p2 || maximum([abs(u[1]),abs(u[2]),abs(u[3]),abs(u[4])]) > barrier
affect_stop!(integrator) = terminate!(integrator)

function condition_hits_PSS_q1(u, t, integrator) # Event when event_f(u,t) == 0
    u[1]
end
function condition_hits_PSS_p2(u, t, integrator) # Event when event_f(u,t) == 0
    u[4]
end

# function affect_update_iterator_q1!(integrator)
#     integrator.u[5] = integrator.u[5] + 1
# end

function affect_update_iterator_q1!(integrator)
    # integrator.u[5] = integrator.u[5] + 1
    q2 = integrator.u[2]
    p1 = integrator.u[3]
    p2 = integrator.u[4]
    P1 = P1_poly(q2, p2, H)
    if !isempty(P1) &&  abs(p1 - P1) < .5
        integrator.u[5] = integrator.u[5] + 1
    end
end

function affect_update_iterator_p2!(integrator)
    integrator.u[6] = integrator.u[6] + 1
end
callback_max_hits_bf = DiscreteCallback(condition_max_hits_back_forward, affect_stop!)
callback_hits_PSS_q1_forward = ContinuousCallback(condition_hits_PSS_q1, nothing,  affect_update_iterator_q1!, rootfind = true)
callback_hits_PSS_q1_backward = ContinuousCallback(condition_hits_PSS_q1,  affect_update_iterator_q1!, nothing, rootfind = true)

callback_hits_PSS_p2_forward = ContinuousCallback(condition_hits_PSS_p2, affect_update_iterator_p2!, rootfind = false)
callback_hits_PSS_p2_backward = ContinuousCallback(condition_hits_PSS_p2, affect_update_iterator_p2!, rootfind = false)

cb_forward = CallbackSet(callback_hits_PSS_q1_forward, callback_hits_PSS_p2_forward, callback_max_hits_bf)
cb_backward = CallbackSet(callback_hits_PSS_q1_backward, callback_hits_PSS_p2_backward, callback_max_hits_bf)

function escape_exit_num_dimer(mesh_list, t_end, H)
    Q2 = mesh_list[1]; P2 = mesh_list[2];
    P1 = P1_poly(Q2, P2, H)
    if ~isempty(P1)
        u0=zeros(6)
        u0[1]=0
        u0[2]=Q2
        u0[3]=P1
        u0[4]=P2
        u0[5]=0
        u0[6]=0

       #constructor for ODE
        # prob = HamiltonianProblem{true}(Hamiltonian_Dimer, q0, p0, (0., t_end));
        prob=ODEProblem(Dimer_Eq_of_M_hits,u0,(0., t_end) )
        sol_f = solve(prob, RK4(), maxiters = 1e20, reltol = 1e-8, abstol = 1e-10, callback = cb_forward, save_start = false, save_end = true, save_everystep = true)

        q1 = sol_f[1,end]
        q2 = sol_f[2,end]
        p1 = sol_f[3,end]
        p2 = sol_f[4,end]

        q = [q1, q2]
        p = [p1, p2]

        forward_hits = sol_f[5,end]

        # println(sol_f.t)

        diff = abs(Hamiltonian_Dimer(q, p, 0) - H)
        if length(sol_f[1,:]) > 1
            sign_sol = sign(sol_f[1,2])
        else
            sign_sol = 0
        end
        if maximum([abs(q1),abs(q2),abs(p1),abs(p2)]) > barrier
            return forward_hits
        else
            return -1
        end
            # # return forward_hits
        #     q0, p0 = [zeros(3) for i in 1:2]
        #     q0[1] = 0
        #     q0[2] = Q2
        #     q0[3] = 0
        #     p0[1] = P1
        #     p0[2] = P2
        #     p0[3] = 0
        #     prob_b = HamiltonianProblem{true}(Hamiltonian_Dimer_Backwards, q0, p0, (0., t_end));
        #     #solve ode , save_everystep=false is important to prevent sol to include all points, not just event points
        #     # sol_b=solve(prob_b,  RK4(), maxiters = 1e20,callback=cb_backward,save_start=false,save_end=true,save_everystep=false)
        #     # sol_b=solve(prob_b,  RK4(), maxiters = 1e20,reltol=1e-8,abstol=1e-10, callback=cb_backward,save_start=false,save_end=true,save_everystep=false)
        #     sol_b=solve(prob_b, Vern9(),maxiters=1e20, reltol=1e-13,abstol=1e-15,callback=cb_backward,save_start=false,save_end=true,save_everystep=false)

        #     uf = zeros(6)
        #     uf[1] = sol_b[1,end]
        #     uf[2] = sol_b[2,end]
        #     uf[3] = sol_b[3,end]
        #     uf[4] = sol_b[4,end]
        #     uf[5] = sol_b[5,end]
        #     uf[6] = sol_b[6,end]
        #     q1 = uf[1]
        #     q2 = uf[2]
        #     p1 = uf[4]
        #     p2 = uf[5]
        #     q = [q1, q2]
        #     p = [p1, p2]

        #     backward_hits = uf[3]
        #     q = [q1, q2]
        #     p = [p1, p2]
        #     diff = abs(Hamiltonian_Dimer(q, p, 0) - H)
        #     if diff < 1e-10 && maximum([abs(q1),abs(q2),abs(p1),abs(p2)]) > barrier
        #         return forward_hits+ backward_hits# sol_f.t, sol_b.t
        #     else
        #         return -1
        #     end
        # else
        #     return -1
        # end
        # if dH<1e-5
        #     if is_it(uf,H,tol_dist)
        #         return sol.u[end][5]
        #     else
        #         return -1
        #     end
        # end
    else
        return -1
    end

end
end


text:
using RobotDynamics: orientation

# Set up model
model = Body{UnitQuaternion{Float64}}()
@test size(model) == (13, 6, 13)
x, u = rand(model)
t, dt = 0, 0.1
z = RD.KnotPoint(x, u, t, dt)
q = orientation(model, x)
ir, iq, iv, iω, iu = RobotDynamics.gen_inds(model)

# Test analytical continuous-time Jacobian
RobotDynamics.velocity_frame(::Body) = :world
function AD_jacobian!(
    ∇f::Matrix,
    model::RD.AbstractModel,
    z::RD.AbstractKnotPoint{n,m},
) where {n,m}
    ix, iu = SVector{n}(1:n), SVector{m}(n+1:n+m)
    t = z.t
    f_aug(z) = RD.dynamics(model, z[ix], z[iu], t)
    s = z.z
    ForwardDiff.jacobian!(∇f, f_aug, s)
end
F = zeros(13, 19)
AD_jacobian!(F, model, z)

F2 = zero(F)
y2 = zeros(13)
RD.jacobian!(RD.StaticReturn(), RD.ForwardAD(), model, F2, y2, z)
@test F2 ≈ F

ω = RD.angular_velocity(model, x)
ForwardDiff.jacobian(
    q -> Rotations.kinematics(UnitQuaternion(q, false), ω),
    Rotations.params(q),
)
Rotations.kinematics(q, ω)

D = RobotDynamics.DynamicsJacobian(13, 6)
RD.jacobian!(RD.StaticReturn(), RD.ForwardAD(), model, D, y2, z)
@test D ≈ F

# @btime AD_jacobian!($F, $model, $z)
# @btime jacobian!($F, $model, $z)

# Test in body frame
RobotDynamics.velocity_frame(::Body) = :body
F3 = zeros(13, 19)
AD_jacobian!(F3, model, z)
@test !(F3 ≈ F)
RD.jacobian!(RD.StaticReturn(), RD.ForwardAD(), model, D, y2, z)
@test D ≈ F3

dmodel = RD.DiscretizedDynamics{RD.RK4}(model)
RD.jacobian!(RD.StaticReturn(), RD.ForwardAD(), dmodel, F, y2, z)

RobotDynamics.velocity_frame(::Body) = :world
RD.jacobian!(RD.StaticReturn(), RD.ForwardAD(), dmodel, F2, y2, z)
@test !(F ≈ F2)
@test !(F ≈ D)
@test !(F2 ≈ D)

# Test state error functions
x0 = rand(model)[1]
dx = RobotDynamics.state_diff(model, x, x0)
@test dx[ir] ≈ (x-x0)[ir]
@test dx[iv.-1] ≈ (x-x0)[iv]
@test dx[iω.-1] ≈ (x-x0)[iω]
@test dx[4:6] ≈ orientation(model, x) ⊖ orientation(model, x0)

G0 = Rotations.∇differential(q)
G = zeros(RD.state_dim(model), RobotDynamics.errstate_dim(model))
@test size(G) == (13, 12)
RD.state_diff_jacobian!(model, G, z)
@test G ≈ cat(I(3), G0, I(6), dims = (1, 2))

b = @SVector rand(13)
∇G0 = Rotations.∇²differential(q, b[SA[4, 5, 6, 7]])
∇G = zeros(RD.errstate_dim(model), RD.errstate_dim(model))
@test size(∇G) == (12, 12)
RobotDynamics.∇²differential!(model, ∇G, x, b)
@test ∇G ≈ cat(zeros(3, 3), ∇G0, zeros(6, 6), dims = (1, 2))

# @btime RobotDynamics.state_diff_jacobian!($G, $model, $z)
# @btime RobotDynamics.∇²differential!($∇G, $model, $x, $b)


text:
using Hotwire
using Polynomials
using CurveFit
using Test

@testset "Hotwire.jl" begin

    # Resistance
    r = Resistor(1e3, 0.0, 298.15)

    @test resistance(r) == 1e3
    @test resistance(r, 298.15) == 1e3

    r = Resistor(1e3, 0.01, 298.15)
    @test resistance(r, 298.15) == 1e3
    @test resistance(r, 299.15) == 1.01e3

    @test temperature(r, 1010.0) == 299.15

    # Thermistor
    r = Thermistor(1e3, 0.0, 298.15)
    @test resistance(r) == 1e3
    @test resistance(r, 303.15) == 1e3

    r = Thermistor(1e3, 3300.0, 298.15)
    @test temperature(r, resistance(r, 50.0+273.15)) ≈ 50.0+273.15
    @test temperature(r, resistance(r, 100.0+273.15)) ≈ 100.0+273.15

    @test 100.0 ≈ resistance(r, temperature(r, 100.0))

    # CTASensor

    r = Resistor(1.0, 0.01, 300.0)
    w = CTASensor(r, 2.0)
    @test r == resistor(w)
    @test optemperature(w) ≈ 400.0
    @test overtemp(w) ≈ 100.0
    @test overheat_ratio(w) ≈ 1.0

    w1 = Wire(1.0, 2.0, 300.0, 0.01)
    @test w == w1

    r = Thermistor(1e3, 3300, 25.0+273.15)
    Rw = 100.0
    w = CTASensor(r, Rw)
    Tw = temperature(r, Rw)
    a = Rw/1e3 - 1.0
    @test r == resistor(w)
    @test optemperature(w) ≈ Tw
    @test overtemp(w) ≈ Tw - 298.15
    @test overheat_ratio(w) ≈ a

    # CCASensor - not much to do
    r = Resistor(1.0, 0.01, 300.0)
    w = CCASensor(r, 1.0)
    @test resistor(w) == r
    @test current(w) == 1.0

    # tempcorr.jl
    # tempcorr(CTASensor)
    Too = 20.0+273.15
    r = Resistor(3.5, 0.0036, Too)
    w = CTASensor(r, 5.5)
    @test 1.0 ≈ tempcorr(w, 1.0, Too) # No temperature correction
    @test tempcorr(w, 1.0, 25.0+273.15) > 1.0 # Increased temperature -> lower anemometer output
    @test tempcorr(w, 1.0, 15.0+273.15) < 1.0 # Decreased temperature -> higher anemometer output
    E1 = 1.0
    Too1 = 25.0+273.15
    E = tempcorr(w, E1, Too1)
    Tw = optemperature(w)
    @test E/E1 ≈ sqrt( (Tw - Too) / (Tw - Too1) )

    r = Thermistor(1e3, 3300, Too)
    w = CTASensor(r, 10.0)
    @test 1.0 ≈ tempcorr(w, 1.0, 20.0+273.15) # No temperature correction
    @test tempcorr(w, 1.0, 25.0+273.15) > 1.0 # Increased temperature -> lower anemometer output
    @test tempcorr(w, 1.0, 15.0+273.15) < 1.0 # Decreased temperature -> higher anemometer output
    E1 = 1.0
    Too1 = 25.0+273.15
    E = tempcorr(w, E1, Too1)
    Tw = optemperature(w)
    @test E/E1 ≈ sqrt( (Tw - Too) / (Tw - Too1) )

    @test 1.0 ≈ tempcorr(w, 1.0, Too, Too)
    @test tempcorr(w, 1.0, 25.0+273.15, Too) > 1.0
    @test tempcorr(w, 1.0, 15.0+273.15, Too) < 1.0
    E = tempcorr(w, E1, Too1, Too)
    Tw = optemperature(w)
    @test E/E1 ≈ sqrt( (Tw - Too) / (Tw - Too1) )

    # tempcorr(CCASensor)
    To = 20.0+273.15
    r = Resistor(3.5, 0.0036, To)
    Two = 100.0+273.15

    Rwo = resistance(r, Two)
    I = 0.2
    Eo = Rwo * I
    hA = Eo*I / (Two - To)
    g = 1.0

    w = CCASensor(r, I, g)

    @test tempcorr(w, Eo, To) ≈ Eo
    @test tempcorr(w, Eo, To, To) ≈ Eo

    Te = 30.0+273.15
    E = Eo
    Eo = tempcorr(w, E, Te)
    Two = temperature(r, Eo/(g*I))
    Tw = temperature(r, E/(g*I))
    Qo = Eo/g * I
    Q = E/g * I
    hAo = Qo / (Two - To)
    hA = Q / (Tw - Te)
    @test Eo / E ≈ (Two - To) / (Tw - Te)
    @test hA ≈ hAo

    Te = 30.0+273.15
    E = Eo
    Eo = tempcorr(w, E, Te, temperature(w.R))
    Two = temperature(r, Eo/(g*I))
    Tw = temperature(r, E/(g*I))
    Qo = Eo/g * I
    Q = E/g * I
    hAo = Qo / (Two - To)
    hA = Q / (Tw - Te)
    @test Eo / E ≈ (Two - To) / (Tw - Te)
    @test hA ≈ hAo

    # Testing root finding equation
    @test Hotwire.funroot(x->(x^2-1.0), 0.5, 0.6) ≈ 1.0 atol=1e-4
    @test Hotwire.funroot(x->sin(x), π*1.1, π*1.02) ≈ π atol=1e-4

    # Testing tempcorr for CCASensor{Thermistor}
    To = 20.0+273.15
    r = Thermistor(1e3, 3300, To)
    g = 1.0
    I = 0.2
    w = CCASensor(r, I, g)
    @test tempcorr(w, 1.0, To) ≈ 1.0

    for Te in [15.0, 30.0] .+ 273.15
        for E in [0.5, 1.0, 2.0]
            Eo = tempcorr(w, E, Te)
            Two = temperature(r, Eo/(g*I))
            Tw = temperature(r, E/(g*I))
            @test Eo/E ≈ (Two-To) / (Tw - Te) atol=1e-6
            Qo = Eo/g * I
            Q = E/g * I
            hAo = Qo / (Two - To)
            hA = Q / (Tw - Te)
            @test hA ≈ hAo atol=1e-5
        end
    end

    for Te in [15.0, 30.0] .+ 273.15
        for E in [0.5, 1.0, 2.0]
            Eo = tempcorr(w, E, Te, To)
            Two = temperature(r, Eo/(g*I))
            Tw = temperature(r, E/(g*I))
            @test Eo/E ≈ (Two-To) / (Tw - Te) atol=1e-6
            Qo = Eo/g * I
            Q = E/g * I
            hAo = Qo / (Two - To)
            hA = Q / (Tw - Te)
            @test hA ≈ hAo atol=1e-5
        end
    end

    # Testing calibration curve
    Rw = Thermistor(5e3, 3500, 25.0+273.15);
    cta = CTASensor(Rw, 100.0, 20/120);

    E = [1.333, 1.385, 1.423, 1.48, 1.521, 1.551, 1.575, 1.624]
    T0 = 20.0+273.15
    coefs = [2527.0, -3916.0, -1548.0, 6110.0, -3945.0, 815.0]
    fit = Polynomial(coefs)
    U = fit.(E)

    cal = calibr_curve(cta, U, E, T0, 5, extrapolate=false)

    @test cal.T0 == T0
    @test cal.fit.coeffs ≈ coefs rtol=1e-4
    @test all(cal.(cta, E) .≈ U)

    U = 1.0:0.5:10.0
    A = 2.0
    B = 0.7
    n = 0.42
    E = sqrt.(A .+ B .* U .^ n)
    fitk = KingFit(E, U)
    efitk = ExtrapolateFit(fitk, sqrt(A), E[1])

    cal = CalibrCurve(E, collect(U), efitk, T0)

    @test efitk.A ≈ A
    @test efitk.B ≈ B
    @test efitk.n ≈ n

    @test all(cal.(cta, E) .≈ U)

    @test cal(cta, sqrt(2.0)) ≈ 0.0

    # Test Probe3d stuff
    k² = (0.0, 0.0, 0.0)
    h² = (1.0, 1.0, 1.0)

    U = collect(1.0:0.5:10.0)
    T0 = 20.0+273.15
    E1 = sqrt.(2.0 .+ 0.7 .* U .^ 0.42)
    E2 = sqrt.(2.2 .+ 0.65 .* U .^ 0.42)
    E3 = sqrt.(1.95 .+ 0.75 .* U .^ 0.41)
    cal = (CalibrCurve(E1, U, KingFit(E1, U), T0),
           CalibrCurve(E2, U, KingFit(E2, U), T0),
           CalibrCurve(E3, U, KingFit(E3, U), T0))

    wires = (Wire(5.0, 5.0*1.8, T0),
             Wire(5.2, 5.2*1.8, T0),
             Wire(5.3, 5.3*1.8, T0))

    probe = Probe3d(wires, cal, k², h²)
    for i in 1:length(U)
        Uⁱ = probe(E1[i], E2[i], E3[i], T0)
        @test Uⁱ[1] ≈ U[i] rtol=1e-4
    end

    k² = (0.03, 0.04, 0.05)
    h² = (1.05, 1.08, 1.02)

    probe = Probe3d(wires, cal, k², h²)
    for i in 1:length(U)
        Uⁱ = probe(E1[i], E2[i], E3[i], T0)
        @test Uⁱ[1] ≈ U[i] rtol=1e-4
    end

end


text:
@doc raw"""
Axisymmetric tokamak equilibrium in (R,Z,ϕ) coordinates with covariant
components of the vector potential given by
```math
A (R, Z, \phi) = \frac{B_0}{2} \, \bigg( R_0 \, \frac{Z}{R} , \, - R_0 \, \ln \bigg( \frac{R}{R_0} \bigg) , \, - \frac{r^2}{q_0} \bigg)^T ,
```
resulting in the magnetic field with covariant components
```math
B (R, Z, \phi) = \frac{B_0}{q_0} \, \bigg( - \frac{Z}{R} , \, \frac{R - R_0}{R} , \, - q_0 R_0 \bigg)^T ,
```
where $r = \sqrt{ (R - R_0)^2 + Z^2 }$.

Parameters:
 * `R₀`: position of magnetic axis
 * `B₀`: B-field at magnetic axis
 * `q₀`: safety factor at magnetic axis
"""
struct AxisymmetricTokamakCylindrical{T <: Number} <: AnalyticEquilibrium
    name::String
    R₀::T
    B₀::T
    q₀::T

    function AxisymmetricTokamakCylindrical{T}(R₀::T, B₀::T, q₀::T) where T <: Number
        new("AxisymmetricTokamakCylindricalEquilibrium", R₀, B₀, q₀)
    end
end

AxisymmetricTokamakCylindrical(R₀::T=1.0, B₀::T=1.0, q₀::T=2.0) where T <: Number = AxisymmetricTokamakCylindrical{T}(R₀, B₀, q₀)

function Base.show(io::IO, equ::AxisymmetricTokamakCylindrical)
    print(io, "Axisymmetric Tokamak Equilibrium in (R,Z,ϕ) Coordinates with\n")
    print(io, "  R₀ = ", equ.R₀, "\n")
    print(io, "  B₀ = ", equ.B₀, "\n")
    print(io, "  q₀ = ", equ.q₀)
end

function X(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    R(x,equ) * cos(ϕ(x,equ))
end

function Y(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    R(x,equ) * sin(ϕ(x,equ))
end

function Z(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    x[2]
end

function R(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    x[1]
end

function r(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    sqrt(r²(x, equ))
end

function r²(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    (R(x,equ) - equ.R₀)^2 + Z(x,equ)^2
end

function θ(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    atan2(Z(x,equ), R(x,equ) - equ.R₀)
end

function ϕ(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    x[3]
end

function J(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    R(x,equ)
end

function periodicity(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    p = zero(x)
    p[3] = 2π
    return p
end

function A₁(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    + equ.B₀ * equ.R₀ * Z(x,equ) / R(x,equ) / 2
end

function A₂(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    - equ.B₀ * equ.R₀ * log(R(x,equ) / equ.R₀) / 2
end

function A₃(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    - equ.B₀ * r²(x,equ) / equ.q₀ / 2
end

function g₁₁(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    one(T)
end

function g₂₂(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    one(T)
end

function g₃₃(x::AbstractArray{T,1}, equ::AxisymmetricTokamakCylindrical) where {T <: Number}
    R(x, equ)^2
end

macro axisymmetric_tokamak_equilibrium_cylindrical(R₀, B₀, q₀)
    generate_equilibrium_code(AxisymmetricTokamakCylindrical(R₀, B₀, q₀); output=false)
end


text:
# This is an example of how to include a dependency on code
# outside of an action.

using Maker
using Base.Test

x = [ 3.14159   6.28319   9.42478
      12.5664   15.708    18.8496
      21.9911   25.1327   28.2743 ]
writecsv("in1.csv", x)

COUNT = 0

f() = readcsv("in1.csv")

Maker.file("in1.csv")

Maker.variable("x", ["in1.csv", "ffun"]) do
    global COUNT += 1
    f()
end

Maker.variable("ffun") do t
    Maker.funhash(f)
end

Maker.variable("y", "x") do
    global COUNT += 1
    3x
end

make("y")
@test COUNT == 2
make("y")
@test COUNT == 2
# redefine f
f() = readcsv("in1.csv")
make("y")
@test COUNT == 4

rm("in1.csv")


text:
import Pkg
Pkg.add("Coverage")
using Coverage
Codecov.submit(Codecov.process_folder())


text:
module CUDNN

using CUDA

include("../libcudnn-5/libcudnn.jl")
include("../libcudnn-5/libcudnn_types.jl")

@windows? (
  begin
    const libcudnn = Libdl.find_library(["cudnn64_5"])
  end : begin
    const libcudnn = Libdl.find_library(["libcudnn"])
  end)
isempty(libcudnn) && throw("CUDNN library cannot be found.")

function checkstatus(status)
  if status != CUDNN_STATUS_SUCCESS
    Base.show_backtrace(STDOUT, backtrace())
    throw(bytestring(cudnnGetErrorString(status)))
  end
end

datatype(::Type{Float32}) = CUDNN_DATA_FLOAT
datatype(::Type{Float64}) = CUDNN_DATA_DOUBLE
datatype(::Type{Float16}) = CUDNN_DATA_HALF

include("handle.jl")
include("tensor.jl")
include("activation.jl")
include("convolution.jl")
include("filter.jl")
include("softmax.jl")

end


text:
X = collect(Float64, reshape(1:40, 10, 4))
x = rand(10) * 10

D = DataFrame(A=rand(10), B=collect(1:10), C=[hex(x) for x in 11:20])
D_NA = deepcopy(D)
D_NA[1, :A] = NA

@testset "Array" begin
    XX = deepcopy(X)
    mu = center!(XX, obsdim=1)
    @test sum(abs.(mean(XX, 1))) == 0
    @test all(std(XX, 1) .== std(X, 1))
    @test all(mu .== vec(mean(X, 1)))

    XX = deepcopy(X)
    mu = center!(XX, ObsDim.First())
    @test sum(abs.(mean(XX, 1))) == 0
    @test all(std(XX, 1) .== std(X, 1))
    @test all(mu .== vec(mean(X, 1)))

    XX = deepcopy(X)
    mu = center!(XX, ObsDim.Last())
    @test sum(abs.(mean(XX, 2))) == 0
    @test all(std(XX, 2) .== std(X, 2))
    @test all(mu .== vec(mean(X, 2)))

    XX = deepcopy(X)
    mu = center!(XX)
    @test sum(abs.(mean(XX, 2))) == 0
    @test all(std(XX, 2) .== std(X, 2))
    @test all(mu .== vec(mean(X, 2)))

    XX = deepcopy(X)
    mu = vec(mean(X, 1))
    center!(XX, mu, obsdim=1)
    @test sum(abs.(mean(XX, 1))) == 0
    @test all(std(XX, 1) .== std(X, 1))

    XX = deepcopy(X)
    mu = vec(mean(X, 1))
    center!(XX, mu, ObsDim.First())
    @test sum(abs.(mean(XX, 1))) == 0
    @test all(std(XX, 1) .== std(X, 1))

    XX = deepcopy(X)
    mu = vec(mean(XX, 2))
    center!(XX, mu, obsdim=2)
    @test sum(abs.(mean(XX, 2))) == 0
    @test all(std(XX, 2) .== std(X, 2))

    XX = deepcopy(X)
    mu = vec(mean(XX, 2))
    center!(XX, mu, ObsDim.Last())
    @test sum(abs.(mean(XX, 2))) == 0
    @test all(std(XX, 2) .== std(X, 2))

    XX = deepcopy(X)
    mu = vec(mean(X[:,[1,3]], 1))
    center!(XX, mu, obsdim=1, operate_on=[1, 3])
    @test sum(abs.(mean(XX[:,[1,3]], 1))) == 0
    @test all(XX[:,2] .== X[:,2])
    @test all(std(XX, 1) .== std(X, 1))

    XX = deepcopy(X)
    mu = vec(mean(X[[1,3],:], 2))
    center!(XX, mu, obsdim=2, operate_on=[1, 3])
    @test sum(abs.(mean(XX[[1,3],:], 2))) == 0
    @test all(XX[2,:] .== X[2,:])
    @test all(std(XX, 2) .== std(X, 2))
    println()

    xx = deepcopy(x)
    center!(xx)
    @test mean(xx) <= 10e-10

    xx = deepcopy(x)
    mu = mean(xx)
    center!(xx, mu)
    @test mean(xx) <= 10e-10

    xx = deepcopy(x)
    mu = ones(xx)
    center!(xx, mu)
    @test mean(xx) - mean(x) ≈ -1

    xx = deepcopy(x)
    mu = ones(xx)
    center!(xx, mu)
    @test mean(xx) - mean(x) ≈ -1
end

@testset "DataFrame" begin
    # Center DataFrame
    DD = deepcopy(D)
    center!(DD)
    @test abs.(mean(DD[:A])) <= 10e-10
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D)
    center!(DD, operate_on=[:B])
    @test all(DD[:A] .== D[:A])
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D)
    mu = center!(DD, operate_on=[:A, :B])
    @test abs.(mean(DD[:A])) <= 10e-10
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[:C] .== D[:C])
    @test all(mu .== [mean(D[:A]), mean(D[:B])])

    DD = deepcopy(D_NA)
    mu = center!(DD, [:A, :B, :C])
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[2:end, :A] .== D[2:end, :A])
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D_NA)
    mu = [0.0, mean(DD[:B]), 0.0]
    mu = center!(DD, mu, [:A, :B, :C])
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[2:end, :A] .== D[2:end, :A])
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D_NA)
    mu = 0.0
    mu = center!(DD, mu, :A)
    @test all(DD[2:end, :A] .== D[2:end, :A])
    @test all(DD[:B] .== D[:B])
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D)
    mu =  [mean(D[:A]), mean(D[:B])]
    @test all(center!(DD, mu, operate_on=[:A, :B]) .== mu)
    @test abs.(mean(DD[:A])) <= 10e-10
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[:C] .== D[:C])

    DD = deepcopy(D_NA)
    center!(DD)
    @test all(DD[2:end, :A] .== D[2:end, :A])
    @test abs.(mean(DD[:B])) <= 10e-10
    @test all(DD[:C] .== D[:C])
    @test isna(DD[1, :A])
end


text:
#ДАНО: Робот - в произвольной клетке ограниченного прямоугольного поля
#РЕЗУЛЬТАТ: Робот - в исходном положении, и все клетки поля промакированы

function kray(r,side) #функция перемещает робота к краю
    while isborder(r,side) == false
        move!(r,side)
    end
end

function moverobot(r,side)
    while isborder(r,side) == false
        move!(r,side)
        putmarker!(r)
    end
end
function up(r)
    if isborder(r,Nord) == false
        move!(r,Nord)
        putmarker!(r)
    end
end

function mpole(r)
    kray(r,Sud)
    kray(r,Ost) #ведет робота в угол
    putmarker!(r)
    while isborder(r,Nord) == false
        moverobot(r,West)
        up(r)
        moverobot(r,Ost)
        up(r)
    end
    if isborder(r,Ost) == false
        moverobot(r,Ost)
    elseif isborder(r,West) == false
        moverobot(r,West)
    end
end

max_stars_repo_name: gardend/mirea-progs

max_issues_repo_name: gardend/mirea-progs

max_forks_repo_name: gardend/mirea-progs


text:
<extend src="layout.jl">
    <block name="content">
        <include src="post-item.jl"/>

            There are no posts.

        <post-item
            for-each=paginator.current.data
            @post=item
            class="post-info" />

            More

    </block>
</extend>


text:
@testset "686.repeated-string-match.jl" begin
    @test repeated_string_match("abcd", "cdabcdab") == 3
    @test repeated_string_match("a", "aa") == 2
    @test repeated_string_match("a", "a") == 1
    @test repeated_string_match("abc", "xyzs") == -1
end


text:
acceptance(v::AbstractArray{Bool}) = mean(v)

function acceptance(v::AbstractArray)
  accepted = 1
  n = length(v)

  for i in 2:n
    if v[i] != v[i-1]
      accepted += 1
    end
  end

  return accepted/n
end

acceptance(v::AbstractArray, region) = mapslices(acceptance, v, region)

function acceptance(s::UnivariateParameterNState; key::Symbol=:accept, diagnostics::Bool=true)
  if diagnostics
    return acceptance(convert(Array{Bool, 2}, s.diagnosticvalues[findfirst(s.diagnostickeys, key), :]))
  else
    return acceptance(s.value)
  end
end

acceptance(s::VariableNState{Univariate}) = acceptance(s.value)

function acceptance(s::MultivariateParameterNState; key::Symbol=:accept, diagnostics::Bool=true)
  if diagnostics
    return acceptance(convert(Array{Bool, 2}, s.diagnosticvalues[findfirst(s.diagnostickeys, key), :]))
  else
    return acceptance(Any[s.value[:, i] for i in 1:s.n])
  end
end

acceptance(s::VariableNState{Multivariate}) = acceptance(Any[s.value[:, i] for i in 1:s.n])


text:
function initialize()
  println("begin.")
  global t = t_start
  global t_output = output_interval
  global output_count = 1

  global lbox = x_max - x_min
  global sph = Array{Ptype,1}(undef, Npart)

  make_shock_tube()

  initialize_exact_riemann_solver()

  global p_initial = 0
  global E_initial = 0
  for i = 1:length(sph)
    p_initial += sph[i].p
    E_initial += sph[i].U + 0.5 * sph[i].p^2 / sph[i].m
  end

  Nthreads = Threads.nthreads()
  println("Number of threads: $Nthreads")

  global hstfilename = @sprintf("%s_%s_%s_%s_%s_%s_%s.hst", outputfile, formulation, gradient, time_integrator, kernel, volume_element, time_dependent_viscosity)
  open(hstfilename, "w") do io
    println(io, "### time dt p_error E_error L1_error")
  end
end

function make_shock_tube()
  mass_left = rho_left * (x_center - x_min)
  mass_right = rho_right * (x_max - x_center)
  mass_total = mass_left + mass_right

  particle_mass = mass_total / Npart

  Nleft = ceil(mass_left / particle_mass)

  interval_left = particle_mass / rho_left
  interval_right = particle_mass / rho_right

  for i = 1:Npart
    if i < Nleft
      x = x_min + interval_left * (i - 0.5)
      rho = rho_left
      v = v_left
      P = P_left
    else
      x = x_center + interval_right * ((i - Nleft) - 0.5)
      rho = rho_right
      v = v_right
      P = P_right
    end
    m = particle_mass
    sph[i] = Ptype(x, v, m, rho, P)
  end

  global active_particle = 1:Npart
end

function initialize_exact_riemann_solver()
  global P_s = P_star()
  global v_s = v_star(P_s)

  step_riemann = lbox / Nsample_riemann
  x_riemann = x_min:step_riemann:x_max
  rho_riemann_array = map(x -> rho_riemann(x, t_end, v_s, P_s, x_center), x_riemann)
  v_riemann_array = map(x -> v_riemann(x, t_end, v_s, P_s, x_center), x_riemann)
  P_riemann_array = map(x -> P_riemann(x, t_end, v_s, P_s, x_center), x_riemann)

  # y-axis limits
  global rho_max = maximum(rho_riemann_array)
  global rho_min = minimum(rho_riemann_array)
  global v_max = maximum(v_riemann_array)
  global v_min = minimum(v_riemann_array)
  global P_max = maximum(P_riemann_array)
  global P_min = minimum(P_riemann_array)
  rho_max = rho_max + 0.1 * (rho_max - rho_min)
  rho_min = rho_min - 0.1 * (rho_max - rho_min)
  v_max = v_max + 0.1 * (v_max - v_min)
  v_min = v_min - 0.1 * (v_max - v_min)
  P_max = P_max + 0.1 * (P_max - P_min)
  P_min = P_min - 0.1 * (P_max - P_min)
end


text:
const onex, oney, onexy = CartesianIndices((0:1, 0:1))[2:end]

function sps(diags::AbstractMatrix{Bool}, d=Matrix{Int}(undef, size(diags) .+ 1))
    d[:, 1] = 0:size(diags, 1)
    d[1, :] = 0:size(diags, 2)
    for i in CartesianIndices(diags)
        @inbounds d[i + onexy] = 1 + (diags[i] ? d[i] : min(d[i+onex], d[i+oney]))
    end
    d
end

@assert Int === Int64
const Distance = Int64

b2d(b::Bool)::Distance = (2 - b) << 32
d2b(d::Distance)::Bool = 2 - (d >> 32)
i2d(i::Int, j::Int=0)::Distance = (i << 32) + j

d2i(d::Distance, ::RoundingMode{:Down})::Int = d >> 32
d2i(d::Distance, ::RoundingMode{:Up})::Int = (d >> 32) + (mod(d, 2^32) > 0)

disorder(rng::StableRNG, d::Distance)::Distance =
    d2i(d,  RoundDown) == 2 ? d : i2d(1, rand(rng, 1:255))

function sps(dd::AbstractMatrix{Distance}, d=Matrix{Distance}(undef, size(dd) .+ 1))
    d[:, 1] .= i2d.(0:size(dd, 1))
    d[1, :] .= i2d.(0:size(dd, 2))
    for i in CartesianIndices(dd)
        @inbounds d[i + onexy] = min(dd[i] + d[i], i2d(1) + min(d[i+onex], d[i+oney]))
    end
    d
end

struct Grid
    dd::Matrix{Distance}
    add::Matrix{Distance}

    Grid(::UndefInitializer, n::Int) =
        new(Matrix{Distance}(undef, n-1, n-1), Matrix{Distance}(undef, n-1, n-1))

    function Grid(diags::AbstractMatrix{Bool}, antidiags::AbstractMatrix{Bool})
        checksquare(diags) == checksquare(antidiags) || throw(DimensionMismatch(
            "size(diags) $(size(diags)) != size(antidiags) $(size(antidiags))"))
        new(b2d.(diags), b2d.(antidiags))
    end
end

Grid(diags::AbstractMatrix{Bool}, n::Int=checksquare(diags)+1) =
    Grid(view(diags, onexy:onexy*(n-1)), view(diags, n-1:-1:1, 1:n-1))

Base.getproperty(g::Grid, s::Symbol) =
    s === :n ? size(g.dd, 1) + 1 :
    s == :diags ? d2b.(g.dd) :
    s == :antidiags ? d2b.(g.add) :
    getfield(g, s)

chessboard(n::Int) = Grid(trues(n-1, n-1), trues(n-1, n-1))
manhattan(n::Int) = Grid(falses(n-1, n-1), falses(n-1, n-1))

function disorder(rng::StableRNG, g::Grid)::Grid
    disordered = Grid(undef, g.n)
    disordered.dd .= disorder.(rng, g.dd)
    disordered.add .= disorder.(rng, g.add)
    disordered
end

isplanar(g::Grid)::Bool = !any(zip(g.dd, g.add)) do (d, ad); d2b(d) && d2b(ad); end

const Vertex = CartesianIndex{2}

Base.clamp(v::Vertex, g::Grid) = Vertex(clamp.(v.I, 1, g.n))

vertices(n::Int) = CartesianIndices((n, n))
vertices(g::Grid) = vertices(g.n)

function sps(g::Grid, v::Vertex, m=g.n)
    br = clamp(v+onexy*m, g)
    tl = clamp(v-onexy*m, g)
    bl = clamp(v+onex*m-oney*m, g)
    tr = clamp(v-onex*m+oney*m, g)

    d = OffsetArrays.OffsetArray(fill(i2d(Int(typemax(Int32))), size(tl:br)), tl:br)

    sps(view(g.dd, v:br-onexy), view(d, v:br))
    sps(view(g.dd, v-onexy:-onexy:tl), view(d, v:-onexy:tl))
    sps(view(g.add, v-oney:onex-oney:bl-onex), view(d, v:onex-oney:bl))
    sps(view(g.add, v-onex:oney-onex:tr-oney), view(d, v:oney-onex:tr))

    d
end

distance(g::Grid, u::Vertex, v::Vertex) = sps(g, u, maximum(abs.(u.I .- v.I)))[v]

eccentricity(g::Grid, v::Vertex, dv=sps(g, v))::Int = d2i(maximum(dv), RoundDown)

euclidean_eccentricity(n::Int, v::Vertex)::Float64 =
    maximum([sqrt(sum((u - v).I.^2)) for u in onexy:onexy*(n-1):onexy*n])

HG(g::Grid, v::Vertex, dv=sps(g, v))::Float64 =
    (eccentricity(g, v, dv) - euclidean_eccentricity(g.n, v)) / log2(g.n^2)
# nb. denominator ought to be sqrt(g.n); gets corrected by draw_hg() in figures.jl

geodesics(g::Grid, u::Vertex, v::Vertex,
          du=sps(g, u, sum(abs.(u.I .- v.I))),
          dv=sps(g, v, sum(abs.(u.I .- v.I))))::BitMatrix =
    [du[i] .+ dv[i] == du[v] for i in intersect(CartesianIndices(du), CartesianIndices(dv))]

NG(g::Grid, u::Vertex, v::Vertex, du=sps(g, u), dv=sps(g, v, d2i(du[v], RoundDown)))::Int =
    sum(geodesics(g, u, v, du, dv))

circle_points(g::Grid, v::Vertex, r::Int, dv=sps(g, v, r)) =
    findall(d->d2i(d, RoundDown) == r, dv)

function midpoints(g::Grid, u::Vertex, v::Vertex, du=sps(g, u),
                   dv=sps(g, v, div(d2i(du[v], RoundDown), 2, RoundUp)))
    duv = d2i(du[v], RoundDown)
    half_duv = div(duv, 2)
    pts = intersect!(circle_points(g, u, half_duv, du),
                     circle_points(g, v, duv - half_duv, dv))
    if duv % 2 == 1
        union!(pts, intersect!(circle_points(g, u, half_duv + 1, du),
                               circle_points(g, v, duv - half_duv - 1, dv)))
    end
    pts
end

function two_circle_points(g::Grid, u::Vertex, v::Vertex, r::Int,
                           du=sps(g, u, r+1), dv=sps(g, v, r+1); strict=true)
    ur = circle_points(g, u, r, du)
    vr = circle_points(g, v, r, dv)
    pts = intersect(ur, vr)
    (!isempty(pts) || strict) && return pts
    union!(
        intersect!(ur, union!(circle_points(g, v, r+1, dv), circle_points(g, v, r-1, dv))),
        intersect!(union!(circle_points(g, u, r+1, du), circle_points(g, u, r-1, du)), vr))
end

TG(g::Grid, d::Int, dmid::Int)::Float64 = (dmid / d - sqrt(3) / 2)

function TG(rng::StableRNG, g::Grid, u::Vertex, v::Vertex, d::Int,
            du=sps(g, u, d+1), dv=sps(g, v, d+1))
    w = rand(rng, two_circle_points(g, u, v, d, du, dv, strict=false))
    @assert w != v

    m = rand(rng, midpoints(g, u, v, du, dv))
    dmid = d2i(distance(g, m, w), RoundDown)

    TG(g, d, dmid)
end

function graphl(g::Grid, k::Int, xhat::Vertex, yhat::Vertex, u::Vertex, du=sps(g, u, k))::Int
    v = u + k*xhat
    for i in 1:k
        d2i(du[v + i*yhat], RoundDown) > k && return i-1
    end
    k
end

function euclideanl(k::Int)::Int
    k2 = k^2
    for i in 1:k
        d = sqrt(k2+i^2)
        @assert abs(d - k) != abs(d - k - 1)
        abs(d - k) >= abs(d - k - 1) && return i-1
    end
    k
end

midpoints(d::Int) = unique!([(div(d, 2, r1), div(d, 2, r2))
                             for (r1, r2) in ((RoundUp, RoundDown), (RoundDown, RoundUp))])

function AG(rng::StableRNG, g::Grid, d::Int, u::Vertex, du=sps(g, u, d+1))::Float64
    v = u + d*rand(rng, (onex, -onex, oney, -oney))
    TG(rng, g, u, v, d, du)
end

function diag_arcs(diags::AbstractMatrix{Bool}, reps::Int, r::Int)::BitMatrix
    step = Int(checksquare(diags) / reps)
    plane = BitMatrix(undef, size(diags))
    for i in onexy:onexy * step:onexy * step * reps
        box = i:i+onexy*(step-1)
        plane[box] .= sps(view(diags, box))[1:step,1:step] .<= r
    end
    plane
end

function grid_circs(g::Grid, reps::Int, r::Int)::BitMatrix
    step = Int((g.n-1) / reps)
    plane = falses(g.n, g.n)
    d = fill(i2d(Int(typemax(Int32))), (g.n, g.n))
    for i in onexy:onexy * step:onexy * step * reps
        box = i:i+onexy*(step-1)
        center = i+onexy*div(step, 2)
        plane[box] .= sps(g, center, r, d)[box] .<= i2d(r)
    end
    plane
end

function wave(diags::AbstractMatrix{Bool})
    d = sps(diags)
    n = checksquare(d)
    [d[i[1], i[2]] .< i[3] for i in CartesianIndices((n, n, n))]
end

max_stars_repo_name: moshelooks/eugrid

max_issues_repo_name: moshelooks/eugrid

max_forks_repo_name: moshelooks/eugrid


text:
# # Optimizing Scene Parameters using Optim.jl
#
# In this tutorial we will explore the exact same problem as demonstrated
# in [Inverse Lighting Tutorial](@ref) but this time we will use the
# Optimization Package [Optim.jl](https://julianlsolvers.github.io/Optim.jl/stable/).
# I would recommend going through a few of the
# [tutorials on Optim](https://julianlsolvers.github.io/Optim.jl/stable/#user/minimization/#_top)
# before starting this one.
#
# If you have already read the previous tutorial, you can safely skip to
# [Writing the Optimization Loop using Optim](@ref). The part previous to this
# is same as the previous tutorial.

using RayTracer, Images, Zygote, Flux, Statistics, Optim

# ## Script for setting up the Scene

screen_size = (w = 64, h = 64)

scene = load_obj("./tree.obj")

cam = Camera(
    Vec3(0.0f0, 6.0f0, -10.0f0),
    Vec3(0.0f0, 2.0f0,  0.0f0),
    Vec3(0.0f0, 1.0f0,  0.0f0),
    45.0f0,
    0.5f0,
    screen_size...
)

origin, direction = get_primary_rays(cam)

function render(light, scene)
    packed_image = raytrace(origin, direction, scene, light, origin, 2)
    array_image = reshape(hcat(packed_image.x, packed_image.y, packed_image.z),
                          (screen_size.w, screen_size.h, 3, 1))
    return array_image
end

showimg(img) = colorview(RGB, permutedims(img[:,:,:,1], (3,2,1)))

light_gt = PointLight(
    Vec3(1.0f0, 1.0f0, 1.0f0),
    20000.0f0,
    Vec3(1.0f0, 10.0f0, -50.0f0)
)

target_img = render(light_gt, scene)

showimg(zeroonenorm(render(light_gt, scene)))

light_guess = PointLight(
    Vec3(1.0f0, 1.0f0, 1.0f0),
    1.0f0,
    Vec3(-1.0f0, -10.0f0, -50.0f0)
)

showimg(zeroonenorm(render(light_guess, scene)))

# ## Writing the Optimization Loop using Optim
#
# Since, there is no direct support of Optim (unlike for Flux) in RayTracer
# the interface might seem a bit ugly. This is mainly due to the way the
# two optimization packages work. Flux allows inplace operation and ideally
# even RayTracer prefers that. But Optim requires us to give the parameters
# as an `AbstractArray`.
#
# Firstly, we shall extract the parameters, using the [`RayTracer.get_params`](@ref)
# function, we want to optimize.
initial_parameters = RayTracer.get_params(light_guess)[end-3:end]

# Since the input to the `loss_function` is an abstract array we need to
# convert it into a form that the RayTracer understands. For this we
# shall use the [`RayTracer.set_params!`](@ref) function which will modify the
# parameters inplace.
#
# In this function we simply compute the loss values and print it for our
# reference
function loss_function(θ::AbstractArray)
    light_optim = deepcopy(light_guess)
    RayTracer.set_params!(light_optim.intensity, θ[1:1])
    RayTracer.set_params!(light_optim.position, θ[2:end])
    loss = sum((render(light_optim, scene) .- target_img) .^ 2)
    @show loss
    return loss
end

# RayTracer uses Zygote's Reverse Mode AD for computing the derivatives.
# However, the default in Optim is ForwardDiff. Hence, we need to override
# that by giving our own gradient function.
function ∇loss_function!(G, θ::AbstractArray)
    light_optim = deepcopy(light_guess)
    RayTracer.set_params!(light_optim.intensity, θ[1:1])
    RayTracer.set_params!(light_optim.position, θ[2:end])
    gs = gradient(light_optim) do L
        sum((render(L, scene) .- target_img) .^ 2)
    end
    G .= RayTracer.get_params(gs[1])[end-3:end]
end

# Now we simply call the `optimize` function with `LBFGS` optimizer.
res = optimize(loss_function, ∇loss_function!, initial_parameters, LBFGS())

@show res.minimizer

# It might be interesting to note that convergence using LBFGS was much faster (only
# 252 iterations) compared to ADAM (401 iterations).

# If we generate a `gif` for the optimization process it will look similar to this
# ```@raw html
# <p align="center">
#      <img width=300 height=300 src="../../assets/inv_lighting_optim.gif">
# </p>
# ```


text:
"""
    UniqueIterator(itr, by=identity)

A lazy iterator over unique elements that have not been produced before.
"""
struct UniqueIterator{T,By}
    itr::T
    by::By
    UniqueIterator(itr::T, by::By=identity) where {T,By} = new{T,By}(itr, by)
end

function Base.iterate(itr::UniqueIterator)
    t = iterate(itr.itr)
    t === nothing && return nothing
    el, state = t
    # we store the hash instead of the object to allow GC to free it
    (el, (state, Set{UInt}(hash(itr.by(el)))))
end

function Base.iterate(itr::UniqueIterator, (state, cache))
    local el, state, h
    while true
        t = iterate(itr.itr, state)
        t === nothing && return nothing
        el, state = t
        h = hash(itr.by(el))
        h ∉ cache && break
    end
    push!(cache, h)
    return el, (state, cache)
end
Base.IteratorEltype(::Type{UniqueIterator{T,F}}) where {T,F} = Base.IteratorEltype(T)
Base.IteratorSize(::Type{UniqueIterator{T,F}}) where {T,F} = Base.SizeUnknown()
Base.eltype(::Type{UniqueIterator{T,F}}) where {T,F} = eltype(T)

Base.IteratorEltype(::Type{Base.Generator{T,F}}) where {T <: UniqueIterator, F} = IteratorEltype(T)
Base.eltype(::Type{Base.Generator{T,F}}) where {T <: UniqueIterator, F} = eltype(T)

iunique(itr; by=identity) = UniqueIterator(itr, by)
iunique(itr, itr2, itrs...; by=identity) = UniqueIterator(Flatten{eltype(itr)}(itr, itr2, itrs...), by)

"""
Flatten with inferrable `eltype`. Requires the given iterators to have the same `eltype`.
"""
struct Flatten{Eltype}
    it
end
Flatten{T}(itrs...) where T = Flatten{T}(itrs)
Base.IteratorEltype(::Type{<:Flatten}) = Base.HasEltype()
Base.IteratorSize(::Type{<:Flatten}) = Base.SizeUnknown()
Base.eltype(::Type{Flatten{T}}) where T = T
function Base.iterate(f::Flatten{T}, state=()) where T
    if state !== ()
        y = iterate(Base.tail(state)...)
        y !== nothing && return (y[1], (state[1], state[2], y[2]))
    end
    x = (state === () ? iterate(f.it) : iterate(f.it, state[1]))
    x === nothing && return nothing
    y = iterate(x[1])
    while y === nothing
         x = iterate(f.it, x[2])
         x === nothing && return nothing
         y = iterate(x[1])
    end
    return convert(T, y[1]), (x[2], x[1], y[2])
end

struct Shuffle{T}
    itr::T
    cacheSize::Int
    Shuffle(itr::T, cacheSize::Int=min(100, Base.IteratorSize(T) isa Base.HasShape || Base.IteratorSize(T) isa Base.HasLength ? length(s.itr)÷5 : 5)) where T = new{T}(itr, cacheSize)
end
Shuffle(v::AbstractArray) = shuffle!(v) # we already have the memory, so shuffle it directly

Base.IteratorEltype(::Shuffle{T}) where T = Base.IteratorEltype(T)
Base.IteratorSize(::Shuffle{T}) where T = Base.IteratorSize(T)
Base.eltype(::Shuffle{T}) where T = eltype(T)
Base.size(s::Shuffle{T}) where T = size(s.itr)
Base.length(s::Shuffle{T}) where T = length(s.itr)

function Base.iterate(s::Shuffle{T}) where T
    els = eltype(T)[]
    sizehint!(els, s.cacheSize)
    it = iterate(s.itr)
    while it !== nothing && length(els) < s.cacheSize
        el, state = it
        push!(els, el)
        it = iterate(s.itr, state)
    end
    shuffle!(els)
    iterate(s, (els, it))
end

function Base.iterate(s::Shuffle, (els, it))
    iszero(length(els)) && return nothing
    shuffle!(reverse!(els))
    ret = popfirst!(els)
    if it !== nothing
        el, state = it
        push!(els, el)
        it = iterate(s.itr, state)
    end
    return (ret, (els, it))
end


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

# # Density Fitting
#
# Density fitting is an extremely useful tool to reduce the computational scaling of many quantum chemical methods.
# Density fitting works by approximating the four-index electron repulsion integral (ERI) tensors from Hartree-Fock
# theory, $g_{\mu\nu\lambda\sigma} = (\mu\nu|\lambda\sigma)$, by
#
# $$(\mu\nu|\lambda\sigma) \approx \widetilde{(\mu\nu|P)}[J^{-1}]_{PQ}\widetilde{(Q|\lambda\sigma)}$$
#
# where the Coulomb metric $[J]_{PQ}$ and the three-index integral $\widetilde{(Q|\lambda\sigma)}$ are defined as
#
# \begin{align}
# [J]_{PQ} &= \int P(r_1)\frac{1}{r_{12}}Q(r_2){\rm d}^3r_1{\rm d}^3r_2\\
# \widetilde{(Q|\lambda\sigma)} &= \int Q(r_1)\frac{1}{r_{12}}\lambda(r_2)\sigma(r_2){\rm d}^3r_1{\rm d}^3r_2
# \end{align}
#
# To simplify the density fitting notation, the inverse Coulomb metric is typically folded into the three-index tensor:
#
# \begin{align}
# (P|\lambda\sigma) &= [J^{-\frac{1}{2}}]_{PQ}\widetilde{(Q|\lambda\sigma)}\\
# g_{\mu\nu\lambda\sigma} &\approx (\mu\nu|P)(P|\lambda\sigma)
# \end{align}
#
# These transformed three-index tensors can then be used to compute various quantities, including the four-index ERIs,
# as well as Coulomb (J) and exchange (K) matrices, and therefore the Fock matrix (F).  Before we go any further, let's
# see how to generate these transformed tensors using  Psi4.
#
# First, let's import  Psi4 and set up some global options, as well as
# define a molecule and initial wavefunction:

# +
# ==> Psi4 & NumPy options, Geometry Definition <==
using PyCall: pyimport
psi4 = pyimport("psi4")
np   = pyimport("numpy") # used only to cast to Psi4 arrays
using TensorOperations: @tensor

function psi4view(psi4matrix)
   # Assumes Float64 type, C ordering
   if !hasproperty(psi4matrix,:__array_interface__)
      throw(ArgumentError("Input matrix cannot be accessed. Cannot be an rvalue"))
   end
   array_interface = psi4matrix.__array_interface__
   array_interface["data"][2] == false   || @warn "Not writable"
   array_interface["strides"] == nothing || @warn "Different ordering than C"
   array_interface["typestr"] == "<f8"   || @warn "Not little-endian Float64 eltype"
   ptr   = array_interface["data"][1]
   shape = reverse(array_interface["shape"])
   ndims = length(shape)
   permutedims(unsafe_wrap(Array{Float64,ndims}, Ptr{Float64}(ptr), shape),reverse(1:ndims))
end

# Set Psi4 memory & output options
psi4.set_memory(Int(2e9))
psi4.core.set_output_file("output.dat", false)

# Geometry specification
mol = psi4.geometry("""
O
H 1 0.96
H 1 0.96 2 104.5
symmetry c1
""")

# Psi4 options
psi4.set_options(Dict("basis"         => "cc-pvdz",
                      "scf_type"      => "df",
                      "e_convergence" => 1e-10,
                      "d_convergence" => 1e-10))

wfn = psi4.core.Wavefunction.build(mol, psi4.core.get_global_option("basis"))
# -

# ## Building the Auxiliary Basis Set
#
# One component of our density-fitted tensors $g_{\mu\nu\lambda\sigma} \approx (\mu\nu|P)(P|\lambda\sigma)$ which
# is unique from their exact, canonical counterparts $(\mu\nu|\lambda\sigma)$ is the additional "auxiliary" index, $P$.
# This index corresponds to inserting a resolution of the identity, which is expanded in an auxiliary basis set $\{P\}$.
# In order to build our density-fitted integrals, we first need to generate this auxiliary basis set.  Fortunately,
# we can do this with the `psi4.core.BasisSet` object:
# ~~~julia
# # Build auxiliary basis set
# aux = psi4.core.BasisSet.build(mol, "DF_BASIS_SCF", "", "JKFIT", "aug-cc-pVDZ")
# ~~~
#
# There are special fitting basis sets that are optimal for a given orbital basis. As we will be building J and K
# objects we want the `JKFIT` basis associated with the orbital `aug-cc-pVDZ` basis. This basis is straightfowardly
# named `aug-cc-pVDZ-jkfit`.

# Build auxiliary basis set
aux = psi4.core.BasisSet.build(mol, "DF_BASIS_SCF", "", "JKFIT", "aug-cc-pVDZ")

# ## Building Density-Fitted ERIs
# Now, we can use our orbital and auxiliary basis sets to construct the `Qpq` object with the inverted metric. As the
# tensors are very similar to full ERI's we can use the same computer for both with the aid of a "zero basis". If we
# think carefully about the $\widetilde{(Q|\lambda\sigma)}$ and $(\mu\nu|\lambda\sigma)$ we should note that on the
# right and left hand sides the two gaussian basis functions are contracted to a single density.
#
# Specifically, for $\widetilde{(Q|\lambda\sigma)}$ the right hand side is a single basis function without being
# multiplied by another, so we can "trick" the MintsHelper object into computing these quanties if we have a "basis
# set" which effectively does not act on another. This is, effectively, what a "zero basis" does.
#
# The $[J^{-\frac{1}{2}}]_{PQ}$ object can be built in a similar way where we use the Psi4 Matrix's built in `power`
# function to raise this to the $-\frac{1}{2}$ power. The call `Matrix.power(-0.5, 1.e-14)` will invert the Matrix to
# the $-\frac{1}{2}$ while guarding against values smaller than 1.e-14. Recall that machine epsilon is ~1.e-16, when
# these small values are taken to a negative fractional power they could become very large and dominate the resulting
# matrix even though they are effectively noise before the inversion.
#
# ~~~julia
# orb = wfn.basisset()
# zero_bas = psi4.core.BasisSet.zero_ao_basis_set()
#
# # Build MintsHelper Instance
# mints = psi4.core.MintsHelper(orb)
#
# # Build (P|pq) raw 3-index ERIs, dimension (1, Naux, nbf, nbf)
# Ppq = mints.ao_eri(zero_bas, aux, orb, orb)
# Ppq = psi4view(Ppq)
#
# # Build and invert the metric
# metric = mints.ao_eri(zero_bas, aux, zero_bas, aux)
# metric.power(-0.5, 1.e-14)
# metric = psi4view(metric)
#
# # Remove the excess dimensions of Ppq & metric
# Ppq = dropdims(Ppq, dims=1)
# metric = dropdims(metric, dims=(1,3))
#
# # Contract Ppq & Metric to build Qso
# @tensor Qso[Q,p,q] := metric[Q,P] * Ppq[P,p,q]
# ~~~

# + code_folding=[]
# ==> Build Density-Fitted Integrals <==
# Get orbital basis & build zero basis
orb = wfn.basisset()
zero_bas = psi4.core.BasisSet.zero_ao_basis_set()

# Build instance of MintsHelper
mints = psi4.core.MintsHelper(orb)

# Build (P|pq) raw 3-index ERIs, dimension (1, Naux, nbf, nbf)
Ppq = mints.ao_eri(zero_bas, aux, orb, orb)
Ppq = psi4view(Ppq)

# Build & invert Coulomb metric, dimension (1, Naux, 1, Naux)
metric = mints.ao_eri(zero_bas, aux, zero_bas, aux)
metric.power(-0.5, 1.e-14)
metric = psi4view(metric)

# Remove excess dimensions of Ppq, & metric
Ppq = dropdims(Ppq, dims=1)
metric = dropdims(metric, dims=(1,3))

# Build the Qso object
@tensor Qpq[Q,p,q] := metric[Q,P] * Ppq[P,p,q];
# -

# ## Example: Building a Density-Fitted Fock Matrix
# Now that we've obtained our `Qpq` tensors, we may use them to build the Fock matrix.  To do so, since we aren't
# implementing a fully density-fitted RHF program, we'll first need to get a density matrix and one-electron Hamiltonian
# from somewhere. Let's get them from a converged HF wavefunction, so we can check our work later:

# ==> Compute SCF Wavefunction, Density Matrix, & 1-electron H <==
scf_e, scf_wfn = psi4.energy("scf", return_wfn=true)
D = np.asarray(scf_wfn.Da())
H = scf_wfn.H();

# Now that we have our density-fitted integrals and a density matrix, we can build a Fock matrix.  There are several
# different algorithms which we can successfully use to do so; for now, we'll use a simple algorithm and `@tensor`
# to illustrate how to perform contractions with these density fitted tensors and leave a detailed discussion of those
# algorithms/different tensor contraction methods elsewhere.  Recall that the Fock matrix, $F$, is given by
#
# $$F = H + 2J - K,$$
#
# where $H$ is the one-electron Hamiltonian matrix, $J$ is the Coulomb matrix, and $K$ is the exchange matrix.  The
# Coulomb and Exchange matrices have elements guven by
#
# \begin{align}
# J[D_{\lambda\sigma}]_{\mu\nu} &= (\mu\nu|\lambda\sigma)D_{\lambda\sigma}\\
# K[D_{\lambda\sigma}]_{\mu\nu} &= (\mu\lambda|\nu\sigma)D_{\lambda\sigma}.
# \end{align}
#
# When employing conventional 4-index ERI tensors, computing both $J$ and $K$ involves contracting over four unique
# indices, which involves four distinct loops -- one over each unique index in the contraction.  Therefore, the
# scaling of this procedure is $\mathcal{O}(N^4)$, where $N$ is the number of iterations in each loop (one for each
# basis function).  The above expressions can be coded using `@tensor` to handle the tensor contractions:
#
# ~~~julia
# @tensor J[p,q] := I[p,q,r,s] * D[r,s]
# @tensor K[p,q] := I[p,r,q,s] * D[r,s]
# ~~~
#
# for exact ERIs `I_pqrs`.  If we employ density fitting, however, we can reduce this scaling by reducing the number
# of unique indices involved in the contractions.  Substituting in the density-fitted $(P|\lambda\sigma)$ tensors into
# the above expressions, we obtain the following:
#
# \begin{align}
# J[D_{\lambda\sigma}]_{\mu\nu} &= (\mu\nu|P)(P|\lambda\sigma)D_{\lambda\sigma}\\
# K[D_{\lambda\sigma}]_{\mu\nu} &= (\mu\lambda|P)(P|\nu\sigma)D_{\lambda\sigma}.
# \end{align}
#
# Naively, this seems like we have actually *increased* the scaling of our algorithm, because we have added the $P$
# index to the expression, bringing the total to five unique indices, meaning this would scale like .  We've actually
# made our lives easier, however: with three different tensors to contract, we can perform one contraction at a time!
#
# For $J$, this works out to the following two-step procedure:
#
# \begin{align}
# \chi_P &= (P|\lambda\sigma)D_{\lambda\sigma} \\
# J[D_{\lambda\sigma}]_{\mu\nu} &= (\mu\nu|P)\chi_P
# \end{align}
#
#
# In the cell below, using `@tensor` and our `Qpq` tensor, try to construct `J`:

# Two-step build of J with Qpq and D
@tensor X_Q[Q] := Qpq[Q,p,q] * D[p,q]
@tensor J[p,q] := Qpq[Q,p,q] * X_Q[Q];

# Each of the above contractions, first constructing the `X_Q` intermediate and finally the full Coulomb matrix `J`, only involve three unique indices.  Therefore, the Coulomb matrix build above scales as $\mathcal{O}(N_{\rm aux}N^2)$.  Notice that we have distinguished the number of auxiliary ($N_{\rm aux}$) and orbital ($N$) basis functions; this is because auxiliary basis sets are usually around double the size of their corresponding orbital counterparts.
#
# We can play the same intermediate trick for building the Exchange matrix $K$:
#
# \begin{align}
# \zeta_{P\nu\lambda} &= (P|\nu\sigma)D_{\lambda\sigma} \\
# K[D_{\lambda\sigma}]_{\mu\nu} &= (\mu\lambda|P)\zeta_{P\nu\lambda}
# \end{align}
#
# Just like with $J$, try building $K$ in the cell below:

# Two-step build of K with Qpq and D
@tensor Z_Qqr[Q,r,q] := Qpq[Q,r,s] * D[s,q]
@tensor K[p,r]       := Qpq[Q,p,q] * Z_Qqr[Q,r,q];

# Unfortunately, our two-step $K$ build does not incur a reduction in the overall scaling of the algorithm, with each contraction above scaling as $\mathcal{O}(N^3N_{\rm aux})$. The major benefit of density fitting for $K$ builds comes in the form of the small storage overhead of the three-index `Qpq` tensors compared to the full four-index `I_pqrs` tensors.  Even when exploiting the full eight-fold symmetry of the $(\mu\nu|\lambda\sigma)$ integrals, storing `I_pqrs` for a system with 3000 AO basis functions will require 81 TB of space, compared to a mere 216 GB to store the full `Qpq` object when exploiting the twofold symmetry of $(P|\lambda\sigma)$.
#
# Now that we've built density-fitted versions of the $J$ and $K$ matrices, let's check our work by comparing a Fock matrix built using our $J$ and $K$ with the fully converged Fock matrix from our original SCF/aug-cc-pVDZ computation.
#
# Below, build F using the one-electron Hamiltonian from the converged SCF wavefuntion and our $J$ and $K$ matrices.  Then, get the converged $F$ from the SCF wavefunction:

# Build F from SCF 1 e- Hamiltonian and our density-fitted J & K
F = H + 2 * J - K
# Get converged Fock matrix from converged SCF wavefunction
scf_F = scf_wfn.Fa();

# Feeling lucky? Execute the next cell to see if you've computed $J$, $K$, and $F$ correctly:

if np.allclose(F, scf_F,atol=1e-4)
    println("Nicely done!! Your density-fitted Fock matrix matches Psi4!")
else
    println("Whoops...something went wrong.  Try again!")
end

# Finally, we can remember the identity of the $D$ matrix for SCF which will be $D_{\lambda\sigma} = C_{\lambda i} C_{\sigma i}$, where $i$ is the occupied index. We can factor our $K$ build once more:
# \begin{align}
# D_{\lambda\sigma} &= C_{\lambda i} C_{\sigma i} \\
# \zeta_{P\nu i} &= (P|\nu\sigma)C_{\sigma i} \\
# K[D_{\lambda\sigma}]_{\mu\nu} &= \zeta_{P\nu i}\zeta_{P\mu i}
# \end{align}
#
# Consider the ratio between the number of basis functions and the size of the occupied index. Why would the above be beneficial?

# ## References
# 1. F. Weigend, Phys. Chem. Chem. Phys. 4, 4285 (2002).
# 2. O. Vahtras, J. Alml ̈of, and M. W. Feyereisen, Chem. Phys. Lett. 213, 514 (1993).
# 3. B. I. Dunlap, J. W. D. Connolly, and J. R. Sabin, J. Chem. Phys. 71, 3396 (1979).
# 4. J. L. Whitten, J. Chem. Phys. 58, 4496 (1973).


text:
using Plots

s_date = 0
e_date = 31

function kwt_alfa(degalfa_p)
    eff = 0.15
    area = 30

    teta_p = deg2rad(170)
    alfa_p = deg2rad(degalfa_p)

    lat = deg2rad(59.3)

    delta_s(n) = deg2rad(-23.44) * cos((2pi / 365) * n)

    omega(t) = deg2rad((15 * t) - 180)

    alfa_s(n, t) = asin(sin(lat) * sin(delta_s(n)) + cos(lat) * cos(delta_s(n)) * cos(omega(t)))

    teta_s(n, t) = pi + (omega(t) / abs(omega(t))) * acos((sin(lat) * sin(alfa_s(n, t)) - sin(delta_s(n))) / (cos(lat) * cos(alfa_s(n, t))))

    I_0 = 1360
    I(n, t) = 1.1 * I_0 * 0.7^((1 / (abs(sin(alfa_s(n, t))) + sin(alfa_s(n, t))) / 2)^0.678)  # underlig abs medelvärde för att få 0 om negativt värde

    I_p(n, t) = I(n, t) * (cos(teta_p - teta_s(n, t)) * cos(alfa_p - alfa_s(n, t)) + sin(alfa_s(n, t)) * sin(alfa_p) * (1 - cos(teta_p - teta_s(n, t))))

    weather_loss = 1 - (1821 / (365 * 12))

    P(x) = (*(((abs(I_p(floor(x / 24), x)) + I_p(floor(x / 24), x)) / 2), weather_loss, eff, area))

    step = 0.1

    x = 24*s_date:step:24*(e_date)+24
    x = filter(x -> x % 12 != 0, x)

    wh = 0

    for i in x
        wh += step * P(i)
    end

    return wh / 1000
end

function maxer()
    values = Float64[]
    for i in 0:1:90
        push!(values, kwt_alfa(i))
    end
    return findmax(values)[2] - 1
end

maxvalues = Int[]
# january
push!(maxvalues, maxer())

#february
s_date = 32
e_date = 59
push!(maxvalues, maxer())

#march
s_date = 60
e_date = 90
push!(maxvalues, maxer())

#april
s_date = 91
e_date = 120
push!(maxvalues, maxer())

#may
s_date = 121
e_date = 151
push!(maxvalues, maxer())

#june
s_date = 152
e_date = 181
push!(maxvalues, maxer())

#july
s_date = 182
e_date = 212
push!(maxvalues, maxer())

#august
s_date = 213
e_date = 243
push!(maxvalues, maxer())

#september
s_date = 244
e_date = 273
push!(maxvalues, maxer())

#october
s_date = 274
e_date = 304
push!(maxvalues, maxer())

#november
s_date = 305
e_date = 334
push!(maxvalues, maxer())

#december
s_date = 335
e_date = 365
push!(maxvalues, maxer())

plot(1:12, maxvalues, label = "", xticks = 1:12)

xlabel!("Optimal elevationsvinkel[grader]")
ylabel!("Månad")

savefig("blub.svg")


text:
"""
    add_moment_constraints(
        ksas::Dict{String,Float64},
        protein_mass_fraction::Float64
    )

A modification that adds enzyme capacity constraints to the problem using a _modified_
version of the MOMENT algorithm. Requires specific activities, `ksas` [mmol product/g
enzyme/h], for each reaction. Proteins are identified by their associated gene IDs. Adds a
variable vector `y` to the problem corresponding to the protein concentration [g enzyme/gDW
cell] of each gene product in the order of `genes(model)`. The total protein concentration
[g protein/gDW cell] is constrained to be less than or equal to the `protein_mass_fraction`.
Reaction flux constraints are changed to the MOMENT constraints (see below) for all
reactions that have a gene reaction rule, otherwise the flux bounds are left unaltered.

See Adadi, Roi, et al. "Prediction of microbial growth rate versus biomass yield by a
metabolic network with kinetic parameters." PLoS computational biology (2012) for more
details of the original algorithm.

Here, a streamlined version of the algorithm is implemented to ensure that the correct units
are used. Specifically, this implementation uses specific activities instead of `kcats`.
Thus, for a reaction that can only proceed forward and is catalyzed by protein `a`, the flux
`x[i]` is bounded by `x[i] <= ksas[i] * y[a]`. If isozymes `a` or `b` catalyse the
reaction, then `x[i] <= ksas[i] * (y[a] + y[b])`. If a reaction is catalyzed by subunits `a`
and `b` then `x[i] <= ksas[i] * min(y[a], y[b])`. These rules are applied recursively in the
model like in the original algorithm. The enzyme capacity constraint is then implemented by
`sum(y) ≤ protein_mass_fraction`. The major benefit of using `ksas` instead of `kcats` is
that active site number and unit issues are prevented.

# Example
```
flux_balance_analysis(
    ...,
    modifications = [ add_moment_constraints(my_kcats, 0.6) ],
)
"""
add_moment_constraints(kcats::Dict{String,Float64}, protein_mass_fraction::Float64) =
    (model, opt_model) -> begin

        lbs, ubs = get_optmodel_bounds(opt_model) # to assign directions
        # get grrs and ignore empty blocks: TODO: fix importing to avoid this ugly conditional see #462
        grrs = Dict(
            rid => reaction_gene_association(model, rid) for
            rid in reactions(model) if !(
                reaction_gene_association(model, rid) == [[""]] ||
                isnothing(reaction_gene_association(model, rid))
            )
        )

        # add protein variables
        y = @variable(opt_model, y[1:n_genes(model)] >= 0)

        # add variables to deal with enzyme subunits
        num_temp = sum(length(values(grr)) for grr in values(grrs)) # number "AND" blocks in grrs
        t = @variable(opt_model, [1:num_temp]) # anonymous variable
        @constraint(opt_model, t .>= 0)

        #=
        Note, not all of t needs to be created, only those with OR GRR rules, however this
        adds a lot of complexity to the code - re-implement this method if efficiency becomes
        an issue.
        =#

        # add capacity constraint
        @constraint(opt_model, sum(y) <= protein_mass_fraction)

        k = 1 # counter
        kstart = 1
        kend = 1
        x = opt_model[:x]
        for (ridx, rid) in enumerate(reactions(model))

            isnothing(get(grrs, rid, nothing)) && continue # only apply MOMENT constraints to reactions with a valid GRR
            grrs[rid] == [[""]] && continue # TODO: remove once #462 is implemented

            # delete original constraints
            delete(opt_model, opt_model[:lbs][ridx])
            delete(opt_model, opt_model[:ubs][ridx])

            #=
            For multi-subunit enzymes, the flux is bounded by the minimum concentration of
            any subunit in the enzyme. If multiple isozymes with subunits exist, then the
            sum of these minima bound the enzyme flux. E.g., suppose that you have a GRR
            [[y, z], [w, u, v]], then x <= sum(min(y, z), min(w, u, v)).

            It is possible to reformulate x <= min(y, z) to preserve convexity:
            x <= t && t <= y && t <= z where t is a subunit variable.

            The remainder of the code implements these ideas.
            =#

            # build up the subunit variables
            kstart = k
            for grr in grrs[rid]
                for gid in grr
                    gidx = first(indexin([gid], genes(model)))
                    @constraint(opt_model, t[k] <= y[gidx])
                end
                k += 1
            end
            kend = k - 1

            # total enzyme concentration is sum of minimums
            isozymes = @expression(opt_model, kcats[rid] * sum(t[kstart:kend]))

            if lbs[ridx] >= 0 && ubs[ridx] > 0 # forward only
                @constraint(opt_model, x[ridx] <= isozymes)
                @constraint(opt_model, 0 <= x[ridx])
            elseif lbs[ridx] < 0 && ubs[ridx] <= 0 # reverse only
                @constraint(opt_model, -isozymes <= x[ridx])
                @constraint(opt_model, x[ridx] <= 0)
            elseif lbs[ridx] == ubs[ridx] == 0 # set to zero
                @constraint(opt_model, x[ridx] == 0)
            else # reversible
                @constraint(opt_model, x[ridx] <= isozymes)
                @constraint(opt_model, -isozymes <= x[ridx])
            end
        end
    end


text:
using JuMP

struct VariableBranch <: AbstractBranch
    lb::Dict{JuMP.VariableRef,Real}
    ub::Dict{JuMP.VariableRef,Real}
end

empty_variable_bound() = Dict{JuMP.VariableRef,Real}()

mutable struct JuMPNode{T<:AbstractBranch} <: AbstractNode
    @abstract_node_fields

    model::Union{Nothing,JuMP.Model}
    auxiliary_data::Dict

    function JuMPNode{T}(
        model = nothing,
        parent = nothing,
        depth = 0,
        bound = -Inf,
        id = -1,
        branch = nothing,
        solution_status = MOI.OPTIMIZE_NOT_CALLED,
        solution = Dict{Any,Real}(),
        auxiliary_data = Dict()) where T<:AbstractBranch
        return new{T}(id, parent, depth, branch, solution_status, bound, solution, model, auxiliary_data)
    end
end

JuMPNode{T}(parent::JuMPNode{T}) where T<:AbstractBranch = JuMPNode{T}(nothing, parent, parent.depth + 1, parent.bound)
JuMPNode(parent::JuMPNode{T}) where T<:AbstractBranch = JuMPNode{T}(parent)

# Initialize BNB tree
function initialize_tree(root::JuMPNode{T})::AbstractTree where T<:AbstractBranch
    tree = AbstractTree()
    push!(tree, root)
    return tree
end
initialize_tree(model::JuMP.Model, T = AbstractBranch)::AbstractTree = initialize_tree(JuMPNode{T}(model))

# Create child node with abstract branch object
function create_child_node(current_node::JuMPNode{T}, branch::AbstractBranch) where T<:AbstractBranch
    node = JuMPNode(current_node)
    node.branch = branch
    return node
end

# Create child node with variable bounds
create_child_node(current_node::JuMPNode{T}, variable::JuMP.VariableRef, lb::Real, ub::Real) where T<:AbstractBranch = create_child_node(current_node, VariableBranch(Dict(variable=>lb), Dict(variable=>ub)))
create_child_node_with_lb(current_node::JuMPNode{T}, variable::JuMP.VariableRef, lb::Real) where T<:AbstractBranch = create_child_node(current_node, VariableBranch(Dict(variable=>lb), empty_variable_bound()))
create_child_node_with_ub(current_node::JuMPNode{T}, variable::JuMP.VariableRef, ub::Real) where T<:AbstractBranch = create_child_node(current_node, VariableBranch(empty_variable_bound(), Dict(variable=>ub)))

function processed!(tree::AbstractTree, node::JuMPNode)
    # Rollback the branch objects
    adjust_branch!(node.auxiliary_data["bounds_changed"])

    Base.push!(tree.processed, node)
end

# basic bounding function
function bound!(node::JuMPNode)
    JuMP.optimize!(node.model)
    node.solution_status = JuMP.termination_status(node.model)

    if node.solution_status in [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
        node.bound = JuMP.objective_value(node.model)
        vrefs = JuMP.all_variables(node.model)
        for v in vrefs
            node.solution[v] = JuMP.value(v)
        end
    elseif node.solution_status in [MOI.INFEASIBLE,MOI.LOCALLY_INFEASIBLE]
        node.bound = -Inf
    elseif node.solution_status == MOI.DUAL_INFEASIBLE
        node.bound = Inf
    else
        @warn "Unexpected node solution status: $(node.solution_status)"
        node.bound = -Inf
    end
end

function apply_changes!(node::JuMPNode)
    branch_objects = AbstractBranch[]
    bounds_changed = AbstractBranch[]

    # Bracktracking tree to the root and collecting branch objects
    pnode = node
    push!(branch_objects, pnode.branch)
    while !isnothing(pnode.parent)
        pnode = pnode.parent
        push!(branch_objects, pnode.branch)
    end
    @assert isnothing(pnode.parent)
    @assert !isnothing(pnode.model)

    # Shallow copy of the root model
    node.model = pnode.model

    # Mark bounds to be changed
    mark_bound_changes!(node, branch_objects)

    # Apply branch objects
    adjust_branch!(branch_objects)
end

function adjust_branch!(branch_objects::Array{AbstractBranch,1})
    for branch in branch_objects
        for (j,v) in branch.lb
            JuMP.set_lower_bound(j, v)
        end
        for (j,v) in branch.ub
            JuMP.set_upper_bound(j, v)
        end
    end
end

function mark_bound_changes!(node::JuMPNode, branch_objects::Array{AbstractBranch,1})
    node.auxiliary_data["bounds_changed"] = AbstractBranch[]
    for branch in branch_objects
        lb_changed = Dict{JuMP.VariableRef,Real}()
        ub_changed = Dict{JuMP.VariableRef,Real}()
        for (j,v) in branch.lb
            lb_changed[j] = JuMP.has_lower_bound(j) ? JuMP.lower_bound(j) : -Inf
        end
        for (j,v) in branch.ub
            ub_changed[j] = JuMP.has_upper_bound(j) ? JuMP.upper_bound(j) : Inf
        end
        Base.push!(node.auxiliary_data["bounds_changed"], VariableBranch(lb_changed, ub_changed))
    end
end


text:
@testset "Particles" begin
    @testset "scalar settling velocity" begin
        T = transportoperator(grd, 100.0)
        @test T isa SparseMatrixCSC
        @test norm(v) / norm(T' * v) > ustrip(upreferred(1u"Myr"))
    end

    @testset "settling velocity as a function" begin
        w(z) = upreferred(100u"m/d" + z * u"m" / 100u"m" * 10u"m/d")
        @testset "no sed remin" begin
            T = transportoperator(grd, w)
            @test T isa SparseMatrixCSC
            @test norm(v) / norm(T' * v) > ustrip(upreferred(1u"Myr"))
        end

        @testset "fractional sed remin" begin
            T = transportoperator(grd, w, fsedremin=0.5)
            @test T isa SparseMatrixCSC
            @test norm(v) / norm(T' * v) < ustrip(upreferred(1u"yr"))
        end
    end

end


text:
using Parameters, GrowableArrays, StatsBase, SpecialFunctions, LinearAlgebra, Statistics, Random

#TODO: change bVAR --> VAR
# function bVAR(y::Array,p::Int64,i::Bool,reps::Int64,burn::Int64,max_try::Int64,prior::Hyperparameter)
#     i == false ? ((Y,X,β,ϵ,Σ,p) = fit_bvar(y,prior)) : ((Y,X,β,ϵ,Σ,p) = fit_bvar(y,Intercept()))
#     return VAR(Y,X,β,ϵ,Σ,p,Intercept())
# end

@with_kw struct Hyperparameter
    λ::AbstractRange{Float64} = 0.1:0.1:1
    τ::AbstractRange{Float64} = 10*(0.1:0.1:1)
    ε::Float64 = 0.0001
    p = 4
    H::Int64 = 10
    reps::Int64 = 10000
    burnin::Int64 = 8000
    max_try::Int64 = 1000
    update::Int64 = 1000
end

function get_prior(y::Array,p::Int64 = 1)
    T,K = size(y)
    μ = mean(y,dims=1)'
    σ = Array{Float64}(undef,0)
    δ = Array{Float64}(undef,0)
    e0 = Array{Float64}(undef,T-p,0)
    for i = 1:K
        ytemp = y[:,i]
        xtemp = [lagmatrix(ytemp,p) ones(length(ytemp)-p)]
        ytemp = ytemp[1+p:end,:]
        btemp=xtemp\ytemp
        etemp=ytemp-xtemp*btemp
        stemp=etemp'*etemp/length(ytemp)
        if abs(btemp[1]) > 1
            btemp[1] = 1
        end
        σ = vcat(σ,stemp)
        δ = vcat(δ,btemp[1])
        e0 = hcat(e0, etemp)
    end
    return μ, σ, δ, e0
end

function lagmatrix(x::Array{F},p::Int64) where F
    sk = 1
    T, K = size(x)
    k    = K*p+1
    idx  = repeat(1:K, p)
    X    = Array{F}(undef,T-p, k)
    # building X (t-1:t-p) allocating data from D matrix - avoid checking bounds
    for j = 1+sk:(sk+K*p)
        for i = 1:(T-p)
            lg = round(Int, ceil((j-sk)/K)) - 1 # create index [0 0 1 1 2 2 ...etc]
            @inbounds X[i, j] = x[i+p-1-lg, idx[j-sk]]
        end
    end
    return X[:,2:end]
end

function lagmatrix(x::Vector{F},p::Int64) where F
    sk = 1
    T = length(x)
    K = 1
    k    = K*p+1
    idx  = repeat(1:K, p)
    X    = Array{F}(undef,T-p, k)
    # building X (t-1:t-p) allocating data from D matrix - avoid checking bounds
    for j = 1+sk:(sk+K*p)
        for i = 1:(T-p)
            lg = round(Int, ceil((j-sk)/K)) - 1 # create index [0 0 1 1 2 2 ...etc]
            @inbounds X[i, j] = x[i+p-1-lg, idx[j-sk]]
        end
    end
    return X[:,2:end]
end

function create_dummies(λ::Float64,τ::Float64,δ::Array,ε::Float64,p::Int64,μ::Array,σ::Array,K::Int64)
    if typeof(σ) == Array{Float64,2}
        σ = σ[:,1]
    end
    if typeof(μ) == Array{Float64,2}
        μ = μ[:,1]
    end
    if λ > 0
        if ε > 0
            yd1 = [diagm(0 => σ.*δ)./λ;
            zeros(K*(p-1),K);
            diagm(0 => σ);
            zeros(1,K)]

            jp = diagm(0 => 1:p)
            xd1 = [hcat(kron(jp,diagm(0 => σ)./λ), zeros((K*p),1));
            zeros(K,(K*p)+1);
            hcat(zeros(1,K*p), ε)]
        else
            yd1 = [diagm(0 => σ.*δ)./λ;
            zeros(K*(p-1),K);
            diagm(0 => σ)]

            jp = diagm(0 => 1:p)
            xd1 = [kron(jp,diagm(0 => σ)./λ);
            zeros(K,(K*p))]
        end
    end
    if τ > 0
        if ε > 0
            yd2 = diagm(0 => δ.*μ)./τ
            xd2 = [kron(ones(1,p),yd2) zeros(K,1)]
        else
            yd2 = diagm(0 => δ.*μ)./τ
            xd2 = [kron(ones(1,p),yd2)]
        end
    end
    y = vcat(yd1,yd2)
    x = vcat(xd1,xd2)
    return y,x
end

function sum_loggamma(K::Int64,v::Int64)
    out = 0
    for i = 1:K
        out = out + logabsgamma((v+1-i)/2)[1]
    end
    return out
end

function mgamln(K::Int64,v::Int64)
    constant = (K*(K-1)/4)*log(pi)
    term2 = sum_loggamma(K,v)
    out = constant + term2
    return out
end

function max_lik_var(y::Array,x::Array,yd::Array,xd::Array)
    T,K = size(y)
    v = size(yd,1)
    y1 = vcat(y,yd)
    x1 = vcat(x,xd)
    #prior moments
    xx0 = xd'*xd
    invxx0 = pinv(xx0)
    b0 = invxx0*xd'*yd
    v0 = size(yd,1)
    e0 = yd-xd*b0
    sigma0 = e0'*e0
    #posterior moments
    xx1 = x1'*x1
    invxx = pinv(xx1)
    b = invxx*x1'*y1
    v1 = v0+T
    e = y1-x1*b
    sigma1 = e'*e

    PP = inv(Matrix(1.0I,T,T)+x*invxx0*x')
    QQ = sigma0

    lngam_ratio = mgamln(K,v0)-mgamln(K,v1)
    py = -(lngam_ratio+(T*K/2)*log(pi))+0.5*K*log(det(PP))+(v0/2)*log(det(QQ))-(v1/2)*log(det(QQ+(y-x*b0)'*PP*(y-x*b0)))
    out = py
    return out
end

function get_opt_lag_prior(y::Array,p::Int64,λ::AbstractRange{Float64},τ::AbstractRange{Float64},δ::Array,ε::Float64,μ::Array,σ::Array)
    K  =  size(y,2)
    outmlik = zeros(length(p),length(λ))
    tableL = zeros(length(p),length(λ))
    tableP = copy(tableL)
    tableD = copy(tableL)
    for i = 1:length(p)
        L = p[i]
        Y = copy(y)
        X = lagmatrix(y,L)
        X = hcat(X, ones(size(X,1)))
        Y = Y[L+1:end,:]
        for j = 1:length(λ)
            yd,xd  =  create_dummies(λ[j],τ[j],δ,ε,L,μ,σ,K)
            # get marginal likelihood
            mlik = max_lik_var(Y,X,yd,xd)
            outmlik[i,j] = mlik
            tableL[i,j] = L
            tableP[i,j] = λ[j]
            tableD[i,j] = τ[j]
        end
    end
    id = outmlik .== maximum(maximum(outmlik))
    optL = tableL[id]
    optP = tableP[id]
    optD = tableD[id]
    return optL, optP, optD
end

function get_coef(β::Array,sigma::Array,ixx::Array,max_try::Int64,K::Int64,p::Int64)
    problem = 0
    vstar = kron(sigma,ixx)
    check = -1
    num_try = 1
    control_stability = 0
    mβ = Array{}
    while check < 0 && num_try < max_try
        mβ = β + (randn(1,K*(K*p+1))*(cholesky(Hermitian(vstar)).U))'

        control_stability = check_stability(mβ,K,p)
        if control_stability == 0
            check = 10
        else
            num_try += 1
        end
    end
    if control_stability > 0
        problem = 1
    end
    return mβ,problem
end

function check_stability(β::Array,K::Int64,p::Int64)
    Kp = K*p
    K_p1 = K*(p-1)
    mCoef = zeros(Kp,Kp)
    mCoef[K+1:Kp,1:K_p1] = Matrix(1.0I,K_p1,K_p1)

    temp = reshape(β,Kp+1,K)
    temp = view(temp,1:Kp,1:K)'
    mCoef[1:K,1:Kp] = temp
    ee = maximum(abs.(eigvals(mCoef)))
    is_stable = ee >= 1
    return is_stable
end

function iwpq(v::Int64,ixpx::Array)
    k = size(ixpx,1)
    z = zeros(v,k)
    mu = zeros(k,1)
    for i = 1:v
        z[i,:]=(cholesky(Hermitian(ixpx)).U'*randn(k,1))'
    end
    out = inv(z'*z)
    return out
end

function get_forecast(K::Int64, H::Int64, T::Int64, p::Int64, Y::SubArray, β::Array, sigma::Array)
    # -------------------------------------------------------------------------
    # get_paths:
    # generates a matrix of simulated paths for Y for given parameters and
    # general set-up (lags, horizon, etc.)
    # -------------------------------------------------------------------------
    # Draw K(0,1) innovations for variance and mean equation:
    csigma = (cholesky(Hermitian(sigma))).U
    uu = randn(H+p,K)
    # Note we only need H*K innovations, but adding an extra p draws makes
    # the indexing in the loop below a bit cleaner.
    #compute forecast
    yhat = zeros(H+p,K)
    yhat[1:p,:] = Y[T-p+1:T,:]
    for fi = p+1:H+p
        xhat= Array{Float64}(undef,0)
        for ji = 1:p
            xhat = vcat(xhat,yhat[fi-ji,:])
        end
        xhat = vcat(xhat,1)
        yhat[fi,:] = xhat'*reshape(β,K*p+1,K) + uu[fi,:]'*csigma
    end
    return yhat
end

function gibbs!(Y::Array,yd::Array,xd::Array,p::Int64,prior::Hyperparameter,mForecast::Array)
    @unpack λ ,τ ,ε ,H ,reps ,burnin ,max_try ,update = prior
    x = lagmatrix(Y,p)
    T = size(x,1)
    K = size(Y,2)
    iT = size(yd,1)
    x = hcat(x, ones(T))
    y = view(Y,p+1:T+p,:)

    y0 = [y; yd]
    x0 = [x; xd]
    # conditional mean of the VAR coefficients
    β0_fix = vec(x0\y0) # ols on the appended data
    β0 = vec(x0\y0) # ols on the appended data
    xx = x0'*x0
    ixx = xx\Matrix(1.0I,size(xx,2),size(xx,2))  # inv(x0'x0) to be used later in the Gibbs sampling algorithm
    σ = Matrix(1.0I,K,K) # starting value for sigma
    igibbs = 1
    jgibbs = 0

    while jgibbs < reps-burnin
        # Display progress:
        if mod(igibbs,update)==0
            println("Replication $igibbs of $reps. Lag $p Prior tightness $λ")
        end
        # #step 1: Sample VAR coefficients
        β1, problem = get_coef(β0_fix,σ,ixx,max_try,K,p)
        if problem == 1
            β1 = β0
        else
            β0 = β1
        end
        #draw covariance
        e = y0-x0*reshape(β1,K*p+1,K)
        scale = e'*e
        v = T+iT
        inv_scale = inv(scale)
        σ = iwpq(v,inv_scale)
        if igibbs > burnin && problem == 0
            jgibbs += 1
            yhat = get_forecast(K, H+1, T, p, y, β1,σ)
            mForecast[jgibbs,:,:] = view(yhat,p+1:H+p,:)
        end
        igibbs += 1
    end
    return mForecast
end

function fit_bvar(y::Array,prior::Hyperparameter)
    @unpack λ ,τ ,ε ,p ,H ,reps ,burnin ,max_try ,update = prior
    T,K = size(y)
    μ, σ, δ, e0  = get_prior(y)
    p_star,λ_star,τ_star = get_opt_lag_prior(y,p,λ,τ,δ,ε,Array(μ),σ)
    yd,xd = create_dummies(λ_star[1],τ_star[1],δ,ε,Int(p_star[1]),Array(μ),σ,K)
    mForecast = zeros(reps-burnin,H,K)
    gibbs!(y,yd,xd,Int(p_star[1]),prior,mForecast)
    return mForecast
end

max_stars_repo_name: lucabrugnolini/VAR

max_issues_repo_name: lucabrugnolini/VAR

max_forks_repo_name: lucabrugnolini/VAR


text:
# String matching functions for Aaron
using StringDistances

# Fairly arbitrarily chosen list of stop words
const stopwords = [
    "the",
    "is",
    "at",
    "which",
    "on"
]

# This function is used as a first pass through the candidate
# string. In particular, we use it to determine which words
# in the candidate string should be removed because of their
# dissimilarity to all words in the target string.
function should_cutword(targ, cand, thresh)
    n = length(cand)

    if targ[1] != cand[1]
        res = true
    # need perfect match for strings of length ≤ 3
    elseif n ≤ 3
        # should allow distance of 1 on pruning
        res = compare(Levenshtein(), targ, cand) != 1.0
    elseif n ≤ 5
        res = compare(Winkler(Jaro()), targ, cand) ≤ 0.85
    else
        res = compare(DamerauLevenshtein(), targ, cand) ≤ thresh
    end
    return res
end

function clean_string(x)
    xlwr = lowercase(x)
    xcln = replace(xlwr, r"[^-a-z]", " ")
    return xcln
end

function cut_stopwords(x_arr, stopwords)
    res = [z for z in x_arr if z ∉ stopwords]
    return res
end

## Sample use
#cut_stopwords(["this", "is", "it"], ["it", "this"])

# Given two arrays, one with target phrase, the other with
# a candidate phrase, this function iterates through the
# words in the candidate array and drops those flagged by
# should_cutword() as being too dissimilar from any word in
# the target array. Note that a prior version of this function
# actually "handled" spell checking by REPLACING words in
# `cand_arr` with their "matched" word in `targ_arr`.

function exclude_words(targ_arr, cand_arr, thresh)
    kept_words = Array{String,1}(0)

    for x in targ_arr
        for y in cand_arr
            if !should_cutword(x, y, thresh)
                if y ∉ kept_words
                    println("Keeping $y")
                    push!(kept_words, y)    # don't keep duplicates
                end
            end
        end
    end
    return kept_words
end

function prune_candidate(targ, cand, thresh, stopwords)
    # iterate over every word in candidate string and if not
    # close enough to word in target string then is dropped
    targ_cln = clean_string(targ)
    cand_cln = clean_string(cand)
    targ_arr = cut_stopwords(split(targ_cln), stopwords)
    cand_arr = split(cand_cln)
    n = length(cand_arr)

    kept_words = exclude_words(targ_arr, cand_arr, thresh)

    return kept_words
end

# Given two strings, `x` and `y`, and two thresholds, this
# function returns a boolean indicating whether or not the two
# strings are approximates matches. The first threshold argument
# controls the similarity level that we want in the first pass
# where irrelevant words are removed. The second threshold is
# used as the criteria for declaring a match in the final
# pruned strings.

#=function is_string_match(x, y, thresh1)
    y_arr = prune_candidate(x, y, thresh1)
    y_cln = join(y_arr, " ")
    x_cln = clean_string(x)
    res = contains(y_cln, x_cln)
    return res
end
=#

function is_anagram(a, b)
    if length(a) ≠ length(b)
        res = false
    else
        a_sorted = sort(collect(a))
        b_sorted = sort(collect(b))
        res = a_sorted == b_sorted
    end
    return res
end

## Testing our function
# is_anagram("abc", "bac")

function string_match(targ, cand, thresh)
    n = length(cand)

    if targ[1] != cand[1]               # first letters much match
        res = true
    elseif n ≤ 3
        if is_anagram(targ, cand)       # anagrams w/ less than 4 letters match
            res = true
        else
            res = dist != 1.0
        end
    else
        similarity = compare(DamerauLevenshtein(), targ, cand)
        warn("dist: $similarity")
        res = similarity ≥ thresh
    end
    return res
end

# Improved is_string_match() function to include only checking
# for same length substrings

function is_string_match(x, y, thresh1, stopwords)
    cand_arr = prune_candidate(x, y, thresh1, stopwords) # create new short candidate array of correctly spelled words
    targ_cln = clean_string(x)
    cand_len = length(cand_arr)
    targ_len = length(split(targ_cln," "))

    res = false
    # Loops through the difference in lengths of the two strings to
    # create comparisons that are of the correct length
    len_diff = cand_len - targ_len
    println("length differences is: $len_diff")
    for i in 0:(len_diff)
        cand_cln = join(cand_arr[i+1:i+targ_len], " ")
        print(cand_cln)
        print("| |")
        println(targ_cln)

        if string_match(targ_cln, cand_cln, thresh1) ###*** needed different compare function to make exact match
            println("Matched $cand_cln to $targ_cln")
            # looking for a match in the correct order
            res = true
            break
        end
    end

    return res
end

# testing our functions
target1 = "pericardial effusion"
candidate1 = "Metastatic non-small cell lung cancer.  Dyspnea. No pulsus on exam. Evaluate for perracardial perracardial effusion/interval change, signs of tamponade."

#x = "pericardial effusion"
#y = "Metastatic non-small cell lung cancer.  Dyspnea. No pulsus on exam. Evaluate for perracardial perracardial effusion/interval change, signs of tamponade."
#thresh1 = 0.75

threshold1 = 0.75

is_string_match(target1, candidate1, threshold1, stopwords)

x1 = "left ventricular left ventricular hypertophy"
y1 = "left ventricular fxn"

is_string_match(x1, y1, 0.75, stopwords)


text:
module JFSM2

using DataFrames

# Functions for compilation

export write_modcfg, compile_fsm, cfg_table, compile_all

include("compil.jl")

end


text:
using CompScienceMeshes, BEAST
using LinearAlgebra

T = tetmeshsphere(1.0,0.25)
X = nedelecc3d(T)
Γ = boundary(T)

X = raviartthomas(Γ)
@show numfunctions(X)

κ,  η  = 1.0, 1.0
κ′, η′ = 2.0κ, η/2.0

T  = Maxwell3D.singlelayer(wavenumber=κ)
T′ = Maxwell3D.singlelayer(wavenumber=κ′)
K  = Maxwell3D.doublelayer(wavenumber=κ)
K′ = Maxwell3D.doublelayer(wavenumber=κ′)

E = Maxwell3D.planewave(direction=ẑ, polarization=x̂, wavenumber=κ)
H = -1/(im*κ*η)*curl(E)

e, h = (n × E) × n, (n × H) × n

@hilbertspace j m
@hilbertspace k l

α, α′ = 1/η, 1/η′
pmchwt = @discretise(
    (η*T+η′*T′)[k,j] -      (K+K′)[k,m] +
         (K+K′)[l,j] + (α*T+α′*T′)[l,m] == -e[k] - h[l],
    j∈X, m∈X, k∈X, l∈X)

# Z = BEAST.sysmatrix(pmchwt)
# b = BEAST.rhs(pmchwt)
# M = BEAST.convert_to_dense(Z)
# u = M \ b
# error()
# W = BEAST.assemble_hide(pmchwt.equation.lhs, pmchwt.test_space_dict, pmchwt.trial_space_dict)

# using BEAST.BlockArrays
# using BEAST.IterativeSolvers

# u = PseudoBlockVector{ComplexF64}(undef, numfunctions.([X,X]))
# fill!(u,0)

# error()
# u, ch = IterativeSolvers.gmres!(u, Z, b, log=true,  maxiter=1000,
#     restart=1000, reltol=1e-5, verbose=true)

u = solve(pmchwt)

Θ, Φ = range(0.0,stop=2π,length=100), 0.0
ffpoints = [point(cos(ϕ)*sin(θ), sin(ϕ)*sin(θ), cos(θ)) for θ in Θ for ϕ in Φ]

# Don't forget the far field comprises two contributions
ffm = potential(MWFarField3D(κ*im), ffpoints, u[m], X)
ffj = potential(MWFarField3D(κ*im), ffpoints, u[j], X)
ff = -η*im*κ*ffj + im*κ*cross.(ffpoints, ffm)

using Plots
plot(xlabel="theta")
plot!(Θ,norm.(ff),label="far field",title="PMCHWT")

#import Plotly
#using LinearAlgebra
#fcrj, _ = facecurrents(u[j],X)
#fcrm, _ = facecurrents(u[m],X)
#Plotly.plot(patch(Γ, norm.(fcrj)))
#Plotly.plot(patch(Γ, norm.(fcrm)))

function nearfield(um,uj,Xm,Xj,κ,η,points,
    Einc=(x->point(0,0,0)),
    Hinc=(x->point(0,0,0)))

    K = BEAST.MWDoubleLayerField3D(wavenumber=κ)
    T = BEAST.MWSingleLayerField3D(wavenumber=κ)

    Em = potential(K, points, um, Xm)
    Ej = potential(T, points, uj, Xj)
    E = -Em + η * Ej + Einc.(points)

    Hm = potential(T, points, um, Xm)
    Hj = potential(K, points, uj, Xj)
    H = 1/η*Hm + Hj + Hinc.(points)

    return E, H
end

Z = range(-1,1,length=100)
Y = range(-1,1,length=100)
nfpoints = [point(0,y,z) for z in Z, y in Y]

import Base.Threads: @spawn
task1 = @spawn nearfield(u[m],u[j],X,X,κ,η,nfpoints,E,H)
task2 = @spawn nearfield(-u[m],-u[j],X,X,κ′,η′,nfpoints)

E_ex, H_ex = fetch(task1)
E_in, H_in = fetch(task2)

E_tot = E_in + E_ex
H_tot = H_in + H_ex

contour(real.(getindex.(E_tot,1)))
contour(real.(getindex.(H_tot,2)))

heatmap(Z, Y, real.(getindex.(E_tot,1)))
#heatmap(Z, Y, real.(getindex.(H_tot,2)))

#plot(real.(getindex.(E_tot[:,51],1)))
#plot!(real.(getindex.(H_tot[:,51],2)))

# Compare the far field and the field far
ffradius = 100.0
E_far, H_far = nearfield(u[m],u[j],X,X,κ,η, ffradius .* ffpoints)
nxE_far = cross.(ffpoints, E_far) * (4π*ffradius) / exp(-im*κ*ffradius)
Et_far = -cross.(ffpoints, nxE_far)

plot()
plot!(Θ, norm.(ff),label="far field")
scatter!(Θ, norm.(Et_far), label="field far")


text:
# Time correlated noise

function timecorr_noise(q, timestep, timedecorr)

	alfa = 1 - timestep / timedecorr

	q = alfa*q + sqrt(1-alfa^2)*randn()

	return q

end

# Resampling function

function resample(wk)

	Ns = length(wk)

	u = cumprod(rand(Ns) .^ (1./collect(Float64, Ns:-1:1)))

	u = u[length(u):-1:1]

	wc = cumsum(wk)

	label = zeros(Int64, Ns)

	k = 1
	for i = 1:Ns
		while wc[k] < u[i]
			k = k + 1
		end
		label[i] = k
	end

	return label

end


text:
import Base: show

function show(io::IO, gr::Graph)
    print(io, "MathGL Graph of size $(width(gr))x$(height(gr))")
end

function image_string(gr::Graph, format)
    name = tempname()"_."string(format)
    write(gr, name)
    tmpfile = open(name, "r")
    str = read(tmpfile)
    close(tmpfile)
    str
end

function show(io::IO, ::MIME"image/png", gr::Graph)
    Base.write(io, image_string(gr, :png))
end

function show(io::IO, ::MIME"image/svg+xml", gr::Graph)
    Base.write(io, image_string(gr, :svg))
end

function show(io::IO, ::MIME"image/jpeg", gr::Graph)
    Base.write(io, image_string(gr, :jpg))
end

function show(io::IO, ::MIME"image/gif", gr::Graph)
    Base.write(io, image_string(gr, :gif))
end

function show(io::IO, ::MIME"image/bmp", gr::Graph)
    Base.write(io, image_string(gr, :bmp))
end

function show(io::IO, ::MIME"application/postscript", gr::Graph)
    Base.write(io, image_string(gr, :eps))
end