text:
##################################
# Crude, i.e., KBS1 implementation
##################################

"""
    forceconfluence!(rws::RewritingSystem, ri, rj,
    [, o::Ordering=ordering(rws)])
Produce (potentially critical) pairs from overlaps of left hand sides of rules
`ri` and `rj`. When failures of local confluence are found, new rules are added
to `rws`.

See [Sims, p. 69].
"""
function forceconfluence!(rws::RewritingSystem, ri, rj, o::Ordering = ordering(rws))
    lhs_i, rhs_i = ri
    lhs_j, rhs_j = rj
    for k in 1:length(lhs_i)
        b = @view lhs_i[end-k+1:end]
        n = longestcommonprefix(b, lhs_j)
        if isone(@view b[n+1:end]) || isone(@view lhs_j[n+1:end])
            a = lhs_i[1:end-k]
            append!(a, rhs_j)
            append!(a, @view b[n+1:end])

            c = rhs_i * @view lhs_j[n+1:end]

            deriverule!(rws, a, c, o)
        end
    end
    return rws
end

##########################
# Naive KBS implementation
##########################

# As of now: default implementation

"""
    forceconfluence!(rws::RewritingSystem, stack, work:kbWork,
        ri, rj[, o::Ordering=ordering(rs)])
Produce (potentially critical) pairs from overlaps of left hand sides of rules
`ri` and `rj`. When failures of local confluence are found, new rules are added
to `rws`.

This version uses `stack` and `work::kbWork` to save allocations and speed-up
the process. See [Sims, p. 77].
"""
function forceconfluence!(
    rws::RewritingSystem{W},
    stack,
    ri,
    rj,
    work::kbWork = kbWork{eltype(W)}(),
    o::Ordering = ordering(rs),
) where {W}
    lhs_i, rhs_i = ri
    lhs_j, rhs_j = rj
    m = min(length(lhs_i), length(lhs_j)) - 1

    for k in 1:m
        if issuffix(@view(lhs_j[1:k]), lhs_i)
            a = store!(work.tmpPair._vWord, @view lhs_i[1:end-k])
            a = append!(a, rhs_j)

            c = store!(work.tmpPair._wWord, rhs_i)
            c = append!(c, @view lhs_j[k+1:end])

            critical, (a, c) = _iscritical(a, c, rws, work)
            if critical
                push!(stack, (a, c))
            end
        end
    end
    deriverule!(rws, stack, work, o)
end

########################################
# KBS using index automata for rewriting
########################################

"""
    forceconfluence!(rws::RewritingSystem, stack, work::kbWork, at::Automaton,
        ri, rj[, o::Ordering=ordering(rws)])
Produce (potentially critical) pairs from overlaps of left hand sides of rules
`ri` and `rj`. When failures of local confluence are found, new rules are added
to `rws`.

This version uses `stack`, `work::kbWork` to save allocations and `at::Automaton`
to speed-up the rewriting process. See [Sims, p. 77].
"""
function forceconfluence!(
    rws::RewritingSystem,
    stack,
    at::Automaton,
    ri,
    rj,
    work::kbWork = kbWork{eltype(W)}(),
    o::Ordering = ordering(rws),
)
    lhs_i, rhs_i = ri
    lhs_j, rhs_j = rj

    m = min(length(lhs_i), length(lhs_j)) - 1

    for k in 1:m
        if issuffix(@view(lhs_j[1:k]), lhs_i)
            a = store!(work.tmpPair._vWord, @view lhs_i[1:end-k])
            a = append!(a, rhs_j)

            c = store!(work.tmpPair._wWord, rhs_i)
            c = append!(c, @view lhs_j[k+1:end])

            critical, (a, c) = _iscritical(a, c, at, work)
            if critical
                push!(stack, (a, c))
            end
        end
    end
    deriverule!(rws, stack, work, at, o)
end


text:
# GMM estimation for a sample from Chi^2(theta)
# compare to two method of moments estimators (see chi2mm.m)
using Econometrics, Random, Distributions, Plots
function main()
n = 30
theta = 3.0
reps = 1000
results = zeros(reps,2)

y = zeros(n) # this is just a place holder to define the moments, below
# define the moment conditions
m = theta -> [(theta.- y) (theta .- 0.5*(y .- theta).^2.0)]
for i = 1:1000
    rand!(Chisq(theta), y) # this replaces the place holder, from above
    thetahat, junk, junk, ms, junk = gmm(m, [3.0], eye(2))
    results[i,1] = sqrt(n)*(thetahat[1]-theta)
    W = inv(cov(ms))
    thetahat, junk, ms, junk = gmm(m, thetahat, W)
    results[i,2] = sqrt(n)*(thetahat[1]-theta)
end
p1 = npdensity(results[:,1])
plot!(p1, title="Inefficient")
p2 = npdensity(results[:,2])
plot!(p2, title="Efficient")
p = plot(p1,p2,layout=(2,1),xlims=(-10,20))
#savefig("Efficient.png")
println("Monte Carlo covariance: (1,1) is inefficient GMM, (2,2) is efficient")
prettyprint(cov(results))
return p
end
main()

ext: jl

max_stars_repo_name: Hiroakiyusheng/Econometrics

max_issues_repo_name: Hiroakiyusheng/Econometrics

max_forks_repo_name: Hiroakiyusheng/Econometrics


text:
#! /usr/bin/env julia

# Copyright 2017 John Hanley.
#
# Permission is hereby granted, free of charge, to any person obtaining a
# copy of this software and associated documentation files (the "Software"),
# to deal in the Software without restriction, including without limitation
# the rights to use, copy, modify, merge, publish, distribute, sublicense,
# and/or sell copies of the Software, and to permit persons to whom the
# Software is furnished to do so, subject to the following conditions:
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
# The software is provided "AS IS", without warranty of any kind, express or
# implied, including but not limited to the warranties of merchantability,
# fitness for a particular purpose and noninfringement. In no event shall
# the authors or copyright holders be liable for any claim, damages or
# other liability, whether in an action of contract, tort or otherwise,
# arising from, out of or in connection with the software or the use or
# other dealings in the software.

# usage:
#     src/store_file_names.jl db file_listing

using SQLite

function create_tables(db)
    create = "create table file_source(pathname  string primary key, src  string not null, pkg string)"
    SQLite.query(db, create)
    tbls = SQLite.tables(db)
    println(tbls)
end

function store(db, listing)
    SQLite.transaction(db)
    ins = SQLite.Stmt(db, "insert into file_source  values (?, 'DEB', ?)")
    SQLite.bind!(ins, 2, pkg_name(listing))
    open(listing) do fin
        for line in eachline(fin)
            fspec = rstrip(line)
            if ! isfile(fspec)
                continue
            end
            SQLite.bind!(ins, 1, fspec)
            SQLite.execute!(ins)
        end
    end
    SQLite.commit(db)
end

# "/a/b/c/def.txt" -> "def"
function pkg_name(listing)
    m = match(r".*/([^/]+)\.txt$", listing)
    return m.captures[1]
end

db_fspec, listing = ARGS
if (!isfile(db_fspec))
    create_tables(SQLite.DB(db_fspec))
end
store(SQLite.DB(db_fspec), listing)


text:
#Pkg.activate(".")
#Pkg.instantiate()

using Documenter
using DocumenterTools: Themes

#Pkg.add(PackageSpec(path="../../QuantumCircuits.jl"))
push!(LOAD_PATH,"../src/")
using QuantumCircuits

makedocs(
    modules = [QuantumCircuits],
    sitename = "QuantumCircuits.jl",
    format = Documenter.HTML(mathengine = Documenter.MathJax(),
                             prettyurls = get(ENV, "CI", nothing) == "true",
                             sidebar_sitename=false
                             ),
    strict = true,
    authors = "Rafał Pracht",
    pages = [
        "Introduction" => "index.md",
        "Quick Start guide" => "quickguide.md",
        "Quantum Gates Library" => [
            "Single-qubit gates s" => "0_Single-qubit-gates.md",
            "Two-qubit gates" => "1_Two-qubit-gates.md",
        ],
        "Novel algorithm to Simulation on NISQ device" => [
            "Problem definition" => "Simulation/problem_definition.md",
            "U4 - Cartan's KAK decomposition" => "Simulation/U4_Cartan_decomposition.md",
            "Algorithm description" => "Simulation/Trotterization.md",
        ],
        "Examples" => [
            "Log-Normal state preparation" => "examples/state_preparation.md",
        ],
        "Library References" => [
            "QCircuits" => [
                "Circuit" => "QCircuits/Circuit.md",
                "Gates" => "QCircuits/Gates.md",
                "Qiskit" => "QCircuits/Qiskit.md",
                "Registers" => "QCircuits/Registers.md",
                "Graph" => "QCircuits/Graph.md",
                "OtherGates" => "QCircuits/OtherGates.md",
                "ComplexGates" => "QCircuits/ComplexGates.md",
                "Instructions" => "QCircuits/Instructions.md",
                "Math" => "QCircuits/Math.md",
                "QBase" => "QCircuits/QBase.md",
            ],
            "Execute" => "Execute.md",
            "QML" => [
                "QML" => "QML/QML.md",
                "Optimization" => "QML/Optimization.md",
                "Gates" => "QML/CircuitOptimization.md",
            ],
        ],
    ],
)
#QuantumCircuits.QCircuits.Gates

deploydocs(
    repo = "github.com/Adgnitio/QuantumCircuits.jl",
    push_preview = true
)


text:
using Printf

function RK2_step!(psi::WfMCTDHF1D, dt::Number; include_kinetic_part::Bool=false,
    include_one_particle_potential_part::Bool=true, freeze_time::Bool=!include_kinetic_part,
    k1::WfMCTDHF1D=psi.m.k1,  k2::WfMCTDHF1D=psi.m.k2 )
    # Usually time is frozen if the kinetic part is not included, because in this case
    # time is propagated by the A (=kinetic) part of the time splitting method.
    m = psi.m
    t = get_time(psi)
    t1 = freeze_time ? t : t+0.5*dt
    t2 = freeze_time ? t : t+dt
    gen_rhs!(k1, psi, include_kinetic_part=include_kinetic_part,
                        include_one_particle_potential_part=include_one_particle_potential_part)
    scale!(k1, 0.5*dt)
    axpy!(k1, psi, 1.0)
    set_time!(k1, t)
    gen_rhs!(k2, k1, include_kinetic_part=include_kinetic_part,
                         include_one_particle_potential_part=include_one_particle_potential_part)
    axpy!(psi, k2, dt)
    set_time!(psi, t2)
    orthonormalize_orbitals!(psi)
end

function RK4_step!(psi::WfMCTDHF1D, dt::Number; include_kinetic_part::Bool=false,
    include_one_particle_potential_part::Bool=true, freeze_time::Bool=!include_kinetic_part,
    k1::WfMCTDHF1D=psi.m.k1,  k2::WfMCTDHF1D=psi.m.k2, k3::WfMCTDHF1D=psi.m.k3,  k4::WfMCTDHF1D=psi.m.k4 )
    # Usually time is frozen if the kinetic part is not included, because in this case
    # time is propagated by the A (=kinetic) part of the time splitting method.
    m = psi.m
    t = get_time(psi)
    t1 = freeze_time ? t : t
    t2 = freeze_time ? t : t+0.5*dt
    t3 = freeze_time ? t : t+0.5*dt
    t4 = freeze_time ? t : t+dt
    #For this storage-efficient implementation, see
    #E. K. Blum: A  Modification of the  Runge-Kutta Fourth-Order Method
    gen_rhs!(k3, psi, include_kinetic_part=include_kinetic_part,
                        include_one_particle_potential_part=include_one_particle_potential_part)
    scale!(k3, dt)
    copy!(k2, k3)
    copy!(k1, psi)
    axpy!(k1, k3, 0.5)
    set_time!(k1, t2)
    gen_rhs!(k3, k1, include_kinetic_part=include_kinetic_part,
                         include_one_particle_potential_part=include_one_particle_potential_part)
    scale!(k3, dt)
    axpy!(k1, k3, 0.5)
    axpy!(k1, k2, -0.5)
    scale!(k2,1.0/6)
    scale!(k3, -0.5)
    set_time!(k1, t3)
    gen_rhs!(k4, k1, include_kinetic_part=include_kinetic_part,
                         include_one_particle_potential_part=include_one_particle_potential_part)
    axpy!(k3, k4, dt)
    axpy!(k1, k3, 1.0)
    axpy!(k2, k3, -1.0)
    scale!(k3, 2.0)
    set_time!(k1, t4)
    gen_rhs!(k4, k1, include_kinetic_part=include_kinetic_part,
                         include_one_particle_potential_part=include_one_particle_potential_part)
    axpy!(k3, k4, dt)
    axpy!(k1, k2, 1.0)
    axpy!(k1, k3, 1.0/6)
    copy!(psi,k1)
    set_time!(psi, t4)
    orthonormalize_orbitals!(psi)
end

function groundstate!(psi::WfMCTDHF1D;  dt::Float64=0.05, max_iter::Int=2000, output_step::Int=20, tol=1e-5,
                      keep_initial_value::Bool=false)
    m = psi.m
    k1 = wave_function(m)
    k2 = wave_function(m)
    psi1 = wave_function(m)

    if !keep_initial_value
        to_frequency_space!(psi)
        nx = get_nx(m.m)
        for k=1:m.N
            u=get_data(psi.o[k].phi,true)
            u[:]=zeros(Complex{Float64}, length(u))
            u[div(k+1,2)]=1
            if k>=3
                u[nx-div(k+1,2)+2]=1
            end
        end
        to_real_space!(psi)
        psi.a[:] = zeros(m.lena)
        for k=1:m.lena #generate singlet state (works only for f=2!!!)
            (i1,i2) = m.slater_indices[k]
            if isodd(i1) && i2==i1+1
                psi.a[k] = 1.0
            end
        end
    end

    orthonormalize_orbitals!(psi)
    psi.a[:] = psi.a[:]/norm(psi.a)
    time0 = time()

    for k=1:max_iter
        if mod(k,output_step)==0
            copy!(psi1, psi)
            imaginary_time_propagate_A!(psi1, 0.25*dt)
            RK2_step!(psi1, -0.5im*dt, k1=k1, k2=k2)
            imaginary_time_propagate_A!(psi1, 0.25*dt)
            normalize!(psi1)
            E1, dev1 = get_energy_expectation_deviation(psi1, tmp=k1)
            err1 = dev1/max(abs(E1),1.0)
        end
        imaginary_time_propagate_A!(psi, 0.5*dt)
        RK2_step!(psi, -1im*dt, k1=k1, k2=k2)
        imaginary_time_propagate_A!(psi, 0.5*dt)
        normalize!(psi)

        if mod(k,output_step)==0
            #E_pot = potential_energy(psi)
            #E_kin = kinetic_energy(psi)
            #E = E_pot + E_kin
            #ctime = time() - time0
            #@printf("step=%4i  E_pot=%14.10f  E_kin=%14.10f  E=%14.10f ctime=%10.2f\n", k, E_pot, E_kin, E, ctime)
            E, dev = get_energy_expectation_deviation(psi, tmp=k1)
            err = dev/max(abs(E),1.0)
            ctime = time() - time0
            @printf("step=%4i  E=%14.10f  err=%12.3e  E1=%14.10f  err1=%12.3e  ctime=%10.2f\n", k, E, err, E1, err1, ctime)
            if err1<err
               @printf("changed step size, old:%24.15e  new:%24.15e\n", dt, 0.5*dt)
               dt = 0.5*dt
               copy!(psi,psi1)
            end
            if (err<tol)||(err1<tol)
                break;
            end
        end

    end

end

function strang_step!(psi::WfMCTDHF1D, dt::Real; include_kinetic_part_in_B::Bool=false,
    include_one_particle_potential_part_in_B::Bool=true)
    if !include_kinetic_part_in_B
        propagate_A!(psi, 0.25*dt)
    end
    if !include_one_particle_potential_part_in_B
        propagate_B!(psi, 0.5*dt)
    end
    if !include_kinetic_part_in_B
        propagate_A!(psi, 0.25*dt)
    end
    RK2_step!(psi, dt, include_kinetic_part=include_kinetic_part_in_B,
        include_one_particle_potential_part=include_one_particle_potential_part_in_B)
    if !include_kinetic_part_in_B
        propagate_A!(psi, 0.25*dt)
    end
    if !include_one_particle_potential_part_in_B
        propagate_B!(psi, 0.5*dt)
    end
    if !include_kinetic_part_in_B
        propagate_A!(psi, 0.25*dt)
    end
end

function run!(psi::WfMCTDHF1D, dt::Real, n::Int, method::Function; output_step::Int=1)
    m = psi.m
    m.k1 = wave_function(m)
    m.k2 = wave_function(m)
    m.k3 = wave_function(m)
    m.k4 = wave_function(m)
    time0 = time()
    set_propagate_time_together_with_A!(m, true)

    orthonormalize_orbitals!(psi)

    for k=1:n
        method(psi, dt)

        if mod(k,output_step)==0
            t = get_time(psi)
            nn = norm(psi)
            E_pot = potential_energy(psi)
            E_kin = kinetic_energy(psi)
            E = E_pot + E_kin
            ctime = time() - time0
            @printf("step=%5i  t=%14.10f  norm=%14.10f  E_pot=%14.10f  E_kin=%14.10f  E=%14.10f  ctime=%10.2f\n", k, t, nn, E_pot, E_kin, E, ctime)
        end
    end

    m.k1 = nothing
    m.k2 = nothing
    m.k3 = nothing
    m.k4 = nothing
end

function local_orders(psi::WfMCTDHF1D, dt::Real, method::Function; method_ref::Function=method,
                      reference_steps=10,
                      rows=8)
    m = psi.m
    m.k1 = wave_function(m)
    m.k2 = wave_function(m)
    m.k3 = wave_function(m)
    m.k4 = wave_function(m)
    tab = Array(Float64, rows, 7)
    set_propagate_time_together_with_A!(m, true)
    wf_save_initial_value = wave_function(m)
    psi_ref = wave_function(m)
    copy!(wf_save_initial_value, psi)

    dt1 = dt
    err_a_old = 0.0
    err_phi_old = 0.0
    err_old = 0.0
    println("             dt       err_a      p      err_phi      p          err      p")
    println("----------------------------------------------------------------------------")
    for row=1:rows
        method(psi, dt1)
        copy!(psi_ref,wf_save_initial_value)
        for k=1:reference_steps
            method_ref(psi_ref, dt1/reference_steps)
        end
        err_phi = -1.0
        for j=1:m.N
            err_phi=max(err_phi, TSSM.distance(psi.o[j].phi, psi_ref.o[j].phi))
        end
        err_a = norm(psi.a-psi_ref.a)
        err = abs(distance(psi, psi_ref))
        if (row==1) then
            @printf("%3i%12.3e%12.3e %19.3e %19.3e\n", row, Float64(dt1), Float64(err_a),Float64(err_phi),Float64(err))
            tab[row,1] = dt1
            tab[row,2] = err_a
            tab[row,3] = 0
            tab[row,4] = err_phi
            tab[row,5] = 0
            tab[row,6] = err
            tab[row,7] = 0
        else
            p_a = log(err_a_old/err_a)/log(2.0);
            p_phi = log(err_phi_old/err_phi)/log(2.0);
            p = log(err_old/err)/log(2.0);
            @printf("%3i%12.3e%12.3e%7.2f %12.3e%7.2f %12.3e%7.2f\n", row, Float64(dt1),
                Float64(err_a), Float64(p_a), Float64(err_phi), Float64(p_phi), Float64(err), Float64(p))
            tab[row,1] = dt1
            tab[row,2] = err_a
            tab[row,3] = p_a
            tab[row,4] = err_phi
            tab[row,5] = p_phi
            tab[row,6] = err
            tab[row,7] = p
        end
        err_a_old = err_a
        err_phi_old = err_phi
        err_old = err
        dt1 = 0.5*dt1
        copy!(psi,wf_save_initial_value)
    end
    m.k1 = nothing
    m.k2 = nothing
    m.k3 = nothing
    m.k4 = nothing
    tab
end

function run!(psi::WfMCTDHF1D, dt::Real, steps::Int, a::Vector{Float64}, b::Vector{Float64};
    method::Symbol=:RK2, iters::Int=0, output_step::Int=1, return_solutions::Bool=false)
    time0 = time()
    s = length(a)
    m = psi.m
    if method==:RK2 || method==:RK4
        m.k1 = wave_function(m)
        m.k2 = wave_function(m)
    end
    if method==:RK4
        m.k3 = wave_function(m)
        m.k4 = wave_function(m)
    end
    if iters>0
        psi0 = wave_function(m)
    end
    if return_solutions
        psi_back = WaveFunction[wave_function(m) for j=1:steps+1]
        copy!(psi_back[1],psi)
    end
    set_propagate_time_together_with_A!(m, true)

    orthonormalize_orbitals!(psi)

    for step=1:steps
        for j = 1:s
            if a[j]!=0.0
                propagate_A!(psi, a[j]*dt)
            end
            if b[j]!=0.0
                if iters>0
                    copy!(psi0, psi)
                end
                if method == :RK2
                    RK2_step!(psi, b[j]*dt)
                elseif method ==:RK4
                    RK4_step!(psi, b[j]*dt)
                end

                for it=1:iters
                    copy!(m.k1, psi0)
                    scale!(m.k1, 0.5)
                    axpy!(m.k1, psi, 0.5)
                    copy!(psi, psi0)
                    set_time!(m.k1, get_time(psi))
                    gen_rhs!(m.k2, m.k1)
                    axpy!(psi, m.k2, b[j]*dt)
                end
            end
            if return_solutions
                copy!(psi_back[step+1], psi)
            end
        end

        if mod(step,output_step)==0
            t = get_time(psi)
            nn = norm(psi)
            E_pot = potential_energy(psi)
            E_kin = kinetic_energy(psi)
            E = E_pot + E_kin
            ctime = time() - time0
            @printf("step=%5i  t=%14.10f  norm=%14.10f  E_pot=%14.10f  E_kin=%14.10f  E=%14.10f  ctime=%10.2f\n",
                     step, t, nn, E_pot, E_kin, E, ctime)
        end
    end

    m.k1 = nothing
    m.k2 = nothing
    m.k3 = nothing
    m.k4 = nothing
    if return_solutions
        return psi_back
    else
        return nothing
    end
end


text:
__precompile__()

module Quaternions

  import Base: +, -, *, /, ^
  import Base: abs, abs2, angle, conj, cos, exp, inv, isfinite, log, real, sin, sqrt
  import Base: convert, promote_rule, float
  import LinearAlgebra: norm, normalize
  import Random: rand, SamplerType, AbstractRNG

  include("Quaternion.jl")
  include("Octonion.jl")
  include("DualQuaternion.jl")

  export Quaternion
  export quat
  export Octonion
  export octo
  export DualQuaternion
  export dualquat
  export angleaxis
  export angle
  export axis
  export linpol
  export normalize
  export normalizea
  export dconj
  export rand
  export quatrand
  export nquatrand
  export octorand
  export dualquatrand
  export ndualquatrand
  export qrotation
  export rotationmatrix
  export slerp
end


text:
# This file contains general algorithms inter alia from chapter 1
# of the book
#
#    Manuel Bronstein, Symbolic Integration I: Transcendental Functions,
#    2nd ed., Springer 2005.
#

using AbstractAlgebra
using Nemo

struct NotImplementedError <: Exception
    msg::String
end

struct AlgorithmFailedError <: Exception
    msg::String
end

abstract type Term end

struct IdTerm <: Term
    arg::RingElement
end

Base.show(io::IO, t::IdTerm) = show(io, t.arg)

struct FunctionTerm <: Term
    op # ::Function
    coeff::RingElement
    arg::RingElement
end

function Base.show(io::IO, t::FunctionTerm)
    if iszero(t.coeff)
        print("0")
        return
    end
    if iszero(t.coeff+one(t.coeff))
        print(io, "-")
    elseif !isone(t.coeff)
        s = string(t.coeff)
        if findnext("+", s, 2)!==nothing || findnext("-", s, 2)!==nothing
            print(io, "(", s, ")")
        else
            print(io, s)
        end
        print(io,*)
    end
    print(io, t.op, "(")
    show(io, t.arg)
    print(io, ")")
end

struct Result
    t::Vector{Term}
 end

 function Base.show(io::IO, r::Result)
     if length(r.t)==0
         print(io, "0")
         return
     end
     show(io, r.t[1])
     for i=2:length(r.t)
         s = string(r.t[i])
         if s[1]!='+' && s[1]!='-'
             print(io, " + ", s)
         else
             print(io, " ", s[1], " ", s[2:end])
         end
     end
 end

isrational(x::T) where T<:Integer = true

isrational(x::T) where T<:Rational = true

isrational(x::fmpq) = true # Nemo rational type

isrational(x::qqbar) = degree(x)==1 && iszero(imag(x)) # Nemo algebraic number type

function isrational(x::P) where P<:PolyElem
    if degree(x)>0
        return false
    else
        return isrational(constant_coefficient(x))
    end
end

function isrational(x::F) where {P<:PolyElem, F<:FracElem{P}}
    isone(denominator(x)) && isrational(numerator(x))
end

rationalize(x::T) where T<:Integer = convert(Rational{BigInt}, x)

rationalize(x::T) where T<:Rational = convert(Rational{BigInt}, x)

rationalize(x::fmpq) = convert(Rational{BigInt}, x) # Nemo rational type

function rationalize(x::qqbar) #Nemo algebraic number type
    (degree(x)==1 && iszero(imag(x))) || error("not rational")
    convert(Rational{BigInt}, fmpq(x))
end

function rationalize(x::P) where P<:PolyElem
    degree(x)<=0 || error("not rational")
    rationalize(constant_coefficient(x))
end

function rationalize(x::F) where  {P<:PolyElem, F<:FracElem{P}}
    isone(denominator(x)) || error("not rational")
    rationalize(numerator(x))
end

rationalize_over_Int(x) = convert(Rational{Int}, rationalize(x))

function common_leading_coeffs(fs::Vector{F}) where F<:FieldElem
    fs
end

function common_leading_coeffs(fs::Vector{F}) where {T<:FieldElement, P<:PolyElem{T}, F<:FracElem{P}}
    l = lcm([denominator(f) for f in fs])
    Fs = [divexact(l*numerator(f), denominator(f)) for f in fs]
    d = maximum([degree(F) for F in Fs])
    l1 = leading_coefficient(l)
    common_leading_coeffs([coeff(F, d)//l1 for F in Fs])
end

function constant_factors(f::PolyElem{T}) where T<:FieldElement
    f0 = parent(f)(common_leading_coeffs(collect(coefficients(f))))
    gcd(f, f0)
end

function rational_roots(f::PolyElem{T}) where T<:FieldElement
    p = map_coefficients(c->fmpq(rationalize(c)), constant_factors(f)) # fmpq needs Nemo
    roots(p)
end

function Nemo.roots(f::PolyElem{qqbar})
    n = degree(f)
    X = gen(parent(f))
    _, xys= PolynomialRing(Nemo.QQBar, vcat(:x,  [Symbol("y$i") for i = 1:n]))
    x = xys[1]
    ys = xys[2:end]

    as = reverse(collect(coefficients(f))[1:end-1])

    G = x^n + sum([ys[i]*x^(n-i) for i=1:n])
    for i=1:n
        G = prod([ G(x, vcat(zeros(Int, i-1), α, ys[i+1:end])...) for α in conjugates(as[i])])
    end

    g = map_coefficients(c->fmpq(c), G(X, zeros(parent(X), n)...))

    rs = roots(g, QQBar)
    [r for r in rs if iszero(f(r))]
end

valuation_infinity(f::F) where {T<:RingElement, P<:PolyElem{T}, F<:FracElem{P}} =
    # See Bronstein's book, Definition 4.3.1, p. 115
    degree(denominator(f)) - degree(numerator(f))

function Remainder(x::FracElem{T}, a::T) where T<:RingElement
    # See Bronstein's book, Section 4.2, p. 115
    b, c = gcdx(a, denominator(x), one(a))
    q, r = divrem(numerator(x)*c, a)
    r
end

function Base.lcm(as::Vector{T}) where T<:RingElement
    m = length(as)
    if m==0
        error("array as must not be empty")
    elseif m==1
        return as[1]
    else
        y = as[1]
        for i=2:m
            y = lcm(y, as[i])
        end
        return y
    end
end

#PolyDivide: Implemented in AbstractAlgebra.jl as divrem

#PolyPseudoDivide: Implemented in AbstractAlgebra.jl as pseudodivrem

#ExtendedEuclidean: implemented in AbstractAlgebra.jl as gcdx
#  Note that the order of output parameters is g, s, t instead of s, t, g.

EuclideanSize(f::T) where T <: Integer = abs(f)
EuclideanSize(f::T) where T <: PolyElem = degree(f)

function Base.gcdx(a::T, b::T, c::T) where T <: RingElem
    # ExtendedEuclidean - diophantine version
    # See Bronstein's book, Section 1.3, p. 13
    R = parent(a)
    (parent(b) == R && parent(c) == R) || error("Incompatible parents")
    g, s, t = gcdx(a, b)
    q, r = divrem(c, g)
    iszero(0) || error("c is not in the ideal generated by a and b")
    s = q*s
    t = q*t
    if !iszero(s) && EuclideanSize(s)>=EuclideanSize(b)
        q, r = divrem(s, b)
        s = r
        t = t + q*a
    end
    s, t
end

function PartialFraction(a::T, d::Vector{T}) where T <: RingElem
    # See Bronstein's book, Section 1.3, p. 15
    n = length(d)
    if n==1
        q, r = divrem(a, d[1])
        return q, [r]
    end
    prod_d2n = prod(d[2:n])
    a0, r = divrem(a, d[1]*prod_d2n)
    a1, t = gcdx(prod_d2n, d[1], r)
    b0, a2n = PartialFraction(t, d[2:n])
    return a0+b0, vcat(a1, a2n)
end

function PartialFraction(a::T, d::Vector{T}, e::Vector{Int}) where T <: RingElem
    # See Bronstein's book, Section 1.3, p. 17
    n = length(d)
    a0, aa = PartialFraction(a, [d[i]^e[i] for i=1:n])
    A = [zeros(parent(a), e[i]) for i=1:n]
    for i=1:n
        for j=e[i]:-1:1
            q, A[i][j] = divrem(aa[i], d[i])
            aa[i] = q
        end
        a0 += aa[i]
    end
    a0, A
end

function SubResultant(A::PolyElem{T}, B::PolyElem{T}) where T <: RingElement
    # See Bronstein's book, Section 1.5, p. 24
    degree(A) >= degree(B) || error("degree(A) must be >= degree(B)")
    T_one = one(leading_coefficient(A)) # 1 of type T
    Rs = [A, B] # indices shifted!
    i = 1
    γ = T[-T_one]
    δ = Int[degree(A)-degree(B)]
    β = T[(-T_one)^(δ[1]+1)]
    r = T[]
    while !iszero(Rs[i+1])
        push!(r, leading_coefficient(Rs[i+1]) )
        Q, R = pseudodivrem(Rs[i-1+1], Rs[i+1])
        q = divexact(R, β[i])
        push!(Rs, q)
        i += 1
        h = δ[i-1]<=1 ? (-r[i-1])^δ[i-1]*γ[i-1]^(1-δ[i-1]) : divexact((-r[i-1])^δ[i-1], γ[i-1]^(δ[i-1]-1))
        push!(γ, h)
        push!(δ, degree(Rs[i-1+1]) - degree(Rs[i+1]) )
        push!(β, -r[i-1]*γ[i]^δ[i])
    end
    k = i - 1
    zero_poly = zero(A)
    if degree(Rs[k+1])>0
        return zero_poly, vcat(Rs[1:k+1], zero_poly)
    end
    if degree(Rs[k-1+1])==1
        return constant_coefficient(Rs[k+1]), vcat(Rs[1:k+1], zero_poly)
    end
    s = 1
    c_num = T_one
    c_den = T_one
    for j=1:k-1
        if isodd(degree(Rs[j-1+1])) && isodd(degree(Rs[j+1]))
            s = -s
        end
        d = -(1+δ[j])*degree(Rs[j+1]) + degree(Rs[j-1+1]) - degree(Rs[j+1+1])
        c_num *= β[j]^degree(Rs[j+1])
        if d>=0
            c_num *= r[j]^d
        else
            c_den *= r[j]^(-d)
        end
    end
    constant_coefficient(divexact(s*c_num*Rs[k+1]^degree(Rs[k-1+1]), c_den)), vcat(Rs[1:k+1], zero_poly)
end

function Squarefree_Musser(A::PolyElem{T}) where T <: RingElement
    # See Bronstein's book, Section 1.7, p. 29
    c = content(A)
    S = divexact(A, c)
    Sminus = gcd(S, derivative(S))
    Sstar = divexact(S, Sminus)
    As = PolyElem{T}[]
    k = 1
    while degree(Sminus)>0
        Y = gcd(Sstar, Sminus)
        push!(As, divexact(Sstar, Y))
        Sstar = Y
        Sminus = divexact(Sminus, Y)
        k += 1
    end
    push!(As, Sstar)
    As[1] = (c*Sminus)*As[1]
    As
end

function Squarefree(A::PolyElem{T}) where T <: RingElement
    # See Bronstein's book, Section 1.7, p. 32
    c = content(A)
    S = divexact(A, c)
    Sprime = derivative(S)
    Sminus = gcd(S, Sprime)
    Sstar = divexact(S, Sminus)
    As = PolyElem{T}[]
    k = 1
    Y = divexact(Sprime, Sminus)
    Z = Y - derivative(Sstar)
    while !iszero(Z)
        push!(As, gcd(Sstar, Z))
        Sstar = divexact(Sstar, As[k])
        Y = divexact(Z, As[k])
        Z = Y - derivative(Sstar)
        k += 1
    end
    push!(As, Sstar)
    As[1] = c*As[1]
    As
end


text:
function write_stl(dataset::AbstractStaticVTKData, filepath_no_ext::String)
    polydata = PyVTK(trianulate(dataset))
    return write_stl(polydata, filepath_no_ext)
end

function write_stl(polydata::PyCall.PyObject, filepath_no_ext::String)
    writer = vtk.vtkSTLWriter()
    writer[:SetFileName](string(filepath_no_ext, ".stl"))
    writer[:SetInputData](polydata)
    writer[:Update]()
    return string(filepath_no_ext, ".stl")
end

function read_stl(filepath::String)
    if splitext(filepath)[2] == ".stl"
        reader = vtk.vtkSTLReader()
        reader[:SetFileName](filepath)
        reader[:Update]()

        polydata = reader[:GetOutput]()
        return extract_simple_block(polydata)
    else
        throw("The file input is not an stl file.")
    end
end


text:
function tag(x::Pair, attr::Pair...)
    t, b = x
    attrstr = join(" $n=\"$p\"" for (n, p) in attr)
    return "<$t$attrstr>$b</$t>"
end

colnames = split(",X,Y,Z", ',')

header = join(tag(:th => txt) for txt in colnames) * "\n"
rows   = collect(tag(:tr => join(tag(:td => i, :style => "font-weight: bold;") * join(tag(:td => rand(1000:9999)) for j in 1:3))) for i in 1:6)
body   = "\n" * join(rows, '\n') * "\n"
table  = tag(:table => string('\n', header, body, '\n'), :style => "width: 60%")
println(table)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
push!(LOAD_PATH, abspath(@__DIR__,".."))

using LinearAlgebra,Statistics,StatsBase,Distributions
using Plots,NamedColors ; theme(:dark) #; plotlyjs();
using SparseArrays
using SpikingRNNs; const global S = SpikingRNNs
using BenchmarkTools
using ProgressMeter
import FileIO: save # save raster png

using InvertedIndices

using Random; Random.seed!(0)

##

# let's make a network with strong E-E connections that explodes

# define E population of conductance based EIF neurons
const dt = 0.1E-3
const Ne = 500

myτe = 20E-3 # seconds
τrefr= 5E-3 # refractoriness
vth_e = 20.   # mV
myτe_ker = 15E-3 # synaptic kernel
vthexp = -50.0 # actual threshold for spike-generation
eifslope = 2.0
v_rest_e = -70.0
v_rev_e = 0.0
v_leak_e = v_rest_e
v_reset_e = v_rest_e
Cap = 300.0 #capacitance mF

nt_e,ps_e = let sker = S.SKExp(myτe_ker)
  sgen = S.SpikeGenEIF(vthexp,eifslope)
  nt = S.NTLIFConductance(sker,sgen,myτe,Cap,
    vth_e,v_reset_e,v_rest_e,τrefr,v_rev_e)
  ps = S.PSLIFConductance(nt,Ne)
  (nt,ps)
end

## define E to E connections
const jee = 4.0
const pee = 0.33333
wmat = S.sparse_constant_wmat(Ne,Ne,pee,jee;no_autapses=true)
conn_e_e = S.ConnGeneralIF2(wmat)

## define external input of Poisson neurons
const ext_rate = 80.0
const je_in = 60.0
nt_in,ps_in = let spkgen = S.SGPoisson(ext_rate)
  sker = S.SKExp(myτe_ker)
  nt = S.NTInputConductance(spkgen,sker,v_rev_e)
  ps = S.PSInputPoissonConductance(nt,je_in,Ne)
  (nt,ps)
end

## Connek !
pop_e = S.Population(ps_e,(conn_e_e,ps_e),(S.FakeConnection(),ps_in))
myntw = S.RecurrentNetwork(dt,pop_e)

## run for a very short time, and check stability
const Ttot = 0.5
const krec = 1
rec_state_e = S.RecStateNow(ps_e,krec,dt,Ttot;idx_save=collect(1:100))
rec_spikes = S.RecSpikes(ps_e,1000.0,Ttot;idx_save=collect(1:100))

## Run
times = (0:myntw.dt:Ttot)
nt = length(times)
# clean up
S.reset!.([rec_spikes,rec_state_e])
S.reset!(ps_e)
S.reset!(conn_e_e)
# initial conditions
ps_e.state_now .= rand(Uniform(v_reset_e,-50),Ne)

@time begin
  @showprogress 1.0 "network simulation " for (k,t) in enumerate(times)
    rec_spikes(t,k,myntw)
    rec_state_e(t,k,myntw)
    S.dynamics_step!(t,myntw)
  end
end

S.add_fake_spikes!(vth_e,rec_state_e,rec_spikes)
## plot a raster
therast1 = S.raster_png(1E-3,rec_spikes;
  spike_height=3,Ttot=Ttot,Nneurons=length(rec_spikes.idx_save));
save("/tmp/tmp.png",therast1)

## plot the state for one neuron

_ = let neu = 1,
  plt=plot()
  plot!(plt,rec_state_e.times,rec_state_e.state_now[neu,:];linewidth=2,leg=false)
end

# explodes quite badly!

## now, let's add global inhibitory stabilization (aka cheating fake stabilization)
ps_istab = let v_rev_i = -70.0
  Aloc = 500.0
  Aglo = 1000.0
  τloc = 50E-3
  τglo = 5E-3
  S.PSConductanceInputInhibitionStabilization(v_rev_i,Aglo,Aloc,τglo,τloc,Ne)
end

## connect with stabilizer
pop_e_stab = S.Population(ps_e,
  (conn_e_e,ps_e),
  (S.FakeConnection(),ps_in),
  (S.FakeConnection(),ps_istab))

myntw_stab = S.RecurrentNetwork(dt,pop_e_stab)

const Ttot = 5.0
const krec = 1
rec_state_e_stab = S.RecStateNow(ps_e,krec,dt,Ttot;idx_save=collect(1:100),
  t_warmup=4.0)
rec_spikes_stab = S.RecSpikes(ps_e,1000.0,Ttot;idx_save=collect(1:100),
  t_warmup=4.0)

## Run
times = (0:myntw.dt:Ttot)
nt = length(times)
# clean up
S.reset!.([rec_spikes_stab,rec_state_e_stab])
S.reset!.([ps_e,ps_in,ps_istab])
S.reset!(conn_e_e)
# initial conditions
ps_e.state_now .= rand(Uniform(v_reset_e,-60),Ne)

@time begin
  @showprogress 1.0 "network simulation " for (k,t) in enumerate(times)
    rec_spikes_stab(t,k,myntw_stab)
    rec_state_e_stab(t,k,myntw_stab)
    S.dynamics_step!(t,myntw_stab)
  end
end

S.add_fake_spikes!(vth_e,rec_state_e_stab,rec_spikes_stab)
## plot a raster
therast2 = S.raster_png(1E-3,rec_spikes_stab;spike_height=3,
  Nneurons=100,Ttot=Ttot)
save("/tmp/tmp2.png",therast2)

_ = let neu = 1,
  plt=plot()
  plot!(plt,rec_state_e_stab.times,rec_state_e_stab.state_now[neu,:];linewidth=1.5,leg=false)
end

# wow, so dynamical stability!  Much asynchronous irregular!


text:
export TreeGradient, TreeGradientNode

mutable struct TreeGradientNode{T}
   ∇::T
   isinitialized::Bool
end

mutable struct TreeGradient{T}
    ∇::T        # overall gradient (sum of all gradients in t)
    n::Int      # total number of samples
    ninit::Int  # number of initialized samples
    t::DataStructures.BalancedTree23{Int,TreeGradientNode{T},Base.Order.ForwardOrdering}
end

function TreeGradient(∇::T, n::Integer; isempty::Bool=false) where T
    t = DataStructures.BalancedTree23{Int,TreeGradientNode{T},Base.Order.ForwardOrdering}(Base.Order.ForwardOrdering())
    insert!(t, 1, TreeGradientNode(∇, !isempty), false)
    ninit = isempty ? 0 : n
    TreeGradient(copy(∇), n, ninit, t)
end

Base.size(tg::TreeGradient) = size(tg.∇)

function Base.iterate(t::DataStructures.BalancedTree23, state=DataStructures.beginloc(t))
    if state < 3
        return nothing
    end
    (t.data[state].k, t.data[state].d), DataStructures.nextloc0(t, state)
end

Base.iterate(tg::TreeGradient) = iterate(tg.t)
Base.iterate(tg::TreeGradient, state) = iterate(tg.t, state)

function Base.collect(t::DataStructures.BalancedTree23{Tk,Td}) where {Tk,Td}
    rv = Vector{Tuple{Tk,Td}}()
    for v in t
        push!(rv, v)
    end
    rv
end

Base.collect(tg::TreeGradient) = collect(tg.t)

initialized_fraction(tg::TreeGradient) = tg.ninit / tg.n

"""

Return the number of samples that make up the sub-gradient at `tg.t.data[i]`.
"""
function subgradient_nsamples(tg::TreeGradient, i::Integer)
    first_sample = tg.t.data[i].k
    next = DataStructures.nextloc0(tg.t, i)
    final_sample = next == 2 ? tg.n : tg.t.data[next].k-1
    final_sample - first_sample + 1
end

"""

Insert the gradient with respect to samples `i` through `j`, with value `∇i`, into `t`. If
`isempty` is `false`, the inserted gradient is counted towards the number of initialized samples.
"""
function Base.insert!(tg::TreeGradient{T}, i::Integer, j::Integer, ∇i::T; isempty::Bool=false) where T
    0 < i <= j <= tg.n || throw(ArgumentError("invalid (i, j): $((i, j))"))
    size(∇i) == size(tg) || throw(DimensionMismatch("∇i has dimensions $(size(∇i)), but existing data has dimensions $(size(tg))"))
    prev, _ = DataStructures.findkey(tg.t, i)
    prev_key = tg.t.data[prev].k
    prev_val = tg.t.data[prev].d
    k, _ = DataStructures.findkey(tg.t, j)
    next = DataStructures.nextloc0(tg.t, k)
    next_key = next == 2 ? tg.n+1 : tg.t.data[next].k

    # take out the previous sub-gradient
    if prev_val.isinitialized
        tg.∇ .-= prev_val.∇
        tg.ninit -= subgradient_nsamples(tg, prev)
    end

    # add the new sub-gradient
    if !isempty
        tg.∇ .+= ∇i
        tg.ninit += j - i + 1
    end

    # delete any intermediate nodes
    k = DataStructures.nextloc0(tg.t, prev)
    key = k == 2 ? tg.n+1 : tg.t.data[k].k
    ∇r1, ∇r2 = ∇i, ∇i # store pointers to data allocated to deleted nodes, so that we can recycle the memory
    while key < next_key
        node = tg.t.data[k]
        ∇r1 = ∇r2
        ∇r2 = node.d.∇
        if node.d.isinitialized
            tg.∇ .-= node.d.∇
            tg.ninit -= subgradient_nsamples(tg, k)
        end
        delete!(tg.t, k)

        k = DataStructures.nextloc0(tg.t, prev)
        key = k == 2 ? tg.n+1 : tg.t.data[k].k
    end

    # either zero out the node to the left and insert a new node, or, if i matches the start of an
    # existing node, update that node in-place (to avoid mutating the tree unnecessarily)
    if prev_key < i
        prev_val.∇ .= 0
        prev_val.isinitialized = false

        # if a node was deleted, we can re-use the memory allocated for it
        if !(∇r2 === ∇i)
            ∇r2 .= ∇i
        else
            ∇r2 = copy(∇i)
        end
        insert!(tg.t, i, TreeGradientNode(∇r2, !isempty), false)
    else
        prev_val.∇ .= ∇i
        prev_val.isinitialized = !isempty
    end
    if j < next_key - 1

        # if at least 2 nodes were deleted, we can re-use memory also here
        if !(∇r1 === ∇i)
            ∇r1 .= 0
        else
            ∇r1 = zero(∇i)
        end
        insert!(tg.t, j+1, TreeGradientNode(∇r1, false), false)
    end
    DataStructures.findkeyless(tg.t, i)
end


text:
using Oceananigans: fields

"""
    RungeKutta3TimeStepper{FT, TG} <: AbstractTimeStepper

Holds parameters and tendency fields for a low storage, third-order Runge-Kutta-Wray
time-stepping scheme described by Le and Moin (1991).
"""
struct RungeKutta3TimeStepper{FT, TG, TI} <: AbstractTimeStepper
                 γ¹ :: FT
                 γ² :: FT
                 γ³ :: FT
                 ζ² :: FT
                 ζ³ :: FT
                 Gⁿ :: TG
                 G⁻ :: TG
    implicit_solver :: TI
end

"""
    RungeKutta3TimeStepper(arch, grid, tracers,
                           Gⁿ = TendencyFields(arch, grid, tracers),
                           G⁻ = TendencyFields(arch, grid, tracers))

Return an `RungeKutta3TimeStepper` object with tendency fields on `arch` and
`grid`. The tendency fields can be specified via optional kwargs.
"""
function RungeKutta3TimeStepper(arch, grid, tracers;
                                implicit_solver::TI = nothing,
                                Gⁿ::TG = TendencyFields(arch, grid, tracers),
                                G⁻ = TendencyFields(arch, grid, tracers)) where {TI, TG}

    !isnothing(implicit_solver) &&
        @warn("Implicit-explicit time-stepping with RungeKutta3TimeStepper is not tested. " *
              "\n implicit_solver: $(typeof(implicit_solver))")

    γ¹ = 8 // 15
    γ² = 5 // 12
    γ³ = 3 // 4

    ζ² = -17 // 60
    ζ³ = -5 // 12

    FT = eltype(grid)

    return RungeKutta3TimeStepper{FT, TG, TI}(γ¹, γ², γ³, ζ², ζ³, Gⁿ, G⁻, implicit_solver)
end

#####
##### Time steppping
#####

"""
    time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt; euler=false)

Step forward `model` one time step `Δt` with a 3rd-order Runge-Kutta method.
The 3rd-order Runge-Kutta method takes three intermediate substep stages to
achieve a single timestep. A pressure correction step is applied at each intermediate
stage.
"""
function time_step!(model::AbstractModel{<:RungeKutta3TimeStepper}, Δt)
    Δt == 0 && @warn "Δt == 0 may cause model blowup!"

    # Be paranoid and update state at iteration 0, in case run! is not used:
    model.clock.iteration == 0 && update_state!(model)

    γ¹ = model.timestepper.γ¹
    γ² = model.timestepper.γ²
    γ³ = model.timestepper.γ³

    ζ² = model.timestepper.ζ²
    ζ³ = model.timestepper.ζ³

    first_stage_Δt  = γ¹ * Δt
    second_stage_Δt = (γ² + ζ²) * Δt
    third_stage_Δt  = (γ³ + ζ³) * Δt

    #
    # First stage
    #

    calculate_tendencies!(model)

    correct_immersed_tendencies!(model, Δt, γ¹, 0)

    rk3_substep!(model, Δt, γ¹, nothing)

    calculate_pressure_correction!(model, first_stage_Δt)
    pressure_correct_velocities!(model, first_stage_Δt)

    tick!(model.clock, first_stage_Δt; stage=true)
    store_tendencies!(model)
    update_state!(model)
    update_particle_properties!(model, first_stage_Δt)

    #
    # Second stage
    #

    calculate_tendencies!(model)

    correct_immersed_tendencies!(model, Δt, γ², ζ²)

    rk3_substep!(model, Δt, γ², ζ²)

    calculate_pressure_correction!(model, second_stage_Δt)
    pressure_correct_velocities!(model, second_stage_Δt)

    tick!(model.clock, second_stage_Δt; stage=true)
    store_tendencies!(model)
    update_state!(model)
    update_particle_properties!(model, second_stage_Δt)

    #
    # Third stage
    #

    calculate_tendencies!(model)

    correct_immersed_tendencies!(model, Δt, γ³, ζ³)

    rk3_substep!(model, Δt, γ³, ζ³)

    calculate_pressure_correction!(model, third_stage_Δt)
    pressure_correct_velocities!(model, third_stage_Δt)

    tick!(model.clock, third_stage_Δt)
    update_state!(model)
    update_particle_properties!(model, third_stage_Δt)

    return nothing
end

#####
##### Time stepping in each substep
#####

stage_Δt(Δt, γⁿ, ζⁿ) = Δt * (γⁿ + ζⁿ)
stage_Δt(Δt, γⁿ, ::Nothing) = Δt * γⁿ

function rk3_substep!(model, Δt, γⁿ, ζⁿ)

    workgroup, worksize = work_layout(model.grid, :xyz)

    barrier = Event(device(model.architecture))

    substep_field_kernel! = rk3_substep_field!(device(model.architecture), workgroup, worksize)

    model_fields = prognostic_fields(model)

    events = []

    for (i, field) in enumerate(model_fields)

        field_event = substep_field_kernel!(field, Δt, γⁿ, ζⁿ,
                                            model.timestepper.Gⁿ[i],
                                            model.timestepper.G⁻[i],
                                            dependencies=barrier)

        # TODO: function tracer_index(model, field_index) = field_index - 3, etc...
        tracer_index = i - 3 # assumption

        implicit_step!(field,
                       model.timestepper.implicit_solver,
                       model.clock,
                       stage_Δt(Δt, γⁿ, ζⁿ),
                       model.closure,
                       tracer_index,
                       model.diffusivity_fields,
                       model.tracers,
                       dependencies = field_event)

        push!(events, field_event)
    end

    wait(device(model.architecture), MultiEvent(Tuple(events)))

    return nothing
end

"""
Time step fields via the 3rd-order Runge-Kutta method

    `U^{m+1} = U^m + Δt (γⁿ G^{m} + ζⁿ G^{m-1})`,

where `m` denotes the substage.
"""

"""
Time step velocity fields with a 3rd-order Runge-Kutta method.
"""
@kernel function rk3_substep_field!(U, Δt, γⁿ, ζⁿ, Gⁿ, G⁻)
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        U[i, j, k] += Δt * (γⁿ * Gⁿ[i, j, k] + ζⁿ * G⁻[i, j, k])
    end
end

"""
Time step velocity fields with a 3rd-order Runge-Kutta method.
"""
@kernel function rk3_substep_field!(U, Δt, γ¹, ::Nothing, G¹, G⁰)
    i, j, k = @index(Global, NTuple)

    @inbounds begin
        U[i, j, k] += Δt * γ¹ * G¹[i, j, k]
    end
end


text:
"""
    AllDerivatives{N}

A singleton type that represents all `m`-th derivatives where `0 ≤ m < N`.
"""
struct AllDerivatives{N}
    function AllDerivatives{N}() where N
        N isa Int || throw(ArgumentError("Type parameter N in AllDerivatives{N} must be an Int"))
        N ≥ 1 || throw(DomainError(N, "Type parameter N in AllDerivatives{N} must be positive"))
        new()
    end
end

"""
    AllDerivatives(N)

A shortcut for `AllDerivatives{N}()`, representing all `m`-th derivatives where `0 ≤ m < N`.
"""
AllDerivatives(N::Integer) = AllDerivatives{Int(N)}()

"""
    Derivative{N}

A singleton type that represents the `N`-th derivative.
"""
struct Derivative{N}
    function Derivative{N}() where N
        N isa Int || throw(ArgumentError("Type parameter N in Derivative{N} must be an Int"))
        N ≥ 0 || throw(DomainError(N, "Type parameter N in Derivative{N} must be positive"))
        new()
    end
end

"""
    Derivative(N)

A shortcut for `Derivative{N}()`, representing the `N`-th derivative.
"""
Derivative(N::Integer) = Derivative{Int(N)}()

"""
    NoDerivative

An alias for `Derivative{0}`. (not exported)
"""
const NoDerivative = Derivative{0}

"""
    NoDerivUnion

An alias for `Union{Derivative{0}, AllDerivatives{1}}`. (not exported)
"""
const NoDerivUnion = Union{Derivative{0}, AllDerivatives{1}}

# AllDerivatives/Derivative act as scalars for broadcasting
Base.Broadcast.broadcastable(x::AllDerivatives) = Ref(x)
Base.Broadcast.broadcastable(x::Derivative) = Ref(x)


text:
using Distributions, HCubature, Printf, LinearAlgebra, Random, ProgressMeter

function test_inside_logarithm()

    rg = MersenneTwister(1)

    N = 5 # number of observations

    y = randn(rg, N)
    f = randn(rg, N)
    A = Diagonal(rand(rg, N))
    b = rand(rg, N)
    S = Diagonal(rand(rg, N)) + 0.2*I
    K = 0.1*randn(rg, N, N) ; K= K*K' + 1.1*I

    # Exact
    a1 = logpdf(MvNormal(A*f + b, S), y)
    a2 = logpdf(MvNormal(zeros(N), K), f)

    # My implementation

    b1 = -0.5*N*log(2π) - 0.5*logdet(S)

    c1 = -0.5*(y .- A*f .- b)'*(S\( y.- A*f .- b))

    b2 = -0.5*N*log(2π)

    c2 = -0.5*logdet(K) - 0.5*f'*(K\f)

    # a1, b1+c1, a2, b2+c
    a1+a2, b1+c1+b2+c2
end


text:
using SafeTestsets
@safetestset "My placeholder tests" begin include("placeholder_tests.jl") end


text:
using MarketData
using Test

@testset "remote" begin
    @testset "FRED" begin
        ta = fred()
        @test ta |> timestamp |> length > 100
    end

    @testset "Yahoo" begin
        t = Dates.now() - Year(2)
        opt = YahooOpt(period1 = t)
        ta = yahoo(:AAPL, opt)
        @test ta |> timestamp |> length > 100
    end
end


text:
using Test
using MPI
using CLIMA.MPIStateArrays
using CLIMA.Mesh.BrickMesh
using Pkg

const ArrayType = CLIMA.array_type()
CLIMA.init()
const comm = MPI.COMM_WORLD

function main()

  crank = MPI.Comm_rank(comm)
  csize = MPI.Comm_size(comm)

  @assert csize == 3

  if crank == 0
    numreal = 4
    numghost = 3

    nabrtorank = [1, 2]

    sendelems = [1, 2, 3, 4, 1, 4]
    nabrtorecv = [1:2, 3:3]
    nabrtosend = [1:4, 5:6]

    vmaprecv = [37, 38, 39, 40, 42, 43, 44, 45, 46, 49, 52, 53, 54,
                57, 60, 61, 62, 63]
    vmapsend = [3, 6, 9, 10, 11, 12, 19, 22, 25, 34, 35, 36,
                1, 2, 3, 28, 31, 34]

    nabrtovmaprecv = [1:13, 14:18]
    nabrtovmapsend = [1:12, 13:18]

    expectedghostdata = [1001, 1002, 1003, 1004, 1006, 1007, 1008, 1009,
                         1010, 1013, 1016, 1017, 1018,
                         2003, 2006, 2007, 2008, 2009]
  elseif crank == 1
    numreal = 2
    numghost = 4

    nabrtorank = [0]

    sendelems = [1, 2]
    nabrtorecv = [1:4]
    nabrtosend = [1:2]

    vmaprecv = [21, 24, 27, 28, 29, 30, 37, 40, 43, 52, 53, 54]
    vmapsend = [1, 2, 3, 4, 6, 7, 8, 9, 10, 13, 16, 17, 18]

    nabrtovmaprecv = [1:length(vmaprecv)]
    nabrtovmapsend = [1:length(vmapsend)]

    expectedghostdata = [ 3,  6,  9,
                         10, 11, 12,
                         19, 22, 25,
                         34, 35, 36]
  elseif crank == 2
    numreal = 1
    numghost = 2

    nabrtorank = [0]

    sendelems = [1]
    nabrtorecv = [1:2]
    nabrtosend = [1:1]

    vmaprecv = [10, 11, 12, 19, 22, 25]
    vmapsend = [3, 6, 7, 8, 9]

    nabrtovmaprecv = [1:length(vmaprecv)]
    nabrtovmapsend = [1:length(vmapsend)]

    expectedghostdata = [ 1,  2,  3,
                         28, 31, 34]
  end

  numelem = numreal+numghost

  realelems = 1:numreal
  ghostelems = numreal .+ (1:numghost)

  weights = Array{Int64}(undef, (0, 0, 0))

  A = MPIStateArray{Int64}(comm, ArrayType, 9, 2, numelem, realelems, ghostelems,
                           ArrayType(vmaprecv), ArrayType(vmapsend), nabrtorank,
                           nabrtovmaprecv, nabrtovmapsend, ArrayType(weights),
                           555)

  Q = Array(A.data)
  Q .= -1
  shift = 100
  Q[:, 1, realelems] .= reshape((crank * 1000)          .+ (1:9*numreal), 9, numreal)
  Q[:, 2, realelems] .= reshape((crank * 1000) .+ shift .+ (1:9*numreal), 9, numreal)
  copyto!(A.data, Q)

  MPIStateArrays.start_ghost_exchange!(A)
  MPIStateArrays.finish_ghost_exchange!(A)

  Q = Array(A.data)
  @test all(         expectedghostdata .== Q[:, 1, :][:][vmaprecv])
  @test all(shift .+ expectedghostdata .== Q[:, 2, :][:][vmaprecv])

end

main()

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
# p37.jl - 2D "wave tank" with Neumann BCs for |y|=1

# x variable in [-A,A], Fourier:
A = 3; Nx = 50; dx = 2*A/Nx; x = -A+dx*(1:Nx);
D2x = (pi/A)^2*toeplitz([-1/(3*(dx/A)^2)-1/6;
        @. .5*(-1).^(2:Nx)/sin((pi*dx/A)*(1:Nx-1)/2)^2]);

# y variable in [-1,1], Chebyshev:
Ny = 15; (Dy,y) = cheb(Ny); D2y = Dy^2;
BC = -Dy[[1,Ny+1],[1,Ny+1]]\Dy[[1,Ny+1],2:Ny];

# Grid and initial data:
vv = @. exp(-8*((x'+1.5)^2+y^2));
dt = 5/(Nx+Ny^2);
vvold = @. exp(-8*((x'+dt+1.5).^2+y.^2));

# Time-stepping by leap frog formula:
plotgap = round(Int,2/dt); dt = 2/plotgap;
for n = 0:2*plotgap
    t = n*dt;
    if rem(n+.5,plotgap)<1
        figure(n/plotgap+1); clf(); surf(x,y,vv);
        gca()[:view_init](60,-100);
        xlim(-A,A); ylim(-1,1); zlim(-0.15,1);
        text3D(-2.5,1,.5,"t = $(round(t))",fontsize=18);
        zticks([]);
    end
    vvnew = 2*vv - vvold + dt^2*(vv*D2x +D2y*vv);
    vvold = vv; vv = vvnew;
    vv[[1,Ny+1],:] = BC*vv[2:Ny,:];       # Neumann BCs for |y|=1
end


text:
using fastPHASE
using Test

@test 1+1 == 2


text:
function inverse(Side::HorizonSide)
    return HorizonSide((mod(Int(Side)+2, 4)))
end

global flag = true

function putmarkers!(r::Robot, side::HorizonSide)
    global flag
    while !isborder(r,side)
        move!(r,side)
        flag = !flag
        if flag
            putmarker!(r)
        end
    end
end

function mark_chess!(r::Robot)
    global flag
    mvs = []
    go_to_corner!(r, mvs)
    print(mvs)
    S = Ost
    mark_line(r,S)
    while !isborder(r, Nord)
        move!(r, Nord)
        flag = !flag
        S = inverse(S)
        mark_line(r, S)
    end
    go_to_corner_Sud_West!(r)
    return_back!(r, mvs)
end

function mark_line(r, S)
    global flag
    if flag
        putmarker!(r)
    end
    putmarkers!(r,S)
end

function go_to_corner!(r::Robot, mvs::AbstractArray)
    while isborder(r, Sud) == false || isborder(r, West) == false
        moveIfPossible!(r, mvs, Sud)
        moveIfPossible!(r, mvs, West)
    end
end

function moveIfPossible!(r::Robot, mvs::AbstractArray, side::HorizonSide)
    if !isborder(r, side)
        move!(r, side)
        push!(mvs, side)
    end
end

function return_back!(r::Robot, mvs::AbstractArray)
    while length(mvs) > 0
        last_move = pop!(mvs)
        move!(r, inverse(last_move))
    end
end

function go_to_corner_Sud_West!(r::Robot)
    while !isborder(r, Sud) || !isborder(r, West)
        while !isborder(r, Sud)
            move!(r, Sud)
        end
        while !isborder(r, West)
            move!(r, West)
        end
    end
end


text:
"""
"""
abstract type SingleFieldFESpace <: FESpace end

"""
"""
function get_cell_dof_basis(f::SingleFieldFESpace)
  @abstractmethod
end

"""
"""
function num_dirichlet_dofs(f::SingleFieldFESpace)
  @abstractmethod
end

"""
"""
function zero_dirichlet_values(f::SingleFieldFESpace)
  @abstractmethod
end

"""
"""
function num_dirichlet_tags(f::SingleFieldFESpace)
  @abstractmethod
end

"""
"""
function get_dirichlet_dof_tag(f::SingleFieldFESpace)
  @abstractmethod
end

"""
"""
function scatter_free_and_dirichlet_values(f::SingleFieldFESpace,free_values,dirichlet_values)
  @abstractmethod
end

"""
"""
function gather_free_and_dirichlet_values!(free_values, dirichlet_values,fs::SingleFieldFESpace,cell_vals)
  @abstractmethod
end

"""
"""
function test_single_field_fe_space(f::SingleFieldFESpace,pred=(==))
  fe_basis = get_cell_basis(f)
  @test isa(fe_basis,CellField)
  @test is_basis(fe_basis)
  test_fe_space(f)
  cell_dofs = get_cell_dofs(f)
  dirichlet_values = zero_dirichlet_values(f)
  @test length(dirichlet_values) == num_dirichlet_dofs(f)
  free_values = zero_free_values(f)
  cell_vals = scatter_free_and_dirichlet_values(f,free_values,dirichlet_values)
  fv, dv = gather_free_and_dirichlet_values(f,cell_vals)
  @test pred(fv,free_values)
  @test pred(dv,dirichlet_values)
  gather_free_and_dirichlet_values!(fv,dv,f,cell_vals)
  @test pred(fv,free_values)
  @test pred(dv,dirichlet_values)
  fv, dv = gather_free_and_dirichlet_values!(fv,dv,f,cell_vals)
  @test pred(fv,free_values)
  @test pred(dv,dirichlet_values)
  fe_function = FEFunction(f,free_values,dirichlet_values)
  @test isa(fe_function, SingleFieldFEFunction)
  test_fe_function(fe_function)
  ddof_to_tag = get_dirichlet_dof_tag(f)
  @test length(ddof_to_tag) == num_dirichlet_dofs(f)
  if length(get_dirichlet_dof_tag(f)) != 0
    @test maximum(get_dirichlet_dof_tag(f)) <= num_dirichlet_tags(f)
  end
  cell_dof_basis = get_cell_dof_basis(f)
  @test isa(CellField(f,get_cell_dofs(f)),CellField)
end

function test_single_field_fe_space(f,matvecdata,matdata,vecdata,pred=(==))
  test_single_field_fe_space(f,pred)
  test_fe_space(f,matvecdata,matdata,vecdata)
end

"""
"""
function CellData.get_cell_map(fs::SingleFieldFESpace)
  fe_basis = get_cell_basis(fs)
  get_cell_map(fe_basis)
end

function CellData.get_cell_axes(f::SingleFieldFESpace)
  get_cell_axes(get_cell_basis(f))
end

"""
    FEFunction(
      fs::SingleFieldFESpace, free_values::AbstractVector, dirichlet_values::AbstractVector)

The resulting FEFunction will be in the space if and only if `dirichlet_values`
are the ones provided by `get_dirichlet_values(fs)`
"""
function FEFunction(
  fs::SingleFieldFESpace, free_values::AbstractVector, dirichlet_values::AbstractVector)
  cell_vals = scatter_free_and_dirichlet_values(fs,free_values,dirichlet_values)
  cell_field = CellField(fs,cell_vals)
  SingleFieldFEFunction(cell_field,cell_vals,free_values,dirichlet_values,fs)
end

function CellData.CellField(fs::SingleFieldFESpace,cell_vals)
  _default_cell_field(fs,cell_vals)
end

function _default_cell_field(fs,cell_vals)
  cell_basis = get_cell_basis(fs)
  cell_field = lincomb(cell_basis,cell_vals)
  cell_field
end

function FEFunction(fe::SingleFieldFESpace, free_values)
  diri_values = get_dirichlet_values(fe)
  FEFunction(fe,free_values,diri_values)
end

"""
"""
function get_dirichlet_values(f::SingleFieldFESpace)
  zero_dirichlet_values(f)
end

"""
"""
function gather_free_and_dirichlet_values(f::SingleFieldFESpace,cell_vals)
  free_values = zero_free_values(f)
  dirichlet_values = zero_dirichlet_values(f)
  gather_free_and_dirichlet_values!(free_values,dirichlet_values,f,cell_vals)
end

"""
"""
function gather_dirichlet_values(f::SingleFieldFESpace,cell_vals)
  dirichlet_values = zero_dirichlet_values(f)
  gather_dirichlet_values!(dirichlet_values,f,cell_vals)
  dirichlet_values
end

"""
"""
function gather_free_values(f::SingleFieldFESpace,cell_vals)
  free_values = zero_free_values(f)
  gather_free_values!(free_values,f,cell_vals)
  free_values
end

"""
"""
function gather_dirichlet_values!(dirichlet_values,f::SingleFieldFESpace,cell_vals)
  free_values = zero_free_values(f)
  gather_free_and_dirichlet_values!(free_values,dirichlet_values,f,cell_vals)
  dirichlet_values
end

"""
"""
function gather_free_values!(free_values,f::SingleFieldFESpace,cell_vals)
    dirichlet_values = zero_dirichlet_values(f)
    gather_free_and_dirichlet_values!(free_values,dirichlet_values,f,cell_vals)
    free_values
end

@deprecate(
  interpolate(fs::SingleFieldFESpace, object),
  interpolate(object, fs::SingleFieldFESpace)
)

@deprecate(
  interpolate!(free_values,fs::SingleFieldFESpace, object),
  interpolate!(object, free_values,fs::SingleFieldFESpace)
)

@deprecate(
  interpolate_everywhere(fs::SingleFieldFESpace, object),
  interpolate_everywhere(object, fs::SingleFieldFESpace)
)

@deprecate(
  interpolate_everywhere!(free_values,dirichlet_values,fs::SingleFieldFESpace, object),
  interpolate_everywhere!(object, free_values,dirichlet_values,fs::SingleFieldFESpace)
)

@deprecate(
  interpolate_dirichlet(fs::SingleFieldFESpace, object),
  interpolate_dirichlet(object, fs::SingleFieldFESpace)
)

@deprecate(
  interpolate_dirichlet!(free_values,dirichlet_values,fs::SingleFieldFESpace, object),
  interpolate_dirichlet!(object, free_values,dirichlet_values,fs::SingleFieldFESpace)
)

"""
The resulting FE function is in the space (in particular it fulfills Dirichlet BCs
even in the case that the given cell field does not fulfill them)
"""
function interpolate(object, fs::SingleFieldFESpace)
  free_values = zero_free_values(fs)
  interpolate!(object, free_values,fs)
end

"""
"""
function interpolate!(object, free_values,fs::SingleFieldFESpace)
    cell_vals = _cell_vals(fs,object)
    gather_free_values!(free_values,fs,cell_vals)
    FEFunction(fs,free_values)
end

function _cell_vals(fs::SingleFieldFESpace,object)
  cdb = get_cell_dof_basis(fs)
  cm = get_cell_map(fs)
  cf = convert_to_cell_field(object,cm,RefStyle(cdb))
  cell_vals = evaluate(cdb,cf)
end

"""
like interpolate, but also compute new degrees of freedom for the dirichlet component.
The resulting FEFunction does not necessary belongs to the underlying space
"""
function interpolate_everywhere(object, fs::SingleFieldFESpace)
  free_values = zero_free_values(fs)
  dirichlet_values = zero_dirichlet_values(fs)
  interpolate_everywhere!(object, free_values,dirichlet_values,fs)
end

"""
"""
function interpolate_everywhere!(object, free_values,dirichlet_values,fs::SingleFieldFESpace)
  cell_vals = _cell_vals(fs,object)
  gather_free_and_dirichlet_values!(free_values,dirichlet_values,fs,cell_vals)
  FEFunction(fs,free_values,dirichlet_values)
end

"""
"""
function interpolate_dirichlet(object, fs::SingleFieldFESpace)
  free_values = zero_free_values(fs)
  dirichlet_values = zero_dirichlet_values(fs)
  interpolate_dirichlet!(object, free_values,dirichlet_values,fs)
end

"""
"""
function interpolate_dirichlet!(object, free_values,dirichlet_values,fs::SingleFieldFESpace)
  cell_vals = _cell_vals(fs,object)
  gather_dirichlet_values!(dirichlet_values,fs,cell_vals)
  fill!(free_values,zero(eltype(free_values)))
  FEFunction(fs,free_values, dirichlet_values)
end

"""
"""
function compute_dirichlet_values_for_tags(f::SingleFieldFESpace,tag_to_object)
  dirichlet_values = zero_dirichlet_values(f)
  compute_dirichlet_values_for_tags!(dirichlet_values,f,tag_to_object)
end

function compute_dirichlet_values_for_tags!(dirichlet_values,f::SingleFieldFESpace,tag_to_object)
  dirichlet_dof_to_tag = get_dirichlet_dof_tag(f)
  _tag_to_object = _convert_to_collectable(tag_to_object,num_dirichlet_tags(f))
  for (tag, object) in enumerate(_tag_to_object)
    cell_vals = _cell_vals(f,object)
    dv = gather_dirichlet_values(f,cell_vals)
    _fill_dirichlet_values_for_tag!(dirichlet_values,dv,tag,dirichlet_dof_to_tag)
  end
  dirichlet_values
end

function _fill_dirichlet_values_for_tag!(dirichlet_values,dv,tag,dirichlet_dof_to_tag)
  for (dof, v) in enumerate(dv)
    if dirichlet_dof_to_tag[dof] == tag
      dirichlet_values[dof] = v
    end
  end
end

function _convert_to_collectable(object,ntags)
  @assert ntags == length(object) "Incorrect number of dirichlet tags provided"
  object
end

function _convert_to_collectable(object::Function,ntags)
  _convert_to_collectable(fill(object,ntags),ntags)
end

function _convert_to_collectable(object::Number,ntags)
  _convert_to_collectable(fill(object,ntags),ntags)
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 9a052fc4-6c5e-11eb-24e7-3f1da02508bc
begin
	using Test
	using Pkg
	using LinearAlgebra

	Pkg.build("PyCall")
	Pkg.add("SymPy")
	Pkg.add("Images")
	Pkg.add("ImageIO")

	using SymPy
	using Images
	using ImageIO

	include("../../lib/KTBC.jl")
	import .KTBC.CreateKTBC
end

# ╔═╡ 437d7d1e-6c5e-11eb-0e0b-91d58c495478
md"""
## \#1
"""

# ╔═╡ 575d4bac-6c5e-11eb-060a-1dabfe5f167f
begin
	Aₜ = [
		-1  1 0
		 0 -1 1
		-1  0 1
	]

	Aₜᵀ = Aₜ'

	@test Aₜᵀ * Aₜ == [
		2 -1 -1
	   -1  2 -1
	   -1 -1  2
	]

	Aₛ = [
	   -1  1  0 0
	   -1  0  1 0
	   -1  0  0 1
		0 -1  0 1
		0  0 -1 1
	]

	Aₛᵀ = Aₛ'

	@test Aₛᵀ * Aₛ == [
		 3   -1  -1  -1
		-1    2   0  -1
		-1    0   2  -1
		-1   -1  -1   3
	]
end

# ╔═╡ 0119c440-6c5f-11eb-2c3d-f1319bccd42d
md"""
## \#2
"""

# ╔═╡ 03a129ce-6c5f-11eb-1883-ab1716ee30a3
begin
	C = sympy.Symbol("C")

	uₜ = ones(3, 1)

	@test convert(Array{Int}, Aₜ * (C * uₜ))[:] == [
		0
		0
		0
	]

	vₜ = [
		 1
		 1
		-1
	]

	@test convert(Array{Int}, Aₜᵀ * (C * vₜ)) == [
		0
		0
		0
	]
end

# ╔═╡ 77c8d8b0-6c5f-11eb-2623-61783b1d556b
md"""
## \#3
"""

# ╔═╡ 6a9bc3b2-6c61-11eb-11e8-9338d319794b
begin
	uₛ = ones(4, 1)

	@test convert(Array{Int}, Aₛ * (C * uₛ))[:] == [
		0
		0
		0
		0
		0
	]

	vₛ₁ = [
		1
		1
	   -2
		1
		1
	]

	@test convert(Array{Int}, Aₛᵀ * (C * vₛ₁)) == [
		0
		0
		0
		0
	]

	vₛ₂ = [
	   -2
		1
		1
	   -2
		1
	]

	@test convert(Array{Int}, Aₛᵀ * (C * vₛ₂)) == [
		0
		0
		0
		0
	]
end

# ╔═╡ 90d10c28-6c62-11eb-1f2e-b346125b3932
md"""
## \#8
"""

# ╔═╡ 9da4bf3a-6c6a-11eb-349a-153eec621c58
begin
	c₁, c₂, c₃, c₄, c₅, c₆ = symbols("c₁ c₂ c₃ c₄ c₅ c₆")

	K = [
		c₁+c₂+c₄       -c₁             -c₂           -c₄
	      -c₁        c₁+c₃+c₅          -c₃           -c₅
	      -c₂          -c₃           c₂+c₃+c₆        -c₆
	      -c₄          -c₅             -c₆         c₄+c₅+c₆
	]

	D = [
		c₁+c₂+c₄      0         0        0
	       0       c₁+c₃+c₅     0        0
	       0          0      c₂+c₃+c₆    0
	       0          0         0     c₄+c₅+c₆
	]

	W = [
		  0      c₁     c₂      c₄
	      c₁     0      c₃      c₅
	      c₂     c₃     0       c₆
	      c₄     c₅     c₆      0
	]

	@test K == D - W
end

# ╔═╡ 0a72fc88-6c63-11eb-0d6f-c790283ecd24
begin
	a₁ = [
		1
	   -1
		0
		0
	]

	K₁ = [
		 c₁ -c₁ 0 0
		-c₁  c₁ 0 0
		 0   0  0 0
		 0   0  0 0
	]

	@test c₁ * a₁ * a₁' == K₁

	a₂ = [
		1
		0
	   -1
		0
	]

	K₂ = [
		c₂ 0 -c₂ 0
		0  0  0  0
	   -c₂ 0  c₂ 0
		0  0  0  0
	]

	@test c₂ * a₂ * a₂' == K₂

	K₃ = [
		0   0   0  0
		0   c₃ -c₃ 0
		0  -c₃  c₃ 0
		0   0   0  0
	]

	a₃ = [
		0
		1
	   -1
	    0
	]

	c₃ * a₃ * a₃'

	@test c₃ * a₃ * a₃' == K₃

	K₄ = [
	    c₄  0  0  -c₄
		0   0  0   0
		0   0  0   0
	   -c₄  0  0   c₄
	]

	a₄ = [
		1
		0
	    0
	   -1
	]

	@test K₄ == c₄ * a₄ * a₄'

	K₅ = [
	    0  0  0  0
		0  c₅ 0 -c₅
		0  0  0  0
		0 -c₅ 0  c₅
	]

    a₅ = [
		0
		1
	    0
	   -1
	]

	@test K₅ == c₅ * a₅ * a₅'

	K₆ = [
         0  0  0   0
		 0  0  0   0
		 0  0  c₆ -c₆
		 0  0 -c₆  c₆
	]

	a₆ = [
		0
		0
	    1
	   -1
	]

	@test K₆ == c₆ * a₆ * a₆'

	@test K == K₁ + K₂ + K₃ + K₄ + K₅ + K₆
end

# ╔═╡ 23e4f7aa-6c6a-11eb-3fe1-8daa7a059e68
begin
	A₆ = [a₁ a₂ a₃ a₄ a₅ a₆]'

	@test A₆ == [
		 1  -1   0   0
 		 1   0  -1   0
 		 0   1  -1   0
 		 1   0   0  -1
 		 0   1   0  -1
 		 0   0   1  -1
	]

	C₆ = diagm(0 => [c₁, c₂, c₃, c₄, c₅, c₆])

	@test K == A₆' * C₆ * A₆
end

# ╔═╡ 03d87872-6c6c-11eb-1043-59a80e58642f
md"""
## \#9
"""

# ╔═╡ 13269750-6c6c-11eb-3262-2937ef0d3f5e
begin
	A₅ = [
		-1  1  0  0 0
		 0 -1  1  0 0
		 0  0 -1  1 0
		 0  0  0 -1 1
	]

	# 1) AᵀA == B₅
	B₅ = CreateKTBC(5)[3]

	A₅ᵀA₅ = A₅' * A₅

	@test B₅ == A₅ᵀA₅

	# 2) T₄ == (T₄ᵀT₄)_reduced

	A₄ = A₅[:,1:4]

	A₄ᵀA₄ = A₄' * A₄

	T₄ = convert(Array{Int}, CreateKTBC(4)[2])

	@test T₄ == A₄ᵀA₄

	# 3) inv(T₄)
	@test inv(T₄) == [
		4.0  3.0  2.0  1.0
 		3.0  3.0  2.0  1.0
 		2.0  2.0  2.0  1.0
 		1.0  1.0  1.0  1.0
	]

	# 4) eigenvalues of T₄
	values, vectors = eigen(T₄)

	round.(values, digits=2) == [
		0.12
		1.0
		2.35
		3.53
	]

	# 5) det(T₄)
	@test det(T₄) == 1.0
end

# ╔═╡ 59da7064-6c70-11eb-096b-4f532330b19c
md"""
## \#12
"""

# ╔═╡ 61271e4e-6c70-11eb-1dba-91041e4280fb
begin
	# Example
	dim = 5

	onesMatrix = ones(dim - 1, dim - 1)
	eyeMatrix = I(dim - 1)

	Kᵣ = diagm(0 => [dim for i=1:dim - 1]) - onesMatrix

	Kᵣinv = (1/dim) * (onesMatrix + eyeMatrix)

	round.(Int, Kᵣ * Kᵣinv) == eyeMatrix

	@test round(Int, det(Kᵣ)) == 125 # Positive definite
end

# ╔═╡ caa0ca62-6c77-11eb-1fc9-a79b32e5eeb3
md"""
```math
KK^{-1} =
\begin{bmatrix}
   n-1 & -1 & \dots & -1 \\
   -1 & n-1 & \dots & -1 \\
   \vdots & \vdots & \ddots & \vdots \\
   -1 & -1 & \dots & n-1 \\
\end{bmatrix} \frac{1}{n}
\begin{bmatrix}
   2 & 1 & \dots & 1 \\
   1 & 2 & \dots & 1 \\
   \vdots & \vdots & \ddots & \vdots \\
   1 & 1 & \dots & 2 \\
\end{bmatrix} =
```
##### K and inv K is (n - 1) * (n - 1)

```math
\frac{1}{n}
\begin{bmatrix}
   n & 0 & \dots & 0 \\
   0 & n & \dots & 0 \\
   \vdots & \vdots & \ddots & \vdots \\
   0 & 0 & \dots & n \\
\end{bmatrix} = I
```

##### Positive definite by det test

###### For n > 1
```math
det(n - 1) = n - 1 > 0
```

###### For 2 * 2, n > 2

```math
det\begin{pmatrix}
   n-1 & -1 \\
   -1 & n-1 \\
\end{pmatrix} = (n - 1)^2 - 1 = n^2 - 2n > 0
```

###### And etc

"""

# ╔═╡ d426ccea-6c7d-11eb-1d42-adde535f324c
md"""
## \#17
"""

# ╔═╡ 7e8ad1e6-6c7b-11eb-053b-01566d56c7f0
load("./17.png")

# ╔═╡ e8e78e24-6c7d-11eb-2a5f-7d9a9665ca06
begin
	A₉ = [
		1 1 0 0 0 0 0 0 0
		0 1 1 0 0 0 0 0 0
		0 0 1 0 0 1 0 0 0
		0 0 0 0 0 1 0 0 1
		0 0 0 0 0 0 0 1 1
		0 0 0 0 0 0 1 1 0
		0 0 0 1 0 0 1 0 0
		1 0 0 1 0 0 0 0 0
		0 0 0 1 1 0 0 0 0
		0 0 0 0 1 1 0 0 0
		0 1 0 0 1 0 0 0 0
		0 0 0 0 1 0 0 1 0
	]

	A₉ᵀA₉ = A₉' * A₉

	n, m = size(A₉ᵀA₉)

	# Non-zero values
	@test n * m - size(filter(iszero, A₉ᵀA₉))[1] == 33

	D₉ = diagm(0 => diag(A₉ᵀA₉))

	W₉ = A₉ᵀA₉ - D₉

	@test W₉ == [
		 0  1  0  1  0  0  0  0  0
		 1  0  1  0  1  0  0  0  0
		 0  1  0  0  0  1  0  0  0
		 1  0  0  0  1  0  1  0  0
		 0  1  0  1  0  1  0  1  0
		 0  0  1  0  1  0  0  0  1
		 0  0  0  1  0  0  0  1  0
		 0  0  0  0  1  0  1  0  1
		 0  0  0  0  0  1  0  1  0
	]

	@test D₉ == [
		 2  0  0  0  0  0  0  0  0
		 0  3  0  0  0  0  0  0  0
		 0  0  2  0  0  0  0  0  0
		 0  0  0  3  0  0  0  0  0
		 0  0  0  0  4  0  0  0  0
		 0  0  0  0  0  3  0  0  0
		 0  0  0  0  0  0  2  0  0
		 0  0  0  0  0  0  0  3  0
		 0  0  0  0  0  0  0  0  2
	]

	# node 5 has 4 edges connected to it
	@test D₉[5, 5] == 4

	# nodes 2, 4, 6, 8 connection in the row
	@test W₉[5:5, :] == [0  1  0  1  0  1  0  1  0]
end

# ╔═╡ Cell order:
# ╠═9a052fc4-6c5e-11eb-24e7-3f1da02508bc
# ╠═437d7d1e-6c5e-11eb-0e0b-91d58c495478
# ╠═575d4bac-6c5e-11eb-060a-1dabfe5f167f
# ╠═0119c440-6c5f-11eb-2c3d-f1319bccd42d
# ╠═03a129ce-6c5f-11eb-1883-ab1716ee30a3
# ╠═77c8d8b0-6c5f-11eb-2623-61783b1d556b
# ╠═6a9bc3b2-6c61-11eb-11e8-9338d319794b
# ╠═90d10c28-6c62-11eb-1f2e-b346125b3932
# ╠═9da4bf3a-6c6a-11eb-349a-153eec621c58
# ╠═0a72fc88-6c63-11eb-0d6f-c790283ecd24
# ╠═23e4f7aa-6c6a-11eb-3fe1-8daa7a059e68
# ╠═03d87872-6c6c-11eb-1043-59a80e58642f
# ╠═13269750-6c6c-11eb-3262-2937ef0d3f5e
# ╟─59da7064-6c70-11eb-096b-4f532330b19c
# ╠═61271e4e-6c70-11eb-1dba-91041e4280fb
# ╟─caa0ca62-6c77-11eb-1fc9-a79b32e5eeb3
# ╠═d426ccea-6c7d-11eb-1d42-adde535f324c
# ╠═7e8ad1e6-6c7b-11eb-053b-01566d56c7f0
# ╠═e8e78e24-6c7d-11eb-2a5f-7d9a9665ca06


text:
@testset "Memory" begin
    game = Kuhn()
    sol = DeepCFRSolver(game; buffer_size=50, traversals=10)
    train!(sol, 100)

    @test length(sol.Mπ) == length(sol.Mπ.I) == length(sol.Mπ.t) == length(sol.Mπ.σ)
    @test length(sol.Mπ) == sol.Mπ.capacity == 50

    Mv1 = first(sol.Mv)
    Mv2 = first(sol.Mv)
    @test length(Mv1) == length(Mv1.I) == length(Mv1.t) == length(Mv1.r)
    @test length(Mv1) == Mv1.capacity == 50
    @test length(Mv2) == length(Mv2.I) == length(Mv2.t) == length(Mv2.r)
    @test length(Mv2) == Mv2.capacity == 50

    empty!(sol.Mπ)
    @test length(sol.Mπ) == 0
    empty!(Mv1)
    @test length(Mv1) == 0
end


text:
using SCIP
using QuadraticToBinary

SCIPopt = SCIP.Optimizer()
MOI.set(SCIPopt, MOI.Silent(), true)
MOI.set(SCIPopt, MOI.TimeLimitSec(), 5)
SCIPopt_CACHED = MOIU.CachingOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()), SCIPopt)

push!(solvers, (opt = SCIPopt, mode = BilevelJuMP.SOS1Mode()))
push!(solvers_sos, (opt = SCIPopt, mode = BilevelJuMP.SOS1Mode()))
push!(solvers_unit, (opt = SCIPopt, mode = BilevelJuMP.SOS1Mode()))
push!(solvers_indicator, (opt = SCIPopt, mode = BilevelJuMP.IndicatorMode()))

push!(solvers_fa, (opt = SCIPopt, mode = BilevelJuMP.FortunyAmatMcCarlMode(with_slack = false)))
push!(solvers_fa, (opt = SCIPopt, mode = BilevelJuMP.FortunyAmatMcCarlMode(with_slack = true)))
push!(solvers_fa2, (opt = SCIPopt, mode = BilevelJuMP.FortunyAmatMcCarlMode(primal_big_M = 100, dual_big_M = 100)))

push!(solvers_cached, (opt = SCIPopt_CACHED, mode = BilevelJuMP.SOS1Mode()))

push!(solvers_bin_exp, (
    opt = QuadraticToBinary.Optimizer{Float64}(SCIPopt_CACHED),
    mode = BilevelJuMP.StrongDualityMode(1e-9, inequality = true)))
MOI.set(solvers_bin_exp[end].opt, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)

push!(solvers_bin_exp, (
    opt = QuadraticToBinary.Optimizer{Float64}(SCIPopt_CACHED),
    mode = BilevelJuMP.StrongDualityMode(inequality = false)))
MOI.set(solvers_bin_exp[end].opt, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)

QTB_SCIPopt_BRIDGED = MOI.Bridges.full_bridge_optimizer(
    QuadraticToBinary.Optimizer{Float64}(SCIPopt_CACHED, lb=-100,ub=100),
    Float64)
MOI.Bridges.add_bridge(QTB_SCIPopt_BRIDGED, MOI.Bridges.Constraint.SOCtoNonConvexQuadBridge{Float64})
MOI.set(QTB_SCIPopt_BRIDGED, QuadraticToBinary.GlobalVariablePrecision(), 1e-5)

push!(solvers_sos_quad_bin, (
    opt = QTB_SCIPopt_BRIDGED,
    mode = BilevelJuMP.SOS1Mode()))

push!(solvers_fa_quad_bin, (
    opt = QTB_SCIPopt_BRIDGED,
    mode = BilevelJuMP.FortunyAmatMcCarlMode(primal_big_M = 100, dual_big_M = 100)))

QTB_SCIPopt = QuadraticToBinary.Optimizer{Float64}(SCIPopt_CACHED, lb=-5, ub=5)
QTB_SCIPopt_B = MOI.Bridges.Constraint.SOCtoNonConvexQuad{Float64}(QTB_SCIPopt)
push!(solvers_fa_quad_bin_mixed, (
    opt = QTB_SCIPopt_B,
    mode = BilevelJuMP.MixedMode(default=
        BilevelJuMP.FortunyAmatMcCarlMode(primal_big_M = 10, dual_big_M = 10))))
MOI.set(QTB_SCIPopt_B, QuadraticToBinary.GlobalVariablePrecision(), 1e-3)


text:
using Minesweeper

# play in REPL
# enter coordinates as: row_index column_index
play()

# pass custom game to play
game = Game(dims=(20,10), n_mines=35)
play(game)


text:
#=============================#
#== put needed modules here ==#
#=============================#
import PyCall
ptl = PyCall.pyimport("qcportal")

import JuliaChem

import JSON

#================================#
#== JuliaChem execution script ==#
#================================#
function script()
  #== get molecule information from QCArchive ==#
  client = ptl.FractalClient()
  s22_database = client.get_collection("ReactionDataset", "S22")

  mol_index = s22_database.get_index()[1]
  mol_object_full = s22_database.get_molecules(first) #doesn't seem to work
  mol_object = get(mol_object_full.molecule, 0)

  #== create input system ==#
  molecule = JSON.parse(mol_object.json())
  display(molecule)

  driver = "energy"

  model = Dict(
    "method" => "RHF",
    "basis" => "STO-3G"
  )

  keywords = Dict(
    "scf" => Dict(
      "niter" => 50,
      "ndiis" => 3,
      "dele" => 1E-8,
      "rmsd" => 1E-6,
      "prec" => "Float64",
      "direct" => false,
      "debug" => false
    )
  )

  #== generate basis set ==#
  basis = JuliaChem.JCBasis.run(molecule, model)

  shell = basis[2]
  display(shell)
end

JuliaChem.initialize()
script()
JuliaChem.finalize()


text:
module MappedCartesianGeometryTests

using Test
using Gridap
import Gridap: ∇

# using LinearAlgebra
# LinearAlgebra.cond(ka)
###

ufun(x) = x[1]
ufun_grad(x) = VectorValue(1.0,)
∇(::typeof(ufun)) = ufun_grad
bfun(x) = 0.0

n = 10
ℓ = 1.0
h = 1.0 / n
β = 5
geomap = ( x -> x[1] <= (ℓ-h) ? x[1] : (x[1]-ℓ+h)*β+ℓ-h)
# Construct the discrete model
model = CartesianDiscreteModel(partition=(n,),map=geomap)
# model = CartesianDiscreteModel(partition=(n,))

order = 1
# Construct the FEspace
diritag = "boundary"
fespace = H1ConformingFESpace(Float64,model,order,diritag)

# Define test and trial
V = TestFESpace(fespace)
U = TrialFESpace(fespace,ufun)

# Define integration mesh and quadrature
trian = Triangulation(model)
points(Grid(model))
quad = CellQuadrature(trian,degree=order*2)
# Define forms
a(v,u) = inner(∇(v), ∇(u))
b(v) = inner(v,bfun)

# Define Assembler
assem = SparseMatrixAssembler(V,U)

# Define the FEOperator
op = LinearFEOperator(a,b,V,U,assem,trian,quad)

# Define the FESolver
ls = BackslashSolver()
solver = LinearFESolver(ls)

# Solve!
uh = solve(solver,op)

# Define exact solution and error
u = CellField(trian,ufun)
e = u - uh

# Define norms to measure the error
l2(u) = inner(u,u)
h1(u) = a(u,u) + l2(u)

# Compute errors
el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))

@test el2 < 1.e-8

writevtk(trian,"onedres",cellfields=["uh"=>uh])

###
F = Gridap.CartesianGeometry

model1 = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0),partition=(1,1))
points(model1.cgrid)
model1.cgrid
# trian1 = Triangulation(model1)
# points(trian1.grid)
# phi1 = CellGeomap(trian1)

gmap(x) = 2*x
model2 = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0),partition=(1,1),map=gmap)
points(model2.cgrid)
trian = Triangulation(model2)
# points(trian2.grid)
phi = CellGeomap(trian)

##

# CellPoints(trian1)

# CellPoints(trian2)

##

quad = CellQuadrature(trian,degree=2)
z = coordinates(quad)
##
self = trian1
Gridap.CellValuesGallery.CellVectorFromLocalToGlobal(cells(self.grid),points(self.grid))

self = trian2
cells(trian1.grid)

##
# phi1.val.cellvalues
# typeof(phi1.val.cellvalues[2]) == typeof(phi2.val.cellvalues[2])
#
# phi1.val.cellvalues[2].lid_to_gid
# phi1.val.cellvalues[2].gid_to_val
# phi1.val.cellvalues[2].cv

  # gid_to_val::V
  # cv::CachedArray{T,1,Array{T,1}}
# phi2.val.cellvalues[2].lid_to_gid
# phi2.val.cellvalues[2].gid_to_val[1]
# phi2.val.cellvalues[2].cv
##

# getindex(phi2.val,1)
# vals = Gridap.CellMapApply._getvalues(1,phi2.val.cellvalues...)
# self = phi1
# vals = Gridap.CellMapApply._getvalues(1,self.cellvalues...)
#
#
#
evaluate(phi1[1],z[1])
evaluate(phi2[1],z[1])

# evaluateI
# pval = collect(phi.val)
# pval[1].inputs
#
# isa(phi,Gridap.CellFieldsOperations.IndexCellFieldLikeAndGradient)
# typeof(phi)
#
# ##
# phi
# uphys = fun ∘ phi
#
# celldofs = fesp.cell_eqclass
# nfdofs = fesp.dim_to_nface_eqclass
# E = dof_type(T)
# free_dofs = zeros(E, num_free_dofs(fesp))
# diri_dofs = zeros(E, num_diri_dofs(fesp))
# aux = zeros(E, length(dofb))
# _interpolate_values_kernel!(
#
# # F._compute_points(cg)
#
# cg[1]
#
# # Here we have the points
# model.ugrid.points

# model.gridgraph.data

end # module


text:
import Base: ==
import Base: hash

struct Node
    state::State
    g::Float64
    h::Float64
    came_from::Union{Node,Nothing}
    action::Union{Action,Nothing}
end

function ==(node1::Node, node2::Node)
    node1.state == node2.state
end

function hash(node::Node)
    hash(node.state)
end

function astar(planner::Planner, num_expansions::Int64, start_state::State)
    closed_set = Set{Node}()
    open_list = MutableBinaryMinHeap{Tuple{Float64,Float64,Int64}}()
    open_dict = Dict{Node,Int64}()
    open_dict_reverse = Dict{Int64,Node}()

    reached_goal = false
    count = 0
    best_node = nothing

    h = getHeuristic(planner, start_state)
    g = 0.0
    f = h
    updated_start_node = Node(start_state, g, h, nothing, nothing)
    handle = push!(open_list, (f, h, count))
    count += 1
    open_dict[updated_start_node] = handle
    open_dict_reverse[handle] = updated_start_node

    for exp = 1:num_expansions
        triplet, handle = top_with_handle(open_list)
        f, h, _ = triplet
        node = open_dict_reverse[handle]
        pop!(open_list)
        delete!(open_dict, node)
        delete!(open_dict_reverse, handle)
        push!(closed_set, node)

        if checkGoal(planner, node.state)
            reached_goal = true
            best_node = node
            break
        end

        for action in getActions(planner, node.state)
            neighbor_state, cost = getSuccessors(planner, node.state, action)
            neighbor_node = Node(neighbor_state, 0, 0, nothing, nothing)
            if neighbor_node in closed_set
                continue
            end

            tentative_g = node.g + cost
            if !haskey(open_dict, neighbor_node)
                h = getHeuristic(planner, neighbor_state)
                f = tentative_g + h
                updated_neighbor_node = Node(neighbor_state, tentative_g, h, node, action)
                handle = push!(open_list, (f, h, count))
                count += 1
                open_dict[updated_neighbor_node] = handle
                open_dict_reverse[handle] = updated_neighbor_node
            else
                handle = open_dict[neighbor_node]
                old_neighbor_node = open_dict_reverse[handle]
                if tentative_g < old_neighbor_node.g
                    h = old_neighbor_node.h
                    f = tentative_g + h
                    updated_neighbor_node =
                        Node(old_neighbor_node.state, tentative_g, h, node, action)
                    delete!(open_dict, neighbor_node)
                    update!(open_list, handle, (f, h, count))
                    count += 1
                    open_dict[updated_neighbor_node] = handle
                    open_dict_reverse[handle] = updated_neighbor_node
                end
            end
        end
        if length(open_list) == 0
            best_node = node
            break
        end
    end

    if !reached_goal && length(open_list) > 0
        _, handle = top_with_handle(open_list)
        best_node = open_dict_reverse[handle]
    end

    info = Dict()
    info["best_node_f"] = best_node.g + best_node.h
    info["reached_goal"] = reached_goal
    info["closed"] = closed_set

    best_action, path = getBestAction(updated_start_node, best_node)
    info["path"] = path

    return best_action, info
end

function getBestAction(start_node::Node, best_node::Node)
    if start_node == best_node
        println("Should not reach this!")
        return nothing, [start_node.state, start_node.state]
    end

    node = best_node.came_from
    action = best_node.action
    path = [best_node.state]
    while true
        if !isnothing(node.came_from)
            # not the start node
            push!(path, node.state)
            next_node = node.came_from
            action = node.action
            if isnothing(next_node.came_from)
                # next node is start node
                break
            end
            node = next_node
        else
            # start node
            break
        end
    end
    push!(path, start_node.state)
    reverse!(path)
    return action, path
end


text:
sqrt 2 = 1.4142135623730951
e = 2.7182818284590455
 pi = 3.141592653340542


text:
module Farey

export farey

farey(x) = sort(unique(n//d for d=1:x for n=0:d))

end # module


text:
module GenieAutoReload

@eval Main __revise_mode__ = :eval
@eval Main using Revise

using Revise
using Genie, Genie.Router, Genie.WebChannels
using Distributed, Logging

export autoreload

Genie.config.websockets_server = true

const WEBCHANNEL_NAME = "autoreload"
const GENIE_AUTORELOAD = true
const WATCHED_EXTENSIONS = String["jl", "html", "md", "js", "css"]
const SCRIPT_URI = "$(Genie.config.base_path)js/plugins/autoreload.js"

function collect_watched_files(files::Vector{String} = String[], extensions::Vector{String} = WATCHED_EXTENSIONS) :: Vector{String}
  result = String[]

  for f in files
    try
      push!(result, Genie.Util.walk_dir(f, only_extensions = extensions)...)
    catch ex
      @error ex
    end
  end

  result
end

function watch(files::Vector{String} = String[], extensions::Vector{String} = WATCHED_EXTENSIONS;
                delay::Int = 0)
  @info "Watching $files"

  Revise.revise()

  @async entr(collect_watched_files(files, extensions); all = true, postpone = true) do
    @info "Reloading!"

    if delay > 0
      @info "Waiting $delay seconds"
      sleep(delay)
    end

    try
      Genie.WebChannels.broadcast(WEBCHANNEL_NAME, "autoreload:full")
      WebChannels.unsubscribe_disconnected_clients(WEBCHANNEL_NAME)
    catch ex
      @warn ex
    end
  end
end

function assets_js() :: String
  """
  function autoreload_subscribe() {
    Genie.WebChannels.sendMessageTo("autoreload", "subscribe");
    console.log("Autoreloading ready");
  }

  setTimeout(autoreload_subscribe, 2000);

  Genie.WebChannels.messageHandlers.push(function(event) {
    if ( event.data == "autoreload:full" ) {
      location.reload(true);
    }
  });
  """
end

function assets_script() :: String
  """
  <script>
  $(assets_js())
  </script>
  """
end

function assets(; devonly = true) :: String
  if (devonly && Genie.Configuration.isdev()) || !devonly
    """<script src="$SCRIPT_URI"></script>"""
  else
    ""
  end
end

function autoreload(files::Vector{String} = String[], extensions::Vector{String} = WATCHED_EXTENSIONS;
                    devonly::Bool = true, delay::Int = 0)
  if devonly && !Genie.Configuration.isdev()
    @warn "AutoReload configured for dev environment only. Skipping."
    return nothing
  end

  route(SCRIPT_URI) do
    assets_js() |> Genie.Renderer.Js.js
  end

  channel("/$(WEBCHANNEL_NAME)/subscribe") do
    WebChannels.subscribe(params(:WS_CLIENT), WEBCHANNEL_NAME)
  end

  GenieAutoReload.watch(files, extensions, delay = delay)
end

function autoreload(files...; extensions::Vector{String} = WATCHED_EXTENSIONS, devonly = true, delay::Int = 0)
  autoreload([files...], [extensions...]; devonly = devonly, delay = delay)
end

end # module


text:
# A pseudo-parser which extracts information from Julia code

using LNR
import JuliaParser.Lexer
include("streams.jl")

# Legacy definitions, should be removed eventually
const identifier_inner = r"[^,⊶⩃↣⩄≗≦⥥∓≨⊘×≥≫⥘⪚⋼⭂⥩⩾⥓∝⤞⇼⥍⭈∽⩁⥎⟑⟉∊∷≳≠⩧⫏⇿⬵≬⩗⥭⦾⥖→∪⫑⪀⩠⥢⤌⋝⊕⪪≈⪏≤⨤⪿⟰≼⫂≹⪣⋴≧∸≐⭋∨⨳⭁∺⋥⟽⊷⟱≡\]⤅⪃⩋⩊⋣⋎⥗⨮⬻⪻≢∙⪕⩓⫺∧⧻⨭⊵≓⥬⥛⋿⭃⫒⫕⩡⬺⧷⥄⊱⨰⊇≊⨬≖>⤕⬴⟿⋘⪇≯⋕⤏⟶⥚⥜⨼∥⪠⥝⬷∘⊴⪈⤔⪍⫄?⊰⪌⋩≟⋜⫀\)⫎⩦⋏⫷⊋⪱⤀⩯⤘⫌⩱≜↓⋗↑≛⋌⪢⫖⋖⩰⊏⊗⪡⋆⟈⤂⥆⧁⊻⤋⤖⩹↦⪳⩸⥅∔⨺⋐≶⟵\}⪙⪧⇺%≭≕⥔⥐⊆⋸⅋⋒≃≝≿⇴⩌⋠⇽≰/⫙⊠⪼⇔\[⟾+≩⊟⨶⥰⪉≎≷⩣⭄&⨲⧣⩭≑⊐⫗⩬⩢⬽⪯⪓⪒≪∈⪘⬿⫸⇹⊅⨥⨩≚⋹⊃⊂⪞⋺⨹⋦∦≮⋧⋛⋾⊁≉￪≔±\{⩒⩑⋫￩⥤⨽⬲⪄⫓⪑∩⧡⩮⪟⪛⋽⪦⇒≁⪝⬳⩝⩳≴⪰⟻≣⦼⩷⇶⋳⪺⪜⩕⥦∛≽⋑⤓⟼⩏≲⊲≸⟺⇷⟹∌⩪⊞⫉⨴⪤⪸⥡⩔⭊⪆⩲⫈⥒⫋⬶⫁⪵∗⫊⩖≙⩐≍⨫⦸⋚⊄⫐⥇⥣⪲↔⪷⨈⧺⭌⨨≄⤟^≵⋭⋊⟷⩅∤⫆⊽\(⬸⤒⪾⩞⥫⥙⋙⨱⬹<⊎⤊⤁⇏≺⋵⥏⩴⋶⪂⥕⪨⋇⊊⫅⊖⪶⋬≻⋍⋓⩍≱⇻⩵↮⋋⪖⨢↠⤎⊈⊮⋪⊓⪔\⨧⩜⥞⫇⪫⬾⋷⤃⧥⫃⨷⥈⤄⩼⋤⥠⬼⤠⩛≂↚⥧|∍⨻⊙⨪∋⪋⋲⤍.\"⊑⩟⇎*:￬⭉⤉⥯⬱⇾⋡÷⥟⥋∉⬰≞≾⫍⨵⩚⩫≅⩿⪎⪴⊒⪽≀⫹⤇⋅⩀⊡⤆∜⤈⨣↛⊩⫔⦷⩺≋\-≇⋨⊜~⫛≌⥉√⋢⊛⤗⋟⧶≏⊔⪗⋞ ⩎⊳∾⥨￫⩘⥌⪹⪩⩻=⨸⪊⨇⧤⇸⊉⥑⥮⭀⧀⊚⊬≒\$⊀⋻⦿⭇⥊≆←⤐≘⋉⊼⥪⧴⪅⩽⪬⪁⋄⤑⨦⩶⇵⪥⊍⫘⩂⪐⟒⪭⪮⤝∻\"\n]"
const identifier = Regex("(?![!0-9])$(identifier_inner.pattern)+")
const identifier_start = Regex("^$(identifier.pattern)")

@defonce type Token{T} end

token(T) = Token{T}()

Base.symbol{T}(t::Token{T}) = T

Lexer.peekchar(r::LineNumberingReader) =
  eof(r)? Lexer.EOF : LNR.peekchar(r)

function lexcomment(ts)
  Lexer.readchar(ts)
  if !eof(ts.io) && Lexer.peekchar(ts) == '='
    try Lexer.skip_multiline_comment(ts, 1) end
    return token(:multicomment)
  else
    Lexer.skip_to_eol(ts)
    return token(:comment)
  end
end

function skipstring(stream::IO, multi=false)
  while !eof(stream)
    if startswith(stream, "\\\"")
    elseif startswith(stream, multi ? "\"\"\"" : "\"")
      break
    else
      read(stream, Char)
    end
  end
end

# TODO: handle string splicing
function lexstring(stream::IO)
  multi = startswith(stream, "\"\"\"")
  multi || startswith(stream, '"')
  skipstring(stream, multi)
  return token(multi ? :multistring : :string)
end

isidentifier(x::Symbol) = !(x in Lexer.syntactic_ops) && !(x in (:(:),))
isidentifier(x) = false

function qualifiedname(ts, name = nexttoken(ts))
  n = [name]
  pos = 0
  while true
    pos = position(ts.io)
    Lexer.next_token(ts) == :(.) || break
    t = Lexer.next_token(ts)
    isidentifier(t) || break
    push!(n, t)
  end
  seek(ts.io, pos)
  return length(n) > 1 ? n : n[1]
end

next_token′(ts) = try Lexer.next_token(ts) catch e return :error end

function nexttoken(ts)
  Lexer.skipws(ts)
  c = Lexer.peekchar(ts)

  t = c == '#' ? lexcomment(ts) :
      c == '"' ? lexstring(ts.io) :
      next_token′(ts)

  isidentifier(t) && (t = qualifiedname(ts, t))

  ts.lasttoken = t
  return t
end

function peektoken(ts)
  t = last = ts.lasttoken
  LNR.withstream(ts.io) do
    t = nexttoken(ts)
  end
  ts.lasttoken = last
  return t
end

# Scope parsing

immutable Scope
  kind::Symbol
  name::UTF8String
end

Scope(kind) = Scope(kind, "")
Scope(kind::Symbol, name::AbstractString) = Scope(kind, convert(UTF8String, name))
Scope(kind, name) = Scope(symbol(kind), string(name))

==(a::Scope, b::Scope) = a.kind == b.kind && a.name == b.name

const blockopeners = Set(map(symbol, ["begin", "function", "type", "immutable",
                                      "let", "macro", "for", "while",
                                      "quote", "if", "else", "elseif",
                                      "try", "finally", "catch", "do",
                                      "module"]))

const blockclosers = Set(map(symbol, ["end", "else", "elseif", "catch", "finally"]))

function nextscope!(scopes, ts)
  lasttoken = ts.lasttoken
  t = nexttoken(ts)
  if t in (:module, :baremodule) && isa(peektoken(ts), Symbol)
    push!(scopes, Scope(:module, nexttoken(ts)))
  elseif t in blockopeners
    push!(scopes, Scope(:block, t))
  elseif t in ('(', '[', '{')
    push!(scopes, Scope(:array, t))
  elseif last(scopes).kind in (:array, :call) && t in (')', ']', '}')
    pop!(scopes)
  elseif t == symbol("end")
    last(scopes).kind in (:block, :module) && lasttoken ≠ :(:) && pop!(scopes)
  elseif t == :using
    push!(scopes, Scope(:using))
  elseif t == '\n'
    last(scopes).kind == :using && pop!(scopes)
  elseif isidentifier(t) || isa(t, Vector{Symbol})
    if peektoken(ts) == '('
      nexttoken(ts)
      push!(scopes, Scope(:call, isa(t, Vector{Symbol}) ? join(collect(t), ".") : t ))
    end
  end
  return t
end

# API functions

Optional(T) = VERSION > v"0.4-" ? Union{T, Void} : Union(T, Void)

function scopes(code::LineNumberingReader, cur::Optional(Cursor) = nothing)
  ts = Lexer.TokenStream(code)
  scs = [Scope(:toplevel)]
  while cur == nothing || cursor(code) < cur
    Lexer.skipws(ts) == true && continue
    nextscope!(scs, ts) == Lexer.EOF && break
  end
  t = ts.lasttoken
  if isa(t, Token) && symbol(t) ≠ :whitespace && (cur == nothing ||
                                                  cur < cursor(code))
    push!(scs, Scope(t))
  elseif last(scs).kind == :call && !(cur == nothing ||
                                      cur >= cursor(code))
    # (only in a :call scope if past the last bracket)
    pop!(scs)
  end
  return scs
end

toLNR(s::AbstractString) = LineNumberingReader(s)
toLNR(r::LineNumberingReader) = r
toCursor(i::Integer) = cursor(i, 1)
toCursor(c::Cursor) = c
toCursor(::Void) = nothing

scopes(code, cur=nothing) = scopes(toLNR(code), toCursor(cur))

scope(code, cur=nothing) = last(scopes(code, cur))

function tokens(code::LineNumberingReader, cur::Optional(Cursor) = nothing)
  ts = Lexer.TokenStream(code)
  words = Set{UTF8String}()
  while true
    Lexer.skipws(ts); start = cursor(code)
    t = nexttoken(ts)
    t == Lexer.EOF && break
    (cur != nothing && start ≤ cur ≤ cursor(code)) && continue
    isa(t, Symbol) && push!(words, string(t))
    isa(t, Vector{Symbol}) && push!(words, join(t, "."))
  end
  return words
end

tokens(code, cur=nothing) = tokens(toLNR(code), toCursor(cur))


text:
@testset "The Graph Manifold" begin
    M1 = Euclidean(3)
    v1 = [1., 0., 0.]
    v2 = [0., 2., 0.]
    d1 = [0., 1., 1.]
    d2 = [1., 0., 1.]

    x1 = RnPoint(v1)
    x2 = RnPoint(v2)
    ξ1 = RnTVector(d1)
    ξ2 = RnTVector(d2)

    A = [0 1;1 0]
    M = Graph(M1,A)
    x = GraphVertexPoint([x1,x2])
    y = GraphVertexPoint([x2,x1])
    ξ = GraphVertexTVector([ξ1,ξ2])
    η = GraphVertexTVector([ξ2,ξ1])

    @test getValue(x) == [x1,x2]
    @test getValue(y) == [x2,x1]
    @test getValue(ξ) == [ξ1,ξ2]
    @test getValue(η) == [ξ2,ξ1]

    @test distance(M,x,y) == sqrt(2)*distance(M1,x1,x2)
    @test dot(M,x,ξ,η) == dot(M1,x1,ξ1,ξ2) + dot(M1,x2,ξ2,ξ1)
    @test exp(M,x,ξ) == GraphVertexPoint([exp(M1,x1,ξ1),exp(M1,x2,ξ2)])
    @test log(M,x,y) == GraphVertexTVector([log(M1,x1,x2),log(M1,x2,x1)])

    @test manifoldDimension(M) == 6
    @test manifoldDimension(x) == 6

    @test norm(M,x,ξ) == sqrt(dot(M,x,ξ,ξ))

    # parallelTransport
    # random
    # validate

    @test typicalDistance(M) ≈ sqrt(6) atol=5*10^(-16)

    @test zeroTVector(M,x) == GraphVertexTVector([
        zeroTVector(M1,x1), zeroTVector(M1,x2)
    ])

    @test "$M" == "The manifold on vertices and edges of a graph of The 3-dimensional Euclidean space of (vertex manifold) dimension 6."
    @test "$x" == "GraphVertexV[Rn([1.0, 0.0, 0.0]), Rn([0.0, 2.0, 0.0])]"
    @test "$ξ" == "GraphVertexT[RnT([0.0, 1.0, 1.0]), RnT([1.0, 0.0, 1.0])]"
end


text:
using Test, YaoBlocks, YaoArrayRegister
using YaoBase

@testset "Yao/#186" begin
    @test getiparams(phase(0.1)) == 0.1
    @test getiparams(Val(2) * phase(0.1)) == ()
    @test_throws NotImplementedError setiparams!(rot(X, 0.5), :nothing)
    @test_throws NotImplementedError setiparams(rot(X, 0.5), :nothing)
end

@testset "block to matrix conversion" begin
    for each in [X, Y, Z, H]
        Matrix{ComplexF64}(each) == Matrix{ComplexF64}(mat(each))
    end
    @test eltype(mat(chain(X))) == ComplexF64
    @test eltype(mat(chain(X, Rx(0.5)))) == ComplexF64
end

@testset "apply lambda" begin
    r = rand_state(3)
    @test apply!(copy(r), put(1 => X)) ≈ apply!(copy(r), put(3, 1 => X))
    r2 = copy(r)
    @test apply(r, put(1 => X)) ≈ apply!(copy(r), put(3, 1 => X))
    @test r2.state == r.state
    f(x::Float32) = x
    @test_throws ErrorException apply!(r, f)
end

@testset "push tests" begin
    # copy return itself by default
    @test copy(X) === X

    # block type can be used as traits
    @test nqubits(X) == 1

    @test isunitary(XGate)
    @test isreflexive(XGate)
    @test ishermitian(XGate)
    @test setiparams(Rx(0.3), 0.5) == Rx(0.5)
    @test setiparams(+, Rx(0.3), 0.5) == Rx(0.8)
    @test YaoBlocks.parameters_range(chain(Z, shift(0.3), phase(0.2), Rx(0.5), time_evolve(X, 0.5), Ry(0.5))) == [(0.0, 2π), (0.0, 2π), (0.0, 2π), (-Inf, Inf), (0.0, 2π)]
end


text:
@testset "209.minimum-size-subarray-sum.jl" begin
    @test min_subarray_len(7, [2, 3, 1, 2, 4, 3]) == 2
    @test min_subarray_len(5, [2, 3, 1, 2, 4, 3]) == 2
    @test min_subarray_len(3, [2, 3, 1, 2, 4, 3]) == 1
    @test min_subarray_len(15, [2, 3, 1, 2, 4, 3]) == 6
end


text:
using StatsBase
using Test

@testset "ECDF" begin
    x = randn(10000000)
    fnecdf = ecdf(x)
    y = [-1.96, -1.644854, -1.281552, -0.6744898, 0, 0.6744898, 1.281552, 1.644854, 1.96]
    @test isapprox(fnecdf(y), [0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975], atol=1e-3)
    @test isapprox(fnecdf(1.96), 0.975, atol=1e-3)
    @test fnecdf(y) ≈ map(fnecdf, y)
    @test extrema(fnecdf) == (minimum(fnecdf), maximum(fnecdf)) == extrema(x)
    fnecdf = ecdf([0.5])
    @test fnecdf([zeros(5000); ones(5000)]) == [zeros(5000); ones(5000)]
    @test extrema(fnecdf) == (minimum(fnecdf), maximum(fnecdf)) == (0.5, 0.5)
end

@testset "Weighted ECDF" begin
    x = randn(10000000)
    w1 = rand(10000000)
    w2 = weights(w1)
    fnecdf = ecdf(x, weights=w1)
    fnecdfalt = ecdf(x, weights=w2)
    @test fnecdf.sorted_values == fnecdfalt.sorted_values
    @test fnecdf.weights == fnecdfalt.weights
    @test fnecdf.weights != w1  #  check that w wasn't accidently modified in place
    @test fnecdfalt.weights != w2
    y = [-1.96, -1.644854, -1.281552, -0.6744898, 0, 0.6744898, 1.281552, 1.644854, 1.96]
    @test isapprox(fnecdf(y), [0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975], atol=1e-3)
    @test isapprox(fnecdf(1.96), 0.975, atol=1e-3)
    @test fnecdf(y) ≈ map(fnecdf, y)
    @test extrema(fnecdf) == (minimum(fnecdf), maximum(fnecdf)) == extrema(x)
    fnecdf = ecdf([1.0, 0.5], weights=weights([3, 1]))
    @test fnecdf(0.75) == 0.25
    @test extrema(fnecdf) == (minimum(fnecdf), maximum(fnecdf)) == (0.5, 1.0)
    @test_throws ArgumentError ecdf(rand(8), weights=weights(rand(10)))
    #  Check frequency weights
    v = randn(100)
    r = rand(1:100, 100)
    vv = vcat(fill.(v, r)...)  #  repeat elements of v according to r
    fw = fweights(r)
    frecdf1 = ecdf(v, weights=fw)
    frecdf2 = ecdf(vv)
    @test frecdf1(y) ≈ frecdf2(y)
    #  Check probability weights
    a = randn(100)
    b = rand(100)
    b̃ = abs(10randn()) * b
    bw1 = pweights(b)
    bw2 = pweights(b̃)
    precdf1 = ecdf(a, weights=bw1)
    precdf2 = ecdf(a, weights=bw2)
    @test precdf1(y) ≈ precdf2(y)
end


text:
# This file was generated, do not modify it. # hide
using PyPlot
figure(figsize=(8, 6))
x = range(-2, 2, length=500)
for α in 1:5
    plot(x, sinc.(α .* x))
end
savefig(joinpath(@OUTPUT, "sinc.svg")) # hide

max_stars_repo_name: mcognetta/franklindocs

max_issues_repo_name: mcognetta/franklindocs

max_forks_repo_name: mcognetta/franklindocs


text:
using DataDeps

register(DataDep(
  "thyroid-disease-thyroid-0387",
  "http://archive.ics.uci.edu/ml/datasets/Thyroid+Disease",
  "http://archive.ics.uci.edu/ml/machine-learning-databases/thyroid-disease/thyroid0387.data",
  "4a20d104d4607afd7b1a477af5aa4381a2027f36cc7027f44bc2bd5a50c6724d",
  post_fetch_method=(path -> begin
    # The ids are appended to each row inside []. Replace with comma separation.
    write(path, replace(replace(read(path, String), "[" => ","), "]" => ""))

    UCIData.process_dataset(path,
        target_index=30,
        feature_indices=1:29,
        id_indices=31,
        categoric_indices=[2:17; 19; 21; 23; 25; 27; 29],
    )
  end),
))

max_stars_repo_name: JackDunnNZ/uci-data

max_issues_repo_name: JackDunnNZ/uci-data

max_forks_repo_name: JackDunnNZ/uci-data