text:
using LinearAlgebra
using HDF5
using Random
#Aの計算
"""
Compute <phi (theta_bra) | phi(theta_ket)>
"""
function overlap(vc::VariationalQuantumCircuit, state0::QulacsQuantumState,
    thetas_left::Vector{Float64}, thetas_right::Vector{Float64})

    circ_tmp = copy(vc)

    # Compute state_left
    update_circuit_param!(circ_tmp, thetas_left)
    state_left = copy(state0)

    update_quantum_state!(circ_tmp, state_left)

    #println("state_left_overlap=", state_left)

    # Compute state_right
    update_circuit_param!(circ_tmp, thetas_right)
    state_right = copy(state0)
    update_quantum_state!(circ_tmp, state_right)
    #println("state_right_overlap=", state_right)

    res = inner_product(state_left, state_right)
    res
end

function compute_B(vc::VariationalQuantumCircuit, state0::QulacsQuantumState, delta_theta=1e-8;
    comm=MPI_COMM_WORLD)
    if is_mpi_on && comm === nothing
        error("comm must be given when mpi is one!")
    end

    num_thetas = num_theta(vc)
    thetas = get_thetas(vc)
    upper_triangle_size = (num_thetas^2 + num_thetas) ÷ 2

    B = zeros(Complex{Float64}, num_thetas, num_thetas)
    local_start, local_size = distribute(upper_triangle_size, mpisize(comm), mpirank(comm))
    idx = 1
    for i in 1:num_thetas
        thetas_i = copy(thetas)
        thetas_i2 = copy(thetas)
        thetas_i[i] += delta_theta
        thetas_i2[i] -= delta_theta
        for j in i:num_thetas
            if local_start <= idx < local_start + local_size
                thetas_j = copy(thetas)
                thetas_j[j] += delta_theta
                thetas_j2 = copy(thetas)
                thetas_j2[j] -= delta_theta
                B[i,j] -= real(
                 ((overlap(vc, state0, thetas_i, thetas)
                - overlap(vc, state0, thetas_i2, thetas))/2*delta_theta)
                *((overlap(vc, state0, thetas, thetas_j)
                - overlap(vc, state0, thetas, thetas_j2))/2*delta_theta)
            )
            B[j, i] = B[i, j]
            end
            idx += 1
        end
    end
    if comm === nothing
      return B
    else
      return Allreduce(B, MPI.SUM, comm)
    end
end

function diff_state(vc::VariationalQuantumCircuit, state0::QulacsQuantumState,thetas, thetas_delta)
    circ_tmp = copy(vc)
    update_circuit_param!(circ_tmp, thetas)
    state_tmp = copy(state0)
    update_quantum_state!(circ_tmp, state_tmp)

    update_circuit_param!(circ_tmp, thetas_delta)
    state_tmp_delta = copy(state0)
    update_quantum_state!(circ_tmp, state_tmp_delta)
    #diff
    res = (get_vector(state_tmp_delta) - get_vector(state_tmp))
    res
end

function compute_B2x(vc::VariationalQuantumCircuit, state0::QulacsQuantumState, delta_theta=1e-8;
    comm=MPI_COMM_WORLD)
    if is_mpi_on && comm === nothing
        error("comm must be given when mpi is one!")
    end

    num_thetas = num_theta(vc)
    thetas = get_thetas(vc)
    upper_triangle_size = (num_thetas^2 + num_thetas) ÷ 2

    B = zeros(Complex{Float64}, num_thetas, num_thetas)
    local_start, local_size = distribute(upper_triangle_size, mpisize(comm), mpirank(comm))
    idx = 1
    for i in 1:num_thetas
        thetas_i = copy(thetas)
        thetas_i2 = copy(thetas)
        thetas_i[i] += delta_theta
        thetas_i2[i] -= delta_theta
        for j in i:num_thetas
            if local_start <= idx < local_start + local_size
                thetas_j = copy(thetas)
                thetas_j[j] += delta_theta
                thetas_j2 = copy(thetas)
                thetas_j2[j] -= delta_theta
                #computeδ|φ>/δθ
                B_right_dif = diff_state(vc, state0, thetas, thetas_i)/delta_theta
                #B_right = inner_product(state, B_right_dif)
                println("B_right_dif_$i$j=",B_right_dif)
                B_left_dif = B_right_dif'
                println("B_left_dif_$i$j=", B_left_dif)
                B[i,j] += real(i*j)
                #println("B[$i $j]=", B[i,j])
            B[j, i] = B[i, j]
            end
            idx += 1
        end
    end
    if comm === nothing
      return B
    else
      return Allreduce(B, MPI.SUM, comm)
    end
end

function compute_B2(vc::VariationalQuantumCircuit, state0::QulacsQuantumState, delta_theta=1e-8;
    comm=MPI_COMM_WORLD)
    if is_mpi_on && comm === nothing
        error("comm must be given when mpi is one!")
    end

    num_thetas = num_theta(vc)
    thetas = get_thetas(vc)
    upper_triangle_size = (num_thetas^2 + num_thetas) ÷ 2

    B = zeros(Complex{Float64}, num_thetas, num_thetas)
    local_start, local_size = distribute(upper_triangle_size, mpisize(comm), mpirank(comm))
    idx = 1
    for i in 1:num_thetas
        thetas_i = copy(thetas)
        thetas_i2 = copy(thetas)
        thetas_i[i] += delta_theta
        thetas_i2[i] -= delta_theta
        for j in i:num_thetas
            if local_start <= idx < local_start + local_size
                thetas_j = copy(thetas)
                thetas_j[j] += delta_theta
                thetas_j2 = copy(thetas)
                thetas_j2[j] -= delta_theta

                B_right = ((overlap(vc, state0, thetas, thetas_j)
                - overlap(vc, state0, thetas, thetas))/delta_theta)
                println("B_right_$i$j=",B_right)
                B_left = B_right'
                println("B_left_$i$j=", B_left)
                B[i,j] += real(B_left*B_right)
                println("B[$i $j]=", B[i,j])
            B[j, i] = B[i, j]
            end
            idx += 1
        end
    end
    if comm === nothing
      return B
    else
      return Allreduce(B, MPI.SUM, comm)
    end
end

function compute_F(vc::VariationalQuantumCircuit, state0::QulacsQuantumState, delta_theta=1e-8;
    comm=MPI_COMM_WORLD)
    if is_mpi_on && comm === nothing
        error("comm must be given when mpi is one!")
    end

    num_thetas = num_theta(vc)
    thetas = get_thetas(vc)
    upper_triangle_size = (num_thetas^2 + num_thetas) ÷ 2

    F = zeros(Complex{Float64}, num_thetas, num_thetas)
    local_start, local_size = distribute(upper_triangle_size, mpisize(comm), mpirank(comm))
    idx = 1
    for i in 1:num_thetas
        thetas_i = copy(thetas)
        thetas_i2 = copy(thetas)
        thetas_i[i] += delta_theta
        thetas_i2[i] -= delta_theta
        for j in i:num_thetas
            if local_start <= idx < local_start + local_size
                thetas_j = copy(thetas)
                thetas_j[j] += delta_theta
                thetas_j2 = copy(thetas)
                thetas_j2[j] -= delta_theta
                F[i, j] = real(
                          overlap(vc, state0, thetas_i, thetas_j)
                        - overlap(vc, state0, thetas_i, thetas, )
                        - overlap(vc, state0, thetas,   thetas_j)
                        + overlap(vc, state0, thetas,   thetas, )
                    )/delta_theta^2
                F[i,j] -= real(
                 ((overlap(vc, state0, thetas_i, thetas)
                - overlap(vc, state0, thetas_i2, thetas))/2*delta_theta)
                *((overlap(vc, state0, thetas, thetas_j)
                - overlap(vc, state0, thetas, thetas_j2))/2*delta_theta)
            )
            #    F[i,j] -= real(
            #     ((overlap(vc, state0, thetas_i, thetas)
            #    - overlap(vc, state0, thetas, thetas))/delta_theta)
            #    *((overlap(vc, state0, thetas, thetas_j)
            #    - overlap(vc, state0, thetas, thetas))/delta_theta)
            #    )
            #    F[i,j] -= real(
            #        overlap(vc, state0, thetas_i, thetas) * overlap(vc, state0, thetas, thetas_j)
            #       -overlap(vc, state0, thetas_i, thetas) * overlap(vc, state0, thetas, thetas)
            #       -overlap(vc, state0, thetas, thetas) * overlap(vc, state0, thetas, thetas_j)
            #       +overlap(vc, state0, thetas, thetas) * overlap(vc, state0, thetas, thetas)
            #       )/delta_theta^2
            F[j, i] = F[i, j]
            end
            idx += 1
            #if MPI_rank == 0
                #println("timing in compute_A: $(1e-9*(t2-t1))")
            #end
        end
    end
    if comm === nothing
      return F
    else
      return Allreduce(F, MPI.SUM, comm)
    end
end

function compute_B3(vc::VariationalQuantumCircuit, state0::QulacsQuantumState, delta_theta=1e-8;
    comm=MPI_COMM_WORLD)
    if is_mpi_on && comm === nothing
        error("comm must be given when mpi is one!")
    end

    num_thetas = num_theta(vc)
    thetas = get_thetas(vc)
    upper_triangle_size = (num_thetas^2 + num_thetas) ÷ 2

    B = zeros(Complex{Float64}, num_thetas, num_thetas)
    local_start, local_size = distribute(upper_triangle_size, mpisize(comm), mpirank(comm))
    idx = 1
    for i in 1:num_thetas
        thetas_i = copy(thetas)
        thetas_i[i] += delta_theta
        for j in i:num_thetas
            if local_start <= idx < local_start + local_size
                thetas_j = copy(thetas)
                thetas_j[j] += delta_theta
                B[i, j] =(
                  (overlap(vc, state0, thetas_i, thetas)
                  - overlap(vc, state0, thetas, thetas)) *
                  (overlap(vc, state0, thetas, thetas_j)
                  - overlap(vc, state0, thetas, thetas))
              /delta_theta^2)
            B[j, i] = B[i, j]
            end
            idx += 1
        end
    end
    if comm === nothing
      return B
    else
      return Allreduce(B, MPI.SUM, comm)
    end
end

function compute_A(vc::VariationalQuantumCircuit, state0::QulacsQuantumState, delta_theta=1e-8;
    comm=MPI_COMM_WORLD)
    if is_mpi_on && comm === nothing
        error("comm must be given when mpi is one!")
    end

    num_thetas = num_theta(vc)
    thetas = get_thetas(vc)
    upper_triangle_size = (num_thetas^2 + num_thetas) ÷ 2

    A = zeros(Complex{Float64}, num_thetas, num_thetas)
    local_start, local_size = distribute(upper_triangle_size, mpisize(comm), mpirank(comm))
    idx = 1
    for i in 1:num_thetas
        thetas_i = copy(thetas)
        thetas_i[i] += delta_theta
        for j in i:num_thetas
            if local_start <= idx < local_start + local_size
                thetas_j = copy(thetas)
                thetas_j[j] += delta_theta
                A[i, j] = real(
                          overlap(vc, state0, thetas_i, thetas_j)
                        - overlap(vc, state0, thetas_i, thetas, )
                        - overlap(vc, state0, thetas,   thetas_j)
                        + overlap(vc, state0, thetas,   thetas, )
                    )/delta_theta^2
                A[j, i] = A[i, j]
            end
            idx += 1
            #if MPI_rank == 0
                #println("timing in compute_A: $(1e-9*(t2-t1))")
            #end
        end
    end
    if comm === nothing
      return A
    else
      return Allreduce(A, MPI.SUM, comm)
    end
end

#Cの計算
"""
Compute <phi (theta_bra) |H| phi(theta_ket)>
"""
function transition(op::OFQubitOperator, vc::VariationalQuantumCircuit, state0::QulacsQuantumState,
    thetas_left::Vector{Float64}, thetas_right::Vector{Float64})

    circ_tmp = copy(vc)

    # Compute state_left
    update_circuit_param!(circ_tmp, thetas_left)
    state_left = copy(state0)
    update_quantum_state!(circ_tmp, state_left)

    # Compute state_right
    update_circuit_param!(circ_tmp, thetas_right)
    state_right = copy(state0)
    update_quantum_state!(circ_tmp, state_right)

    # Compute <state_right|H|state_right>
    res = get_transition_amplitude(op, state_left, state_right)
    res
end

function compute_C(op::OFQubitOperator, vc::VariationalQuantumCircuit,
    state0::QulacsQuantumState, delta_theta=1e-8;comm=MPI_COMM_WORLD)
    if is_mpi_on && comm === nothing
        error("comm must be given whn mpi is one!")
    end
    num_thetas = num_theta(vc)
    thetas = get_thetas(vc)

    C = zeros(Complex{Float64}, num_thetas)
    local_start, local_size = distribute(num_thetas, mpisize(comm), mpirank(comm))
    for i in local_start:local_start+local_size-1
        thetas_i = copy(thetas)
        thetas_i[i] += delta_theta

        thetas_i2 = copy(thetas)
        thetas_i2[i] -= delta_theta

        C[i] = -real(
            transition(op, vc, state0, thetas, thetas_i)
            -transition(op, vc, state0, thetas, thetas_i2)
        )/(2*delta_theta)
    end
    if comm === nothing
      return C
    else
      return Allreduce(C, MPI.SUM, comm)
    end
end

function compute_next_thetas_vqs(op::OFQubitOperator, vc::VariationalQuantumCircuit,
    state0::QulacsQuantumState, dtau; delta_theta=1e-8, comm=MPI_COMM_WORLD, verbose=false)
    thetas_dot = compute_thetadot(op, vc, state0, delta_theta, comm=comm, verbose=verbose)
    thetas_dot = convert(Vector{Float64}, thetas_dot)
    return get_thetas(vc) .+ dtau * thetas_dot
end

"Improved VQS"

function compute_next_thetas_direct(op::OFQubitOperator, vc::VariationalQuantumCircuit,
    state0::QulacsQuantumState, dtau;
    comm=MPI_COMM_WORLD, maxiter=100, gtol=1e-5,verbose=true
    )

    state_tau = copy(state0)
    update_quantum_state!(vc, state_tau)
    Etau = get_expectation_value(op, state_tau)

    function cost(thetas::Vector{Float64})
        vc_ = copy(vc)
        state_ = copy(state0)
        update_circuit_param!(vc_, thetas)
        update_quantum_state!(vc_, state_)
        cost_term1 = dtau * get_transition_amplitude(op, state_, state_tau)
        cost_term2 = (dtau * Etau + 1) * inner_product(state_, state_tau)
        real(cost_term1 - cost_term2)
    end

    cost_history = Float64[] #コスト関数の箱
    thetas_init = get_thetas(vc)
    init_theta_list = thetas_init

    function callback(x)
        push!(cost_history, cost(x))
        if verbose && rank == 0
            println("iter ", length(cost_history), " ", cost_history[end])
        end
    end

    if comm !== nothing
        MPI.Bcast!(thetas_init, 0, comm)
    end

    scipy_opt = pyimport("scipy.optimize")
    method = "BFGS"
    options = Dict("disp" => verbose, "maxiter" => maxiter, "gtol" => gtol)
    opt = scipy_opt.minimize(cost, init_theta_list, method=method,
        callback=callback, jac=generate_numerical_grad(cost),
        options=options)
    return opt["x"]
end

function _expval(vc, state0, thetas, op)
    state = copy(state0)
    vc_ = copy(vc)
    update_circuit_param!(vc_, thetas)
    update_quantum_state!(vc_, state)
    return get_expectation_value(op, state)
end

function compute_next_thetas_safe(op::OFQubitOperator, vc::VariationalQuantumCircuit,
    state0::QulacsQuantumState, dtau, tau::Float64;
    algorithm="direct", comm=MPI_COMM_WORLD, maxiter=100, gtol=1e-5, delta_theta=1e-8,
    verbose=true, max_recursion=10)

    Etau = _expval(vc, state0, get_thetas(vc), op)
    if algorithm == "direct"
        thetas_next = compute_next_thetas_direct(op, vc, state0, dtau,
            comm=comm, maxiter=maxiter, gtol=gtol,verbose=verbose)
    else
        thetas_next = compute_next_thetas_vqs(op, vc, state0, dtau,
            comm=comm, verbose=verbose, delta_theta=delta_theta)
    end
    Etau_next = _expval(vc, state0, thetas_next, op)

    if verbose
       println("dtau: $(dtau), Etau: $(Etau) -> $(Etau_next)")
    end

    if max_recursion == 0 || Etau_next <= Etau
        if verbose && mpirank(comm) == 0
            println("next (recusive) tau point= ", tau+dtau)
            println("Etau=")
        end
        return thetas_next
    end

    # Bad case: we need to decrease dtau
    if verbose
       println("Falling back to recursive model with max_recursion = $(max_recursion)")
    end
    thetas_next1 = compute_next_thetas_safe(op, vc, state0, 0.5*dtau, tau,
        comm=comm, maxiter=maxiter, gtol=gtol,verbose=verbose, max_recursion=max_recursion-1)
    vc_ = copy(vc)
    update_circuit_param!(vc_, thetas_next1)
    return compute_next_thetas_safe(op, vc_, state0, 0.5*dtau, tau+0.5*dtau,
        comm=comm, maxiter=maxiter, gtol=gtol,verbose=verbose, max_recursion=max_recursion-1)
end

function compute_next_thetas_unsafe(op::OFQubitOperator, vc::VariationalQuantumCircuit,
    state0::QulacsQuantumState, dtau;
    algorithm="direct", comm=MPI_COMM_WORLD, maxiter=100, gtol=1e-5, delta_theta=1e-8,
    verbose=true
    )
    if algorithm == "direct"
        thetas_next = compute_next_thetas_direct(op, vc, state0, dtau,
            comm=comm, maxiter=maxiter, gtol=gtol, verbose=verbose)
    else
        thetas_next = compute_next_thetas_vqs(op, vc, state0, dtau,
            comm=comm, verbose=verbose, delta_theta=delta_theta)
    end
    return thetas_next
end

#theta(tau)の微分の計算
"""
Compute thetadot = A^(-1) C
"""

function compute_thetadot(op::OFQubitOperator, vc::VariationalQuantumCircuit,
    state0::QulacsQuantumState,delta_theta=1e-8; comm=MPI_COMM_WORLD, verbose=false)
    if is_mpi_on && comm === nothing
        error("comm must be given when mpi is one!")
    end
    num_thetas = num_theta(vc)
    t1 = time_ns()
    A = compute_A(vc, state0, delta_theta; comm=comm)
    C = compute_C(op, vc, state0, delta_theta)
    t2 = time_ns()

    # 今並列計算を実行している。各プロセスで違う乱数がよばれている可能性や変なバグなどがあるかもしれないのでroot=0の値を書くプロセスにブロードキャストして合わせる必要がある。
    if comm !== nothing
        MPI.Bcast!(A, 0, comm)
        MPI.Bcast!(C, 0, comm)
    end

    t3 = time_ns()
    if verbose && mpirank(comm) == 0
        println("timing in compute_thetadot: $(1e-9*(t2-t1)) sec for A $(1e-9*(t3-t2)) sec for C")
        println("condition number of A =", LinearAlgebra.cond(A))
        U, s, Vt = LinearAlgebra.svd(A)
        println("maximum of singular value of A=", findmax(s)[1])
        println("minimum of singular value of A=", findmin(s)[1])
        for i in eachindex(s)
            println("singular values of A = ", s[i])
        end
    end

    #thetadot, r = LinearAlgebra.LAPACK.gelsy!(A, C, 1e-5)
    thetadot = fit_svd(C, A, 1e-5)
    #thetadot = tikhonov(C, A, 1e-3)
    thetadot
end

"""
Perform imaginary-time evolution.

ham_op:
    Hamiltonian
vc:
    Variational circuit. The current value of variational parameters are
    used as the initial value of the imaginary-time evolution.
state0:
    The initial state to which the Variational circuit is applied to
taus:
    list of imaginary times in ascending order
    The first element must be 0.0.
return:
    list of variational parameters at the given imaginary times.
"""
function imag_time_evolve(ham_op::OFQubitOperator, vc::VariationalQuantumCircuit, state0::QulacsQuantumState,
    taus::Vector{Float64}, delta_theta=1e-8;
    comm=MPI_COMM_WORLD, verbose=false, algorithm::String="direct", tol_dE_dtau=1e-5, recursive=true
    )::Tuple{Vector{Vector{Float64}}, Vector{Float64}}
    if is_mpi_on && comm === nothing
        error("comm must be given when mpi is one!")
    end
    if taus[1] != 0.0
        error("The first element of taus must be 0!")
    end
    thetas_tau = [copy(get_thetas(vc))]
    log_norm_tau = zeros(Float64, length(taus))  #エルミートの期待値だから必ず実数

    stop_imag_time_evol = false
    for i in 1:length(taus)-1
        if verbose && mpirank(comm) == 0
          println("imag_time_evolve: starting step $(i)...")
        end

        #compute expectation value
        Etau = _expval(vc, state0, thetas_tau[end], ham_op)
        if verbose && mpirank(comm) == 0
            println("tau_ite, Etau_ite=, ", taus[i], " ", Etau)
            #println("tau_ite=", taus[i])
        end

        dtau = taus[i+1] - taus[i]
        if dtau < 0.0
            error("taus must be in strictly asecnding order!")
        end
        if recursive
            # Compute theta

            if stop_imag_time_evol
                if verbose && mpirank(comm) == 0
                    println("Skipping imag_time_evolve")
                end
                thetas_next_ = copy(thetas_tau[end])
            else
                vc_ = copy(vc)
                update_circuit_param!(vc_, thetas_tau[i])
                thetas_next_ = compute_next_thetas_safe(ham_op, vc_, state0, dtau, taus[i],
                    algorithm=algorithm, comm=comm, verbose=verbose, delta_theta=delta_theta)
                Etau_next = _expval(vc_, state0, thetas_next_, ham_op)
                if abs(Etau - Etau_next) / dtau < tol_dE_dtau
                    if verbose && mpirank(comm) == 0
                        println("Energy converged!")
                    end
                    stop_imag_time_evol = true
                end
            end
            push!(thetas_tau, thetas_next_)
            # Compute norm
            log_norm_tau[i+1] = log_norm_tau[i] - Etau * (taus[i+1] - taus[i])

        else
            vc_ = copy(vc)
            update_circuit_param!(vc_, thetas_tau[i])
            thetas_next_ = compute_next_thetas_unsafe(ham_op, vc_, state0, dtau,
                algorithm=algorithm, comm=comm, delta_theta=delta_theta, verbose=verbose)
            Etau_next = _expval(vc_, state0, thetas_next_, ham_op)
            push!(thetas_tau, thetas_next_)

            log_norm_tau[i+1] = log_norm_tau[i] - Etau * (taus[i+1] - taus[i])
        end
    end
    thetas_tau, log_norm_tau
end

function _create_quantum_state(c, state0)
    state0_ = copy(state0)
    update_quantum_state!(c, state0_)
    state0_
end

function _create_quantum_state(c, theta::Vector{Float64}, state0::QuantumState)
    c_ = copy(c)
    state0_ = copy(state0)
    update_circuit_param!(c_, theta)
    _create_quantum_state(c_, state0_)
end

"""
Calculate green function based on imaginary-time evolution.

ham_op:
    Hamiltonian
c_op:
    annihilation operator
cdagg_op:
    creation operator
vc:
    Variational circuit. The current value of variational parameters are
    used as the initial value of the imaginary-time evolution.
state_gs:
    The state is the ground state of the hamiltonian
state0_ex:
    The excited state to which the creation operator  is applied to
taus:
    list of imaginary times in ascending order
    The first element must be 0.0. The last element must be beta.
return:
    The list of Green function at each tau
"""

function compute_gtau(
    ham_op::OFQubitOperator,
    left_op::OFQubitOperator,
    right_op::OFQubitOperator,
    vc_ex::VariationalQuantumCircuit,
    state_gs::QulacsQuantumState,
    state0_ex::QulacsQuantumState,
    taus::Vector{Float64}, delta_theta=1e-8;
    comm=MPI_COMM_WORLD, verbose=false, algorithm::String="direct",recursive=true
    )
    if is_mpi_on && comm === nothing
        error("comm must be given when mpi is one!")
    end

    if taus[1] != 0.0
        error("The first element of taus must be 0!")
    end

    if !all(taus[2:end] .> taus[1:end-1])
       error("taus must in strictly asecnding order!")
    end

    if verbose && mpirank(comm) == 0
        println("Applying to an operator to the ket...")
    end
    circuit_right_ex = copy(vc_ex)
    right_squared_norm = apply_qubit_op!(right_op, state_gs,
       circuit_right_ex, state0_ex,
       minimizer=mk_scipy_minimize(
           options = Dict("disp" => verbose),
           verbose=verbose)
        )
    state_right_ex = copy(state0_ex)
    update_quantum_state!(circuit_right_ex, state_right_ex)
    if verbose
        println("Successfully applied an operator to the ket!")
    end

    # exp(-tau H)c^{dag}_j|g.s>
    if algorithm == "vqs"
        thetas_tau_right, log_norm_tau_right = imag_time_evolve(
            ham_op, circuit_right_ex, state0_ex, taus, delta_theta,
            comm=comm, verbose=verbose, algorithm="vqs", recursive=recursive)
    elseif algorithm == "direct"
        thetas_tau_right, log_norm_tau_right = imag_time_evolve(
            ham_op, circuit_right_ex, state0_ex, taus, delta_theta,
            comm=comm, verbose=verbose, algorithm="direct", recursive=recursive)
    end
    #thetas_tau_right, log_norm_tau_right = imag_time_evolve(ham_op, circuit_right_ex, state0_ex, taus, delta_theta, comm=comm)

    Gfunc_ij_list = Complex{Float64}[]
    E_gs = get_expectation_value(ham_op, state_gs)

    norm = Complex{Float64}[]

    for t in eachindex(taus)
        state_right = _create_quantum_state(vc_ex, thetas_tau_right[t], state0_ex)
        state_left = copy(state_gs)
        # Divide the qubit operator of c_i into its real and imaginary parts.
        op_re, op_im = divide_real_imag(left_op)
        g_re = get_transition_amplitude(op_re, state_left, state_right)
        g_im = get_transition_amplitude(op_im, state_left, state_right)
        #基底エネルギーを基準にする。
        push!(norm, exp(log_norm_tau_right[t] + E_gs * taus[t]))
        push!(Gfunc_ij_list, -(g_re + im * g_im) * right_squared_norm * exp(log_norm_tau_right[t] + E_gs *  taus[t]))
    end
    Gfunc_ij_list, norm
end

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
## reverse complement
rcm(M) = reverse(reverse(M, dims=1), dims=2)

### loop for forward and reverse motif scanning
## no doubt there are efficiency savings available here
function scanmotif(seq, mot)
    n = length(seq)
    m = size(mot, 2)
    fscores = zeros(Float64, n - m + 1)
    rscores = zeros(Float64, n - m + 1)
    rot = rcm(mot)
    for i = 1:(n - m + 1)
        for j = 1:m
            if seq[i + j - 1] == DNA_A
                fscores[i] += mot[1, j]
                rscores[i] += rot[1, j]
            elseif seq[i + j - 1] == DNA_C
                fscores[i] += mot[2, j]
                rscores[i] += rot[2, j]
            elseif seq[i + j - 1] == DNA_G
                fscores[i] += mot[3, j]
                rscores[i] += rot[3, j]
            elseif seq[i + j - 1] == DNA_T
                fscores[i] += mot[4, j]
                rscores[i] += rot[4, j]
            end
        end
    end
    fscores, rscores
end

### motif scann stats
function scanmotstats(mot, seq::T, thr=5) where{T}

    fs, rs = scanmotif(seq, mot.pbg)

    ## neg ecdf
    nec = ecdf([-fs ; -rs])

    fsi = findall(fs .> thr)
    rsi = findall(rs .> thr)

    res = DataFrame(Motif=String[], start=Int[], stop=Int[], score=Float64[], strand=String[], emp_p=Float64[], prmax=Float64[], seq=Vector{T}())

    maxscore = sum(maximum(mot.pbg, dims=1))

    n = size(mot.pbg, 2)
    for f in fsi
       push!(res, (mot.name, f, f + n - 1, fs[f], "+", nec(-fs[f]), fs[f]/maxscore, seq[f:f+n-1]))
    end

    for r in rsi
       push!(res, (mot.name, r, r + n - 1, rs[r], "-", nec(-rs[r]), rs[r]/maxscore, reverse_complement(seq[r:r+n-1])))
    end
    res
end

function consensus(mot)
    bases = (DNA_A, DNA_C, DNA_G, DNA_T)
    LongDNASeq(mapreduce(c -> bases[argmax(c)], vcat, eachcol(mot.pwm)))
end

function motifind(motname, mots)
    ms = findall(m -> occursin(motname, m.name), mots)
    if isempty(ms)
       error("$motname not found")
    elseif length(ms) > 1
        error("Multiple matches for $motname\n$([mots[m].name for m in ms])")
    else
       first(ms)
    end

end

function consensus(motname, mots)
    ms = filter(m -> occursin(motname, m.name), mots)
    if isempty(ms)
       error("$motname not found")
    elseif length(ms) > 1
        error("Multiple matches for $motname\n$([m.name for m in ms])")
    else
        consensus(first(ms))
    end
end

ext: jl


text:
# colorFunctions.jl
include("ManagePic.jl")

# auxiliary functions

"""
Get maximum, middle and minimum values of each color.

# Arguments
- `rgb::Array{Array{Float64,2},1}`: r, g, b matrices list.
"""
function getValues(rgb)

    r, g, b = rgb
    Cmax = max.(r, g, b)
    Cmin = min.(r, g, b)
    sum = r .+ g .+ b
    middleValue = sum .- Cmax .- Cmin
    return [Cmax, middleValue, Cmin]
end
"""
Change contrast.

# Arguments
- `rgb::Array{Array{Float64,2},1}`: r, g, b matrices list.
- `parameter::Float64`: in range [-100, 100], negative values create image with lower contrast.
"""
function changeContrast(rgb, parameter)
    if parameter >= -100 && parameter <= 100
        factor = (102 * (parameter + 100)) ./ (100 * (102 - parameter))
        r, g, b = rgb
        arraySize = size(r)
        newRed = zeros(Float64, arraySize[1], arraySize[2])
        newGreen = zeros(Float64, arraySize[1], arraySize[2])
        newBlue = zeros(Float64, arraySize[1], arraySize[2])

        newRed = factor .* (r .- 0.5) .+ 0.5
        newGreen = factor .* (g .- 0.5) .+ 0.5
        newBlue = factor .* (b .- 0.5) .+ 0.5
    else
        println("Error")
    end
    return [newRed, newGreen, newBlue]
end

"""
Change brightness.

# Arguments
- `rgb::Array{Array{Float64,2},1}`: r, g, b matrices list.
- `parameter::Float64`: in range [-100, 100]; negative value create image with lower brightness.
"""
function changeBrightness(rgb, value) # wartosc [-100,100]
    if value <= 100 && value >= -100 && (typeof(value) == Float64 || typeof(value) == Int64)
        point = value / 100
        r, g, b = rgb
        newR = min.(max.(0.5* point .+ r, 0), 1)
        newG =  min.(max.(0.5* point .+ g, 0), 1)
        newB = min.(max.(0.5*point .+ b, 0), 1)
            return [newR, newG,newB]
    else
        return ("Error")
    end
end

"""
Change lightness.

# Arguments
- `rgb::Array{Array{Float64,2},1}`: r, g, b matrices list.
- `parameter::Float64`: in range [-100, 100]; negative value create image with lower lightness.
"""
function changeLightness(rgb, value) # wartosc [-100,100]
    if value <= 100 && value >= -100 && (typeof(value) == Float64 || typeof(value) == Int64)
        point = value / 100

        h, s, l = rgb2hsl(rgb)
        newL = l .+ l .* point
        tryplet = h, s, newL
        newR, newG, newB = hsl2rgb(tryplet)
        return [newR, newG,newB]

    else
        return ("Error")
    end
end

"""
Change saturation.

#Arguments
- `rgb::Array{Array{Float64,2},1}`: r, g, b matrices list.
- `parameter::Float64`: in range [-100, 100], negative value create image with lower saturation.
"""
function changeSaturation(rgb, value)
    if value <= 100 && value >= -100 && (typeof(value) == Float64 || typeof(value) == Int64)
        point = value / 100
        h, s, l = rgb2hsl(rgb)
        if value >= 0
            complementS = 1 .- s
            newS = complementS .* point .+ s
        else
            newS = s .+ s .* point
        end
        tryplet = h, newS, l
        r, g, b = hsl2rgb(tryplet)
        return [r,g,b]
    else
        return ("Error")
    end
end

##################################################################################
# FILTERS
##################################################################################

# GRAYSCALES

"""
Count average value of each pixel.

# Arguments
- `rgb::Array{Array{Float64,2},1}`: r, g, b matrices list.
"""
function average(rgb) # as a average value of pixel

    r, g, b = rgb
    average = (r .+ g .+ b) ./ 3
    return [average, average, average]
end

"""
Get Lightness from HSL and set as a grayscale.

# Arguments
- `rgb::Array{Array{Float64,2},1}`: r, g, b matrices list.
"""
function grayscaleLightness(rgb) # as HSL lightness

    l = lightness(rgb)
    return [l, l, l]
end

"""
Gets grayscale by using some ratio.

# Arguments
- `rgb::Array{Array{Float64,2},1}`: r, g, b matrices list.
"""
function grayscaleLuminosity(rgb) # calculated by ratio

    r, g, b = rgb
    luminosity = 0.25 .* r .+ 0.75 .* g .+ 0.07 .* b
    return [luminosity, luminosity, luminosity]
end

# OTHER

"""
Create a negative of a image.
- `rgb::Array{Array{Float64,2},1}`: r, g, b matrices list.
"""
function negative(rgb)
    r, g, b = rgb
    newRed = 1 .- r
    newGreen = 1 .- g
    newBlue = 1 .- b
    return [newRed, newGreen, newBlue]
end

# Only
function onlyRed(rgb)
    r, g, b = rgb
    arraySize = size(g)
    g = zeros(Float64, arraySize[1], arraySize[2])
    b = zeros(Float64, arraySize[1], arraySize[2])
    return [r, g, b]
end

function onlyGreen(rgb)
    r, g, b = rgb
    arraySize = size(g)
    r = zeros(Float64, arraySize[1], arraySize[2])
    b = zeros(Float64, arraySize[1], arraySize[2])
    return [r, g, b]
end

function onlyBlue(rgb)
    r, g, b = rgb
    arraySize = size(g)
    g = zeros(Float64, arraySize[1], arraySize[2])
    r = zeros(Float64, arraySize[1], arraySize[2])
    return [r, g, b]
end

# Without
function withoutRed(rgb)
    r, g, b = rgb
    arraySize = size(g)
    r = zeros(Float64, arraySize[1], arraySize[2])
    return [r, g, b]
end

function withoutGreen(rgb)
    r, g, b = rgb
    arraySize = size(g)
    g = zeros(Float64, arraySize[1], arraySize[2])
    return [r, g, b]
end

function withoutBlue(rgb)
    r, g, b = rgb
    arraySize = size(b)
    b = zeros(Float64, arraySize[1], arraySize[2])
    return [r, g, b]
end

# AsGray
function redAsAGrayscale(rgb)
    red, green, blue = rgb
    average = (red .+ green .+ blue) ./ 3
    return [average, green, blue]
end

function greenAsAGrayscale(rgb)
    red, green, blue = rgb
    average = (red .+ green .+ blue) ./ 3
    return [red, average, blue]
end

function blueAsAGrayscale(rgb)
    red, green, blue = rgb
    average = (red .+ green .+ blue) ./ 3
    return [red, green, average]
end

function onlyRedAndGrayscale(rgb)
    red, green, blue = rgb
    average = (red .+ green .+ blue) ./ 3
    return [red, average, average]
end

function onlyGreenAndGrayscale(rgb)
    red, green, blue = rgb
    average = (red .+ green .+ blue) ./ 3
    return [average, green, average]
end

function onlyBlueAndGrayscale(rgb)
    red, green, blue = rgb
    average = (red .+ green .+ blue) ./ 3
    return [average, average, blue]
end

# Double rainbow? what is tHAT? :o

function changeColors(rgb, point) # point is in range [0,360]
    if point <= 360 && point >= 0
        h, s, l = rgb2hsl(rgb)
        Hue = h .+ point
        trypletHSL = Hue, s, l
        newRGB = hsl2rgb(trypletHSL)
        return newRGB
    end
end

max_stars_repo_name: o-Mateo-o/ImageEditorApp

max_issues_repo_name: o-Mateo-o/ImageEditorApp

max_forks_repo_name: o-Mateo-o/ImageEditorApp


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

using Markdown
using InteractiveUtils

# ╔═╡ Cell order:


text:
@testset "SLIC" begin
  Z = [ones(10,10) 2ones(10,10); 3ones(10,10) 4ones(10,10)]
  𝒮 = georef((Z=Z,))
  p = partition(𝒮, SLIC(4, 1.0))
  @test length(p) == 4
  @test all(nelements.(p) .== 100)
  p1, p2, p3, p4 = p
  @test mean(coordinates(centroid(p1, ind)) for ind in 1:nelements(p1)) == [5.0,5.0]
  @test mean(coordinates(centroid(p2, ind)) for ind in 1:nelements(p2)) == [15.0,5.0]
  @test mean(coordinates(centroid(p3, ind)) for ind in 1:nelements(p3)) == [5.0,15.0]
  @test mean(coordinates(centroid(p4, ind)) for ind in 1:nelements(p4)) == [15.0,15.0]

  C = cluster(𝒮, SLIC(4, 1.0))
  @test C.cluster == vec(Z')

  𝒮 = georef((z=[√(i^2+j^2) for i in 1:100, j in 1:100],))
  p = partition(𝒮, SLIC(50, 0.001))
  @test 50 ≤ length(p) ≤ 60

  if visualtests
    @test_reference "data/slic.png" plot(p)
  end

  # test SLIC with heterogeneous data
  Z = (a=rand(10), b=1:10, x=rand(10), y=rand(10))
  𝒮 = georef(Z, (:x, :y))
  C = cluster(𝒮, SLIC(2, 1.0))
  @test domain(C) == domain(𝒮)
  @test Set(C.cluster) ⊆ Set(1:2)

  # test SLIC for orphaned points
  a = [0,0,0,0,0,0,0,0,0,0]
  x = [0.4993029939801461, 0.14954882636793432, 0.23118957975519616, 0.6816610871344635, 0.6665309965318731, 0.691522274292691, 0.012495903053589608, 0.9831177095525963, 0.4445263730141056, 0.2175871587746574]
  y = [0.32721108209880256, 0.11427387079564899, 0.826401075107011, 0.6164294766961782, 0.6562529361193601, 0.43388375115444644, 0.7624847842129086, 0.1516623758764959, 0.07641616063237144, 0.8669098569279463]
  Z = (a=a, x=x, y=y)
  𝒮 = georef(Z, (:x, :y))
  C = cluster(𝒮, SLIC(2, 1.0))
  @test Set(C.cluster) ⊆ Set(1:2)

  # test SLIC with weights in attribute columns
  z1 = [√((i-0)^2+(j-0)^2) for i in 1:100, j in 1:100]
  z2 = [√((i-100)^2+(j-100)^2) for i in 1:100, j in 1:100]
  𝒮 = georef((z1=z1, z2=z2))
  w1 = Dict(:z1 => 10, :z2 => 0.1)
  w2 = Dict(:z1 => 0.1, :z2 => 10)
  p1 = partition(𝒮, SLIC(50, 0.001, weights=w1))
  p2 = partition(𝒮, SLIC(50, 0.001, weights=w2))
  @test 50 ≤ length(p1) ≤ 60
  @test 50 ≤ length(p2) ≤ 60

  if visualtests
    @test_reference "data/slic-w1.png" plot(p1)
    @test_reference "data/slic-w2.png" plot(p2)
  end

  # test GeoClustering.slic_srecursion function
  k = 20
  l = [10.0, 100.0, 1000.0]
  s = GeoClustering.slic_srecursion(k, l)
  @test s[1] == 10/3 && s[2] == 100/3 && s[3] == 1000/3

  # the following test deals with the case where the bounding box
  # of the data has very different sides, one of which is too small
  # we want to make sure that the initialization of centroids always
  # returns a non-empty set
  k = 1
  m = 0.000001
  x = LinRange(550350.6224548942, 552307.2106300013, 1200)
  y = LinRange(9.35909841165263e6, 9.36050447440832e6, 1200)
  z = LinRange(-44.90690201082941, 351.4007207008662, 1200)
  Z = (x = x, y = y, z = z, a = rand(1200))
  𝒮 = georef(Z, (:x, :y, :z))
  s = GeoClustering.slic_spacing(𝒮, SLIC(k, m))
  bbox = boundingbox(𝒮)
  lo, up = coordinates.(extrema(bbox))
  ranges = [(l+sᵢ/2):sᵢ:u for (l, sᵢ, u) in zip(lo, s, up)]
  @test !isempty(Iterators.product(ranges...))
  c = GeoClustering.slic_initialization(𝒮, s)
  @test !isempty(c)

  # visual SLIC test for the μCT image
  k = 45
  m = 0.55
  μCT = load(joinpath(datadir,"muCT.tif"))
  𝒮 = georef((μCT = Float64.(μCT),))
  C = cluster(𝒮, SLIC(45, 0.55))

  if visualtests
    @test_reference "data/slic-muCT.png" plot(C)
  end
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct ConnectionMonitorSource <: SwaggerModel
    resourceId::Any # spec type: Union{ Nothing, String } # spec name: resourceId
    port::Any # spec type: Union{ Nothing, Int32 } # spec name: port

    function ConnectionMonitorSource(;resourceId=nothing, port=nothing)
        o = new()
        validate_property(ConnectionMonitorSource, Symbol("resourceId"), resourceId)
        setfield!(o, Symbol("resourceId"), resourceId)
        validate_property(ConnectionMonitorSource, Symbol("port"), port)
        setfield!(o, Symbol("port"), port)
        o
    end
end # type ConnectionMonitorSource

const _property_map_ConnectionMonitorSource = Dict{Symbol,Symbol}(Symbol("resourceId")=>Symbol("resourceId"), Symbol("port")=>Symbol("port"))
const _property_types_ConnectionMonitorSource = Dict{Symbol,String}(Symbol("resourceId")=>"String", Symbol("port")=>"Int32")
Base.propertynames(::Type{ ConnectionMonitorSource }) = collect(keys(_property_map_ConnectionMonitorSource))
Swagger.property_type(::Type{ ConnectionMonitorSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ConnectionMonitorSource[name]))}
Swagger.field_name(::Type{ ConnectionMonitorSource }, property_name::Symbol) =  _property_map_ConnectionMonitorSource[property_name]

function check_required(o::ConnectionMonitorSource)
    (getproperty(o, Symbol("resourceId")) === nothing) && (return false)
    true
end

function validate_property(::Type{ ConnectionMonitorSource }, name::Symbol, val)
end


text:
@testset "serialization" begin
    function test_handler(session, request)
        obs1 = Observable(Float16(2.0))
        obs2 = Observable(DOM.div("Data: ", obs1, dataTestId="hey"))
        return DOM.div(obs2)
    end
    xx = "hey"; some_js = js"var"; x = [1f0, 2f0]
    @test string(js"console.log($xx); $x; $((2, 4)); $(some_js) hello = 1;") == "console.log(\"hey\"); [1.0,2.0]; [2,4]; var hello = 1;"
    test_throw() = sprint(io->JSServe.serialize_readable(io, JSServe.Asset("file.dun_exist")))
    Test.@test_throws ErrorException("Unrecognized asset media type: dun_exist") test_throw()
    testsession(test_handler) do app
        hey = query_testid("hey")
        @test evaljs(app, js"$(hey).innerText") == "Data: Float16(2.0)"
        float16_obs = children(children(app.dom)[1][])[2]
        float16_obs[] = Float16(77)
        @test evaljs(app, js"$(hey).innerText") == "Data: Float16(77.0)"
    end
end

@testset "http" begin
    @test_throws ErrorException("Invalid sessionid: lol") JSServe.request_to_sessionid((target="lol",))
    @test JSServe.request_to_sessionid((target="lol",), throw=false) == nothing
end

@testset "hyperscript" begin
    function handler(session, request)
        the_script = DOM.script("window.testglobal = 42")
        s1 = Hyperscript.Style(css("p", fontWeight="bold"), css("span", color="red"))
        the_style = DOM.style(Hyperscript.styles(s1))
        return DOM.div(:hello, the_style, the_script, dataTestId="hello")
    end
    testsession(handler) do app
        @test evaljs(app, js"window.testglobal")  == 42
        hello_div = query_testid("hello")
        @test evaljs(app, js"$(hello_div).innerText")  == "hello"
        @test evaljs(app, js"$(hello_div).children.length") == 3
        @test evaljs(app, js"$(hello_div).children[0].tagName") == "P"
        @test evaljs(app, js"$(hello_div).children[1].tagName") == "STYLE"
        @test evaljs(app, js"$(hello_div).children[2].tagName") == "SCRIPT"
    end
end


text:
@testset "Cholesky" begin
    import DiffLinearAlgebra: level2partition, level3partition, chol_unblocked_rev,
        chol_blocked_rev

    let rng = MersenneTwister(123456), N = 5
        A = randn(rng, N, N)
        r, d, B2, c = level2partition(A, 4, false)
        R, D, B3, C = level3partition(A, 4, 4, false)
        @test all(r .== R.')
        @test all(d .== D)
        @test B2[1] == B3[1]
        @test all(c .== C)

        # Check that level2partition with 'U' is consistent with 'L'.
        rᵀ, dᵀ, B2ᵀ, cᵀ = level2partition(transpose(A), 4, true)
        @test r == rᵀ
        @test d == dᵀ
        @test B2.' == B2ᵀ
        @test c == cᵀ

        # Check that level3partition with 'U' is consistent with 'L'.
        R, D, B3, C = level3partition(A, 2, 4, false)
        Rᵀ, Dᵀ, B3ᵀ, Cᵀ = level3partition(transpose(A), 2, 4, true)
        @test transpose(R) == Rᵀ
        @test transpose(D) == Dᵀ
        @test transpose(B3) == B3ᵀ
        @test transpose(C) == Cᵀ
    end

    let rng = MersenneTwister(123456), N = 10
        A, Ā = full.(LowerTriangular.(randn.(rng, [N, N], [N, N])))
        B, B̄ = transpose.([A, Ā])
        @test chol_unblocked_rev(Ā, A, false) ≈ chol_blocked_rev(Ā, A, 1, false)
        @test chol_unblocked_rev(Ā, A, false) ≈ chol_blocked_rev(Ā, A, 3, false)
        @test chol_unblocked_rev(Ā, A, false) ≈ chol_blocked_rev(Ā, A, 5, false)
        @test chol_unblocked_rev(Ā, A, false) ≈ chol_blocked_rev(Ā, A, 10, false)
        @test chol_unblocked_rev(Ā, A, false) ≈ transpose(chol_unblocked_rev(B̄, B, true))

        @test chol_unblocked_rev(B̄, B, true) ≈ chol_blocked_rev(B̄, B, 1, true)
        @test chol_unblocked_rev(B̄, B, true) ≈ chol_blocked_rev(B̄, B, 5, true)
        @test chol_unblocked_rev(B̄, B, true) ≈ chol_blocked_rev(B̄, B, 10, true)
    end

    # Check sensitivities for lower-triangular version.
    let P = 15, rng = MersenneTwister(123456), N = 10

        Σ_gen = ()->(A = randn(rng, P, P); A'A + 1e-3I)
        S_gen, H_gen = ()->Symmetric(Σ_gen()), ()->Hermitian(Σ_gen())
        U_gen = ()->chol(S_gen())

        _chol = Σ->chol(Symmetric(Σ))
        _∇chol = (Y, Ȳ, Σ)->∇(chol, Val{1}, (), Y, Ȳ, Symmetric(Σ))

        @test check_errs(N, _chol, _∇chol, U_gen, Σ_gen, Σ_gen)
        @test check_errs(N, _chol, _∇chol, U_gen, H_gen, H_gen)
        @test check_errs(N, _chol, _∇chol, U_gen, S_gen, S_gen)
    end
end


text:
using ComponentArrays
using DiffEqFlux, DifferentialEquations, Plots
using UnPack

"""
    To test whether `ComponentArray` (exported from `ComponentArrays`) is compatible with AD (auto-diff) systems.

Borrowed from DiffEqFlux.jl documentation:
https://diffeqflux.sciml.ai/stable/examples/neural_ode_flux/
"""
function main()
    # u0 = Float32[2.; 0.]
    u0 = ComponentArray(a=2.0, b=0.0)
    datasize = 30
    tspan = (0.0f0,1.5f0)

    function trueODEfunc(du,u,p,t)
        true_A = [-0.1 2.0; -2.0 -0.1]
        @unpack a, b = u
        tmp = ([a^3, b^3]'true_A)'
        du.a, du.b = tmp
    end
    t = range(tspan[1],tspan[2],length=datasize)
    prob = ODEProblem(trueODEfunc,u0,tspan)
    ode_data = Array(solve(prob,Tsit5(),saveat=t))

    dudt2 = Chain(x -> x.^3,
                  Dense(2,50,tanh),
                  Dense(50,2))
    p,re = Flux.destructure(dudt2) # use this p as the initial condition!
    dudt(u,p,t) = re(p)(u) # need to restructure for backprop!
    prob = ODEProblem(dudt,u0,tspan)

    function predict_n_ode()
        Array(solve(prob,Tsit5(),u0=u0,p=p,saveat=t))
    end

    function loss_n_ode()
        pred = predict_n_ode()
        loss = sum(abs2,ode_data .- pred)
        loss
    end

    loss_n_ode() # n_ode.p stores the initial parameters of the neural ODE

    cb = function (;doplot=false) # callback function to observe training
        pred = predict_n_ode()
        display(sum(abs2,ode_data .- pred))
        # plot current prediction against data
        pl = scatter(t,ode_data[1,:],label="data")
        scatter!(pl,t,pred[1,:],label="prediction")
        display(plot(pl))
        return false
    end

    # Display the ODE with the initial parameter values.
    cb()

    data = Iterators.repeated((), 1000)
    Flux.train!(loss_n_ode, Flux.params(u0,p), data, ADAM(0.05), cb = cb)
end


text:
"""
    is_closed_downwards(X::AbstractMatrix)

Interpret each row of the `N x D` matrix `X` as a vector-valued observation.
"""
is_closed_downwards(X::AbstractMatrix{<:Real}) = true
function is_closed_downwards(X::AbstractMatrix{<:Union{<:Real, Missing}})
    for d in 2:size(X, 2), n in 1:size(X, 1)
        X[n, d-1] isa Missing && !(X[n, d] isa Missing) && return false
    end
    return true
end


text:
include("../source/Ledger.jl")
using .Ledger
using Dates

dummy_transaction = Transaction(1000, "sender", 2000, "recepient", 10.20, "thanks", "mysignature")

first_block = Block(42, DateTime(0), Transaction[dummy_transaction, dummy_transaction], "previous_dude_made_some_string", 69)

println(first_block)

new_blockchain = Blockchain("dude")

old_blockchin = Blockchain("yomama", Dates.now(), Dates.now(), Block[first_block, first_block], 5, 3, 42)

max_stars_repo_name: marchunter/litcoin

max_issues_repo_name: marchunter/litcoin

max_forks_repo_name: marchunter/litcoin


text:
using StatsBase: Histogram
using AbstractPlotting, StatsBase
import Distributions
using KernelDensity

using Random: seed!
using GeometryBasics: FRect2D

seed!(0)

@testset "histogram" begin
    v = randn(1000)
    h = fit(Histogram, v)
    fig, ax, plt = plot(h)

    @test plt isa BarPlot
    x = h.edges[1]
    @test plt[1][] ≈ Point{2, Float32}.(x[1:end-1] .+ step(x)/2, h.weights)

    v = (randn(1000), randn(1000))
    h = fit(Histogram, v, nbins = 30)
    fig, ax, plt = plot(h)
    @test plt isa Heatmap
    x = h.edges[1]
    y = h.edges[2]
    @test plt[1][] ≈ x[1:end-1] .+ step(x)/2
    @test plt[2][] ≈ y[1:end-1] .+ step(y)/2
    @test plt[3][] ≈ h.weights

    fig, ax, plt = surface(h)
    @test plt isa Surface
    x = h.edges[1]
    y = h.edges[2]
    @test plt[1][] ≈ x[1:end-1] .+ step(x)/2
    @test plt[2][] ≈ y[1:end-1] .+ step(y)/2
    @test plt[3][] ≈ h.weights

    v = (randn(1000), randn(1000), randn(1000))
    edges = ntuple(_ -> -3:0.3:3, 3)
    h = fit(Histogram, v, edges)
    fig, ax, plt = plot(h)
    @test plt isa Volume
    x = h.edges[1]
    y = h.edges[2]
    z = h.edges[3]
    @test plt[1][] ≈ x[1:end-1] .+ step(x)/2
    @test plt[2][] ≈ y[1:end-1] .+ step(y)/2
    @test plt[3][] ≈ z[1:end-1] .+ step(z)/2
    @test plt[4][] == h.weights
end

@testset "density" begin
    v = randn(1000)
    d = kde(v)
    fig, ax, p1 = plot(d)
    @test p1 isa Lines
    fig, ax, p2 = lines(d.x, d.density)
    @test p1[1][] == p2[1][]

    x = randn(1000)
    y = randn(1000)
    d = kde((x, y))
    fig, ax, p1 = plot(d)
    @test p1 isa Heatmap
    fig, ax, p2 = heatmap(d.x, d.y, d.density)
    @test p1[1][] == p2[1][]
    @test p1[2][] == p2[2][]
    @test p1[3][] == p2[3][]

    fig, ax, p1 = surface(d)
    @test p1 isa Surface
    fig, ax, p2 = surface(d.x, d.y, d.density)
    @test p1[1][] == p2[1][]
    @test p1[2][] == p2[2][]
    @test p1[3][] == p2[3][]
end

@testset "distribution" begin
    d = Distributions.Normal()
    rg = AbstractPlotting.support(d)
    @test minimum(rg) ≈ -3.7190164854556866
    @test maximum(rg) ≈ 3.719016485455714
    fix, ax, plt = plot(d)
    @test plt isa Lines
    @test !AbstractPlotting.isdiscrete(d)
    @test first(plt[1][][1]) ≈ minimum(rg) rtol = 1f-6
    @test first(plt[1][][end]) ≈ maximum(rg) rtol = 1f-6

    for (x, pd) in plt[1][]
        @test pd ≈ Distributions.pdf(d, x) rtol = 1f-6
    end

    d = Distributions.Poisson()
    rg = AbstractPlotting.support(d)
    @test rg == 0:6
    fig, ax, p = plot(d)
    @test p isa ScatterLines
    plt = p.plots[1]
    @test AbstractPlotting.isdiscrete(d)

    @test first.(plt[1][]) == 0:6
    @test last.(plt[1][]) ≈ Distributions.pdf.(d, first.(plt[1][]))
end

@testset "qqplot" begin
    v = randn(1000)
    q = Distributions.qqbuild(fit(Distributions.Normal, v), v)
    fig, ax, p = qqnorm(v)

    @test length(p.plots) == 2
    plt = p.plots[1]
    @test plt isa Scatter
    @test first.(plt[1][]) ≈ q.qx rtol = 1e-6
    @test last.(plt[1][]) ≈ q.qy rtol = 1e-6

    plt = p.plots[2]
    @test plt isa LineSegments
    @test first.(plt[1][]) ≈ [extrema(q.qx)...] rtol = 1e-6
    @test last.(plt[1][]) ≈ [extrema(q.qx)...] rtol = 1e-6

    fig, ax, p = qqnorm(v, qqline = nothing)
    @test length(p.plots) == 1
    plt = p.plots[1]
    @test plt isa Scatter
    @test first.(plt[1][]) ≈ q.qx rtol = 1e-6
    @test last.(plt[1][]) ≈ q.qy rtol = 1e-6

    fig, ax, p = qqnorm(v, qqline = :fit)
    plt = p.plots[2]
    itc, slp = hcat(fill!(similar(q.qx), 1), q.qx) \ q.qy
    xs = [extrema(q.qx)...]
    ys = slp .* xs .+ itc
    @test first.(plt[1][]) ≈ xs rtol = 1e-6
    @test last.(plt[1][]) ≈ ys rtol = 1e-6

    fig, ax, p = qqnorm(v, qqline = :quantile)
    plt = p.plots[2]
    xs = [extrema(q.qx)...]
    quantx, quanty = quantile(q.qx, [0.25, 0.75]), quantile(q.qy, [0.25, 0.75])
    slp = diff(quanty) ./ diff(quantx)
    ys = quanty .+ slp .* (xs .- quantx)
    @test first.(plt[1][]) ≈ xs rtol = 1e-6
    @test last.(plt[1][]) ≈ ys rtol = 1e-6
end

@testset "crossbar" begin
    fig, ax, p = crossbar(1, 3, 2, 4)
    @test p isa CrossBar
    @test p.plots[1] isa Poly
    @test p.plots[1][1][] == [FRect2D(Float32[0.6, 2.0], Float32[0.8, 2.0]),]
    @test p.plots[2] isa LineSegments
    @test p.plots[2][1][] == Point{2,Float32}[Float32[0.6, 3.0], Float32[1.4, 3.0]]

    fig, ax, p = crossbar(1, 3, 2, 4; show_notch = true, notchmin = 2.5, notchmax = 3.5);
    @test p isa CrossBar
    @test p.plots[1] isa Poly
    @test p.plots[1][1][][1] isa AbstractPlotting.AbstractMesh
    poly = Point{2,Float32}[[0.6, 2.0], [1.4, 2.0], [1.4, 2.5], [1.2, 3.0], [1.4, 3.5],
                            [1.4, 4.0], [0.6, 4.0], [0.6, 3.5], [0.8, 3.0], [0.6, 2.5]]
    @test map(Point2f0, p.plots[1][1][][1].position) == poly
    @test p.plots[2] isa LineSegments
    @test p.plots[2][1][] == Point{2,Float32}[Float32[0.8, 3.0], Float32[1.2, 3.0]]
end

@testset "boxplot" begin
    a = repeat(1:5, inner = 20)
    b = 1:100
    fig, ax, p = boxplot(a, b)
    plts = p.plots
    @test length(plts) == 3
    @test plts[1] isa Scatter
    @test isempty(plts[1][1][])

    # test categorical
    a = repeat(["a", "b", "c", "d", "e"], inner = 20)
    b = 1:100
    fig, ax, p = boxplot(a, b; whiskerwidth = 1.0)
    plts = p.plots
    @test length(plts) == 3
    @test plts[1] isa Scatter
    @test isempty(plts[1][1][])

    @test plts[2] isa LineSegments
    pts = Point{2, Float32}[
        [1.0, 5.75], [1.0, 1.0], [0.6, 1.0], [1.4, 1.0], [1.0, 15.25],
        [1.0, 20.0], [1.4, 20.0], [0.6, 20.0], [2.0, 25.75], [2.0, 21.0],
        [1.6, 21.0], [2.4, 21.0], [2.0, 35.25], [2.0, 40.0], [2.4, 40.0],
        [1.6, 40.0], [3.0, 45.75], [3.0, 41.0], [2.6, 41.0], [3.4, 41.0],
        [3.0, 55.25], [3.0, 60.0], [3.4, 60.0], [2.6, 60.0], [4.0, 65.75],
        [4.0, 61.0], [3.6, 61.0], [4.4, 61.0], [4.0, 75.25], [4.0, 80.0],
        [4.4, 80.0], [3.6, 80.0], [5.0, 85.75], [5.0, 81.0], [4.6, 81.0],
        [5.4, 81.0], [5.0, 95.25], [5.0, 100.0], [5.4, 100.0], [4.6, 100.0]
    ]
    @test plts[2][1][] == pts

    @test plts[3] isa CrossBar
    @test plts[3].plots[1] isa Poly

    poly = [
        FRect2D(Float32[0.6, 5.75], Float32[0.8, 9.5]),
        FRect2D(Float32[1.6, 25.75], Float32[0.8, 9.5]),
        FRect2D(Float32[2.6, 45.75], Float32[0.8, 9.5]),
        FRect2D(Float32[3.6, 65.75], Float32[0.8, 9.5]),
        FRect2D(Float32[4.6, 85.75], Float32[0.8, 9.5]),
    ]

    @test plts[3].plots[1][1][] == poly

    #notch
    fig, ax, p = boxplot(a, b, show_notch=true)
    plts = p.plots

    @test length(plts) == 3

    pts = Point{2,Float32}[
        [1.0, 5.75], [1.0, 1.0], [1.0, 1.0], [1.0, 1.0], [1.0, 15.25],
        [1.0, 20.0], [1.0, 20.0], [1.0, 20.0], [2.0, 25.75], [2.0, 21.0],
        [2.0, 21.0], [2.0, 21.0], [2.0, 35.25], [2.0, 40.0], [2.0, 40.0],
        [2.0, 40.0], [3.0, 45.75], [3.0, 41.0], [3.0, 41.0], [3.0, 41.0],
        [3.0, 55.25], [3.0, 60.0], [3.0, 60.0], [3.0, 60.0], [4.0, 65.75],
        [4.0, 61.0], [4.0, 61.0], [4.0, 61.0], [4.0, 75.25], [4.0, 80.0],
        [4.0, 80.0], [4.0, 80.0], [5.0, 85.75], [5.0, 81.0], [5.0, 81.0],
        [5.0, 81.0], [5.0, 95.25], [5.0, 100.0], [5.0, 100.0], [5.0, 100.0],
    ]

    @test plts[2] isa LineSegments
    @test plts[2][1][] == pts

    @test plts[3] isa CrossBar
    @test plts[3].plots[1] isa Poly

    notch_boxes = Vector{Point{2,Float32}}[map(Point2f0, [[0.6, 5.75], [1.4, 5.75], [1.4, 7.14366], [1.2, 10.5], [1.4, 13.8563], [1.4, 15.25], [0.6, 15.25], [0.6, 13.8563], [0.8, 10.5], [0.6, 7.14366]]),
                                           map(Point2f0, [[1.6, 25.75], [2.4, 25.75], [2.4, 27.1437], [2.2, 30.5], [2.4, 33.8563], [2.4, 35.25], [1.6, 35.25], [1.6, 33.8563], [1.8, 30.5], [1.6, 27.1437]]),
                                           map(Point2f0, [[2.6, 45.75], [3.4, 45.75], [3.4, 47.1437], [3.2, 50.5], [3.4, 53.8563], [3.4, 55.25], [2.6, 55.25], [2.6, 53.8563], [2.8, 50.5], [2.6, 47.1437]]),
                                           map(Point2f0, [[3.6, 65.75], [4.4, 65.75], [4.4, 67.1437], [4.2, 70.5], [4.4, 73.8563], [4.4, 75.25], [3.6, 75.25], [3.6, 73.8563], [3.8, 70.5], [3.6, 67.1437]]),
                                           map(Point2f0, [[4.6, 85.75], [5.4, 85.75], [5.4, 87.1437], [5.2, 90.5], [5.4, 93.8563], [5.4, 95.25], [4.6, 95.25], [4.6, 93.8563], [4.8, 90.5], [4.6, 87.1437]])]
    meshes = plts[3].plots[1][1][]
    @testset for (i, mesh) in enumerate(meshes)
        @test mesh isa AbstractPlotting.AbstractMesh
        vertices = map(Point2f0, mesh.position)
        @test vertices ≈ notch_boxes[i]
    end
end

@testset "violin" begin
    x = repeat(1:4, 250)
    y = x .+ randn.()
    fig, ax, p = violin(x, y, side = :left, color = :blue)
    @test p isa Violin
    @test p.plots[1] isa Poly
    @test p.plots[1][:color][] == :blue
    @test p.plots[2] isa LineSegments
    @test p.plots[2][:color][] == :white
    @test p.plots[2][:visible][] == :false

    # test categorical
    x = repeat(["a", "b", "c", "d"], 250)
    fig2, ax2, p2 = violin(x, y, side = :left, color = :blue)
    @test p2 isa Violin
    @test p2.plots[1] isa Poly
    @test p2.plots[1][:color][] == :blue
    @test p2.plots[2] isa LineSegments
    @test p2.plots[2][:color][] == :white
    @test p2.plots[2][:visible][] == :false
end


text:
# Bogumił Kamiński, 2019-03-25

using LinearAlgebra

function uppersum(M)
    n = size(M, 1)
    s = 0
    for i in 2:n, j in 1:i
        s += M[j, i]
    end
    s
end

@time m = randn(10_000, 10_000);
@time m = randn(10_000, 10_000);

m

@time sum(m)
@time sum(m)

@time sum(UpperTriangular(m))
@time sum(UpperTriangular(m))

@time uppersum(m)
@time uppersum(m)

max_stars_repo_path: 3. Uppersum/uppersum_1.jl

max_issues_repo_path: 3. Uppersum/uppersum_1.jl

max_forks_repo_path: 3. Uppersum/uppersum_1.jl


text:
using Rayuela
using Test

# Common utility functions
include("common.jl")

# IO - fvecs and ivecs read/write
include("xvecs.jl")

# Chain quantization
include("chainq.jl")


text:
# Copyright (c) 2020: Matthew Wilhelm & Matthew Stuber.
# This work is licensed under the Creative Commons Attribution-NonCommercial-
# ShareAlike 4.0 International License. To view a copy of this license, visit
# http://creativecommons.org/licenses/by-nc-sa/4.0/ or send a letter to Creative
# Commons, PO Box 1866, Mountain View, CA 94042, USA.
#############################################################################
# Dynamic Bounds Base
# The underlying abstraction layer for DynamicBounds.
# See https://github.com/PSORLab/DynamicBoundsBase.jl
#############################################################################
# src/attribute_handlers.jl
# Defines API for accessing attributes of the integrators and problems.
#############################################################################

const AnyDEAttribute = Union{AbstractIntegratorAttribute,
                             AbstractRelaxProblemAttribute}

"""
$(FUNCTIONNAME)

Return a `Bool` indicating whether `::AbstractDERelaxIntegrator` supports
the `::Union{AbstractDERelaxIntegrator, AbstractDERelaxProblem}`.
"""
function supports end

supports(::AbstractDERelaxIntegrator, ::AbstractIntegratorAttribute) = false
supports(::AbstractDERelaxProblem, ::AbstractRelaxProblemAttribute) = false

"""
$(FUNCTIONNAME)

Returns the value of attribute set for the `::Union{AbstractDERelaxIntegrator, AbstractDERelaxProblem}`. May be vector
valued.
"""
function get end
function get(integrator::AbstractDERelaxIntegrator, attr::Vector{AbstractIntegratorAttribute}, idxs::Vector)
	[get(integrator, i) for i in idxs]
end
function get(problem::AbstractDERelaxProblem, attr::Vector{AbstractRelaxProblemAttribute}, idxs::Vector)
	[get(problem, i) for i in idxs]
end
function get(m::AbstractDERelaxIntegrator, attr::AbstractIntegratorAttribute)
    throw(ArgumentError("AbstractDERelaxIntegrator of type $(typeof(m)) does not support accessing the attribute $attr via get"))
end
function get(m::AbstractDERelaxProblem, attr::AbstractRelaxProblemAttribute)
    throw(ArgumentError("AbstractDERelaxProblem of type $(typeof(m)) does not support accessing the attribute $attr via get"))
end

"""
$(FUNCTIONNAME)

An version of `get` which retreives all data asssociated with a particular integrator attribute.
"""
function getall end
function getall(m::AbstractDERelaxIntegrator, attr::AbstractIntegratorAttribute, args...)
    throw(ArgumentError("AbstractDERelaxIntegrator of type $(typeof(m)) does not support accessing the attribute $attr via getall"))
end
function getall(m::AbstractDERelaxProblem, attr::AbstractRelaxProblemAttribute, args...)
    throw(ArgumentError("AbstractDERelaxProblem of type $(typeof(m)) does not support accessing the attribute $attr via getall"))
end

"""
$(FUNCTIONNAME)

An in-place version of `get!` which retreives all data asssociated with a particular integrator attribute.
"""
function getall! end
function getall!(output, m::AbstractDERelaxIntegrator, attr::AbstractIntegratorAttribute, args...)
    throw(ArgumentError("AbstractDERelaxIntegrator of type $(typeof(m)) does not support accessing the attribute $attr via getall!"))
end
function getall!(output, m::AbstractDERelaxProblem, attr::AbstractRelaxProblemAttribute, args...)
    throw(ArgumentError("AbstractDERelaxProblem of type $(typeof(m)) does not support accessing the attribute $attr via getall!"))
end

"""
$(FUNCTIONNAME)

Assigns a `value` to the `attr` attribute of the integrator, problem, or relaxation.
"""
function set! end

function set!(m::AbstractDERelaxIntegrator, attr::AbstractIntegratorAttribute, idxs::Vector, vec_of_val::Vector)
    if length(idxs) !== length(vec_of_val)
        throw(DimensionMismatch("Number of indices ($(length(idxs))) does " *
                                "not match the number of values " *
                                "($(length(vec_of_val))) set to `$attr`."))
    end
	for (i, v) in enumerate(vec_of_val)
		set!(m, attr, @inbounds idxs[i], v)
	end
end

function set!(m::AbstractDERelaxProblem, attr::AbstractRelaxProblemAttribute, idxs::Vector, vec_of_val::Vector)
    if length(idxs) !== length(vec_of_val)
        throw(DimensionMismatch("Number of indices ($(length(idxs))) does " *
                                "not match the number of values " *
                                "($(length(vec_of_val))) set to `$attr`."))
    end
	for (i, v) in enumerate(vec_of_val)
		set!(m, attr, @inbounds idxs[i], v)
	end
end

function set!(m::AbstractDERelaxIntegrator, attr::AbstractIntegratorAttribute, args...)
    throw_set_error_fallback(m, attr, args...)
end

function set!(m::AbstractDERelaxProblem, attr::AbstractRelaxProblemAttribute, args...)
    throw_set_error_fallback(m, attr, args...)
end

"""
$(FUNCTIONNAME)

An in-place version of `set!` which sets all data asssociated with a particular
integrator attribute.
"""
function setall!(m::AbstractDERelaxIntegrator, attr::AbstractIntegratorAttribute, args...)
	throw_set_error_fallback(m, attr, args...)
end

"""
UnsupportedError <: Exception

Abstract type for error thrown when an element is not supported by the integrator.
"""
abstract type UnsupportedError <: Exception end

"""
$(TYPEDEF)

An error indicating that the attribute `attr` is unsupported.
"""
struct UnsupportedRelaxAttribute{AttrType<:AnyDEAttribute} <: UnsupportedError
    attr::AttrType
    message::String
end
UnsupportedRelaxAttribute(attr::AnyDEAttribute) = UnsupportedRelaxAttribute(attr, "")
UnsupportedRelaxAttribute(err::UnsupportedRelaxAttribute) = "Attribute $(err.attr)"

"""
NotAllowedError <: Exception

Abstract type for error thrown when an operation is supported but cannot be
applied in the current state of the integrator
.
"""
abstract type NotAllowedError <: Exception end

"""
$(TYPEDEF)

An error indicating that the attribute `attr` is supported but cannot be set
for some reason given in the error string.
"""
struct SetRelaxAttributeNotAllowed{AttrType<:AnyDEAttribute} <: NotAllowedError
    attr::AttrType
	message::String
end
SetRelaxAttributeNotAllowed(attr::AnyDEAttribute) = SetRelaxAttributeNotAllowed(attr, "")
operation_name(err::SetRelaxAttributeNotAllowed) = "Setting attribute $(err.attr)"
message(err::SetRelaxAttributeNotAllowed) = err.message

for (mtype, attr_type) in ((AbstractDERelaxIntegrator, AbstractIntegratorAttribute),
	                       (AbstractDERelaxProblem, AbstractRelaxProblemAttribute))
	@eval function throw_set_error_fallback(storage::$mtype,
	                                  attr::$attr_type,
	                                  value;
	                                  error_if_supported = SetRelaxAttributeNotAllowed(attr))
	    if supports(storage, attr)
	        throw(error_if_supported)
	    else
	        throw(UnsupportedRelaxAttribute(attr))
	    end
	end
	@eval function throw_set_error_fallback(storage::$mtype,
	                                  attr::$attr_type,
	                                  index, value;
	                                  error_if_supported = SetRelaxAttributeNotAllowed(attr))
	    if supports(storage, attr, typeof(index))
	        throw(error_if_supported)
	    else
	        throw(UnsupportedRelaxAttribute(attr))
	    end
	end
end


text:
abstract type Properties end

include("filters/filters.jl")
include("drivers/drivers.jl")

Base.cconvert(::Type{API.hid_t}, obj::Properties) = obj
Base.unsafe_convert(::Type{API.hid_t}, obj::Properties) = obj.id

function Base.close(obj::Properties)
    if obj.id != -1
        if isvalid(obj)
            API.h5p_close(obj)
        end
        obj.id = -1
    end
    nothing
end

Base.isvalid(obj::Properties) = obj.id != -1 && API.h5i_is_valid(obj)

# By default, properties objects are only initialized lazily
function init!(prop::P) where {P<:Properties}
    if !isvalid(prop)
        prop.id = API.h5p_create(classid(P))
    end
    return prop
end

function (::Type{P})(;kwargs...) where {P <: Properties}
    obj = P(API.H5P_DEFAULT)
    for (k, v) in kwargs
        setproperty!(obj, k, v)
    end
    return obj
end

function Base.getproperty(p::P, name::Symbol) where {P <: Properties}
    name === :id    ? getfield(p, :id) :
    class_getproperty(P, p, name)
end

function Base.setproperty!(p::P, name::Symbol, val) where {P<:Properties}
    if name === :id
        return setfield!(p, :id, API.hid_t(val))
    end
    init!(p)
    class_setproperty!(P, p, name, val)
end

Base.propertynames(p::P) where {P<:Properties} = (all_propertynames(P)..., :id)
all_propertynames(::Type{P}) where {P<:Properties} = (class_propertynames(P)..., all_propertynames(superclass(P))...,)

# defaults: refer to super class
class_getproperty(::Type{P}, props, name) where {P<:Properties} =
    class_getproperty(superclass(P), props, name)
class_setproperty!(::Type{P}, p, name, val) where {P<:Properties} =
    class_setproperty!(superclass(P), p, name, val)
class_propertynames(::Type{P}) where {P<:Properties} = ()

"""
    @propertyclass P classid

Define a new subtype of `P <: Properties` corresponding to a HDF5 property list
with class identifier `classid`.

Once defined, the following interfaces can be defined:

    superclass(::Type{P})

This should return the type from which `P` inherits. If not defined, it will
inherit from `GenericProperties`.

    class_propertynames(::Type{P})

This should return a `Tuple` of `Symbol`s, being the names of the properties
associated with `P`.

    class_getproperty(::Type{P}, p::Properties, name::Symbol)

If `name` is an associated property of type `P`, this should return the value of
the propery, otherwise call `class_getproperty(superclass(P), p, name)`.

    class_setproperty!(::Type{P}, p::Properties, name::Symbol, val)

If `name` is an associated property of type `P`, this should set the value of
the propery, otherwise call `class_setproperty!(superclass(P), p, name, val)`.
"""
macro propertyclass(name, classid)
    expr = quote
        Core.@__doc__ mutable struct $name <: Properties
            id::API.hid_t
            function $name(id)
                obj = new(id)
                finalizer(close, obj)
                obj
            end
        end
        classid(::Type{$name}) = $classid
    end
    return esc(expr)
end

@propertyclass GenericProperties API.H5P_DEFAULT

superclass(::Type{P}) where {P<:Properties} = GenericProperties

class_getproperty(::Type{GenericProperties}, props, name) =
    error("$(typeof(props)) has no property $name")
class_setproperty!(::Type{GenericProperties}, props, name, val) =
    error("$(typeof(props)) has no property $name")
all_propertynames(::Type{GenericProperties}) = ()

# for initializing multiple Properties from a set of keyword arguments
"""
    setproperties!(props::Properties...; kwargs...)

For each `(key, value)` pair in `kwargs`, set the corresponding properties in
each `Properties` object in `props`. Returns a `Dict` of any pairs which didn't
match properties in `props`.
"""
function setproperties!(props::Properties...; kwargs...)
    filter(kwargs) do (k,v)
        found = false
        for prop in props
            if k in all_propertynames(typeof(prop))
                setproperty!(prop,k,v)
                found = true
            end
        end
        return !found
    end
end

###
### Convenience macros for defining getter/setter functions
###

"""
    @tuple_property(name)
"""
macro tuple_property(property)
    get_property = Symbol(:get_,property)
    set_property! = Symbol(:set_,property,:!)
    api_get_property = :(API.$(Symbol(:h5p_get_,property)))
    api_set_property = :(API.$(Symbol(:h5p_set_,property)))
    quote
        function $(esc(get_property))(p::Properties)
            return $api_get_property(p)
        end
        function $(esc(set_property!))(p::Properties, val::Tuple)
            return $api_set_property(p, val...)
        end
    end
end

"""
    @enum_property(name, sym1 => enumvalue1, sym2 => enumvalue2, ...)

Wrap property getter/setter API functions that use enum values to use symbol instead.
"""
macro enum_property(property, pairs...)
    get_property = Symbol(:get_,property)
    set_property! = Symbol(:set_,property,:!)
    api_get_property = :(API.$(Symbol(:h5p_get_,property)))
    api_set_property = :(API.$(Symbol(:h5p_set_,property)))

    get_expr = :(error("Unknown $property value $enum"))
    set_expr = :(throw(ArgumentError("Invalid $property $val")))

    for pair in reverse(pairs)
        @assert pair isa Expr && pair.head == :call && pair.args[1] == :(=>)
        _, val, enum = pair.args
        get_expr = :(enum == $enum ? $val : $get_expr)
        set_expr = :(val == $val ? $enum : $set_expr)
    end
    quote
        function $(esc(get_property))(p::Properties)
            property = $(QuoteNode(property))
            enum = $api_get_property(p)
            return $get_expr
        end
        function $(esc(set_property!))(p::Properties, val)
            property = $(QuoteNode(property))
            enum = $set_expr
            return $api_set_property(p, enum)
        end
        function $(esc(set_property!))(p::Properties, enum::Integer)
            # deprecate?
            return $api_set_property(p, enum)
        end
    end
end

"""
    @bool_property(name)

Wrap property getter/setter API functions that use `0`/`1` to use `Bool` values
"""
macro bool_property(property)
    get_property = Symbol(:get_,property)
    set_property! = Symbol(:set_,property,:!)
    api_get_property = :(API.$(Symbol(:h5p_get_,property)))
    api_set_property = :(API.$(Symbol(:h5p_set_,property)))
    quote
        function $(esc(get_property))(p::Properties)
            return $api_get_property(p) != 0
        end
        function $(esc(set_property!))(p::Properties, val)
            return $api_set_property(p, val)
        end
    end
end

###
### Define Properties types
###

"""
    ObjectCreateProperties(;kws...)

Properties used when creating a new object. Available options:

- `obj_track_times :: Bool`: governs the recording of times associated with an
  object. If set to `true`, time data will be recorded. See
  $(h5doc("H5P_SET_OBJ_TRACK_TIMES")).

"""
@propertyclass ObjectCreateProperties API.H5P_OBJECT_CREATE

@bool_property(obj_track_times)

class_propertynames(::Type{ObjectCreateProperties}) = (
    :obj_track_times,
    :track_times,
    )
function class_getproperty(::Type{ObjectCreateProperties}, p::Properties, name::Symbol)
    name === :obj_track_times ? get_obj_track_times(p) :
    # deprecated
    name === :track_times ? (depwarn("`track_times` property is deprecated, use `obj_track_times` instead",:track_times); get_obj_track_times(p)) :
    class_getproperty(superclass(ObjectCreateProperties), p, name)
end
function class_setproperty!(::Type{ObjectCreateProperties}, p::Properties, name::Symbol, val)
    name === :obj_track_times ? set_obj_track_times!(p, val) :
    # deprecated
    name === :track_times ? (depwarn("`track_times=$val` keyword option is deprecated, use `obj_track_times=$val` instead",:track_times); set_obj_track_times!(p, val)) :
    class_setproperty!(superclass(ObjectCreateProperties), p, name, val)
end

"""
    GroupCreateProperties(;kws...)

Properties used when creating a new `Group`. Inherits from
[`ObjectCreateProperties`](@ref), with additional options:

- `local_heap_size_hint :: Integer`: the anticipated maximum local heap size in
  bytes. See $(h5doc("H5P_SET_LOCAL_HEAP_SIZE_HINT")).

"""
@propertyclass GroupCreateProperties API.H5P_GROUP_CREATE
superclass(::Type{GroupCreateProperties}) = ObjectCreateProperties

class_propertynames(::Type{GroupCreateProperties}) = (
    :local_heap_size_hint,
    )
function class_getproperty(::Type{GroupCreateProperties}, p::Properties, name::Symbol)
    name === :local_heap_size_hint ? API.h5p_get_local_heap_size_hint(p) :
    class_getproperty(superclass(GroupCreateProperties), p, name)
end
function class_setproperty!(::Type{GroupCreateProperties}, p::Properties, name::Symbol, val)
    name === :local_heap_size_hint ? API.h5p_set_local_heap_size_hint(p, val) :
    class_setproperty!(superclass(GroupCreateProperties), p, name, val)
end

"""
    FileCreateProperties(;kws...)

Properties used when creating a new `File`. Inherits from
[`ObjectCreateProperties`](@ref),  with additional properties:

- `userblock :: Integer`: user block size in bytes. The default user block size
  is 0; it may be set to any power of 2 equal to 512 or greater (512, 1024,
  2048, etc.). See $(h5doc("H5P_SET_USERBLOCK")).

"""
@propertyclass FileCreateProperties API.H5P_FILE_CREATE
superclass(::Type{FileCreateProperties}) = ObjectCreateProperties

class_propertynames(::Type{FileCreateProperties}) = (
    :userblock,
    )
function class_getproperty(::Type{FileCreateProperties}, p::Properties, name::Symbol)
    name === :userblock   ? API.h5p_get_userblock(p) :
    class_getproperty(superclass(FileCreateProperties), p, name)
end
function class_setproperty!(::Type{FileCreateProperties}, p::Properties, name::Symbol, val)
    name === :userblock   ? API.h5p_set_userblock(p, val) :
    class_setproperty!(superclass(FileCreateProperties), p, name, val)
end

"""
    DatatypeCreateProperties(;kws...)

"""
@propertyclass DatatypeCreateProperties API.H5P_DATATYPE_CREATE
superclass(::Type{DatatypeCreateProperties}) = ObjectCreateProperties

"""
    DatasetCreateProperties(;kws...)

Properties used when creating a new `Dataset`. Inherits from
[`ObjectCreateProperties`](@ref), with additional properties:

- `alloc_time`: the timing for the allocation of storage space for a dataset's
  raw data; one of:
   - `:default`

   - `:early`: allocate all space when the dataset is created

   - `:incremental`: Allocate space incrementally, as data is  written to the
     dataset

   - `:late`: Allocate all space when data is first written to the dataset.

  See $(h5doc("H5P_SET_ALLOC_TIME")).

- `chunk`: a tuple containing the size of the chunks to store each dimension.
  See $(h5doc("H5P_SET_CHUNK")) (note that this uses Julia's column-major
  ordering).

- `external`: A tuple of `(name,offset,size)`, See $(h5doc("H5P_SET_EXTERNAL")).

- `filters` (only valid when `layout=:chunked`): a filter or vector of filters
  that are applied to applied to each chunk of a dataset, see [Filters](@ref).
  When accessed, will return a [`Filters.FilterPipeline`](@ref) object that can
  be modified in-place.

- `layout`: the type of storage used to store the raw data for a dataset. Can be
  one of:

   - `:compact`: Store raw data in the dataset object header in file. This
     should only be used for datasets with small amounts of raw data.

   - `:contiguous`: Store raw data separately from the object header in one
     large chunk in the file.

   - `:chunked`: Store raw data separately from the object header as chunks of
     data in separate locations in the file.

   - `:virtual`:  Draw raw data from multiple datasets in different files.

  See $(h5doc("H5P_SET_LAYOUT")).

The following options are shortcuts for the various filters, and are set-only.
They will be appended to the filter pipeline in the order in which they appear

- `blosc = true | level`: set the [`Filters.BloscFilter`](@ref) compression
  filter; argument can be either `true`, or the compression level.

- `deflate = true | level`: set the [`Filters.Deflate`](@ref) compression
  filter; argument can be either `true`, or the compression level.

- `fletcher32 = true`: set the [`Filters.Fletcher32`](@ref) checksum filter.

- `shuffle = true`: set the [`Filters.Shuffle`](@ref) filter.

"""
@propertyclass DatasetCreateProperties API.H5P_DATASET_CREATE
superclass(::Type{DatasetCreateProperties}) = ObjectCreateProperties

@enum_property(alloc_time,
    :default     => API.H5D_ALLOC_TIME_DEFAULT,
    :early       => API.H5D_ALLOC_TIME_EARLY,
    :incremental => API.H5D_ALLOC_TIME_INCR,
    :late        => API.H5D_ALLOC_TIME_LATE)

# reverse indices
function get_chunk(p::Properties)
    dims, N = API.h5p_get_chunk(p)
    ntuple(i -> Int(dims[N-i+1]), N)
end
set_chunk!(p::Properties, dims) = API.h5p_set_chunk(p, length(dims), API.hsize_t[reverse(dims)...])

@enum_property(layout,
    :compact    => API.H5D_COMPACT,
    :contiguous => API.H5D_CONTIGUOUS,
    :chunked    => API.H5D_CHUNKED,
    :virtual    => API.H5D_VIRTUAL)

# filters getters/setters
get_filters(p::Properties) = Filters.FilterPipeline(p)
set_filters!(p::Properties, val::Filters.Filter) = push!(empty!(Filters.FilterPipeline(p)), val)
set_filters!(p::Properties, vals::Union{Tuple, AbstractVector}) = append!(empty!(Filters.FilterPipeline(p)), vals)

# convenience
set_deflate!(p::Properties, val::Bool) = val && push!(Filters.FilterPipeline(p), Filters.Deflate())
set_deflate!(p::Properties, level::Integer) = push!(Filters.FilterPipeline(p), Filters.Deflate(level=level))
set_shuffle!(p::Properties, val::Bool) = val && push!(Filters.FilterPipeline(p), Filters.Shuffle())
set_fletcher32!(p::Properties, val::Bool) = val && push!(Filters.FilterPipeline(p), Filters.Fletcher32())

class_propertynames(::Type{DatasetCreateProperties}) = (
    :alloc_time,
    :chunk,
    :external,
    :filters,
    :layout,
    # convenience
    :blosc,
    :deflate,
    :fletcher32,
    :shuffle,
    # deprecated
    :compress,
    )

function class_getproperty(::Type{DatasetCreateProperties}, p::Properties, name::Symbol)
    name === :alloc_time  ? get_alloc_time(p) :
    name === :chunk       ? get_chunk(p) :
    #name === :external    ? API.h5p_get_external(p) :
    name === :filters     ? get_filters(p) :
    name === :layout      ? get_layout(p) :
    class_getproperty(superclass(DatasetCreateProperties), p, name)
end
function class_setproperty!(::Type{DatasetCreateProperties}, p::Properties, name::Symbol, val)
    name === :alloc_time  ? set_alloc_time!(p, val) :
    name === :chunk       ? set_chunk!(p, val) :
    name === :external    ? API.h5p_set_external(p, val...) :
    name === :filters     ? set_filters!(p, val) :
    name === :layout      ? set_layout!(p, val) :
    # set-only for convenience
    name === :blosc       ? set_blosc!(p, val) :
    name === :deflate     ? set_deflate!(p, val) :
    name === :fletcher32  ? set_fletcher32!(p, val) :
    name === :shuffle     ? set_shuffle!(p, val) :
    # deprecated
    name === :compress    ? (depwarn("`compress=$val` keyword option is deprecated, use `deflate=$val` instead",:class_setproperty!); set_deflate!(p, val)) :
    class_setproperty!(superclass(DatasetCreateProperties), p, name, val)
end

"""
    StringCreateProperties(;kws...)

"""
@propertyclass StringCreateProperties API.H5P_STRING_CREATE

@enum_property(char_encoding,
    :ascii => API.H5T_CSET_ASCII,
    :utf8  => API.H5T_CSET_UTF8)

class_propertynames(::Type{StringCreateProperties}) = (
    :char_encoding,
    )
function class_getproperty(::Type{StringCreateProperties}, p::Properties, name::Symbol)
    name === :char_encoding ? get_char_encoding(p) :
    class_getproperty(superclass(StringCreateProperties), p, name)
end
function class_setproperty!(::Type{StringCreateProperties}, p::Properties, name::Symbol, val)
    name === :char_encoding ? set_char_encoding!(p, val) :
    class_setproperty!(superclass(StringCreateProperties), p, name, val)
end

"""
    LinkCreateProperties(;kws...)

Properties used when creating links.

- `char_encoding`: the character enconding, either `:ascii` or `:utf8`.

- `create_intermediate_group :: Bool`: if `true`, will create missing
  intermediate groups.
"""
@propertyclass LinkCreateProperties API.H5P_LINK_CREATE
superclass(::Type{LinkCreateProperties}) = StringCreateProperties

@bool_property(create_intermediate_group)

class_propertynames(::Type{LinkCreateProperties}) = (
    :create_intermediate_group,
    )
function class_getproperty(::Type{LinkCreateProperties}, p::Properties, name::Symbol)
    name === :create_intermediate_group ? get_create_intermediate_group(p) :
    class_getproperty(superclass(LinkCreateProperties), p, name)
end
function class_setproperty!(::Type{LinkCreateProperties}, p::Properties, name::Symbol, val)
    name === :create_intermediate_group ? set_create_intermediate_group!(p, val) :
    class_setproperty!(superclass(LinkCreateProperties), p, name, val)
end

"""
    AttributeCreateProperties(;kws...)

Properties used when creating attributes.

- `char_encoding`: the character enconding, either `:ascii` or `:utf8`.
"""
@propertyclass AttributeCreateProperties API.H5P_ATTRIBUTE_CREATE
superclass(::Type{AttributeCreateProperties}) = StringCreateProperties

"""
    FileAccessProperties(;kws...)

Properties used when accessing files.

- `alignment :: Tuple{Integer, Integer}`: a `(threshold, alignment)` pair: any
  file object greater than or equal in size to threshold bytes will be aligned
  on an address which is a multiple of alignment. Default values are 1, implying
  no alignment.

- `driver`: the file driver used to access the file. See [Drivers](@ref).

- `driver_info` (get only)

- `fclose_degree`: file close degree property. One of:

  - `:weak`
  - `:semi`
  - `:strong`
  - `:default`

- `libver_bounds`: a `(low, high)` pair: `low` sets the earliest possible format
  versions that the library will use when creating objects in the file; `high`
  sets the latest format versions that the library will be allowed to use when
  creating objects in the file. Values can be a `VersionNumber` for the hdf5
  library, `:earliest`, or `:latest` . See $(h5doc("H5P_SET_LIBVER_BOUNDS"))

"""
@propertyclass FileAccessProperties API.H5P_FILE_ACCESS

@tuple_property(alignment)

@enum_property(fclose_degree,
               :weak    => API.H5F_CLOSE_WEAK,
               :semi    => API.H5F_CLOSE_SEMI,
               :strong  => API.H5F_CLOSE_STRONG,
               :default => API.H5F_CLOSE_DEFAULT)

# getter/setter for libver_bounds
libver_bound_to_enum(val::Integer) = val
function libver_bound_to_enum(val::VersionNumber)
    val >= v"1.12"   ? API.H5F_LIBVER_V112 :
    val >= v"1.10"   ? API.H5F_LIBVER_V110 :
    val >= v"1.8"    ? API.H5F_LIBVER_V18 :
    throw(ArgumentError("libver_bound must be >= v\"1.8\"."))
end
function libver_bound_to_enum(val::Symbol)
    val == :earliest ? API.H5F_LIBVER_EARLIEST :
    val == :latest   ? libver_bound_to_enum(libversion) :
    throw(ArgumentError("Invalid libver_bound $val."))
end
function libver_bound_from_enum(enum)
    enum == API.H5F_LIBVER_EARLIEST ? :earliest :
    enum == API.H5F_LIBVER_V18      ? v"1.8" :
    enum == API.H5F_LIBVER_V110     ? v"1.10" :
    enum == API.H5F_LIBVER_V112     ? v"1.12" :
    error("Unknown libver_bound value $enum")
end
function get_libver_bounds(p::Properties)
    low, high = API.h5p_get_libver_bounds(p)
    return libver_bound_from_enum(low), libver_bound_from_enum(high)
end
function set_libver_bounds!(p::Properties, (low, high)::Tuple{Any,Any})
    API.h5p_set_libver_bounds(p, libver_bound_to_enum(low), libver_bound_to_enum(high))
end
function set_libver_bounds!(p::Properties, val)
    API.h5p_set_libver_bounds(p, libver_bound_to_enum(val), libver_bound_to_enum(val))
end

class_propertynames(::Type{FileAccessProperties}) = (
    :alignment,
    :driver,
    :driver_info,
    :fapl_mpio,
    :fclose_degree,
    :libver_bounds,
    )

function class_getproperty(::Type{FileAccessProperties}, p::Properties, name::Symbol)
    name === :alignment     ? get_alignment(p) :
    name === :driver        ? Drivers.get_driver(p) :
    name === :driver_info   ? API.h5p_get_driver_info(p) : # get only
    name === :fclose_degree ? get_fclose_degree(p) :
    name === :libver_bounds ? get_libver_bounds(p) :
    # deprecated
    name === :fapl_mpio     ? (depwarn("The `fapl_mpio` property is deprecated, use `driver=HDF5.Drivers.MPIO(...)` instead.", :fapl_mpio); drv = get_driver(p, MPIO); (drv.comm, drv.info)) :
    class_getproperty(superclass(FileAccessProperties), p, name)
end
function class_setproperty!(::Type{FileAccessProperties}, p::Properties, name::Symbol, val)
    name === :alignment     ? set_alignment!(p, val) :
    name === :driver        ? Drivers.set_driver!(p, val) :
    name === :fclose_degree ? set_fclose_degree!(p, val) :
    name === :libver_bounds ? set_libver_bounds!(p, val) :
    # deprecated
    name === :fapl_mpio     ? (depwarn("The `fapl_mpio` property is deprecated, use `driver=HDF5.Drivers.MPIO(...)` instead.", :fapl_mpio); p.driver = Drivers.MPIO(val...)) :
    class_setproperty!(superclass(FileAccessProperties), p, name, val)
end

@propertyclass LinkAccessProperties API.H5P_LINK_ACCESS

@propertyclass GroupAccessProperties API.H5P_GROUP_ACCESS
superclass(::Type{GroupAccessProperties}) = LinkAccessProperties

@propertyclass DatatypeAccessProperties API.H5P_DATATYPE_ACCESS
superclass(::Type{DatatypeAccessProperties}) = LinkAccessProperties

@propertyclass DatasetAccessProperties API.H5P_DATASET_ACCESS
superclass(::Type{DatasetAccessProperties}) = LinkAccessProperties

@propertyclass AttributeAccessProperties API.H5P_ATTRIBUTE_ACCESS
superclass(::Type{AttributeAccessProperties}) = LinkAccessProperties

"""
    DatasetTransferProperties(;kws...)

Properties used when transferring data to/from datasets

- `dxpl_mpio`: MPI transfer mode:
   - `:independent`: use independent I/O access (default),
   - `:collective`: use collective I/O access.
"""
@propertyclass DatasetTransferProperties API.H5P_DATASET_XFER

@enum_property(dxpl_mpio,
               :independent => API.H5FD_MPIO_INDEPENDENT,
               :collective  => API.H5FD_MPIO_COLLECTIVE)

class_propertynames(::Type{DatasetTransferProperties}) = (
    :dxpl_mpio,
    )
function class_getproperty(::Type{DatasetTransferProperties}, p::Properties, name::Symbol)
    name === :dxpl_mpio  ? get_dxpl_mpio(p) :
    class_getproperty(superclass(DatasetTransferProperties), p, name)
end
function class_setproperty!(::Type{DatasetTransferProperties}, p::Properties, name::Symbol, val)
    name === :dxpl_mpio  ? set_dxpl_mpio!(p, val) :
    class_setproperty!(superclass(DatasetTransferProperties), p, name, val)
end

@propertyclass FileMountProperties API.H5P_FILE_MOUNT
@propertyclass ObjectCopyProperties API.H5P_OBJECT_COPY

const DEFAULT_PROPERTIES = GenericProperties()
# These properties are initialized in __init__()
const ASCII_LINK_PROPERTIES = LinkCreateProperties()
const UTF8_LINK_PROPERTIES = LinkCreateProperties()
const ASCII_ATTRIBUTE_PROPERTIES = AttributeCreateProperties()
const UTF8_ATTRIBUTE_PROPERTIES = AttributeCreateProperties()

_link_properties(::AbstractString) = UTF8_LINK_PROPERTIES
_attr_properties(::AbstractString) = UTF8_ATTRIBUTE_PROPERTIES


text:
using AtomsBase
using Unitful
using UnitfulAtomic
using StaticArrays

function load_data(;num_entries = 2000, file = "LJCluster/curated_lj_cluster.xyz" )
    systems  = Vector{AbstractSystem}(undef, num_entries)
    energies = Vector{Float64}(undef, num_entries)
    forces    = Vector{Vector{T} where T<:SVector{3, <:Float64}}(undef, num_entries)
    bias = @SVector [4.0, 4.0, 4.0]
    box = [[8.0, 0.0, 0.0], [0.0, 8.0, 0.0], [0.0, 0.0, 9.0]] * 1u"Å"
    bcs = [DirichletZero(), DirichletZero(), DirichletZero()]
    open(file, "r") do io
        count = 1
        while !eof(io) && (count <= num_entries)
            line = readline(io)
            num_atoms = parse(Int, line)
            info_line = split(readline(io))
            energies[count] = parse(Float64, info_line[2][8:end])

            atoms = Vector{Atom}(undef, 13)
            force = Vector{SVector{3, Float64}}(undef, 13)
            for i = 1:num_atoms
                line = split(readline(io))
                element = Symbol(line[1])
                position = @SVector [ parse(Float64, line[2]), parse(Float64, line[3]), parse(Float64, line[4]) ]
                position += bias
                atoms[i] = Atom(element, position * 1u"Å")

                force[i] = @SVector [ parse(Float64, line[5]), parse(Float64, line[6]), parse(Float64, line[7])]
            end

            forces[count] = force
            systems[count] = FlexibleSystem(atoms, box, bcs)
            count += 1
        end
        println(count)
    end
    return systems, energies, forces
end


text:
#--------------------------------------------------------------------
# Spacetime Discretization methods in Julia
# Soham 01-2018
#--------------------------------------------------------------------

function getPatchIC{T<:Integer}(fn::Function, s::T, Nx::T, M::T, loc::Int)::Boundary
    if s == 0
        return Boundary(0, projectonPatchBndbyRestriction(fn, Nx, M, loc))
    elseif s == 1
        return Boundary(1, projectonPatchBndbyRestriction(fn, Nx, M, loc))
    else
        error("Unknown direction passed.")
    end
end

function getPatchBnd(patch::Patch, s::Int)::Boundary
    if s == 0
        boundary = Boundary(0, patch.value[:, end])
    elseif s == 1
        boundary = Boundary(1, patch.value[end, :])
    else
        error("Unknown direction passed")
    end
    return boundary
end

function calcPatch(bndx::Boundary, bndy::Boundary, RHS::Array{Float64,2},
                   derivOP::Array{Float64,4}, boundaryOP::Array{Float64,4},
                   loc::Array{Int,1})::Patch
    Nx   = size(bndx.value)[1] - 1
    Ny   = size(bndy.value)[1] - 1
    bval = zeros(Nx+1, Ny+1)
    bval[:, 1] = bndx.value
    bval[1, :] = bndy.value
    return Patch(loc, shapeL2H(shapeH2L(derivOP + boundaryOP) \ shapeH2L(RHS + bval), Nx, Ny))
end

function extractPatchCoeffs(patch::Patch)::Array{Float64,2}
    fnodal  = patch.value
    Nx      = size(fnodal)[1] - 1
    Ny      = size(fnodal)[2] - 1
    invndmx = inv(vandermonde(Nx))
    invndmy = inv(vandermonde(Ny))
    fmodal  = zeros(Nx+1, Ny+1)
    for index in CartesianRange(size(fmodal))
        i = index.I[1]
        j = index.I[2]
        fmodal[i,j] = sum(invndmx[i,m]*fnodal[m,n]*invndmy[j,n] for m in 1:Nx+1, n in 1:Ny+1)
    end
    return fmodal
end

function interpolatePatch(patch::Patch, Nx::Int, Ny::Int)::Patch
    fmodal = extractPatchCoeffs(patch)
    Nxp    = size(fmodal)[1] - 1
    Nyp    = size(fmodal)[2] - 1
    pvndmx = pseudovandermonde(Nxp, chebgrid(Nx))
    pvndmy = pseudovandermonde(Nyp, chebgrid(Ny))
    fnodal = zeros(Nx+1, Ny+1)
    for index in CartesianRange(size(fnodal))
        i = index.I[1]
        j = index.I[2]
        fnodal[i,j] = sum(pvndmx[i,m]*fmodal[m,n]*pvndmy[j,n] for m in 1:Nxp+1, n in 1:Nyp+1)
    end
    return Patch(patch.loc, fnodal)
end


text:
struct ProdMean{T<:NTuple{N,Mean} where N} <: CompositeMean
    means::T
end

ProdMean(means::Mean...) = ProdMean(means)

mean(pm::ProdMean, x::AbstractVector) = prod(mean(m, x) for m in components(pm))

get_param_names(pm::ProdMean) = composite_param_names(components(pm), :pm)

function grad_mean(pm::ProdMean, x::AbstractVector)
    dm = Array{Float64}(undef, num_params(pm))
    means = [mean(m, x) for m in components(pm)]
    v = 1
    for (i, m) in enumerate(components(pm))
        p = prod(m2 for (j, m2) in enumerate(means) if i != j)
        w = v + num_params(m)
        dm[v:(w - 1)] = p * grad_mean(m, x)
        v = w
    end
    dm
end

# Multiplication operators
Base.:*(m1::ProdMean, m2::Mean) = ProdMean(m1.means..., m2)
Base.:*(m1::ProdMean, m2::ProdMean) = ProdMean(m1.means..., m2.means...)
Base.:*(m1::Mean, m2::Mean) = ProdMean(m1, m2)
Base.:*(m1::Mean, m2::ProdMean) = ProdMean(m1, m2.means...)


text:
module KdotP
# ---------------------------------------------------------------------------------------- #

using Crystalline, GellMannMatrices, LinearAlgebra, PrettyTables, Unicode
import Crystalline: irdim

export kdotp, irdim

# ---------------------------------------------------------------------------------------- #
const ATOL_DEFAULT = 1e-10
# ---------------------------------------------------------------------------------------- #

include("types.jl")
include("show.jl")
include("constraints.jl")
include("compare.jl")
include("timereversal.jl")

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

end # module


text:
const h_crit = Float32(1/sqrt(2))

const fc_vals = Tuple(vcat(Float64[1/120, 1/60],  # STS-2, CMG-3T
                       repeat([1/30], 2),         # CMG-40
                       [0.2],                     # Lennartz LE-3D
                       repeat([1.0], 5),          # ...everything...
                       repeat([2.0], 2),          # passive geophones
                       [4.5]))                    # passive industry geophones
const fs_vals = (0.1, 1.0, 2.0, 5.0, 10.0, 20.0, 25.0, 40.0, 50.0, 60.0, 62.5, 80.0, 100.0)

const irregular_units = ("%",
 "%{cloud_cover}",
 "{direction_vector}",
 "Cel",
 "K",
 "{none}",
 "Pa",
 "T",
 "V",
 "W",
 "m",
 "m/m",
 "m/s",
 "m/s2",
 "m3/m3",
 "rad",
 "rad/s",
 "rad/s2",
 "t{SO_2}")

# Acceptable type codes in :misc
const OK = (      0x00, 0x01,
                  0x10, 0x11, 0x12, 0x13, 0x14,
                  0x20, 0x21, 0x22, 0x23, 0x24,
                  0x30, 0x31, 0x32,
                  0x50, 0x51, 0x52, 0x53, 0x54,
                  0x60, 0x61, 0x62, 0x63, 0x64,
                  0x70, 0x71, 0x72,
                  0x80, 0x81,
                  0x90, 0x91, 0x92, 0x93, 0x94,
                  0xa0, 0xa1, 0xa2, 0xa3, 0xa4,
                  0xb0, 0xb1, 0xb2,
                  0xd0, 0xd1, 0xd2, 0xd3, 0xd4,
                  0xe0, 0xe1, 0xe2, 0xe3, 0xe4,
                  0xf0, 0xf1, 0xf2 )

const evtypes = (       "not_existing",
                        "not_reported",
                        "anthropogenic_event",
                        "collapse",
                        "cavity_collapse",
                        "mine_collapse",
                        "building_collapse",
                        "explosion",
                        "accidental_explosion",
                        "chemical_explosion",
                        "controlled_explosion",
                        "experimental_explosion",
                        "industrial_explosion",
                        "mining_explosion",
                        "quarry_blast",
                        "road_cut",
                        "blasting_levee",
                        "nuclear_explosion",
                        "induced_or_triggered_event",
                        "rock_burst",
                        "reservoir_loading",
                        "fluid_injection",
                        "fluid_extraction",
                        "crash",
                        "plane_crash",
                        "train_crash",
                        "boat_crash",
                        "other_event",
                        "atmospheric_event",
                        "sonic_boom",
                        "sonic_blast",
                        "acoustic_noise",
                        "thunder",
                        "avalanche",
                        "snow_avalanche",
                        "debris_avalanche",
                        "hydroacoustic_event",
                        "ice_quake",
                        "slide",
                        "landslide",
                        "rockslide",
                        "meteorite",
                        "volcanic_eruption" )
const phase_list = ("P",
                    "PKIKKIKP",
                    "PKIKKIKS",
                    "PKIKPPKIKP",
                    "PKPPKP",
                    "PKiKP",
                    "PP",
                    "PS",
                    "PcP",
                    "S",
                    "SKIKKIKP",
                    "SKIKKIKS",
                    "SKIKSSKIKS",
                    "SKKS",
                    "SKS",
                    "SKiKP",
                    "SP",
                    "SS",
                    "ScS",
                    "pP",
                    "pPKiKP",
                    "pS",
                    "pSKS",
                    "sP",
                    "sPKiKP",
                    "sS",
                    "sSKS")

const pol_list = ('U', 'D', '-', '+', '_', ' ')
const loc_types = ("HYPOCENTER", "CENTROID", "AMPLITUDE", "MACROSEISMIC", "RUPTURE_START", "RUPTURE_END")
const loc_methods = ("HYPOELLIPSE", "HypoDD", "Velest", "centroid")

const hln = ('H','L','N')
const iclist = ('A','B','D','F','G','I','J','K','M','O','P','Q','R','S','T','U','V','W','Z')
const oid = ('O','I','D')
const zne = ('Z','N','E')
const nvc = ('A','B','C','1','2','3','U','V','W')
const oidfhu = ('O','I','D','F','H','U')
const icfo = ('I','C','F','O')
const junits = ("rad", "rad/s", "rad/s2")

const geodetic_datum = ("ETRS89", "GRS 80", "JGD2011")


text:
"""
least_squares_function_fit.jl

Author: Philippe Nadon
Date: April 26, 2019

Demonstrates LU decomposition and Least Squares fitting.
This code is not production-ready, and is frankly poorly coded.
The purpose was not to design elegant code,
but rather to demonstrate the capabilities of Julia.
"""

using Pkg; Pkg.add("QuadGK")
using QuadGK
using LinearAlgebra

function main()
	func( t) = 2*t^2 + 4
	bounds = [ 0, 2]
	num_digits = 3

	x = approx_func( func, bounds, num_digits)
	display( x)
end

"""
func: The function to approximate
bounds: the lower, upper bounds within which to approximate
num_digits: The number of digits of precision

returns: the coefficients of a quadratic function fitted onto *func*
"""
function approx_func( func, bounds, num_digits)
	lower, upper = bounds

	coefs = [
		[ ( upper - lower) ((upper^2 - lower^2) / 2) (( upper^3 - lower^3) / 3)];
		[ ((upper^2 - lower^2) / 2) (( upper^3 - lower^3) / 3) (( upper^4 - lower^4) / 4)];
		[ ((upper^3 - lower^3) / 3) (( upper^4 - lower^4) / 4) (( upper^5 - lower^5) / 5)]]

	right_side = [
		quadgk( func, lower, upper, rtol=1e-5);
		quadgk( t -> t * func(t), lower, upper, rtol=1e-5);
		quadgk( t -> t^2 * func(t), lower, upper, rtol=1e-5)]
	right_side = first.(right_side)

	L,U = lu( coefs, Val(false))

	y = solve_lower( L, right_side)
	x = solve_upper( U, y)

	x
end

function solve_lower( A, b)
	x1 = b[1] / A[1,1]
	x2 = (b[2] - A[2,1]*x1) / A[2,2]
	x3 = (b[3] - A[3,1]*x1 - A[3,2]*x2) / A[3,3]

	[ x1; x2; x3]
end

function solve_upper( A, b)
	x3 = b[3] / A[3,3]
	x2 = (b[2] - A[2,3]*x3) / A[2,2]
	x1 = (b[1] - A[1,3]*x3 - A[1,2]*x2) / A[1,1]

	x = [ x1; x2; x3]
	x = map( x -> round(x; sigdigits = 3), x)
end

main()

max_stars_repo_name: pnadon/least-squares-function-fit

max_issues_repo_name: pnadon/least-squares-function-fit

max_forks_repo_name: pnadon/least-squares-function-fit