text:
import MPI

include("types.jl")
include("init.jl")
include("get_sim_params.jl")
include("update.jl")
include("io.jl")
include("get_sim_params.jl")

function main()
    #MPI Initialization
    MPI.Init()
    comm = MPI.COMM_WORLD
    my_id = MPI.Comm_rank(comm)
    nproc = MPI.Comm_size(comm)
    println("Number of MPI processes: ", nproc)

    #temp1_init: temperature init on borders
    temp1_init = 10.0
    #temp2_init: temperature init inside domain
    temp2_init = -10.0
    #diffusion coefficient
    k0 = Float64(1)

    #define rank of root process
    root = 0

    #simulation parameters
    params_int = Array{Int64}(undef, 5)
    params_double = Array{Float64}(undef, 2)
    if (my_id == root)
        params_int, params_double, output_path = get_sim_params(ARGS)
    end

    MPI.Barrier(comm)
    MPI.Bcast!(params_int, 5, 0, comm)
    MPI.Bcast!(params_double, 2, 0, comm)

    size_x    = params_int[1]
    size_y    = params_int[2]
    nx_domains = params_int[3]
    ny_domains = params_int[4]
    maxStep   = params_int[5]
    dt1       = params_double[1]
    epsilon   = params_double[2]

    #Warning message if dimensions and number of processes don't match
    if ((my_id == root) && (nproc != (nx_domains * ny_domains)))
        println("ERROR - Number of processes not equal to number of subdomains")
        println("nprocs: ", nproc, " domains: ",  nx_domains * ny_domains )
    end

    #Various other variables
    size_global_x = size_x + 2
    size_global_y = size_y + 2
    hx = Float64(1.0 / size_global_x)
    hy = Float64(1.0 / size_global_y)
    dt2 = 0.25 * min(hx, hy)^2 / k0 #fraction of CFL condition
    size_total_x = size_x + 2 * nx_domains + 2 #including ghost cells
    size_total_y = size_y + 2 * ny_domains + 2 #including ghost cells

    #Take a right time step for convergence
    if (dt1 >= dt2)
        if (my_id == 0)
            println()
            println("Time step too large ==> Taking convergence criterion")
        end
        dt = dt2
    else
        dt = dt1
    end

    #2D solution including ghost cells
    u0 = twoDArray(Float64, size_total_x, size_total_y)
    u = twoDArray(Float64, size_total_x, size_total_y)

    #Allocate coordinates of processes (start cell, end cell)
    xs = oneDArray(Int, nproc)
    xe = oneDArray(Int, nproc)
    ys = oneDArray(Int, nproc)
    ye = oneDArray(Int, nproc)

    #Size of each physical domain
    xcell = Int(size_x / nx_domains)
    ycell = Int(size_y / ny_domains)

    #allocate flattened (1D) local physical solution (i.e., relative to sub-domain)
    u_local = oneDArray(Float64, xcell * ycell)
    #allocate flattened (1D) global physical solution
    u_global = oneDArray(Float64, size_x * size_y)

    #find processes surrounding mine (if they exist)
    my_neighbors = neighbors(my_id, nproc, nx_domains, ny_domains)

    #compute coordinates of processes for each sub-domain
    process_coordinates!(xs, ys, xe, ye, xcell, ycell, nx_domains, ny_domains, nproc)

    #initialize domain
    init_values(u0, size_total_x, size_total_y, temp1_init, temp2_init)

    #update ghost cells
    updateBound!(u0, size_total_x, size_total_y, my_neighbors, comm,
                my_id, xs, ys, xe, ye, xcell, ycell, nproc)

    #Initialize step, time and convergence boolean
    step = 0
    t = 0.0
    converged = false

    #Starting time
    if (my_id==0)
        time_init = time()
    end

    #Main loop : until convergence
    while (!converged)
        #increment step and time
        step += 1
        t += dt
        #perform one step of the explicit scheme
        local_diff = computeNext!(u0, u, size_total_x, size_total_y, dt, hx, hy,
            my_id, xs, ys, xe, ye, nproc, k0)

        #update ghost cells
        updateBound!(u0, size_total_x, size_total_y, my_neighbors, comm,
            my_id, xs, ys, xe, ye, xcell, ycell, nproc)

        MPI.Barrier(comm)

        #sum local_diff to get global difference
        global_diff = MPI.Allreduce(local_diff, MPI.SUM, comm)
        global_diff = sqrt(global_diff)
        #break if convergence reached or step greater than maxStep
        if ((global_diff <= epsilon) || (step >= maxStep))
            converged = true
        end
    end

    #get ending time
    if (my_id == 0)
        time_final = time()
        #elapsed time
        elapsed_time = time_final - time_init
        println("Elapsed time = ", elapsed_time)
        println("Steps = ", step)
    end

    #find solution on my sub-domain (as a 1-dimensional array)
    i = 1
    for j = ys[my_id+1]:ye[my_id+1]
        u_local[(i-1)*xcell+1:i*xcell] = u0[xs[my_id+1]:xe[my_id+1],j]
        i = i+1
    end

    #gather local solution to global solution (as a 1-dimensional array)
    u_global = MPI.Gather(u_local, root, comm)

    #write to disk
    if (my_id == 0)
        write_to_disk(u_global, nx_domains, ny_domains, xcell, ycell, temp1_init, output_path)
    end

    MPI.Finalize()

end

main()

max_stars_repo_name: ornl-training/julia-basics

max_issues_repo_name: ornl-training/julia-basics

max_forks_repo_name: ornl-training/julia-basics


text:
module LinearSolverTests

using Gridap
using GridapPardiso
using Test
using SparseArrays
using SparseMatricesCSR

tol = 1.0e-13

#####################################################
# SparseMatrixCSC
#####################################################
#
# Matrix from Intel MKL Pardiso examples
#
#        DATA ia /1,4,6,9,12,14/
#        DATA ja
#     1  /1,2,  4,
#     2   1,2,
#     3       3,4,5,
#     4   1,  3,4,
#     5     2,    5/
#        DATA a
#     1  / 1.d0,-1.d0,      -3.d0,
#     2   -2.d0, 5.d0,
#     3                4.d0, 6.d0, 4.d0,
#     4   -4.d0,       2.d0, 7.d0,
#     5          8.d0,            -5.d0/
#####################################################
I_ = [1,1,1,2,2,3,3,3,4,4,4,5,5]
J_ = [1,2,4,1,2,3,4,5,1,3,4,2,5]
V_ = [1,-1,-3,-2,5,4,6,4,-4,2,7,8,-5]
rows = 5
cols = 5

# pardiso!
I = Vector{Int32}(); J = Vector{Int32}(); V = Vector{Float64}()
for (ik, jk, vk) in zip(I_,J_,V_)
    push_coo!(SparseMatrixCSC,I,J,V,ik,jk,vk)
end
finalize_coo!(SparseMatrixCSC,I,J,V,rows, cols)
A = sparse(I,J,V,rows,cols)
b = ones(size(A)[2])
x = similar(b)
ps = PardisoSolver(GridapPardiso.MTYPE_REAL_NON_SYMMETRIC, new_iparm(A), GridapPardiso.MSGLVL_VERBOSE)
ss = symbolic_setup(ps, A)
ns = numerical_setup(ss, A)
solve!(x, ns, b)
@test maximum(abs.(A*x-b)) < tol
test_linear_solver(ps, A, b, x)

if Int == Int64
    # pardiso_64!
    I = Vector{Int64}(); J = Vector{Int64}(); V = Vector{Float64}()
    for (ik, jk, vk) in zip(I_,J_,V_)
        push_coo!(SparseMatrixCSC,I,J,V,ik,jk,vk)
    end
    finalize_coo!(SparseMatrixCSC,I,J,V,rows, cols)
    A = sparse(I,J,V,rows, cols)
    b = ones(size(A)[2])
    x = similar(b)
    ps = PardisoSolver(GridapPardiso.MTYPE_REAL_NON_SYMMETRIC, new_iparm(A), GridapPardiso.MSGLVL_VERBOSE)
    ss = symbolic_setup(ps, A)
    ns = numerical_setup(ss, A)
    solve!(x, ns, b)
    @test maximum(abs.(A*x-b)) < tol
    test_linear_solver(ps, A, b, x)
end

#####################################################
# SparseMatrixCSR
#####################################################
#
# Matrix from Intel MKL Pardiso examples
#
#        DATA ia /1,4,6,9,12,14/
#        DATA ja
#     1  /1,2,  4,
#     2   1,2,
#     3       3,4,5,
#     4   1,  3,4,
#     5     2,    5/
#        DATA a
#     1  / 1.d0,-1.d0,      -3.d0,
#     2   -2.d0, 5.d0,
#     3                4.d0, 6.d0, 4.d0,
#     4   -4.d0,       2.d0, 7.d0,
#     5          8.d0,            -5.d0/
#####################################################
I_ = [1,1,1,2,2,3,3,3,4,4,4,5,5]
J_ = [1,2,4,1,2,3,4,5,1,3,4,2,5]
V_ = [1,-1,-3,-2,5,4,6,4,-4,2,7,8,-5]
rows = 5
cols = 5

# pardiso!

for Bi in (0,1)
    I = Vector{Int32}(); J = Vector{Int32}(); V = Vector{Float64}()
    for (ik, jk, vk) in zip(I_,J_,V_)
        push_coo!(SparseMatrixCSR,I,J,V,ik,jk,vk)
    end
    finalize_coo!(SparseMatrixCSR,I,J,V,rows,cols)
    A = sparsecsr(SparseMatrixCSR{Bi},I,J,V,rows,cols)
    b = ones(size(A)[2])
    x = similar(b)
    ps = PardisoSolver(GridapPardiso.MTYPE_REAL_NON_SYMMETRIC, new_iparm(A), GridapPardiso.MSGLVL_VERBOSE)
    ss = symbolic_setup(ps, A)
    ns = numerical_setup(ss, A)
    solve!(x, ns, b)
    @test maximum(abs.(A*x-b)) < tol
    test_linear_solver(ps, A, b, x)

    if Int == Int64
        # pardiso_64!
        I = Vector{Int64}(); J = Vector{Int64}(); V = Vector{Float64}()
        for (ik, jk, vk) in zip(I_,J_,V_)
            push_coo!(SparseMatrixCSR,I,J,V,ik,jk,vk)
        end
        finalize_coo!(SparseMatrixCSR,I,J,V,rows,cols)
        A = sparsecsr(SparseMatrixCSR{Bi},I,J,V,rows,cols)
        b = ones(size(A)[2])
        x = similar(b)
        ps = PardisoSolver(GridapPardiso.MTYPE_REAL_NON_SYMMETRIC, new_iparm(A), GridapPardiso.MSGLVL_VERBOSE)
        ss = symbolic_setup(ps, A)
        ns = numerical_setup(ss, A)
        solve!(x, ns, b)
        @test maximum(abs.(A*x-b)) < tol
        test_linear_solver(ps, A, b, x)
    end
end

#####################################################
# SymSparseMatrixCSR
#####################################################
#
# Matrix from Intel MKL Pardiso examples
#
# ia = (/ 1, 5, 8, 10, 12, 15, 17, 18, 19 /)
# ja = (/ 1,    3,       6, 7,    &
#            2, 3,    5,          &
#               3,             8, &
#                  4,       7,    &
#                     5, 6, 7,    &
#                        6,    8, &
#                           7,    &
#                              8 /)
# a = (/ 7.d0,        1.d0,             2.d0, 7.d0,        &
#              -4.d0, 8.d0,       2.d0,                    &
#                     1.d0,                         5.d0,  &
#                           7.d0,             9.d0,        &
#                                 5.d0, 1.d0, 5.d0,        &
#                                      -1.d0,       5.d0,  &
#                                            11.d0,        &
#                                                   5.d0 /)
#####################################################
I_ = [1,1,1,1,2,2,2,3,3,4,4,5,5,5,6,6,7,8]
J_ = [1,3,6,7,2,3,5,3,8,4,7,5,6,7,6,8,7,8]
V_ = [7,1,2,7,-4,8,2,1,5,7,9,5,1,5,-1,5,11,5]
rows = 8
cols = 8

# pardiso!

for Bi in (0,1)
    I = Vector{Int32}(); J = Vector{Int32}(); V = Vector{Float64}()
    for (ik, jk, vk) in zip(I_,J_,V_)
        push_coo!(SymSparseMatrixCSR,I,J,V,ik,jk,vk)
    end
    finalize_coo!(SymSparseMatrixCSR,I,J,V,rows,cols)
    A = symsparsecsr(SymSparseMatrixCSR{Bi},I,J,V,rows,cols)
    b = ones(size(A)[2])
    x = similar(b)
    ps = PardisoSolver(GridapPardiso.MTYPE_REAL_SYMMETRIC_INDEFINITE, new_iparm(A), GridapPardiso.MSGLVL_VERBOSE)
    ss = symbolic_setup(ps, A)
    ns = numerical_setup(ss, A)
    solve!(x, ns, b)
    @test maximum(abs.(A*x-b)) < tol
    test_linear_solver(ps, A, b, x)

    if Int == Int64
        # pardiso_64!
        I = Vector{Int64}(); J = Vector{Int64}(); V = Vector{Float64}()
        for (ik, jk, vk) in zip(I_,J_,V_)
            push_coo!(SymSparseMatrixCSR,I,J,V,ik,jk,vk)
        end
        finalize_coo!(SymSparseMatrixCSR,I,J,V,rows,cols)
        A = symsparsecsr(SymSparseMatrixCSR{Bi},I,J,V,rows,cols)
        b = ones(size(A)[2])
        x = similar(b)
        ps = PardisoSolver(GridapPardiso.MTYPE_REAL_SYMMETRIC_INDEFINITE, new_iparm(A), GridapPardiso.MSGLVL_VERBOSE)
        ss = symbolic_setup(ps, A)
        ns = numerical_setup(ss, A)
        solve!(x, ns, b)
        @test maximum(abs.(A*x-b)) < tol
        test_linear_solver(ps, A, b, x)
    end
end

end

ext: jl


text:
# import LinearAlgebra: mul!, *, ldiv!, \, adjoint, transpose
# import BlockArrays: BlockArray

# #
# # constructors
# #

# function BlockArray(L::LowerTriangular, Ps::Vararg{AbstractVector{Int}, 2})
#     return LowerTriangular(BlockArray(L.data, Ps...))
# end
# function BlockArray(U::UpperTriangular, Ps::Vararg{AbstractVector{Int}, 2})
#     return UpperTriangular(BlockArray(U.data, Ps...))
# end

# #
# # adjoint / transpose
# #

# adjoint(X::BlockLowerTriangular) = UpperTriangular(adjoint(X.data))
# transpose(X::BlockLowerTriangular) = UpperTriangular(transpose(X.data))

# adjoint(X::BlockUpperTriangular) = LowerTriangular(adjoint(X.data))
# transpose(X::BlockUpperTriangular) = LowerTriangular(transpose(X.data))

# #
# # mul! and *
# #

# function mul!(y::BlockVector, L::BlockLowerTriangular, x::BlockVector)
#     @assert are_conformal(L, x) && are_conformal(L', y)

#     for r in 1:nblocks(L, 1)
#         fill!(y[Block(r)], 0.0)
#         for c in 1:r-1
#             y[Block(r)] = y[Block(r)] + L[Block(r, c)] * x[Block(c)]
#         end
#         y[Block(r)] = y[Block(r)] + LowerTriangular(L[Block(r, r)]) * x[Block(r)]
#     end
#     return y
# end

# function mul!(Y::BlockMatrix, L::BlockLowerTriangular, X::BlockMatrix)
#     @assert are_conformal(L, X)
#     @assert cumulsizes(L, 1) == cumulsizes(Y, 1)
#     @assert cumulsizes(Y, 2) == cumulsizes(X, 2)

#     P, Q, R = nblocks(L, 1), nblocks(L, 2), nblocks(X, 2)
#     for p in 1:P, r in 1:R
#         fill!(Y[Block(p, r)], 0.0)
#         for q in 1:p-1
#             Y[Block(p, r)] = Y[Block(p, r)] + L[Block(p, q)] * X[Block(q, r)]
#         end
#         Y[Block(p, r)] = Y[Block(p, r)] + LowerTriangular(L[Block(p, p)]) * X[Block(p, r)]
#     end
#     return Y
# end

# function mul!(y::BlockVector, U::BlockUpperTriangular, x::BlockVector)
#     @assert are_conformal(U, x) && are_conformal(U', y)

#     for r in 1:nblocks(U, 1)
#         mul!(y[Block(r)], UpperTriangular(U[Block(r, r)]), x[Block(r)])
#         for c in r+1:nblocks(U, 2)
#             y[Block(r)] = y[Block(r)] + U[Block(r, c)] * x[Block(c)]
#         end
#     end
#     return y
# end

# function mul!(Y::BlockMatrix, U::BlockUpperTriangular, X::BlockMatrix)
#     @assert are_conformal(U, X)
#     @assert cumulsizes(U, 1) == cumulsizes(Y, 1)
#     @assert cumulsizes(Y, 2) == cumulsizes(X, 2)

#     P, Q, R = nblocks(U, 1), nblocks(U, 2), nblocks(X, 2)
#     for p in 1:P, r in 1:R
#         mul!(Y[Block(p, r)], UpperTriangular(U[Block(p, p)]), X[Block(p, r)])
#         for q in p+1:Q
#             Y[Block(p, r)] = Y[Block(p, r)] + U[Block(p, q)] * X[Block(q, r)]
#         end
#     end
#     return Y
# end

# *(A::BlockTriangular, x::BlockVector) = mul!(copy(x), A, x)
# *(A::BlockTriangular, X::BlockMatrix) = mul!(copy(X), A, X)

# @adjoint function *(A::BlockTriangular, x::BlockVector)
#     back(Δ::BlockVector) = (Δ * x', A' * Δ)
#     back(Δ::AbstractVector) = back(BlockArray(Δ, blocksizes(x)))
#     return A * x, back
# end
# @adjoint function *(A::BlockTriangular, B::BlockMatrix)
#     back(Δ::BlockMatrix) = (Δ * B', A' * Δ)
#     back(Δ::AbstractMatrix) = back(BlockArray(Δ, blocksizes(B)))
#     return A * B, back
# end

# #
# # ldiv! and \
# #

# function ldiv!(L::BlockLowerTriangular, x::BlockVector)
#     @assert are_conformal(L, x)
#     for p in 1:nblocks(x, 1)
#         for p′ in 1:p-1
#             x[Block(p)] = x[Block(p)] - L[Block(p, p′)] * x[Block(p′)]
#         end
#         ldiv!(LowerTriangular(L[Block(p, p)]), x[Block(p)])
#     end
#     return x
# end

# function ldiv!(L::BlockLowerTriangular, X::BlockMatrix)
#     @assert are_conformal(L, X)
#     for q in 1:nblocks(X, 2), p in 1:nblocks(X, 1)
#         for p′ in 1:p-1
#             X[Block(p, q)] = X[Block(p, q)] - L[Block(p, p′)] * X[Block(p′, q)]
#         end
#         ldiv!(LowerTriangular(L[Block(p, p)]), X[Block(p, q)])
#     end
#     return X
# end

# function ldiv!(U::BlockUpperTriangular, x::BlockVector)
#     @assert are_conformal(U, x )
#     for p in reverse(1:nblocks(x, 1))
#         for p′ in p+1:nblocks(x, 1)
#             x[Block(p)] = x[Block(p)] - U[Block(p, p′)] * x[Block(p′)]
#         end
#         ldiv!(UpperTriangular(U[Block(p, p)]), x[Block(p)])
#     end
#     return x
# end

# function ldiv!(U::BlockUpperTriangular, X::BlockMatrix)
#     @assert are_conformal(U, X)
#     for q in 1:nblocks(X, 2), p in reverse(1:nblocks(X, 1))
#         for p′ in p+1:nblocks(X, 1)
#             X[Block(p, q)] = X[Block(p, q)] - U[Block(p, p′)] * X[Block(p′, q)]
#         end
#         ldiv!(UpperTriangular(U[Block(p, p)]), X[Block(p, q)])
#     end
#     return X
# end

# \(L::BlockLowerTriangular, x::BlockVector) = ldiv!(L, copy(x))
# \(L::BlockLowerTriangular, X::BlockMatrix) = ldiv!(L, copy(X))
# \(U::BlockUpperTriangular, x::BlockVector) = ldiv!(U, copy(x))
# \(U::BlockUpperTriangular, X::BlockMatrix) = ldiv!(U, copy(X))


text:
using CSV: File

using DataFrames

using Dates

using Query

function query_regions(
    regions::Dict,
    vcf::DataFrame,
    output_dir::String,
    )

    start_time = now()

    for (key, value) in regions

        println(key)

        println(value)

        query = @from i in vcf begin
            @where i.CHROM==key
            @where i.POS==value
            @select {i.CHROM, i.POS, i.REF, i.ALT, i.FORMAT}
            @collect DataFrame

        end

        println("\n")

        println(query)

    end

    end_time = now()

    println("\nTook $(canonicalize(Dates.CompoundPeriod(end_time - start_time))).\n")

end

export query_regions


text:
using Test
using Plots
function test_backwardpass_sqrt!(res::SolverIterResults,solver::Solver,bp::BackwardPass)
    # Get problem sizes
    res = results2
    n,m,N = get_sizes(solver)
    m̄,mm = get_num_controls(solver)
    n̄,nn = get_num_states(solver)

    # Objective
    costfun = solver.obj.cost

    dt = solver.dt

    X = res.X; U = res.U; K = res.K; d = res.d; Su = res.S; s = res.s

    for k = 1:N-1
        res.S[k] = zeros(nn+mm,nn)
    end
    Su = res.S
    # # for now just re-instantiate
    # bp = BackwardPass(nn,mm,N)
    # Qx = bp.Qx; Qu = bp.Qu; Qxx = bp.Qxx; Quu = bp.Quu; Qux = bp.Qux
    # Quu_reg = bp.Quu_reg; Qux_reg = bp.Qux_reg

    # Boundary Conditions
    lxx,lx = taylor_expansion(costfun, X[N][1:n])

    # Initialize expected change in cost-to-go
    Δv = zeros(2)

    # Terminal constraints
    if res isa ConstrainedIterResults
        C = res.C; Iμ = res.Iμ; λ = res.λ
        Cx = res.Cx; Cu = res.Cu
        Iμ_sqrt = sqrt.(Iμ[N])

        Su[N][1:nn,1:nn] = chol_plus(cholesky(lxx).U,Iμ_sqrt*Cx[N])
        s[N] = lx + Cx[N]'*(Iμ[N]*C[N] + λ[N])

        @test isapprox(lxx + Cx[N]'*Iμ[N]*Cx[N],Su[N]'*Su[N])
    else
        Su[N] = cholesky(lxx).U
        s[N] = lx
    end

    # Backward pass
    k = N-1
    while k >= 1

        x = X[k][1:n]
        u = U[k][1:m]
        h = sqrt(dt)

        expansion = taylor_expansion(costfun,x,u)
        lxx,luu,lux,lx,_lu = expansion

        # Compute gradients of the dynamics
        fdx, fdu = res.fdx[k], res.fdu[k]

        # Gradients and Hessians of Taylor Series Expansion of Q
        Qx = dt*lx + fdx'*s[k+1]
        Qu = dt*_lu + fdu'*s[k+1]
        Wxx = chol_plus(cholesky(dt*lxx).U, Su[k+1]*fdx)
        Wuu = chol_plus(cholesky(dt*luu).U, Su[k+1]*fdu)
        Qux = dt*lux + (fdu'*Su[k+1]')*(Su[k+1]*fdx)

        @test isapprox(dt*lxx + fdx'*Su[k+1]'*Su[k+1]*fdx, Wxx'*Wxx)
        @test isapprox(dt*luu + fdu'*Su[k+1]'*Su[k+1]*fdu, Wuu'*Wuu)
        # Constraints
        if res isa ConstrainedIterResults
            Iμ_sqrt = sqrt.(Iμ[k])

            Qx += Cx[k]'*(Iμ[k]*C[k] + λ[k])
            Qu += Cu[k]'*(Iμ[k]*C[k] + λ[k])
            Wxx = chol_plus(Wxx,Iμ_sqrt*Cx[k])
            Wuu = chol_plus(Wuu,Iμ_sqrt*Cu[k])
            Qux += Cu[k]'*Iμ[k]*Cx[k]

            @test isapprox(dt*lxx + fdx'*Su[k+1]'*Su[k+1]*fdx + Cx[k]'*Iμ[k]*Cx[k], Wxx'*Wxx)
            @test isapprox(dt*luu + fdu'*Su[k+1]'*Su[k+1]*fdu + Cu[k]'*Iμ[k]*Cu[k], Wuu'*Wuu)
        end
        #
        if solver.opts.bp_reg_type == :state
            Wuu_reg = chol_plus(Wuu,sqrt(res.ρ[1])*I*fdu)
            Qux_reg = Qux + res.ρ[1]*fdu'*fdx
        elseif solver.opts.bp_reg_type == :control
            Wuu_reg = chol_plus(Wuu,sqrt(res.ρ[1])*Matrix(I,m,m))
            Qux_reg = Qux
        end

        #TODO find better PD check for Wuu_reg
        # # Regularization
        # if !isposdef(Hermitian(Array(Wuu_reg)))  # need to wrap Array since isposdef doesn't work for static arrays
        #     # increase regularization
        #     regularization_update!(res,solver,:increase)
        #
        #     # reset backward pass
        #     k = N-1
        #     Δv[1] = 0.
        #     Δv[2] = 0.
        #     continue
        # end

        # Compute gains
        K[k] = -Wuu_reg\(Wuu_reg'\Qux_reg)
        d[k] = -Wuu_reg\(Wuu_reg'\Qu)

        # Calculate cost-to-go
        s[k] = Qx + (K[k]'*Wuu')*(Wuu*d[k]) + K[k]'*Qu + Qux'*d[k]

        tmp1 = (Wxx')\Qux'
        tmp2 = cholesky(Wuu'*Wuu - tmp1'*tmp1).U
        Su[k][1:nn,1:nn] = Wxx + tmp1*K[k]
        Su[k][nn+1:nn+mm,1:nn] = tmp2*K[k]

        # calculated change is cost-to-go over entire trajectory
        Δv[1] += d[k]'*Qu
        Δv[2] += 0.5*d[k]'*Wuu'*Wuu*d[k]

        k = k - 1;
    end

    # decrease regularization after backward pass
    regularization_update!(res,solver,:decrease)

    return Δv
end

function chol_plus(A,B)
    n1,m = size(A)
    n2 = size(B,1)
    P = zeros(n1+n2,m)
    P[1:n1,:] = A
    P[n1+1:end,:] = B
    return qr(P).R
end

# model_dubins, obj_uncon_dubins = TrajectoryOptimization.Dynamics.dubinscar
#
# # dubins car
# u_min_dubins = [-1; -1]
# u_max_dubins = [1; 1]
# x_min_dubins = [0; -100; -100]
# x_max_dubins = [1.0; 100; 100]
# obj_con_dubins = ConstrainedObjective(obj_uncon_dubins, u_min=u_min_dubins, u_max=u_max_dubins, x_min=x_min_dubins, x_max=x_max_dubins)
#
# # -Constrained objective
#
# model = model_dubins
# obj = obj_con_dubins
# u_max = u_max_dubins
# u_min = u_min_dubins
# Solver options

# Set up quadrotor model, objective, solver
intergrator = :rk4
dt = 0.005
N = 201
tf = 10.0
r_quad = 3.0
model, = TrajectoryOptimization.Dynamics.quadrotor
n = model.n
m = model.m

# -initial state
x0 = zeros(n)
x0[1:3] = [0.; 0.; 0.]
q0 = [1.;0.;0.;0.]
x0[4:7] = q0

# -final state
xf = copy(x0)
xf[1:3] = [0.;40.;0.] # xyz position
xf[4:7] = q0

# -control limits
u_min = 0.0
u_max = 5.0

Q = (1e-1)*Matrix(I,n,n)
Q[4,4] = 1.0
Q[5,5] = 1.0
Q[6,6] = 1.0
Q[7,7] = 1.0
R = (1.0)*Matrix(I,m,m)
Qf = (1000.0)*Matrix(I,n,n)
# obstacles constraint
r_sphere = 3.0
spheres = ((0.,10.,0.,r_sphere),(0.,20.,0.,r_sphere),(0.,30.,0.,r_sphere))
n_spheres = 3

function cI(c,x,u)
    for i = 1:n_spheres
        c[i] = sphere_constraint(x,spheres[i][1],spheres[i][2],spheres[i][3],spheres[i][4]+r_quad)
    end
    c
end

# unit quaternion constraint
function cE(c,x,u)
    c = sqrt(x[4]^2 + x[5]^2 + x[6]^2 + x[7]^2) - 1.0
end

obj_uncon = LQRObjective(Q, R, Qf, tf, x0, xf)
obj = TrajectoryOptimization.ConstrainedObjective(obj_uncon,u_min=u_min,u_max=u_max,cI=cI,cE=cE)

# Solver
solver = Solver(model,obj,integration=intergrator,N=N)
solver.opts.bp_reg_initial = 1.0
U0 = ones(solver.model.m,solver.N)

results1 = init_results(solver,Array{Float64}(undef,0,0),U0)
results2 = init_results(solver,Array{Float64}(undef,0,0),U0)
rollout!(results1,solver)
rollout!(results2,solver)
update_jacobians!(results1, solver)
update_jacobians!(results2, solver)

n,m,N = get_sizes(solver)
m̄,mm = get_num_controls(solver)
n̄,nn = get_num_states(solver)
bp = BackwardPass(nn,mm,solver.N)

v1 = _backwardpass!(results1,solver,bp)
v2 = test_backwardpass_sqrt!(results2,solver,bp)

@test isapprox(v1[1:2], v2[1:2])
@test isapprox(to_array(results1.K),to_array(results2.K))
@test isapprox(to_array(results1.d),to_array(results2.d))

S_sqrt = [zeros(nn,nn) for k = 1:N]
cond_normal = zeros(N)
cond_sqrt = zeros(N)

for k = 1:N
    S_sqrt[k] = results2.S[k]'*results2.S[k]
    cond_normal[k] = cond(results1.S[k])
    cond_sqrt[k] = cond(results2.S[k])
end

@test isapprox(to_array(results1.S),to_array(S_sqrt))

plot(cond_normal,label="bp")
plot(cond_sqrt,title="Condition number",label="sqrt bp")


text:
## -------------  Model

@with_kw mutable struct Model
    # Number of parameters
    K :: Int64 = 12
    # Number of households
    N :: Int64 = 10
    # Number of options
    M :: Int64 = 20
    # How often to simulate each type
    nSim :: Int64 = 200

    # Governs how deviation function shows intermediate results
    showFreq :: Int64 = 100
    currentIter :: Int64 = 0
    # No of classes when computing stats
    nEndowCl :: Int = 6
    nChoiceCl :: Int = 8
end

# No of calibrated params
function n_params(mdl :: Model)
    return mdl.K
end

## ---------  Model parameters

mutable struct ModelParams
    alphaV :: Vector{Float64}
    # betaV :: Vector{Float64}
    # gammaV :: Vector{Float64}
    # deltaV :: Vector{Float64}
end

function show_params(mp :: ModelParams)
    println("Alpha:  ");
    println(round.(mp.alphaV, digits = 3));
    # print("Beta:   "); println(round.(mp.betaV, digits = 3));
    # print("Gamma:  "); println(round.(mp.gammaV, digits = 3));
    # print("Delta:  "); println(round.(mp.deltaV, digits = 3));
end

## --------  Model stats

mutable struct ModelStats
    # Fraction of choices in each [choice class, endowment class, by endowment] bin
    frac_ceM  ::  Array{Float64, 2}
end

function n_tg_moments(ms :: ModelStats)
    return length(ms.frac_ceM)
end

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


text:
module Logdet
    using LinearAlgebra
    using Distributions
    using SparseArrays
    export trlog

    function trlog(A)
        itemax = 80
        num = 80
        return trlog(A,itemax,num)
    end

    function calclog(T)
        e,v = eigen(T)
        valuei = 0.0
        kmax = size(T)[1]
        for l=1:kmax
            valuei +=  log(e[l])*abs(v[1,l])^2
        end
        return valuei
    end

    function trlog(A,itemax,num)
        eps = 1e-4
        return trlog(A,itemax,num,eps)
    end

    function trlog(A,itemax,num,eps)
        println("Estimation of log(det(A))=Tr(log(A)) start------------------------------------------")
        println("Num. of Lanczos steps: $itemax")
        println("Num. of random vectors: $num")

        n = size(A)[1]
        value = 0.0
        oldvalue = 0.0
        eps2= 1e-6
        g = zeros(Float64,2,2)
        y = zeros(Float64,2)
        #value1 = 0
        #println("trace(A) = ",tr(A))
        #trA = tr(A)
        #trvalue = 0
        #trvalue1 = 0
        d = Normal(0,1/sqrt(2))
        #d = Normal(0,1)
        #ene,ve = eigen(A)
        #println(ene)
        #logdet = sum(log.(ene))
#        logdet = tr(log(A))
        #meanvalue = 0.0
        icount = 0
        #fp = open("hi2.dat","w")
        for i=1:num
            #phii = (rand(n)+im*rand(n))
            #phii = (rand(n) .- 0.5).*2 .+ im*(rand(n) .- 0.5).*2
            phii = rand(d,n)+im*rand(d,n)
            #nphi = norm(phii)
            #dotd = dot(phii,phii)
            #println("dot0 ", dotd)
            #phii = phii ./norm(phii)
            #println("dot ", dot(phii,phii))
            #phii = zeros(Float64,n)
            #phii[i] = 1
            #phii = rand([-1.0,1.0],n)
#            phii = rand([-1.0,1.0,im,-im],n)
            nphi = norm(phii)
            phii = phii/norm(phii)
            valuei = Lanczos(A,phii,itemax,eps)
            #T = Lanczos(A,phii,itemax)
            #e,v = eigen(T)

            #valuei = 0
            #valuei1 = 0
            #println(e)
            #exit()
            #for l=1:itemax
                #println(e[l],"\t",abs(v[1,l])^2)
            #    valuei +=  log(e[l])*abs(v[1,l])^2
            #    valuei1 +=  e[l]*abs(v[1,l])^2
            #end
            #println(valuei1,"\t",dot(phii,A*phii))
            #exit()
            #trvalue += dot(phii,A*phii)
            #trvalue1 += dot(phii,phii)
            value += valuei*nphi^2
            #=
            if i > itemax-20
                icount += 1
                meanvalue += value/i
            end
            =#
            if i % div(num,10) == 0
                println("$i-th ", value/i)
            end
            #=
            g[1,1] += 1
            g[1,2] += 1/i
            g[2,1] += 1/i
            g[2,2] += 1/i^2
            y[1] += value/i
            y[2] += (1/i)*value/i
            if i > 20
                ceff = g \ y
                #println(ceff)
                #println(ceff[1]+ceff[2]/i)
            end
            =#

            if i >1
                hi = abs(oldvalue-value/i)/abs(oldvalue)
                if abs(oldvalue-value/i) < eps2
                    println("In the $i-th step, this is converged. Absolute error might be ",abs(oldvalue-value/i))
                    println("Estimation of log(det(A))=Tr(log(A)) end------------------------------------------")
                    return value/i
                end
                ##println(i,"\t",hi,"\t",value/i,"\t",oldvalue-value/i)#,"\t",meanvalue/icount)
                #if i > 20
                    #println(fp,i,"\t",hi,"\t",value/i,"\t",oldvalue-value/i)#,"\t",ceff[1]+ceff[2]/i)
                #end
            end

            oldvalue = value/i
            #value1 += valuei1*nphi^2
            #println("$i-th value = ",value,"\t",logdet,"\t",abs(value-logdet)/abs(logdet))
            #println("$i-th trvalue = ",trvalue,"\t",trA,"\t",abs(trvalue-trA)/abs(trA),"\t",value1)
            #println("$i-th value = ",value/i,"\t",logdet,"\t",abs(value/i-logdet)/abs(logdet))
            #println("$i-th trvalue = ",trvalue/i,"\t",trA,"\t",abs(trvalue/i-trA)/abs(trA),"\t",value1/i)
            #println("$i-th trvalue1 = ",trvalue1/i)
            #println("eigenvalues = ",real.(e))
        end
        #close(fp)
        println("Estimation of log(det(A))=Tr(log(A)) end------------------------------------------")
        #exit()

        return value/num

    end

    function Lanczos(A,v0,itemax)
        eps = 1e-8
        return Lanczos(A,v0,itemax,eps)
    end

    function Lanczos(A,v0,itemax,eps)
        ipt = Int64[1,2]
        n = length(v0)
        ftype = eltype(v0)
        β = zeros(Float64,itemax+1)
        α = zeros(Float64,itemax+1)
        w = zeros(ftype,n,2)
        x = similar(v0)
        w[:,2] = v0[:]/norm(v0)
        T = zeros(ftype,itemax+1,itemax+1)
        trlog = 0.0
        e = 0
        V = 0
        oldvaluei = 0.0
        valuei = 0.0
        hi = 0.0

        @inbounds for k in 1:itemax
            #vj = w[:,ipt[2]]
            #vj-1 = w[:,ipt[1]]
            mul!(x,A,w[:,ipt[2]])
            #x = A*w[:,ipt[2]] #Avj = A*v
            #mul!(w[:,ipt[1]],A,w[:,ipt[2]],1,- β[k]) #mul!(C,A,B,alpha,beta) -> C = A*B*alpha + C*beta
            @. w[:,ipt[1]] = x[:] - β[k]*w[:,ipt[1]] #wj = Avj - betaj*vj-1
            α[k] = real(dot(w[:,ipt[1]],w[:,ipt[2]])) # alphaj = (wj,vj)
            @. w[:,ipt[1]] = w[:,ipt[1]]-α[k]*w[:,ipt[2]] #wj = wj - alphaj*vj
            β[k+1] = norm(w[:,ipt[1]]) #betaj+1 = norm(w[:,ipt[1]])
            @. w[:,ipt[1]] /= β[k+1] #vj = wj/betaj+1
            ii = ipt[1]
            ipt[1] = ipt[2]
            ipt[2] = ii
            T[k,k] = α[k]
            T[k,k+1] =  β[k+1]
            T[k+1,k] =  β[k+1]
            valuei = calclog(T[1:k,1:k])
            #println("$k-th value: ",valuei)

            if k>1
                hi = abs(oldvaluei-valuei)/abs(valuei)
                #println("$k-th ",hi)
                if hi < eps
                    return  valuei
                end
            end
            oldvaluei = valuei

            #println("$k-th alpha: ",abs(α[k]-α[k-1])/abs(α[k]))
            #println("$k-th beta: ",abs(β[k+1]-β[k])/abs(β[k+1]))
            #println("$k-th beta: ",β[k+1])

            #ista = ifelse(k+1-4 > 0,k+1-4,k+1)
            #trlog = 0.0
            #for i=1:k+1
            #    trlog += log(e[i])*abs(V[i,1])^2
            #end
            #println(e[ista:end])
            #println(trlog)
            #e,V =  eigen(T[1:itemax,1:itemax])
            #println(e)
        end

        #T =  eigen(T[1:itemax,1:itemax])
        #println(e)
        #println(V*v0)
        #exit()
        println("Warning!: Lanczos in logdet.jl is not converged. error: ",hi)
        return valuei

        #return T[1:itemax,1:itemax]
    end
end

#=

using LinearAlgebra
using .Logdet

using Random
using Distributions
using SparseArrays
Random.seed!(2234)
d = Normal()

n = 2000
#A = sprand(Float64,n,n,3/n)
A = zeros(Float64,n,n)
t = -1
for i=1:n
    j = i+1
    j += ifelse(j > n,-n,0)
    A[i,j] = -t

    j = i-1
    j += ifelse(j <1,n,0)
    A[i,j] = -t
    A[i,i] = rand()
end

#A = rand(d,n,n) ./n

AdagA = A'*A  + Diagonal(ones(n))

#AdagA = Matrix(AdagA)
e,v = eigen(AdagA)
println(e[end-4:end])

itemax = 50*2
v0 = rand(n)
#Logdet.Lanczos(AdagA,v0,itemax)

trlog = Logdet.trlog(AdagA,itemax,40)

#logdet = log(det(AdagA))
trlog_exact = tr(log(AdagA))
println("trlog_exact AdagA = ", trlog_exact)
println("trlog AdagA = ", trlog)

=#


text:
# ====================================================================
# Kinnosuke.jl : A lightweight toolbox for scientific machine learning
# Tianbai Xiao (c) 2020
# ====================================================================

module Kinnosuke

using OrdinaryDiffEq, Flux, DiffEqFlux, Optim
using FileIO, JLD2, Plots

include("diff.jl")
include("widget.jl")
include("train.jl")

export central_diff, central_diff!, upwind_diff, upwind_diff!
export track, device
export vis_train

end # module


text:
# This file is part of Kpax3. License is MIT.

# the following code has been modified from
# http://thirld.com/blog/2015/05/30/julia-profiling-cheat-sheet/

#=
# do not forget to run julia with --track-allocation=user
# to view results
using ProfileView
f = open("build/profile.bin", "r")
r = deserialize(f);
ProfileView.view(r[1], lidict=r[2])
=#
function benchmark()
  include("src/boot.jl")

  settings = KSettings("test/data/mcmc_6.fasta", "test/data/mcmc_6",
                       α=0.5, θ=-0.25, γ=[0.6; 0.35; 0.05],
                       r=log(0.001)/log(0.95), maxclust=1, maxunit=1,
                       verbose=false, verbosestep=100000, popsize=50,
                       xrate=0.9, mrate=0.05, T=1000000, burnin=500000,
                       tstep=1, op=[0.6; 0.3; 0.1], λs1=1.0, λs2=1.0,
                       parawm=5.0)

  x = AminoAcidData(settings)

  # Run once, to force compilation.
  println("-- First run --")
  Random.seed!(20150326)
  @time kpax3mcmc(x, partition, settings)

  # Run a second time, with profiling.
  println("-- Second run --")
  Random.seed!(20150326)
  Profile.init(delay=0.01)
  Profile.clear()
  Profile.clear_malloc_data()
  @profile @time kpax3mcmc(x, partition, settings)

  # Write profile results
  r = Profile.retrieve()
  f = open("build/profile.bin", "w")
  serialize(f, r)
  close(f)

  nothing
end


text:
"""
    markov_transition(μ::Function, Σ::Function, minp::AbstractFloat, statevectors::AbstractVector...)

Returns the sparse, Markov Transition Matrix (that is, `P[i,j] = Pr(i|j)`) where all
elements are `> minp`. The state-space is defined as the Cartesian Product of the `statevectors`.
The functions `μ` and `Σ` should take a tuple from `Base.product(statevectors...)` and return EITHER
  - the mean and **STANDARD DEVIATION** of `Distributions.Normal`, or
  - the mean and **VARIANCE MATRIX** of `Distributions.MvNormal`
"""
function markov_transition(μ::Function, Σ::Function, minp::AbstractFloat, statevectors::AbstractVector{T}...) where {T<:AbstractFloat}

  0. <= minp < 1. || throw(DomainError())

  state_prod = Base.product(statevectors...)

  P = zeros(T, length(state_prod), length(state_prod))

  for (j, s2) in enumerate(state_prod)
    for (i, s1) in enumerate(state_prod)
      dist = myDist(μ(s1), Σ(s1))
      P[i,j] = mypdf(dist, s2)
    end
  end

  return sparsify!(P, minp)
end

"""
    markovswitching_transition(μ::Function, Σ::Function, π::Matrix{Float64}, minp::AbstractFloat, statevectors::AbstractVector...)

Returns a sparse, Markov Transition Matrix (that is, `P[i,j] = Pr(i|j)`) where all
elements are `> minp`. The state space is defined as the Cartesian Product of the
`statevectors` with `1:size(π,1)`. The functions `μ(::Int, ::Tuple)` and `Σ(::Int, ::Tuple)`
should take the index of a regime `1:size(π,1)` and a tuple from `Base.product(statevectors...)`
and return EITHER
  - the mean and **STANDARD DEVIATION** of `Distributions.Normal`, or
  - the mean and **VARIANCE MATRIX** of `Distributions.MvNormal`
The markov-switching matrix is NOT transposed and equals `π[i,j] = Pr(j|i)`

"""
function markovswitching_transition(μ::Function, Σ::Function, π::Matrix{Float64}, minp::AbstractFloat, statevectors::AbstractVector{T}...) where {T<:AbstractFloat}

  k = size(π,1)
  k == size(π,2)   || throw(DimensionMismatch())
  all(sum(π, dims=2) .≈ 1.0) || throw(error("each row of π must sum to 1"))

  n = prod(map(length, statevectors))
  regimes = Base.OneTo(k)
  P = zeros(T, k*n, k*n)

  for r1 in regimes
    for r2 in regimes
      P[(r1-1)*n+1:r1*n,  (r2-1)*n+1:r2*n] .= π[r1, r2] .* markov_transition( (s) -> μ(r1, s), (s) -> Σ(r1, s), minp, statevectors... )
    end
  end

  return whichP(P, Val{minp > 0.0})
end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _sunpos(jd::AbstractFloat, radians::Bool)
    # Number of Julian centuries since 1899-12-31T12:00:00
    t = (jd - 2415020) / JULIANCENTURY
    # Sun's mean longitude
    l = (279.696678 + mod(36000.768925*t, 360)) * 3600
    # Allow for ellipticity of the orbit (equation of centre) using the Earth's
    # mean anomaly ME
    me = 358.475844 + mod(35999.049750*t, 360)
    ellcor = (6910.1 - 17.2 * t) * sind(me) + 72.3 * sind(2 * me)
    l += ellcor
    # Allow for the Venus perturbations using the mean anomaly of Venus MV
    mv = 212.603219 + mod(58517.803875*t, 360)
    vencorr = 4.8*cosd(299.1017 + mv - me) +
        5.5*cosd(148.3133 +  2 * mv - 2 * me) +
        2.5*cosd(315.9433 +  2 * mv - 3 * me) +
        1.6*cosd(345.2533 +  3 * mv - 4 * me) +
        1.0*cosd(318.15   +  3 * mv - 5 * me)
    l += vencorr
    # Allow for the Mars perturbations using the mean anomaly of Mars MM.
    mm = 319.529425  +  mod(19139.858500*t, 360)
    marscorr = 2 * cosd(343.8883 - 2 * mm + 2 * me) +
        1.8 * cosd(200.4017 - 2 * mm + me)
    l += marscorr
    # Allow for the Jupiter perturbations using the mean anomaly of Jupiter MJ
    mj = 225.328328 + mod(3034.6920239 * t, 360)
    jupcorr = 7.2*cosd(179.5317 - mj + me) +
        2.6*cosd(263.2167 - mj) +
        2.7*cosd( 87.1450 - 2 * mj + 2 * me) +
        1.6*cosd(109.4933 - 2 * mj + me)
    l += jupcorr
    # Allow for the Moons perturbations using the mean elongation of the Moon
    # from the Sun D
    d = 350.7376814 + mod(445267.11422 * t, 360)
    mooncorr  = 6.5*sind(d)
    l += mooncorr
    # Allow for long period terms
    longterm = 6.4*sind(231.19 + 20.20*t)
    l += longterm
    l  = mod(l + 2592000, 1296000)
    longmed = sec2rad(l)
    # Allow for Aberration
    l -= 20.5
    # Allow for Nutation using the longitude of the Moons mean node OMEGA
    sin_omega, cos_omega = sincos(deg2rad(259.183275 - mod(1934.142008 * t, 360)))
    l -= 17.2 * sin_omega
    # True Obliquity
    oblt = 23.452294 - 0.0130125 * t + 9.2 * cos_omega / 3600
    # Right Ascension and Declination
    l /= 3600
    oblt = deg2rad(oblt)
    sin_oblt, cos_oblt = sincos(oblt)
    sin_l, cos_l = sincos(deg2rad(l))
    ra = mod2pi(atan(sin_l * cos_oblt, cos_l))
    dec = asin(sin_l * sin_oblt)
    if radians
        return ra, dec, longmed, oblt
    else
        return rad2deg(ra), rad2deg(dec), rad2deg(longmed), rad2deg(oblt)
    end
end

"""
    sunpos(jd[, radians=true]) -> ra, dec, elong, obliquity

### Purpose ###

Compute the right ascension and declination of the Sun at a given date.

### Arguments ###

* `jd`: the Julian date of when you want to calculate Sun position.  It can be
  either a scalar or a vector.  Use `jdcnv` to get the Julian date for a given
  date and time.
* `radians` (optional boolean keyword): if set to `true`, all output quantities
  are given in radians.  The default is `false`, so all quantities are given in
  degrees.

### Output ###

The 4-tuple `(ra, dec, elong, obliquity)`:

* `ra`: the right ascension of the Sun at that date
* `dec`: the declination of the Sun at that date
* `elong`: ecliptic longitude of the Sun at that date
* `obliquity`: the obliquity of the ecliptic

All quantities are given in degrees, unless `radians` keyword is set to `true`
(see "Arguments" section).  If `jd` is an array, arrays of the same given as
`jd` are returned.

### Method ###

Uses a truncated version of Newcomb's Sun.  Adapted from the IDL routine SUN_POS
by CD Pike, which was adapted from a FORTRAN routine by B. Emerson (RGO).

### Example ###

(1) Find the apparent right ascension and declination of the Sun on May 1, 1982

```jldoctest
julia> using AstroLib

julia> adstring(sunpos(jdcnv(1982, 5, 1))[1:2], precision=2)
" 02 31 32.614  +14 54 34.92"
```

The Astronomical Almanac gives `02 31 32.58 +14 54 34.9` so the error for this
case is < 0.5".

(2) Plot the apparent right ascension, in hours, and declination of the Sun, in
degrees, for every day in 2016.  Use
[PyPlot.jl](https://github.com/JuliaPlots/Plots.jl/) for plotting.

```julia
using PyPlot
using Dates

days = DateTime(2016):Day(1):DateTime(2016, 12, 31);
ra, declin = sunpos(jdcnv.(days));
plot(days, ra/15); plot(days, declin)
```

### Notes ###

Patrick Wallace (Rutherford Appleton Laboratory, UK) has tested the accuracy of
a C adaptation of the present algorithm and found the following results.  From
1900-2100 `sunpos` gave 7.3 arcsec maximum error, 2.6 arcsec RMS.  Over the
shorter interval 1950-2050 the figures were 6.4 arcsec max, 2.2 arcsec RMS.

The returned `ra` and `dec` are in the given date's equinox.

Code of this function is based on IDL Astronomy User's Library.
"""
sunpos(jd::Real; radians::Bool=false) = _sunpos(float(jd), radians)

function sunpos(jd::AbstractArray{J}; radians::Bool=false) where {J<:Real}
    typej = float(J)
    ra = similar(jd, typej)
    dec = similar(jd, typej)
    longmed = similar(jd, typej)
    oblt = similar(jd, typej)
    for i in eachindex(jd)
        ra[i], dec[i], longmed[i], oblt[i] = sunpos(jd[i], radians=radians)
    end
    return ra, dec, longmed, oblt
end


text:
module TestConversions
    using Base.Test
    using DataTables
    using DataStructures: OrderedDict, SortedDict

    dt = DataTable()
    dt[:A] = 1:5
    dt[:B] = [:A, :B, :C, :D, :E]
    @test isa(convert(Array, dt), Matrix{Any})
    @test convert(Array, dt) == convert(Array, convert(NullableArray, dt))
    @test isa(convert(Array{Any}, dt), Matrix{Any})

    dt = DataTable()
    dt[:A] = 1:5
    dt[:B] = 1.0:5.0
    # Fails on Julia 0.4 since promote_type(Nullable{Int}, Nullable{Float64}) gives Nullable{T}
    if VERSION >= v"0.5.0-dev"
        @test isa(convert(Array, dt), Matrix{Float64})
    end
    @test convert(Array, dt) == convert(Array, convert(NullableArray, dt))
    @test isa(convert(Array{Any}, dt), Matrix{Any})
    @test isa(convert(Array{Float64}, dt), Matrix{Float64})

    dt = DataTable()
    dt[:A] = NullableArray(1.0:5.0)
    dt[:B] = NullableArray(1.0:5.0)
    a = convert(Array, dt)
    aa = convert(Array{Any}, dt)
    ai = convert(Array{Int}, dt)
    @test isa(a, Matrix{Float64})
    @test a == convert(Array, convert(NullableArray, dt))
    @test a == convert(Matrix, dt)
    @test isa(aa, Matrix{Any})
    @test aa == convert(Matrix{Any}, dt)
    @test isa(ai, Matrix{Int})
    @test ai == convert(Matrix{Int}, dt)

    dt[1,1] = Nullable()
    @test_throws ErrorException convert(Array, dt)
    na = convert(NullableArray, dt)
    naa = convert(NullableArray{Any}, dt)
    nai = convert(NullableArray{Int}, dt)
    @test isa(na, NullableMatrix{Float64})
    @test isequal(na, convert(NullableMatrix, dt))
    @test isa(naa, NullableMatrix{Any})
    @test isequal(naa, convert(NullableMatrix{Any}, dt))
    @test isa(nai, NullableMatrix{Int})
    @test isequal(nai, convert(NullableMatrix{Int}, dt))

    a = NullableArray([1.0,2.0])
    b = NullableArray([-0.1,3])
    c = NullableArray([-3.1,7])
    di = Dict("a"=>a, "b"=>b, "c"=>c)

    dt = convert(DataTable,di)
    @test isa(dt,DataTable)
    @test names(dt) == Symbol[x for x in sort(collect(keys(di)))]
    @test isequal(dt[:a], NullableArray(a))
    @test isequal(dt[:b], NullableArray(b))
    @test isequal(dt[:c], NullableArray(c))

    od = OrderedDict("c"=>c, "a"=>a, "b"=>b)
    dt = convert(DataTable,od)
    @test isa(dt, DataTable)
    @test names(dt) == Symbol[x for x in keys(od)]
    @test isequal(dt[:a], NullableArray(a))
    @test isequal(dt[:b], NullableArray(b))
    @test isequal(dt[:c], NullableArray(c))

    sd = SortedDict("c"=>c, "a"=>a, "b"=>b)
    dt = convert(DataTable,sd)
    @test isa(dt, DataTable)
    @test names(dt) == Symbol[x for x in keys(sd)]
    @test isequal(dt[:a], NullableArray(a))
    @test isequal(dt[:b], NullableArray(b))
    @test isequal(dt[:c], NullableArray(c))

    a = [1.0]
    di = Dict("a"=>a, "b"=>b, "c"=>c)
    @test_throws DimensionMismatch convert(DataTable,di)

end


text:
IsingMeasurement = MonteCarlo.IsingMeasurement
AbstractObservable = MonteCarloObservable.AbstractObservable

struct DummyMeasurement <: AbstractMeasurement end

@testset "Type Hierarchy" begin
    @test IsingEnergyMeasurement <: IsingMeasurement
    @test IsingMagnetizationMeasurement <: IsingMeasurement
    @test IsingMeasurement <: AbstractMeasurement
    @test MonteCarlo.DQMCMeasurement <: AbstractMeasurement
end

@testset "Interface" begin
    m = DummyMeasurement()
    dummy_model = DummyModel()
    model = IsingModel(dims=2, L=2)
    mc = MC(model, beta=1.0)

    @test nothing == MonteCarlo.prepare!(m, mc, model)
    @test_throws MethodError MonteCarlo.measure!(m, mc, model, 1)
    @test nothing == MonteCarlo.finish!(m, mc, model)
    @test MonteCarlo.default_measurements(mc, dummy_model) == Dict{Symbol, AbstractMeasurement}()
end

@testset "Checking defaults" begin
    model = IsingModel(dims=2, L=2)
    mc = MC(model, beta=1.0)

    defaults = MonteCarlo.default_measurements(mc, model)
    @test !isempty(defaults)

    @test haskey(defaults, :Magn) && defaults[:Magn] isa IsingMagnetizationMeasurement
    obs = observables(defaults[:Magn])
    @test haskey(obs, "Magnetization per site") && obs["Magnetization per site"] isa AbstractObservable
    @test haskey(obs, "Magnetic susceptibility") && obs["Magnetic susceptibility"] isa AbstractObservable
    @test haskey(obs, "Total magnetization") && obs["Total magnetization"] isa AbstractObservable
    @test haskey(obs, "Total magnetization squared") && obs["Total magnetization squared"] isa AbstractObservable

    @test haskey(defaults, :Energy) && defaults[:Energy] isa IsingEnergyMeasurement
    obs = observables(defaults[:Energy])
    @test haskey(obs, "Energy per site") && obs["Energy per site"] isa AbstractObservable
    @test haskey(obs, "Specific heat") && obs["Specific heat"] isa AbstractObservable
    @test haskey(obs, "Total energy") && obs["Total energy"] isa AbstractObservable
    @test haskey(obs, "Total energy squared") && obs["Total energy squared"] isa AbstractObservable

    @test isempty(mc.thermalization_measurements)
    @test !isempty(mc.measurements)
    @test haskey(mc.measurements, :Magn) && mc.measurements[:Magn] isa IsingMagnetizationMeasurement
    @test haskey(mc.measurements, :Energy) && mc.measurements[:Energy] isa IsingEnergyMeasurement

    model = HubbardModelAttractive(2, 2)
    mc = DQMC(model, beta=1.0)

    @test isempty(mc.measurements)
    @test isempty(mc.thermalization_measurements)
end

@testset "Retrieving Measurements" begin
    model = IsingModel(dims=2, L=2)
    mc = MC(model, beta=1.0)

    @test mc.thermalization_measurements == measurements(mc, :TH)
    @test mc.measurements == measurements(mc)

    obs = observables(mc, :all)
    @test keys(obs[:TH]) == keys(measurements(mc, :TH))
    @test keys(obs[:ME]) == keys(measurements(mc))

    @test haskey(obs[:ME][:Energy], "Total energy")
    @test typeof(obs[:ME][:Energy]["Total energy"]) <: AbstractObservable
    @test haskey(obs[:ME][:Energy], "Total energy squared")
    @test typeof(obs[:ME][:Energy]["Total energy squared"]) <: AbstractObservable
    @test haskey(obs[:ME][:Energy], "Energy per site")
    @test typeof(obs[:ME][:Energy]["Energy per site"]) <: AbstractObservable
    @test haskey(obs[:ME][:Energy], "Specific heat")
    @test typeof(obs[:ME][:Energy]["Specific heat"]) <: AbstractObservable

    @test haskey(obs[:ME][:Magn], "Total magnetization")
    @test typeof(obs[:ME][:Magn]["Total magnetization"]) <: AbstractObservable
    @test haskey(obs[:ME][:Magn], "Total magnetization squared")
    @test typeof(obs[:ME][:Magn]["Total magnetization squared"]) <: AbstractObservable
    @test haskey(obs[:ME][:Magn], "Magnetization per site")
    @test typeof(obs[:ME][:Magn]["Magnetization per site"]) <: AbstractObservable
    @test haskey(obs[:ME][:Magn], "Magnetic susceptibility")
    @test typeof(obs[:ME][:Magn]["Magnetic susceptibility"]) <: AbstractObservable
end

@testset "Interacting with Measurements" begin
    model = IsingModel(dims=2, L=2)
    mc = MC(model, beta=1.0)

    @test length(mc.measurements) == 2
    delete!(mc, AbstractMeasurement)
    @test isempty(mc.measurements)

    @test_throws ErrorException push!(mc, :E => Int64)
    @test_throws MethodError MonteCarlo.unsafe_push!(mc, :E => 1.0)
    @test_throws ErrorException push!(mc, :E => IsingEnergyMeasurement, :bad_stage)
    @test_throws ErrorException MonteCarlo.unsafe_push!(mc, :E => IsingEnergyMeasurement(mc, model), :bad_stage)
    @test_throws ErrorException delete!(mc, :E, :bad_stage)
    @test_throws ErrorException delete!(mc, IsingEnergyMeasurement, :bad_stage)

    push!(mc, :E => IsingEnergyMeasurement)
    @test haskey(mc.measurements, :E) && mc.measurements[:E] isa IsingEnergyMeasurement
    delete!(mc, :E)
    @test !haskey(mc.measurements, :E)

    MonteCarlo.unsafe_push!(mc, :Magn => IsingMagnetizationMeasurement(mc, model))
    @test haskey(mc.measurements, :Magn) && mc.measurements[:Magn] isa IsingMagnetizationMeasurement
    delete!(mc, IsingMagnetizationMeasurement)
    @test !haskey(mc.measurements, :Magn)

    @test isempty(mc.thermalization_measurements)
    push!(mc, :E => IsingEnergyMeasurement, :TH)
    @test haskey(mc.thermalization_measurements, :E) &&
        mc.thermalization_measurements[:E] isa IsingEnergyMeasurement
    delete!(mc, :E, :TH)

    MonteCarlo.unsafe_push!(mc, :E => IsingEnergyMeasurement(mc, model), :TH)
    @test haskey(mc.thermalization_measurements, :E) &&
        mc.thermalization_measurements[:E] isa IsingEnergyMeasurement
    delete!(mc, IsingEnergyMeasurement, :TH)
    @test !haskey(mc.thermalization_measurements, :TH)
end

@testset "Statistics of measurements" begin
    m = IsingModel(dims=2, L=4);
    mc = MC(m, beta=1.0);
    run!(mc, sweeps=100, thermalization=0, verbose=false);

    ms = measurements(mc)[:Energy]
    obs = observables(mc)[:Energy]

    @test mean(ms)["Total energy"] == mean(obs["Total energy"])
    @test var(ms)["Total energy"] == var(obs["Total energy"])
    @test std_error(ms)["Total energy"] == std_error(obs["Total energy"])
    # This wont work because we're not using LightObservables
    # @test tau(ms)["Total energy"] == tau(obs["Total energy"])
end

@testset "Saving and Loading" begin
    model = IsingModel(dims=2, L=2)
    mc = MC(model, beta=1.0)
    run!(mc, thermalization=10, sweeps=10, verbose=false)
    push!(mc, :E => IsingEnergyMeasurement, :TH)

    meas = measurements(mc, :all)
    MonteCarlo.save_measurements("testfile.jld", mc, overwrite=true)
    _meas = load("testfile.jld")
    for (k, v) in meas
        for (k2, v2) in v
            for f in fieldnames(typeof(v2))
                @test getfield(v2, f) == getfield(_meas[k][k2], f)
            end
        end
    end
    rm("testfile.jld")
end

function calc_measured_greens(mc::DQMC, G::Matrix)
    eThalfminus = mc.stack.hopping_matrix_exp
    eThalfplus = mc.stack.hopping_matrix_exp_inv

    eThalfplus * G * eThalfminus
end

@testset "Measured Greens function" begin
    m = HubbardModelAttractive(8, 2, mu=0.5)
    mc = DQMC(m, beta=5.0, safe_mult=1)
    MonteCarlo.build_stack(mc, mc.stack)
    MonteCarlo.propagate(mc)

    # greens(mc) matches expected output
    @test greens(mc).val ≈ calc_measured_greens(mc, mc.stack.greens)

    # wrap greens test
    for k in 0:9
        MonteCarlo.wrap_greens!(mc, mc.stack.greens, mc.stack.current_slice - k, -1)
    end
    # greens(mc) matches expected output
    @test greens(mc).val ≈ calc_measured_greens(mc, mc.stack.greens)
end

@testset "DQMC Measurement constructors" begin
    for m1 in (HubbardModelAttractive(4, 2), HubbardModelRepulsive(4, 2))
        mc = DQMC(m1, beta=1.0, safe_mult=1)

        # Greens
        m = greens_measurement(mc, m1)
        @test m isa MonteCarlo.DQMCMeasurement
        @test m.greens_iterator == Greens()
        @test m.lattice_iterator === nothing
        @test m.kernel == MonteCarlo.greens_kernel
        @test m.observable isa LogBinner{Matrix{Float64}}
        @test m.temp === nothing

        # Occupation
        m = occupation(mc, m1)
        @test m isa MonteCarlo.DQMCMeasurement
        @test m.greens_iterator == Greens()
        @test m.lattice_iterator == EachSiteAndFlavor()
        @test m.kernel == MonteCarlo.occupation_kernel
        @test m.observable isa LogBinner{Vector{Float64}}
        @test m.temp isa Vector{Float64}

        for time in (:equal, :unequal)
            # Charge densities
            if time == :equal
                m = charge_density_correlation(mc, m1)
                @test m.greens_iterator == Greens()
            else
                m = charge_density_susceptibility(mc, m1)
                @test m.greens_iterator == CombinedGreensIterator(mc)
            end
            @test m isa MonteCarlo.DQMCMeasurement
            @test m.lattice_iterator == EachSitePairByDistance()
            @test m.kernel == MonteCarlo.cdc_kernel
            @test m.observable isa LogBinner{Vector{Float64}}
            @test m.temp isa Vector{Float64}

            # Spin densities
            for dir in (:x, :y, :z)
                if time == :equal
                    m = spin_density_correlation(mc, m1, dir)
                    @test m.greens_iterator == Greens()
                else
                    m = spin_density_susceptibility(mc, m1, dir)
                    @test m.greens_iterator == CombinedGreensIterator(mc)
                end
                @test m isa MonteCarlo.DQMCMeasurement
                @test m.lattice_iterator == EachSitePairByDistance()
                @test m.kernel == Core.eval(MonteCarlo, Symbol(:sdc_, dir, :_kernel))
                @test m.observable isa LogBinner{Vector{Float64}}
                @test m.temp isa Vector{Float64}
            end

            # pairings
            if time == :equal
                m = pairing_correlation(mc, m1)
                @test m.greens_iterator == Greens()
            else
                m = pairing_susceptibility(mc, m1)
                @test m.greens_iterator == CombinedGreensIterator(mc)
            end
            @test m isa MonteCarlo.DQMCMeasurement
            @test m.lattice_iterator == EachLocalQuadByDistance(1:5)
            @test m.kernel == MonteCarlo.pc_combined_kernel
            @test m.observable isa LogBinner{Array{Float64, 3}}
            @test m.temp isa Array{Float64, 3}
        end

        # Magnetizations
        for dir in (:x, :y, :z)
            m = magnetization(mc, m1, dir)
            @test m isa MonteCarlo.DQMCMeasurement
            @test m.greens_iterator == Greens()
            @test m.lattice_iterator == EachSite()
            @test m.kernel == Core.eval(MonteCarlo, Symbol(:m, dir, :_kernel))
            @test m.observable isa LogBinner{Vector{Float64}}
            @test m.temp isa Vector{Float64}
        end

        # Current Current susceptibility
        m = current_current_susceptibility(mc, m1, lattice_iterator = EachLocalQuadBySyncedDistance(2:5))
        @test m isa MonteCarlo.DQMCMeasurement
        @test m.greens_iterator == CombinedGreensIterator(mc)
        @test m.lattice_iterator == EachLocalQuadBySyncedDistance(2:5)
        @test m.kernel == MonteCarlo.cc_kernel
        @test m.observable isa LogBinner{Matrix{Float64}}
        @test m.temp isa Matrix{Float64}
    end

    # m =
    # @test m isa MonteCarlo.DQMCMeasurement
    # @test m.greens_iterator ==
    # @test m.lattice_iterator ==
    # @test m.kernel == MonteCarlo.
    # @test m.observable isa LogBinner{Vector{Float64}}
    # @test m.temp isa
end

# @testset "Uniform Fourier" begin
#     A = rand(64, 64)
#     @test uniform_fourier(A, 64) == sum(A) / 64
#     @test uniform_fourier(A, 10) == sum(A) / 10

#     m = HubbardModelAttractive(8, 2)
#     mc = DQMC(m, beta=5.0)
#     @test uniform_fourier(A, mc) == sum(A) / 64

#     mask = MonteCarlo.DistanceMask(MonteCarlo.lattice(m))
#     MonteCarlo.unsafe_push!(mc, :CDC => ChargeDensityCorrelationMeasurement(mc, m, mask=mask))
#     MonteCarlo.unsafe_push!(mc, :SDC => SpinDensityCorrelationMeasurement(mc, m, mask=mask))
#     MonteCarlo.unsafe_push!(mc, :PC => PairingCorrelationMeasurement(mc, m, mask=mask))
#     run!(mc, verbose=false)
#     measured = measurements(mc)

#     @test uniform_fourier(measured[:CDC]) isa MonteCarlo.UniformFourierWrapped
#     @test_throws MethodError uniform_fourier(measured[:SDC])
#     @test uniform_fourier(measured[:SDC], :x) isa MonteCarlo.UniformFourierWrapped
#     @test uniform_fourier(measured[:SDC].y) isa MonteCarlo.UniformFourierWrapped
#     @test uniform_fourier(measured[:PC]) isa MonteCarlo.UniformFourierWrapped

#     @test mean(uniform_fourier(measured[:CDC])) == sum(mean(measured[:CDC])) / 64
#     @test var(uniform_fourier(measured[:SDC], :x)) == sum(var(measured[:SDC].x)) / 64
#     @test std_error(uniform_fourier(measured[:SDC].z)) == sum(std_error(measured[:SDC].z)) / 64
#     @test tau(uniform_fourier(measured[:PC])) == maximum(tau(measured[:PC]))
# end


text:
@testset "Thunk" begin
    @test @thunk(3) isa Thunk

    @testset "show" begin
        rep = repr(Thunk(rand))
        @test occursin(r"Thunk\(.*rand.*\)", rep)
    end

    @testset "Externing" begin
        @test extern(@thunk(3)) == 3
        @test extern(@thunk(@thunk(3))) == 3
    end

    @testset "unthunk" begin
        @test unthunk(@thunk(3)) == 3
        @test unthunk(@thunk(@thunk(3))) isa Thunk
    end

    @testset "calling thunks should call inner function" begin
        @test (@thunk(3))() == 3
        @test (@thunk(@thunk(3)))() isa Thunk
    end

    @testset "erroring thunks should include the source in the backtrack" begin
        expected_line = (@__LINE__) + 2  # for testing it is at right palce
        try
            x = @thunk(error())
            extern(x)
        catch err
            err isa ErrorException || rethrow()
            st = stacktrace(catch_backtrace())
            # Should be 2nd last line, as last line will be the `error` function
            stackframe = st[2]
            @test stackframe.line == expected_line
            @test stackframe.file == Symbol(@__FILE__)
        end
    end

    @testset "Broadcast" begin
        @testset "Array" begin
            was_unthunked = 0
            array_thunk = @thunk begin
                was_unthunked += 1;
                [1.0, 2.0, 3.0]
            end

            was_unthunked = 0
            @test array_thunk .+ fill(10, 3) .+  fill(10, 3) == [21.0, 22.0, 23.0]
            @test was_unthunked == 1

            was_unthunked = 0
            @test array_thunk .+ 10.0 .+ 10.0 == [21.0, 22.0, 23.0]
            @test was_unthunked == 1

        end

        @testset "Scalar" begin
            was_unthunked=0
            scalar_thunk = @thunk begin
                was_unthunked += 1;
                sqrt(4.0)
            end

            was_unthunked = 0
            @test scalar_thunk .+ fill(10, 3) .+  fill(10, 3) == [22.0, 22.0, 22.0]
            @test was_unthunked == 1

            was_unthunked = 0
            @test scalar_thunk .+ 10.0 .+ 10.0 == 22.0
            @test was_unthunked == 1
        end
    end
end


text:
# AbstractGridMesh

export AbstractGridMesh
"""
    abstract type AbstractGridMesh{S,T,D} <: AbstractGridFunction{S,T,D} end
"""
abstract type AbstractGridMesh{S,T,D} <: AbstractGridFunction{S,T,D} end

export bbox
bbox(::AbstractGridMesh) = error("undefined")

export test_AbstractGridMesh
function test_AbstractGridMesh(gf::AbstractGridMesh)
    test_AbstractGridFunction(gf)

    @testset "$(typeof(gf)) <: AbstractGridMesh" begin
        D = ndims(gf)
        gfbbox = bbox(gf)
        @test gfbbox isa BBox{D}

        gfdomain = domain(gf)

        for n in 1:10
            ind = map((a, b) -> rand(a:b), min(gfbbox), max(gfbbox))
            coord = coordinate(gf, ind)
            @test all(min(gfdomain) .≤ coord .≤ max(gfdomain))
            ind′ = index(gf, coord)
            @test ind′ ≈ ind
        end
    end
end


text:
type DropoutDesc
    ptr::Ptr{Void}

    function DropoutDesc(h, droprate::Float64; seed=0)
        p = Ptr{Void}[0]
        cudnnCreateDropoutDescriptor(p)
        desc = new(p[1])
        finalizer(desc, cudnnDestroyDropoutDescriptor)

        p = Cint[0]
        cudnnDropoutGetStatesSize(h, p)
        statessize = p[1]
        states = CuArray{Int8}(Int(statessize))
        cudnnSetDropoutDescriptor(desc, h, droprate, states, statessize, seed)
        desc
    end
end

Base.unsafe_convert(::Type{Ptr{Void}}, desc::DropoutDesc) = desc.ptr

function dropout{T}(x::CuArray{T}, droprate::Float64)
    h = handle(x)
    dropdesc = DropoutDesc(h, droprate)

    xdesc = TensorDesc(x)
    p = Cint[0]
    cudnnDropoutGetReserveSpaceSize(xdesc, p)
    reservesize = p[1]
    reservespace = CuArray{Int8}(Int(reservesize))

    y = similar(x)
    cudnnDropoutForward(h, dropdesc, xdesc, x, xdesc, y, reservespace, reservesize)

    function backward!(gy, gx)
        isvoid(gx) && return
        cudnnDropoutBackward(h, dropdesc, xdesc, dy, xdesc, dx, reservespace, reservesize)
    end
    y, backward!
end


text:
function StatsBase.fit{N}(m::AAModel, img::Array{Float64, N},
                          init_shape::Shape, max_it::Int)
    @assert size(m.s0, 1) == size(init_shape, 1) * 2 "Shape has wrong size!"
    cur_shape = init_shape
    mean_shape = reshape(m.s0, m.np, 2)
    comp_warp = zeros(size(init_shape))
    iter = 1
    while iter <= max_it
        warped = warp_to_mean_shape(m, img, cur_shape)
        error_img = flatten(warped) .- m.A0
        println("Sum of squared errors: $(sum(error_img .^ 2))")
        if iter > 5 || max_it < 10
            delta_qp = m.R * error_img
            A, trans = q_params_to_affine(m, -delta_qp[1:4])
            d_s0_ = reshape(m.s0 - m.S * delta_qp[5:end], m.np, 2)
            d_s0 = d_s0_ * A + repmat(trans, m.np, 1)
            comp_warp = compose_warps(m, cur_shape, d_s0)
        else
            delta_q = m.R[1:4, :] * error_img
            A, trans = q_params_to_affine(m, -delta_q)
            d_s0_ = reshape(m.s0, m.np, 2)
            d_s0 = (d_s0_ * A + repmat(trans, m.np, 1)) - d_s0_
            comp_warp = cur_shape + d_s0
        end
        iter += 1
        cur_shape = comp_warp
    end
    fitted_shape = cur_shape
    # wrp = pa_warp(img, cur_shape, mean_shape, m.trigs)
    warped = warp_to_mean_shape(m, img, cur_shape)
    error_img = flatten(warped) - m.A0
    fitted_app = reshape(m.A0 + m.A * (m.A' * error_img),
                         m.frame.h, m.frame.w, m.nc)
    return fitted_shape, fitted_app
end


text:
function payoff(S1::abstractMatrix, bmPayoff::BermudanPayoff, rfCurve::abstractZeroRateCurve, T1::num2 = maturity(bmPayoff)) where {abstractZeroRateCurve <: AbstractZeroRateCurve, abstractMatrix <: AbstractMatrix{num}} where {num <: Number, num2 <: Number}
    T_opt = bmPayoff.T_ex
    T = maturity(bmPayoff)
    NStep = size(S1, 2) - 1
    index1 = round(Int, T / T1 * NStep) + 1
    S = collect(S1[:, 1:index1])
    S0 = first(S)
    (Nsim, Nstep) = size(S)
    Nstep -= 1
    r = rfCurve.r
    dt = T / Nstep
    # initialize
    exerciseTimes = Nstep .* ones(Nsim)
    phi(Sti::numtype_) where {numtype_ <: Number} = payout(Sti, bmPayoff)
    idx_opt = round.(Int, T_opt ./ dt)
    @views V = phi.(S[:, end]) #payoff
    # Backward Procedure
    for j in reverse(idx_opt)
        @views Tmp = phi.(S[:, j])
        inMoneyIndexes = findall(Tmp .> 0.0)
        if !isempty(inMoneyIndexes)
            #S_I=S[inMoneyIndexes,j];
            S_I = view(S, inMoneyIndexes, j)
            #-- Intrinsic Value
            @views IV = Tmp[inMoneyIndexes]
            #-- Continuation Value
            #- Linear Regression on Quadratic Form
            A = [ones(length(S_I)) S_I S_I .^ 2]
            #@views b=V[inMoneyIndexes].*exp.(-r*dt*(exerciseTimes[inMoneyIndexes].-j));
            @views b = V[inMoneyIndexes] .* exp.(-[integral(r, dt * (exerciseTime - j)) for exerciseTime in exerciseTimes[inMoneyIndexes]])
            #MAT=A'*A;
            LuMat = lu(A' * A)
            Btilde = A' * b
            alpha = LuMat \ Btilde
            #alpha=A\b;
            #Continuation Value
            CV = A * alpha
            #----------
            # Find premature exercise times
            Index = findall(IV .> CV)
            @views exercisePositions = inMoneyIndexes[Index]
            # Update the outputs
            @views V[exercisePositions] = IV[Index]
            @views exerciseTimes[exercisePositions] .= j - 1
        end
    end
    Out = V .* exp.(-[integral(r, dt * exerciseTime) for exerciseTime in exerciseTimes])

    return Out
end


text:
using Pkg
Pkg.instantiate()
Pkg.develop(PackageSpec(path=pwd()))
Pkg.build("EzXML")


text:
export ft_write

function ft_write(ft_handle::Culong,
                  buffer::Array{UInt8,1},
                  bytestowrite::Integer = length(buffer))
  byteswritten = Ref{Cuint}()
  ft_status = ccall((:FT_Write, d2xx),
                     Cuint,
                     (Culong, Ptr{UInt8}, Cuint, Ref{Cuint}),
                     ft_handle, buffer, bytestowrite, byteswritten)
  checkstatus(ft_status)
  return byteswritten[]
end


text:
using Distributions
using Klara
using PGUManifoldMC

include("rv_model.jl")

using RvModelKeplerian

DATADIR = "../../data"
SUBDATADIR = "amsmmala"

dataset = readdlm(joinpath(DATADIR, "example2.csv"), ',', header=false);  # read observational data
obs_times = dataset[:,1]
obs_rv = dataset[:,2]
sigma_obs = dataset[:,3]
set_times(obs_times);     # set data to use for model evaluation
set_obs( obs_rv);
set_sigma_obs(sigma_obs);

include("utils_ex.jl")
param_true = make_param_true_ex2()
param_perturb_scale = make_param_perturb_scale(param_true)
param_init = 0
param_init = param_true.+0.01*param_perturb_scale.*randn(length(param_true))
println("param_init= ",param_init)

sampler = AMSMMALA(
  0.06,
  # update=(sstate, pstate, i, tot) -> mod_update!(sstate, pstate, i, tot, 10),
  update=(sstate, pstate, i, tot) -> rand_exp_decay_update!(sstate, pstate, i, tot, 7.),
  transform=H -> softabs(H, 1000.)
)

nmcmc = 50000
nburnin = 10000
mcrange = BasicMCRange(nsteps=nmcmc, burnin=nburnin)

outopts = Dict{Symbol, Any}(:monitor=>[:value, :logtarget], :diagnostics=>[:accept])

v0 = Dict(:p=>param_init)

p = BasicContMuvParameter(
  :p,
  logtarget=plogtarget,
  diffopts=DiffOptions(mode=:forward, order=2)
)

model = likelihood_model(p, false)
#model = likelihood_model([p], isindexed=false)

# target_accept_rates = [0.1, 0.2, 0.25, 0.3, 0.5, 0.75, 0.8, 0.85, 0.9, 0.93, 0.95, 0.97, 0.99]
target_accept_rates = [0.234]
#target_accept_rates = [0.20, 0.234, 0.25, 0.275]
#target_accept_rates = [0.234]
amsmmala_times = Array(Float64, length(target_accept_rates))
amsmmala_stepsizes = Array(Float64, length(target_accept_rates))
amsmmala_acceptance = Array(Float64, length(target_accept_rates))
amsmmala_esses = Array(Float64, length(target_accept_rates))
amsmmala_iacts = Array(Float64, length(target_accept_rates))
amsmmala_chains = Array(Any, length(target_accept_rates))

for i in 1:length(target_accept_rates)
  target_accept_rate = target_accept_rates[i]
  println("# i= ", i, ": ",target_accept_rate)
  mctuner = PSMMALAMCTuner(
    VanillaMCTuner(verbose=true),
    VanillaMCTuner(verbose=true),
    AcceptanceRateMCTuner(target_accept_rate, verbose=true)
  )

  job = BasicMCJob(
    model,
    sampler,
    mcrange,
    v0,
    tuner=mctuner,
    outopts=outopts
  )

  tic()
  run(job)
  runtime = toc()

  chain = output(job)
  ratio = acceptance(chain)

  println("# mean[",i,"] = ",mean(chain))
  println("# std[",i,"] =",vec(std(chain.value,2)))
  println("# ess[",i,"] = ",ess(chain))

  amsmmala_times[i] = runtime
  amsmmala_stepsizes[i] = job.sstate.tune.totaltune.step
  amsmmala_acceptance[i] = ratio
  amsmmala_esses[i] = minimum(ess(chain))
  amsmmala_iacts[i] = maximum(iact(chain))
  amsmmala_chains[i] = chain
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: pinpoint_sms_voice
using AWS.Compat
using AWS.UUIDs

"""
    CreateConfigurationSet()

Create a new configuration set. After you create the configuration set, you can add one or more event destinations to it.

# Optional Parameters
- `ConfigurationSetName`: The name that you want to give the configuration set.
"""
create_configuration_set(; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("POST", "/v1/sms-voice/configuration-sets"; aws_config=aws_config)
create_configuration_set(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("POST", "/v1/sms-voice/configuration-sets", args; aws_config=aws_config)

"""
    CreateConfigurationSetEventDestination()

Create a new event destination in a configuration set.

# Required Parameters
- `ConfigurationSetName`: ConfigurationSetName

# Optional Parameters
- `EventDestination`:
- `EventDestinationName`: A name that identifies the event destination.
"""
create_configuration_set_event_destination(ConfigurationSetName; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("POST", "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations"; aws_config=aws_config)
create_configuration_set_event_destination(ConfigurationSetName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("POST", "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations", args; aws_config=aws_config)

"""
    DeleteConfigurationSet()

Deletes an existing configuration set.

# Required Parameters
- `ConfigurationSetName`: ConfigurationSetName

"""
delete_configuration_set(ConfigurationSetName; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("DELETE", "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)"; aws_config=aws_config)
delete_configuration_set(ConfigurationSetName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("DELETE", "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)", args; aws_config=aws_config)

"""
    DeleteConfigurationSetEventDestination()

Deletes an event destination in a configuration set.

# Required Parameters
- `ConfigurationSetName`: ConfigurationSetName
- `EventDestinationName`: EventDestinationName

"""
delete_configuration_set_event_destination(ConfigurationSetName, EventDestinationName; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("DELETE", "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations/$(EventDestinationName)"; aws_config=aws_config)
delete_configuration_set_event_destination(ConfigurationSetName, EventDestinationName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("DELETE", "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations/$(EventDestinationName)", args; aws_config=aws_config)

"""
    GetConfigurationSetEventDestinations()

Obtain information about an event destination, including the types of events it reports, the Amazon Resource Name (ARN) of the destination, and the name of the event destination.

# Required Parameters
- `ConfigurationSetName`: ConfigurationSetName

"""
get_configuration_set_event_destinations(ConfigurationSetName; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("GET", "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations"; aws_config=aws_config)
get_configuration_set_event_destinations(ConfigurationSetName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("GET", "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations", args; aws_config=aws_config)

"""
    ListConfigurationSets()

List all of the configuration sets associated with your Amazon Pinpoint account in the current region.

# Optional Parameters
- `NextToken`: A token returned from a previous call to the API that indicates the position in the list of results.
- `PageSize`: Used to specify the number of items that should be returned in the response.
"""
list_configuration_sets(; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("GET", "/v1/sms-voice/configuration-sets"; aws_config=aws_config)
list_configuration_sets(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("GET", "/v1/sms-voice/configuration-sets", args; aws_config=aws_config)

"""
    SendVoiceMessage()

Create a new voice message and send it to a recipient's phone number.

# Optional Parameters
- `CallerId`: The phone number that appears on recipients' devices when they receive the message.
- `ConfigurationSetName`: The name of the configuration set that you want to use to send the message.
- `Content`:
- `DestinationPhoneNumber`: The phone number that you want to send the voice message to.
- `OriginationPhoneNumber`: The phone number that Amazon Pinpoint should use to send the voice message. This isn't necessarily the phone number that appears on recipients' devices when they receive the message, because you can specify a CallerId parameter in the request.
"""
send_voice_message(; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("POST", "/v1/sms-voice/voice/message"; aws_config=aws_config)
send_voice_message(args::AbstractDict{String, Any}; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("POST", "/v1/sms-voice/voice/message", args; aws_config=aws_config)

"""
    UpdateConfigurationSetEventDestination()

Update an event destination in a configuration set. An event destination is a location that you publish information about your voice calls to. For example, you can log an event to an Amazon CloudWatch destination when a call fails.

# Required Parameters
- `ConfigurationSetName`: ConfigurationSetName
- `EventDestinationName`: EventDestinationName

# Optional Parameters
- `EventDestination`:
"""
update_configuration_set_event_destination(ConfigurationSetName, EventDestinationName; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("PUT", "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations/$(EventDestinationName)"; aws_config=aws_config)
update_configuration_set_event_destination(ConfigurationSetName, EventDestinationName, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = pinpoint_sms_voice("PUT", "/v1/sms-voice/configuration-sets/$(ConfigurationSetName)/event-destinations/$(EventDestinationName)", args; aws_config=aws_config)


text:
function create_param_dual_cum_sum(model::ParametricOptimizer)
    return zeros(model.number_of_parameters_in_model)
end

function calculate_dual_of_parameters(model::ParametricOptimizer)
    param_dual_cum_sum = create_param_dual_cum_sum(model)
    update_duals_with_affine_constraint_cache!(param_dual_cum_sum, model)
    update_duals_with_quadratic_constraint_cache!(param_dual_cum_sum, model)
    update_duals_in_affine_objective!(param_dual_cum_sum, model.affine_objective_cache)
    update_duals_in_quadratic_objective!(param_dual_cum_sum, model.quadratic_objective_cache_pc)
    empty_and_feed_duals_to_model(model, param_dual_cum_sum)
    return model
end

function update_duals_with_affine_constraint_cache!(param_dual_cum_sum::Vector{Float64}, model::POI.ParametricOptimizer)
    for S in SUPPORTED_SETS
        affine_constraint_cache_inner = model.affine_constraint_cache[MOI.ScalarAffineFunction{Float64}, S]
        if !isempty(affine_constraint_cache_inner)
            update_duals_with_affine_constraint_cache!(param_dual_cum_sum, model.optimizer, affine_constraint_cache_inner)
        end
    end
    return
end

function update_duals_with_affine_constraint_cache!(param_dual_cum_sum::Vector{Float64}, optimizer::OT, affine_constraint_cache_inner::DD.WithType{F, S}) where {OT,F,S}
    for (ci, param_array) in affine_constraint_cache_inner
        calculate_parameters_in_ci!(param_dual_cum_sum, optimizer, param_array, ci)
    end
    return
end

function update_duals_with_quadratic_constraint_cache!(param_dual_cum_sum::Vector{Float64}, model::POI.ParametricOptimizer)
    for S in SUPPORTED_SETS
        quadratic_constraint_cache_pc_inner = model.quadratic_constraint_cache_pc[MOI.ScalarQuadraticFunction{Float64}, S]
        if !isempty(quadratic_constraint_cache_pc_inner)
            update_duals_with_quadratic_constraint_cache!(param_dual_cum_sum, model, quadratic_constraint_cache_pc_inner)
        end
    end
    return
end

function update_duals_with_quadratic_constraint_cache!(param_dual_cum_sum::Vector{Float64}, model::ParametricOptimizer, quadratic_constraint_cache_pc_inner::DD.WithType{F, S}) where {F,S}
    for (poi_ci, param_array) in quadratic_constraint_cache_pc_inner
        moi_ci = model.quadratic_added_cache[poi_ci]
        calculate_parameters_in_ci!(param_dual_cum_sum, model.optimizer, param_array, moi_ci)
    end
    return
end

function calculate_parameters_in_ci!(param_dual_cum_sum::Vector{Float64}, optimizer::OT, param_array::Vector{MOI.ScalarAffineTerm{T}}, ci::CI) where {OT, CI, T}
    cons_dual = MOI.get(optimizer, MOI.ConstraintDual(), ci)

    for param in param_array
        param_dual_cum_sum[param.variable_index.value - PARAMETER_INDEX_THRESHOLD] += cons_dual*param.coefficient
    end
    return
end

function update_duals_in_affine_objective!(param_dual_cum_sum::Vector{Float64}, affine_objective_cache::Vector{MOI.ScalarAffineTerm{T}}) where T
    for param in affine_objective_cache
        param_dual_cum_sum[param.variable_index.value - PARAMETER_INDEX_THRESHOLD] += param.coefficient
    end
    return
end

function update_duals_in_quadratic_objective!(param_dual_cum_sum::Vector{Float64}, quadratic_objective_cache_pc::Vector{MOI.ScalarAffineTerm{T}}) where T
    for param in quadratic_objective_cache_pc
        param_dual_cum_sum[param.variable_index.value - PARAMETER_INDEX_THRESHOLD] += param.coefficient
    end
    return
end

function empty_and_feed_duals_to_model(model::ParametricOptimizer, param_dual_cum_sum::Vector{Float64})
    empty!(model.dual_value_of_parameters)
    model.dual_value_of_parameters = zeros(model.number_of_parameters_in_model)
    for (vi_val_minus_threshold, param_dual) in enumerate(param_dual_cum_sum)
        model.dual_value_of_parameters[vi_val_minus_threshold] = param_dual
    end
    return
end

struct ParameterDual <: MOI.AbstractVariableAttribute end

function MOI.get(model::ParametricOptimizer, ::ParameterDual, vi_val::Int)
    if !is_additive(model, MOI.ConstraintIndex{MOI.SingleVariable, Parameter}(vi_val))
        error("Cannot calculate the dual of a multiplicative parameter")
    end
    return model.dual_value_of_parameters[vi_val - PARAMETER_INDEX_THRESHOLD]
end

function MOI.get(model::ParametricOptimizer, ::MOI.ConstraintDual, cp::MOI.ConstraintIndex{MOI.SingleVariable,POI.Parameter})
    if !is_additive(model, cp)
        error("Cannot calculate the dual of a multiplicative parameter")
    end
    return model.dual_value_of_parameters[cp.value - PARAMETER_INDEX_THRESHOLD]
end

function MOI.get(model::ParametricOptimizer, attr::MOI.ConstraintDual, ci::MOI.ConstraintIndex)
    return MOI.get(model.optimizer, attr, ci)
end

function is_additive(model::ParametricOptimizer, cp::MOI.ConstraintIndex)
    if cp.value in model.multiplicative_parameters
        return false
    end
    return true
end


text:
using VortexDistributions, FourierGPE, Plots, LaTeXStrings, SpecialFunctions

## Initialize simulation
# units of healing length, chemical potential
L = (100.,100.)
N = (512,512)
sim = Sim(L,N)
@unpack_Sim sim

## homogeneous state
μ = 1.0
g = 0.01
ψ0(x,y,μ,g) = sqrt(μ/g)*sqrt(max(1.0-V(x,y,0.0)/μ,0.0)+im*0.0)
healinglength(x,y,μ,g) = 1/sqrt(g*abs2(ψ0(x,y,μ,g)))
x,y = X
kx,ky = K
ψg = ψ0.(x,y',μ,g)
showpsi(x,y,ψg)

## initial dipole, with periodic phase
ψd = copy(ψg)
d = 12
ξv = healinglength(0.,0.,μ,g)
# can use default ξ = 1
xp,yp = 0.,d/2
xn,yn = 0.,-d/2
pv = PointVortex(xp,yp,1)
nv = PointVortex(xn,yn,-1)
dipole = [nv;pv]

psi = Torus(ψd,x,y) # set field topology for VortexDistributions
vortex!(psi,dipole) # make dipole

## make another dipole
# xp,yp = 6.,d/2
# xn,yn = 6.,-d/2
# pv = PointVortex(xp,yp,1)
# nv = PointVortex(xn,yn,-1)
# dipole = [nv;pv]
# vortex!(psi,dipole) # make dipole

showpsi(x,y,ψd)

## test new methods
kd = 2*pi/d
kL = 2*pi/L[1]
kξ = 2*pi
kmin = 0.1*kL
kmax = kξ

Np = 300
k = log10range(kmin,kmax,Np)

# @time Ek = kespectrum(k,ψd,X,K)
# plot(k,Ek,scale=:log10)

# ## incompressible spectrum
# @time Eki = ikespectrum(k,ψd,X,K)
# plot!(k,Eki,scale=:log10)

## power-law plot in logspace
@time Eki = ikespectrum(k,ψd,X,K)
Eqp = qpespectrum(k,ψd,X,K)
Ekc = ckespectrum(k,ψd,X,K)

n0 = abs2.(ψd[1,1])
E0 = π*n0
p1 = plot(k,Eki/E0,scale=:log10,grid=false,label=L"\epsilon_I(k)",legend=:bottomleft)
plot!(k,Eqp/E0,label=L"\epsilon_Q(k)")
plot!(k,2.5k.^(-3),label=L"k^{-3}")
plot!(k,200*k,label=L"k")
ylims!(5e-4,1e2)
vline!([kL],ls=:dash,label=L"k_L")
vline!([kd],ls=:dash,label=L"k_d")
vline!([kξ],ls=:dash,label=L"k_\xi")
xlabel!(L"k\xi")
ylabel!(L"\epsilon_I(k)")

# Analytic form (pure dipole)
# vort = findvortices(Torus(ψd,x,y)) |> vortex_array
# d = vort[:,2] |> diff
# Λ = 0.8249
# f(x) = x*(besselk(1,x)*besseli(0,x)-besselk(0,x)*besseli(1,x))
# F(x,Λ) = f(x/(2Λ))^2/x
# F(x) = F(x,Λ)
# Ed(k,d) = 2*F(k)*(1-besselj0(k*d))
# plot!(k,Ed.(k,d),label=L"\epsilon_{a}^v(k)")

## two point correlator (approx)
ddk = diff(k); push!(ddk,last(ddk))
gtwo(r)=sum(@. k^2*Eki*besselj0(k*r)*ddk)
Gtwo(r)=gtwo(r)/gtwo(0)

## plot Gtwo
r = LinRange(0,10,100)
plot(r,Gtwo.(r),label=L"G_2(r)")
vline!([d],ls=:dash,label=L"d")

## FULL SPECTRAL DENSITIES

## ke spectrum
kd = 2*pi/d
kL = 2*pi/L[1]
kξ = 2*pi
kmin = 0.1*kL
kmax = kξ

Np = 300
k = log10range(kmin,kmax,Np)

## quantum spectra
# all kinetic
Ek = kespectrum(k,ψd,X,K)
plot(k,Ek,scale=:log10,label="total",legend=:bottom)
# plot!(k,k.^3,label=L"k^3")

## try Eh^v spectrum (incompressible)
using SpecialFunctions

function qikespectrum(k,ψ,X,K)
    x,y = X; kx,ky = K
 	dx,dy = diff(x)[1],diff(y)[1]
	DX,DK = dfftall(X,K)
    k2field = k2(K)
    psi = XField(ψ,X,K,k2field)
    vx,vy = velocity(psi)

    ux = @. abs(ψ)*vx; uy = @. abs(ψ)*vy
    Wi, Wc = helmholtz(ux,uy,psi)
    wix,wiy = Wi
    @. wix *= exp(im*angle(ψ)) # restore phase factors
    @. wiy *= exp(im*angle(ψ))
	cix = autocorrelate(wix,X,K)
	ciy = autocorrelate(wiy,X,K)
    C = cix .+ ciy

    return bessel_projection(k,x,y,C)
end

Eki = qikespectrum(k,ψd,X,K)
plot!(k,Eki,label="incompressible")
plot!(k,1000k.^(-3))

## sound
function skespectrum(k,ψ,X,K)
    x,y = X; kx,ky = K
 	dx,dy = diff(x)[1],diff(y)[1]
	DX,DK = dfftall(X,K)
    k2field = k2(K)
    psi = XField(ψ,X,K,k2field)
    vx,vy = velocity(psi)

    ux = @. abs(ψ)*vx; uy = @. abs(ψ)*vy
    Wi, Wc = helmholtz(ux,uy,psi)
    wcx,wcy = Wc
    @. wcx *= exp(im*angle(ψ)) # restore phase factors
    @. wcy *= exp(im*angle(ψ))
	ccx = autocorrelate(wcx,X,K)
	ccy = autocorrelate(wcy,X,K)
    C = ccx .+ ccy

    return bessel_projection(k,x,y,C)
end

Ekc = ckespectrum(k,ψd,X,K)
plot!(k,Ekc,label="compressible")

## quantum pressure
function qkespectrum(k,ψ,X,K)
    x,y = X; kx,ky = K
 	dx,dy = diff(x)[1],diff(y)[1]
	DX,DK = dfftall(X,K)
    k2field = k2(K)
    psi = XField(abs.(ψ) |> complex,X,K,k2field)
    rnx,rny = gradient(psi)
    @. rnx *= exp(im*angle(ψ)) # restore phase factors
    @. rny *= exp(im*angle(ψ)) # restore phase factors
	cx = autocorrelate(rnx,X,K)
	cy = autocorrelate(rny,X,K)
    C = cx .+ cy

    return bessel_projection(k,x,y,C)
end

Ekq = qkespectrum(k,ψd,X,K)
plot!(k,Ekq,label="quantum pressure")
ylims!(0.1,1e5)
vline!([kL],ls=:dash,label=L"k_L")
vline!([kd],ls=:dash,label=L"k_d")
vline!([kξ],ls=:dash,label=L"k_\xi")


text:
using GridapGmsh
using GridapGmsh: gmsh, GmshDiscreteModel

function add_pin!(gmsh, o, r, t, lc)
    factory = gmsh.model.geo

    # inner and outer circle
    l1 = add_circle!(gmsh, o, r, lc)
    l2 = add_circle!(gmsh, o, r + t, lc)

    s1 = factory.addPlaneSurface([l1])
    s2 = factory.addPlaneSurface([l2, l1]) # l1 is a hole

    return s1, s2
end

function add_circle!(gmsh, o, r, lc)
    factory = gmsh.model.geo
    ox, oy = o

    p1 = factory.addPoint(ox, oy, 0, lc)  # center - origin
    p2 = factory.addPoint(ox + r, oy, 0, lc)  # right
    p3 = factory.addPoint(ox, oy + r, 0, lc)  # up
    p4 = factory.addPoint(ox - r, oy, 0, lc)  # left
    p5 = factory.addPoint(ox, oy - r, 0, lc)  # down

    c1 = factory.addCircleArc(p2, p1, p3)
    c2 = factory.addCircleArc(p3, p1, p4)
    c3 = factory.addCircleArc(p4, p1, p5)
    c4 = factory.addCircleArc(p5, p1, p2)

    l1 = factory.addCurveLoop([c1, c2, c3, c4])

    return l1
end

function add_square!(gmsh, s)
    factory = gmsh.model.geo

    p1 = factory.addPoint(0, 0, 0, lc)
    p2 = factory.addPoint(s, 0, 0, lc)
    p3 = factory.addPoint(s, s, 0, lc)
    p4 = factory.addPoint(0, s, 0, lc)

    l1 = factory.addLine(p1, p2)
    l2 = factory.addLine(p2, p3)
    l3 = factory.addLine(p3, p4)
    l4 = factory.addLine(p4, p1)

    cl1 = factory.addCurveLoop([l1, l2, l3, l4])

    return cl1
end

N = 4

# in cm
p = 1.6        # pitch
ri = 0.5       # internal radius
t = 0.1        # wall thickness
ro = ri + t    # external radius
lc = 0.1

gmsh.initialize()
gmsh.model.add("bwr")

# tags
pinTags = Int32[]
cladTags = Int32[]
gdPinTags = Int32[]

# gd pins positions
GD_pos = [(2, 3), (3, 2)]

for i in 1:N, j in 1:N
    xo = p / 2 + (i - 1) * p
    yo = p / 2 + (j - 1) * p

    pinTag, cladTag = add_pin!(gmsh, (xo, yo), ri, t, lc)
    push!(cladTags, cladTag)
    if (i, j) in GD_pos
        push!(gdPinTags, pinTag)
    else
        push!(pinTags, pinTag)
    end
end

#! TODO: h2oTag = add_reflector!(gmsh, 4p)
s = N * p
factory = gmsh.model.geo

p1 = factory.addPoint(0, 0, 0, lc)
p2 = factory.addPoint(s, 0, 0, lc)
p3 = factory.addPoint(s, s, 0, lc)
p4 = factory.addPoint(0, s, 0, lc)

l1 = factory.addLine(p1, p2)
l2 = factory.addLine(p2, p3)
l3 = factory.addLine(p3, p4)
l4 = factory.addLine(p4, p1)

cl1 = factory.addCurveLoop([l1, l2, l3, l4])

h2oTag = factory.addPlaneSurface(vcat(cl1, pinTags, cladTags, gdPinTags))

pg1 = gmsh.model.geo.addPhysicalGroup(2, pinTags)
pg2 = gmsh.model.geo.addPhysicalGroup(2, cladTags)
pg3 = gmsh.model.geo.addPhysicalGroup(2, gdPinTags)
pg4 = gmsh.model.geo.addPhysicalGroup(2, [h2oTag])
pg5 = gmsh.model.geo.addPhysicalGroup(1, [l1])
pg6 = gmsh.model.geo.addPhysicalGroup(1, [l2])
pg7 = gmsh.model.geo.addPhysicalGroup(1, [l3])
pg8 = gmsh.model.geo.addPhysicalGroup(1, [l4])

gmsh.model.setPhysicalName(2, pg1, "pin")
gmsh.model.setPhysicalName(2, pg2, "cladding")
gmsh.model.setPhysicalName(2, pg3, "pin-gd")
gmsh.model.setPhysicalName(2, pg4, "water")
gmsh.model.setPhysicalName(1, pg5, "bottom")
gmsh.model.setPhysicalName(1, pg6, "right")
gmsh.model.setPhysicalName(1, pg7, "top")
gmsh.model.setPhysicalName(1, pg8, "left")

# removemos puntos de la geometria duplicados (origenes por ejemplo)
gmsh.model.geo.removeAllDuplicates()

gmsh.model.geo.synchronize()

gmsh.model.mesh.generate(2)

gmsh.write("bwr.msh")

if !("-nopopup" in ARGS)
    gmsh.fltk.run()
end

gmsh.finalize()

using Gridap
mshfile = joinpath(@__DIR__,"../bwr.msh")
model = GmshDiscreteModel(mshfile; renumber=true)
Gridap.Io.to_json_file(model, "bwr.json")

jsonfile = joinpath(@__DIR__,"../bwr.json")
model = DiscreteModelFromFile(jsonfile)

# number of azimuthal angles
nφ = 16

# azimuthal spacing
δ = 0.002

# boundary conditions
bcs = BoundaryConditions(top=Reflective, bottom=Reflective, left=Reflective, right=Reflective)

# initialize track generator
tg = TrackGenerator(model, nφ, δ, bcs=bcs)

# perform ray tracing
trace!(tg)

# proceed to segmentation
segmentize!(tg)


text:
using CUDAdrv, CUDAnative
using CuArrays
using TimerOutputs

function backendname(useCUDA)
    return useCUDA ? "GPU" : "CPU"
end

function run()

    to = TimerOutput()

    nx = 1000
    ny = 1000
    N = 1000000

    # Model parameters
    dx = 10.0
    dy = 10.0
    dt = 10.0
    K  = 1.0e-4

    println("Hit c for cuda...")
    c = readline()
    useCUDA = (c == "c")

    h = fill(5.0, N)
    # for i = 1:N
    #     h[i] = 10.0*i/N
    # end
    u = fill(2.0, N)
    v = fill(1.0, N)

    if useCUDA
        h = CuArray(h)
        u = CuArray(u)
        v = CuArray(v)
    end
    result = similar(h)
    state  = (h, u, v)
    params = (dx, dy, dt, K, nx, ny)

    # first separate kernels:
    fill!(result, 0.0)
    println("separate, in: ", h[1:10])
    @timeit to "separate " * backendname(useCUDA) begin
    gridloop(advection, result, state, params)
    gridloop(diffusion, result, state, params)
    end
    println("result:  ", result[1:10])

    # then merged into one:
    fill!(result, 0.0)
    println("fused, in: ", h[1:10])
    @timeit to "fused " * backendname(useCUDA) begin
    gridloop(kernel_fused, result, state, params)
    end
    println("result:  ", result[1:10])

    # then composed:
    fill!(result, 0.0)
    println("composed, in: ", h[1:10])
    @timeit to "composed " * backendname(useCUDA) begin
    gridloop(kernel_composed, result, state, params)
    end
    println("result:  ", result[1:10])

    print_timer(to)
end

@inline function advection((i, j), result, (h, u, v), (dx, dy, dt, K, nx, ny))
    # u*du/dx
    uxm = (h[i,j] - h[i - 1,j]) / dx
    uxp = (h[i + 1,j] - h[i,j]) / dx
    am  = 0.5 * min(u[i,j] + u[i + 1,j], 0.0)
    ap  = 0.5 * max(u[i - 1,j] + u[i,j], 0.0)
    advx = ap * uxm + am * uxp

    # v*du/dy
    uym = (h[i,j] - h[i,j - 1]) / dy
    uyp = (h[i,j + 1] - h[i,j]) / dy
    am  = 0.5 * min(v[i,j] + v[i,j + 1], 0.0)
    ap  = 0.5 * max(v[i,j - 1] + v[i,j], 0.0)
    advy = ap * uym + am * uyp

    # adv = 0.5*( (u[i,j]+u[i-1,j])*(u[i,j]-u[i-1,j])/dx + (v[i,j]+v[i,j-1])*(u[i,j]-u[i,j-1])/dy)
    @inbounds result[i,j] -= dt * (advx + advy)
end

@inline function diffusion((i, j), result, (h, u, v), (dx, dy, dt, K, nx, ny))
    @inbounds result[i,j] += dt * K * ( (h[i + 1,j] + h[i - 1,j] - 2.0 * h[i,j]) / (dx * dx) + (h[i,j + 1] + h[i,j - 1] - 2.0 * h[i,j]) / (dy * dy) )
end

@inline function kernel_fused((i, j), result, (h, u, v), (dx, dy, dt, K, nx, ny))
    result[i,j] = h[i,j]
    # u*du/dx
    uxm = (h[i,j] - h[i - 1,j]) / dx
    uxp = (h[i + 1,j] - h[i,j]) / dx
    am  = 0.5 * min(u[i,j] + u[i + 1,j], 0.0)
    ap  = 0.5 * max(u[i - 1,j] + u[i,j], 0.0)
    advx = ap * uxm + am * uxp

    # v*du/dy
    uym = (h[i,j] - h[i,j - 1]) / dy
    uyp = (h[i,j + 1] - h[i,j]) / dy
    am  = 0.5 * min(v[i,j] + v[i,j + 1], 0.0)
    ap  = 0.5 * max(v[i,j - 1] + v[i,j], 0.0)
    advy = ap * uym + am * uyp

    # adv = 0.5*( (u[i,j]+u[i-1,j])*(u[i,j]-u[i-1,j])/dx + (v[i,j]+v[i,j-1])*(u[i,j]-u[i,j-1])/dy)
    @inbounds result[i,j] -= dt * (advx + advy)
    @inbounds result[i,j] += dt * K * ( (h[i + 1,j] + h[i - 1,j] - 2.0 * h[i,j]) / (dx * dx) + (h[i,j + 1] + h[i,j - 1] - 2.0 * h[i,j]) / (dy * dy) )
end

@inline function kernel_composed((i, j), result, (h, u, v), (dx, dy, dt, K, nx, ny))
    result[i,j] = h[i,j]
    diffusion((i, j), result, (h, u, v), (dx, dy, dt, K, nx, ny))
    advection((i, j), result, (h, u, v), (dx, dy, dt, K, nx, ny))
end

function cuda_wrap_kernel(kernel::Function, result::AbstractArray{Float64}, state, args)
    nx = args[5]
    ny = args[6]

    ix = (blockIdx().x - 1) * blockDim().x  + threadIdx().x
    iy = (blockIdx().y - 1) * blockDim().y  + threadIdx().y
    if (ix > 1 && ix < nx && iy > 1 && iy < ny)
        kernel((ix, iy), result, state, args)
    end
    return nothing
end

function gridloop(kernel::Function, result::CuArray{Float64}, state, args)
    # Grid loop for normal arrays on CPU for CUarrays (CUDA)
    nx = args[5]
    ny = args[6]
    result = reshape(result, (nx, ny))
    h = state[1]
    u = state[2]
    v = state[3]
    h = reshape(h, (nx, ny))
    u = reshape(u, (nx, ny))
    v = reshape(v, (nx, ny))
    state = (h, u, v)

    ths = 256
    bls = Int(ceil(length(result) / ths))
    @cuda threads = ths blocks = bls cuda_wrap_kernel(kernel, result, state, args)
end

function gridloop(kernel::Function, result::Array{Float64}, state, args)
    # Grid loop for normal arrays on CPU for normal arrays
    nx = args[5]
    ny = args[6]
    result = reshape(result, (nx, ny))
    h = state[1]
    u = state[2]
    v = state[3]
    h = reshape(h, (nx, ny))
    u = reshape(u, (nx, ny))
    v = reshape(v, (nx, ny))
    state = (h, u, v)

    for j = 2:ny - 1
        for i = 2:nx - 1
            kernel((i, j), result, state, args)
        end
    end
    result = vec(result)
    h = vec(h)
    u = vec(u)
    v = vec(v)
    end

# function gridloop(kernel::Function, result::AbstractArray{Float64}, state, args)
#     nx = args[5]
#     ny = args[6]
#     result = reshape(result, (nx,ny))
#     h = state[1]
#     u = state[2]
#     v = state[3]
#     h = reshape(h, (nx,ny))
#     u = reshape(u, (nx,ny))
#     v = reshape(v, (nx,ny))
#     state = (h,u,v)
#
#     if result isa CuArray
#         ths = 256
#         bls = Int(ceil(length(result) / ths))
#         @cuda threads=ths blocks=bls cuda_wrap_kernel(kernel, result, args)
#     else
#         for i = 1 : length(result)
#             kernel(i, result, args)
#         end
#     end
#     result = vec(result)
#     h = vec(h)
#     u = vec(u)
#     v = vec(v)
# end

function gridloop(kernel::Function, result::AbstractArray{Float64}, state, args)
    nx = args[5]
    ny = args[6]
    result = reshape(result, (nx, ny))
    h = state[1]
    u = state[2]
    v = state[3]
    h = reshape(h, (nx, ny))
    u = reshape(u, (nx, ny))
    v = reshape(v, (nx, ny))
    state = (h, u, v)

    if result isa CuArray
        ths = 256
        bls = Int(ceil(length(result) / ths))
        @cuda threads = ths blocks = bls cuda_wrap_kernel(kernel, result, state, args)
    else
        for i = 1:length(result)
            kernel(i, result, state, args)
        end
    end
    result = vec(result)
    h = vec(h)
    u = vec(u)
    v = vec(v)
    end