text:
module MTL

import Libdl
using CEnum
using Scratch
using JSON

const libcmt = joinpath(@get_scratch!("cmt"), "lib", "libcmt.dylib")

## source code includes

# low-level autogenerated wrappers
export MtSize, MtlDim, MtlDim3
include("libcmt.jl")

# convert from 1 based indexing to 0 based indexing
Base.convert(::Type{NsRange}, range::UnitRange{T}) where T <: Integer =
	NsRange(first(range), length(range))
# used for byte ranges.
Base.convert(::Type{NsRange}, range::StepRange{T}) where T <: Integer =
	NsRange(first(range)-step(range), length(range)*step(range))

# Basic types
include("error.jl")
include("helpers.jl")

# julia wrappers
include("storage_type.jl")
include("resource.jl")
include("device.jl")
include("compile-opts.jl")
include("library.jl")
include("function.jl")
include("events.jl")
include("fences.jl")
include("heap.jl")
include("buffer.jl")
include("command_queue.jl")
include("command_buf.jl")
include("compute-pipeline.jl")
include("argument_encoder.jl")

include("command_enc.jl")
include("command_enc_blit.jl")
include("command_enc_compute.jl")

include("reflection/reflection_compute.jl")

include("size.jl")

end # module


text:
# %% imports

include("cmaes.jl")
include("CMAES_GaussGrid_decision.jl")
include("save_param.jl")

# %%

function optim()
    tunekeys = ["vthresh", "ae", "b", "c", "de","imin", "imax",
                "m", "std", "apost", "apre"]
    cfg = YAML.load(open("/Users/titou/Documents/PFE/RL-STDP/Julia/CartPole/CMAES/cfg_cmaes.yml"))
    N=length(tunekeys)
    best = Array{Float64}[]
    for i in 1:10
        push!(best,[])
    end
    best_fit = [-1.0 for i in 1:10]
    init_params = [cfg[tunekeys[idx]]*1.0 for idx in 1:N]
    init_genes  = descale(init_params, tunekeys)
    lambda = 4 + floor(Int64, 3*log(N))
    mu = floor(Int64, lambda/2)
    c = CMAES(init_genes, N=N, μ=mu, λ=lambda, τ=sqrt(N), τ_c=N^2, τ_σ=sqrt(N))
    for i in 1:1000
        @time step!(c)
        bestind = argmin(c.F_λ)
        maxfit = -c.F_λ[bestind]
        println(i," ", maxfit, " (current best = $(maximum(best_fit)))")
        if maxfit > best_fit[1]
            best_fit[1] = maxfit
            best[1] = copy(c.offspring[bestind])
            idx = sortperm(best_fit)
            best_fit = best_fit[idx]
            best = best[idx]
        end
        if length(best[1]) > 0 && i%10 == 0
            save_param(best_fit,best,tunekeys)
        end
        if sum(best_fit) > 800
            break
        end
    end
    return best_fit, best
end

# %% run optim

best_fit, best = optim()

# %%

function viz_optim(best)
    tunekeys = ["vthresh", "ae", "b", "c", "de","imin", "imax"
                 "m", "std", "apost", "apre"]
    for i in 1:length(best)
        params = Dict(tunekeys .=> rescale(best[i], tunekeys))
        @show params
    end
end

# %%
viz_optim(best)

ext: jl

max_stars_repo_name: SuReLI/RL-STDP

max_issues_repo_name: SuReLI/RL-STDP

max_forks_repo_name: SuReLI/RL-STDP


text:
function ci_matrix(b::AbstractBasis, orbitals::AbstractMatrix, Vn; Ve=x->exp(-x^2))
    @argcheck length(b) == size(orbitals, 1)
    # index conversion from HF space to CI space
    indx = [ (i,j) for i in axes(orbitals, 2) for j in axes(orbitals, 2) ]
    w = get_weight(b)
    ci = zeros(length(indx), length(indx))

    h1 = one_electron_operator(b, Vn)

    # Electron repulsion tensor to help calculation
    ve = zeros( eltype(orbitals), (length(b), length(b)) )
    for j in axes(ve,2)
        ve[:,j] = Ve.(b[j].-b)
    end

    Threads.@threads for i in axes(ci,1)
        ϕ₁ = @view orbitals[:,indx[i][1]]
        ϕ₂ = @view orbitals[:,indx[i][2]]
        for j in axes(ci,2)
            ψ₁ = @view orbitals[:,indx[j][1]]
            ψ₂ = @view orbitals[:,indx[j][2]]
            # Calculate Coulomb repulsion
            @tullio tmp = w[n] * ϕ₁[n] * ψ₁[n] * ve[n,m] * ϕ₂[m] * ψ₂[m] * w[m]
            # Orbitals are orthogonal
            if indx[i][1] == indx[j][1]
                tmp += ϕ₂' * h1 * ψ₂
            end
            if indx[i][2] == indx[j][2]
                tmp += ϕ₁' * h1 * ψ₁
            end
            ci[i,j] = tmp
        end
    end
    return ci, indx
end

function full_ci(b::AbstractBasis, orbitals::AbstractMatrix, Vn; Ve=x->exp(-x^2))
    # We need to adjust threads so that we don't oversubscribe
    bt = BLAS.get_num_threads()
    nt = Threads.nthreads()
    BLAS.set_num_threads(1)

    @info "Creating CI matrix"
    @time ci, indx = ci_matrix(b, orbitals, Vn; Ve=Ve)

    BLAS.set_num_threads(nt)

    @info "Diagonalizing CI matrix"
    @time ev, ve = eigen(ci)

    BLAS.set_num_threads(bt)

    return Dict(
        "orbitals"=>orbitals,
        "energies"=>ev,
        "states"=>ve,
        "indexes"=>indx
    )
end

function reduced_ci_orbitals(b::BasisLobatto, orbitals::AbstractMatrix; pointlike=false)
    @argcheck length(b) == size(orbitals, 1)
    # Index range, for elements-2
    function _index_gen(l, i)
        if i == 0
            return 1
        else
            return l + _index_gen(l, i-1) - 1
        end
    end
    # rows for each element contributions
    function _row_range(l, i)
        return 2+(i-1)*(l-2):i*(l-2)+1
    end
    nelements = length(b.egvector)
    new_orbitals = zeros(length(b), length(b)-nelements)
    le = length(get_element(b, 1))  # all elements have same ammount of points
    w = get_weight(b)
    if pointlike
        for ne in 1:nelements
            nr = _index_gen(le, ne-1):_index_gen(le, ne)
            tmp = diagm(ones(le))[:,begin:end-1]
            tmp[:,1] = orbitals[nr]
            t = gram_schmit(tmp, w[nr])
            new_orbitals[nr, _row_range(le, ne)] = t[:,begin+1:end]
        end
    else
        for ne in 1:nelements
            nr = _index_gen(le, ne-1):_index_gen(le, ne)
            tmp = zeros(length(b.egvector[ne]), length(b.egvector[ne])-1)
            tmp[:,1] = orbitals[nr]
            for i in 2:size(tmp,2)
                tmp[:, i] = particle_in_box(b.egvector[ne], i-1)
            end
            t = gram_schmit(tmp, w[nr])
            new_orbitals[nr, _row_range(le, ne)] = t[:,begin+1:end]
        end
    end
    new_orbitals[:,1] = orbitals[:,1]
    return gram_schmit(new_orbitals, w)
end

function boundary_ci_orbitals(b::BasisLobatto, orbitals::AbstractMatrix)
    out = zeros(length(b), number_of_elements(b)+2)
    out[:,1] = orbitals[:,1] # HF ground state
    out[1,2] = 1
    nb = 1
    for j in 1:number_of_elements(b)
        le = length( get_element(b,j) )
        nb += le - 1
        out[nb,j+2] = 1
    end
    return gram_schmit(out, get_weight(b))
end

function element_ci_orbitals(b::BasisLobatto, orbitals::AbstractMatrix, i)
    @argcheck length(b) == size(orbitals, 1)
    collumns = get_element_indexes(b, i)
    norbs = zeros(length(b), length(collumns)+1)
    norbs[:,1] = orbitals[:,1]
    for i in axes(collumns,1)
        norbs[collumns[i], i+1] = 1
    end
    return gram_schmit(norbs, get_weight(b))
end

"""
    block_ci(b::BasisLobatto, orbitals::AbstractMatrix, V; ne=2)

Calculate CI on element blocks

# Arguments
- `b::BasisLobatto`            :  Basis for the calculation
- `orbitals::AbstractMatrix`   :  First collumn is taken as initial state
- `V`                          :  Nuclear potential

# Keywords
- `ne=2`  :  number of elements in a block
"""
function block_ci(b::BasisLobatto, orbitals::AbstractMatrix, V; ne=2, nstates=1)
    @argcheck 1 <= ne < number_of_elements(b)
    @argcheck nstates >= 1
    ib = [i for i in 1:ne]
    rci = []
    @info "Calculating CI blocks"
    for i in 0:number_of_elements(b)-ne
        tmp = element_ci_orbitals(b, orbitals, ib .+ i )
        push!(rci, full_ci(b, tmp, V) )
    end

    @info "Combining blocks"
    cie = CIHamilton(b, V)
    cio = CIOverlap(b)

    l = length(rci)*nstates
    HH = zeros(l,l)
    SS = zeros(l,l)

    p = Progress( Int( l*(l-1)/2) )
    for i in axes(rci, 1)
        for j in i:length(rci)
            for ii in 1:nstates
                for jj in 1:nstates
                    HH[ (i-1)*nstates+ii, (j-1)*nstates+jj ] = ci_vector_product(cie, rci[i], rci[j]; i1=ii, i2=jj)
                    SS[ (i-1)*nstates+ii, (j-1)*nstates+jj ] = ci_vector_product(cio, rci[i], rci[j]; i1=ii, i2=jj)
                    next!(p)
                end
            end
        end
    end
    @info "Calculating final state"
    H = Symmetric(HH)
    S = Symmetric(SS)
    S_inv = inv(S)

    e, v = eigen(S_inv*H)
    return Dict("energy"=>e[1:nstates], "state"=>v[:,1:nstates])
end

function ci_vector_product(op, ci1, ci2; i1=1, i2=1)
    states1 = ci1["states"][:,i1]
    states2 = ci2["states"][:,i2]

    orbitals1 = ci1["orbitals"]
    orbitals2 = ci2["orbitals"]

    vector1 = ci1["indexes"]
    vector2 = ci2["indexes"]

    out = Threads.Atomic{Float64}(0)
    Threads.@threads for i in length(vector1):-1:1
        ψ1 = @view orbitals1[:, vector1[i][1] ]
        ψ2 = @view orbitals1[:, vector1[i][2] ]
        tmp = 0.0
        for j in reverse( axes(vector2,1) )
            ϕ1 = @view orbitals2[:, vector2[j][1] ]
            ϕ2 = @view orbitals2[:, vector2[j][2] ]
            tmp += op( ψ1, ψ2, ϕ1, ϕ2 ) * ( states1[i] * states2[j] )
        end
        Threads.atomic_add!(out, tmp)
    end
    return out.value
end

"""
    CIHamilton

Used to calculate energy on CI-vector elements

# Fields
- `w::Vector{Float64}`   : integration weights
- `h1::Matrix{Float64}`  : one electron Hamilton
- `ve::Matrix{Float64}`  : electron-electron repulsion

"""
struct CIHamilton
    w::Vector{Float64}
    h1::Matrix{Float64}
    ve::Matrix{Float64}
    function CIHamilton(b::AbstractBasis, Vn, Ve=x->exp(-x^2))
        h1 = one_electron_operator(b, Vn)
        w = get_weight(b)

        # Electron repulsion tensor
        ve = zeros( (length(b), length(b)) )
        for j in axes(ve,2)
            ve[:,j] = Ve.(b[j].-b)
        end
        new(w, h1, ve)
    end
end

function (cie::CIHamilton)(psi1, psi2, phi1, phi2)
    @tullio tmp = cie.w[n] * psi1[n] * phi1[n] * cie.ve[n,m] * psi2[m] * phi2[m] * cie.w[m]
    @tullio s1 = phi1[n] * psi1[n] * cie.w[n]
    @tullio s2 = phi2[n] * psi2[n] * cie.w[n]
    tmp += ( psi1' * cie.h1 * phi1 ) * s2
    tmp += ( psi2' * cie.h1 * phi2 ) * s1
    return tmp
end

"""
    CIOverlap

Calculate overlap on CI-vector elements
"""
struct CIOverlap
    g::Diagonal{Float64, Vector{Float64}}
    function CIOverlap(b::AbstractBasis)
        new(metric_tensor(b))
    end
end

function (cio::CIOverlap)(psi1, psi2, phi1, phi2)
    return ( psi1' * cio.g * phi1 ) * ( psi2' * cio.g * phi2 )
end

max_stars_repo_name: tjjarvinen/OneDimensionalCorrelation

max_issues_repo_name: tjjarvinen/OneDimensionalCorrelation

max_forks_repo_name: tjjarvinen/OneDimensionalCorrelation


text:
module ReflectionTests
using  Base.Test
using  ExpressionPatterns
using  ExpressionPatterns.Dispatch
using  ExpressionPatterns.Dispatch.Reflection

# pattern preference

Dispatch.set_conflict_warnings(:no)

xpy = :(x+y)
ypx = :(y+x)

@macromethod f(x+y, z)  = [x,y]
@macromethod f(z, x+y)  = [x,y]

@metafunction f(x+y, z) = [x,y]
@metafunction f(z, x+y) = [x,y]

@test @f(x+y,y+x) == [:x, :y]
@test  f(xpy,ypx) == [:x, :y]

@prefer @f(z, x+y) over @f(x+y, z)
@prefer  f(z, x+y) over  f(x+y, z)

@test @f(x+y,y+x) == [:y,:x]
@test  f(xpy,ypx) == [:y,:x]

@prefer @f(x+y, z) over @f(z, x+y)
@prefer  f(x+y, z) over  f(z, x+y)

@test @f(x+y,y+x) == [:x, :y]
@test  f(xpy,ypx) == [:x, :y]

# label preference

@macromethod g(x+y, z)[a] = [x,y]
@macromethod g(z, x+y)[b] = [x,y]

@metafunction g(x+y, z)[a] = [x,y]
@metafunction g(z, x+y)[b] = [x,y]

@test @g(x+y,y+x) == [:x, :y]
@test  g(xpy,ypx) == [:x, :y]

@prefer b over a in @g
@prefer b over a in  g

@test @g(x+y,y+x) == [:y, :x]
@test  g(xpy,ypx) == [:y, :x]

@prefer a over b in @g
@prefer a over b in  g

@test @g(x+y,y+x) == [:x, :y]
@test  g(xpy,ypx) == [:x, :y]

# method removal

@macromethod h(x+y)[plus] = [x,y]
@macromethod h(x)[any]    = [x]

@metafunction h(x+y)[plus] = [x,y]
@metafunction h(x)[any]    = [x]

@test @h(x+y)    == [:x, :y]
@test  h(:(x+y)) == [:x, :y]

@remove plus from @h
@remove plus from  h

@test @h(x+y)    == [:(x+y)]
@test  h(:(x+y)) == [:(x+y)]

@remove @h(x)
@remove  h(x)

@test_throws MetaMethodError try @eval @h(x); catch err throw(err isa LoadError ? err.error : err) end
@test_throws MetaMethodError h(:x)

# which

MM = Dispatch.Applications.MACROMETHODS[ReflectionTests]
MF = Dispatch.Applications.METAFUNCTIONS[ReflectionTests]

@macromethod m(x)   = [x]
@macromethod m(x+y) = [x,y]

@metafunction m(x)   = [x]
@metafunction m(x+y) = [x,y]

const mmac = Symbol("@m")
const mfun = :m

@test MM[mmac].methods[1] == @whichmeta @m(x+y)
@test MM[mmac].methods[2] == @whichmeta @m(x)

@test MF[mfun].methods[1] == @whichmeta m(x+y)
@test MF[mfun].methods[2] == @whichmeta m(x)

# conflicts

@macromethod n(x,x+y)[a] ()
@macromethod n(x+y,x)[b] ()

let stdout = STDOUT
  conflicts_io, = redirect_stdout()
  @metaconflicts @n
  print_with_color(:red,    "[a]<x x + y>")
  print_with_color(:yellow, " | ")
  print_with_color(:red,    "[b]<x + y x>")
  println()

  metaconflicts_result = readline(conflicts_io)
  readline(conflicts_io)
  expected = readline(conflicts_io)

  redirect_stdout(stdout)
  @test metaconflicts_result == expected
end

# Cross-module reflection

module A
using ..Dispatch
export @f, @g

  @macromethod f(x) (x,)
  @macromethod g(x) (x,)

end

module B
using  ExpressionPatterns.Dispatch
using  ExpressionPatterns.Dispatch.Reflection

  @metamodule import ..A.@f

  @macromethod f(x+y) (x,y)
  @macromethod g(x+y) (x,y)

  @macromethod h(x+y, z) (x,y)
  @macromethod h(z, x+y) (x,y)

  @macromethod m(x+y, z)[a] (x,y)
  @macromethod m(z, x+y)[b] (x,y)

end

@test A.@f(x+y) == B.@f(x+y)
@test A.@g(x+y) != B.@g(x+y)

@test B.@h(x+y,y+x) == (:x,:y)

@prefer B.@h(z, x+y) over B.@h(x+y, z)

@test B.@h(x+y,y+x) == (:y,:x)

@test B.@m(x+y,y+x) == (:x,:y)

@prefer b over a in B.@m

@test B.@m(x+y,y+x) == (:y,:x)

@prefer a over b in B.@m

@test B.@m(x+y,y+x) == (:x,:y)

# Errors and warnings

macro test_warning(ex, expected_warning::String)
  @gensym had_color
  quote
    let
      token  = gensym()
      stderr = STDERR
      nerr,  = redirect_stderr()

      $(esc(had_color)) = Base.have_color
      eval($Base, :(have_color = false))
      try
        $(esc(ex))
        println(STDERR)
        println(STDERR, "$token")
        warning = readline(nerr)
        line = (eof(nerr) ? "" : readline(nerr))
        while line != "$token"
          if line != ""
            warning = "$warning\n$line"
          end
          line = readline(nerr)
        end
        @test warning == $expected_warning
      finally
        redirect_stderr(stderr)
        eval($Base, :(have_color = $$had_color))
      end
    end
  end
end

@test_throws ArgumentError ExpressionPatterns.Dispatch.TableManipulation.set_conflict_warnings(:bleh)

ExpressionPatterns.Dispatch.TableManipulation.set_conflict_warnings(:yes)

@macromethod m()[a] ()
@test_warning (@prefer a over a in @m) "WARNING: Using `prefer` with non-conflicting methods. (pattern`()` and pattern`()` in macromethod m)"
@test_warning (@remove b from @m)      "WARNING: Couldn't remove the method [b] from the macromethod m"

ExpressionPatterns.Dispatch.TableManipulation.set_conflict_warnings(:no)

end


text:
using SortingAlgorithms
using Test
using StatsBase
using Random

a = rand(1:10000, 1000)

for alg in [TimSort, HeapSort, RadixSort]
    b = sort(a, alg=alg)
    @test issorted(b)
    ix = sortperm(a, alg=alg)
    b = a[ix]
    @test issorted(b)
    @test a[ix] == b

    b = sort(a, alg=alg, rev=true)
    @test issorted(b, rev=true)
    ix = sortperm(a, alg=alg, rev=true)
    b = a[ix]
    @test issorted(b, rev=true)
    @test a[ix] == b

    b = sort(a, alg=alg, by=x->1/x)
    @test issorted(b, by=x->1/x)
    ix = sortperm(a, alg=alg, by=x->1/x)
    b = a[ix]
    @test issorted(b, by=x->1/x)
    @test a[ix] == b

    c = copy(a)
    permute!(c, ix)
    @test c == b

    invpermute!(c, ix)
    @test c == a

    if alg != RadixSort  # RadixSort does not work with Lt orderings
        c = sort(a, alg=alg, lt=(>))
        @test b == c
    end

    c = sort(a, alg=alg, by=x->1/x)
    @test b == c
end

randnans(n) = reinterpret(Float64,[rand(UInt64)|0x7ff8000000000000 for i=1:n])

function randn_with_nans(n,p)
    v = randn(n)
    x = findall(rand(n).<p)
    v[x] = randnans(length(x))
    return v
end

Random.seed!(0xdeadbeef)

for n in [0:10..., 100, 101, 1000, 1001]
    r = 1:10
    v = rand(1:10,n)
    h = fit(Histogram, v, r)

    for ord in [Base.Order.Forward, Base.Order.Reverse]
        # insertion sort (stable) as reference
        pi = sortperm(v, alg=InsertionSort, order=ord)
        @test isperm(pi)
        si = v[pi]
        @test fit(Histogram, si, r) == h
        @test issorted(si, order=ord)
        @test all(issorted,[pi[si.==x] for x in r])
        c = copy(v)
        permute!(c, pi)
        @test c == si
        invpermute!(c, pi)
        @test c == v

        # stable algorithms
        for alg in [TimSort, RadixSort]
            p = sortperm(v, alg=alg, order=ord)
            @test p == pi
            s = copy(v)
            permute!(s, p)
            @test s == si
            invpermute!(s, p)
            @test s == v
        end

        # unstable algorithms
        for alg in [HeapSort]
            p = sortperm(v, alg=alg, order=ord)
            @test isperm(p)
            @test v[p] == si
            s = copy(v)
            permute!(s, p)
            @test s == si
            invpermute!(s, p)
            @test s == v
        end
    end

    v = randn_with_nans(n,0.1)
    for ord in [Base.Order.Forward, Base.Order.Reverse],
        alg in [TimSort, HeapSort, RadixSort]
        # test float sorting with NaNs
        s = sort(v, alg=alg, order=ord)
        @test issorted(s, order=ord)
        @test reinterpret(UInt64, v[map(isnan, v)]) == reinterpret(UInt64, s[map(isnan, s)])

        # test float permutation with NaNs
        p = sortperm(v, alg=alg, order=ord)
        @test isperm(p)
        vp = v[p]
        @test isequal(vp,s)
        @test reinterpret(UInt64,vp) == reinterpret(UInt64,s)
    end
end


text:
using CSV, DataFrames, Plots

clusfile = "Data/Clustered/Healthy/Healthy-chr1-clusters.tsv"
print(clusfile)
df = CSV.read(clusfile, delim=',', copycols=true)
print(first(df,5))

nclusters = maximum(df[!,"clusterid"])
sort!(df, [:chromosome, :clusterid])

df = df[1:5,:]
M = Array{Int64}(zeros(nrow(df), nrow(df)))
genes = df[!,:gname]
cl = df[:,:clusterid]
for i in 1:length(genes)
	clidsrc = cl[i]
	M[i,i] = cl[i]
	for j in i+1:length(genes)
		# println("$i $j")
		cliddst = cl[j]
		clidsrc == cliddst ? M[i,j] = clidsrc : M[i,j] = 0
	end
end

# print(M)
# myplot = spy(M)
# draw(PNG("Plots/myplot.png", 3inch, 3inch), myplot)
heatmap(M, title="Threshold")
plot(size=(1200,1200))

savefig("Plots/clustering.png")

max_stars_repo_name: josemaz/gene-matrices

max_issues_repo_name: josemaz/gene-matrices

max_forks_repo_name: josemaz/gene-matrices


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 845b7d0a-7ca2-11eb-2683-cd370811bf68
using NiLang

# ╔═╡ 88912bea-7ca2-11eb-1cde-db111d594c20
using Viznet, PlutoUI, Compose

# ╔═╡ f42405a8-7ca2-11eb-133b-df4f1ce9bb19
html"""
download this notebook
"""

# ╔═╡ 66e972c2-7ca2-11eb-0a7a-a9305dd20511
md"# BBMCA - NiLang implementation"

# ╔═╡ c8722cea-7ca5-11eb-3e7d-0f39322ea40a
md"Check [Physics-like models of computation](https://www.sciencedirect.com/science/article/abs/pii/0167278984902525) (Norman Margolus, 1984) for theories about Billiard ball celluar automata (BBMCA)."

# ╔═╡ 845d1458-7ca2-11eb-0e69-11a10b9894a8
@i function load_and_clear!(x, config, i, j)
	m, n ← size(config)
	x ⊻= config[i,j]  # to make it faster, should put `@inbounds` before it
	config[i,j] ⊻= x
	x ⊻= config[i,mod1(j+1, n)] << 1
	config[i,mod1(j+1, n)] ⊻= x >> 1
	x ⊻= config[mod1(i+1, m),j] << 2
	config[mod1(i+1, m),j] ⊻= x >> 2
	x ⊻= config[mod1(i+1, m),mod1(j+1, n)] << 3
	config[mod1(i+1, m),mod1(j+1, n)] ⊻= x >> 3
end

# ╔═╡ 84680136-7ca2-11eb-1df5-ffa0b4b9d126
@i function margolus_rule(y, x)
	# remove reversibility check to make it run faster
    @invcheckoff if x==6
        y ⊻= 9
    elseif x==9
        y ⊻= 6
    elseif x==4
        y ⊻= 2
    elseif x==2
        y ⊻= 4
    elseif x==1
        y ⊻= 8
    elseif x==8
        y ⊻= 1
    else
        y ⊻= x
    end
end

# ╔═╡ 845c1292-7ca2-11eb-3e56-996aa5229b4e
@i function update_bbmca!(config, iseven)
	# computing offsets, and borrow some ancillas from system
	@routine begin
		offset ← 1
		m, n ← size(config)
		if !iseven
			offset += 1
		end
	end
	for j=offset:2:n
		for i=offset:2:m
			x ← 0
			y ← 0
			# load block to `x` and clean up original data
			load_and_clear!(x, config, i, j)
			# compute new config to `y`
			margolus_rule(y, x)
			# clean up `x` with the following observation:
			# applying margolus rule twice restores the configuration
			margolus_rule(x, y)
			# store `y` to block
			(~load_and_clear!)(y, config, i, j)
			# ancillas `x` and `y` are returned to the pool automatically
		end
	end
	# uncompute `offset`
	~@routine
end

# ╔═╡ 34006d60-7ca3-11eb-3c51-c9758394b838
md"# Visualization"

# ╔═╡ 846e21da-7ca2-11eb-05a3-51e22ed04147
function showconfig(ba::AbstractMatrix)
    m, n = size(ba, 1), size(ba, 2)
    lt = Viznet.SquareLattice(n, m)
    brush1 = nodestyle(:square, fill("black"), stroke("#888888"), linewidth(unit(lt)*mm); r=unit(lt)/2.2)
    brush0 = nodestyle(:square, fill("white"), stroke("#888888"), linewidth(unit(lt)*mm); r=unit(lt)/2.2)
    canvas() do
        for i=1:m, j=1:n
            (ba[i, j] == 1 ? brush1 : brush0) >> lt[j,i]
        end
    end
end

# ╔═╡ 846e8170-7ca2-11eb-351e-eb45c629f6a6
@bind btn Clock(0.1)

# ╔═╡ 84730f04-7ca2-11eb-3b13-3fd82a9e2109
# initial configuration
config = let
	x=zeros(Int, 10, 10)
	x[1,1] = 1
	x
end;

# ╔═╡ 84763e9c-7ca2-11eb-356c-c391966cdc98
# parity - Note: BBMCA is a two state CA
bbmca_parity = Ref(true)

# ╔═╡ 847711e6-7ca2-11eb-326a-15f6bfc05347
let
	btn
	# update
	update_bbmca!(config, bbmca_parity[])
	# change parity
	bbmca_parity[] = !(bbmca_parity[])
	# visualize
	Compose.set_default_graphic_size(10cm, 10cm)
	showconfig(config)
end

# ╔═╡ Cell order:
# ╟─f42405a8-7ca2-11eb-133b-df4f1ce9bb19
# ╟─66e972c2-7ca2-11eb-0a7a-a9305dd20511
# ╟─c8722cea-7ca5-11eb-3e7d-0f39322ea40a
# ╠═845b7d0a-7ca2-11eb-2683-cd370811bf68
# ╠═845c1292-7ca2-11eb-3e56-996aa5229b4e
# ╠═845d1458-7ca2-11eb-0e69-11a10b9894a8
# ╠═84680136-7ca2-11eb-1df5-ffa0b4b9d126
# ╟─34006d60-7ca3-11eb-3c51-c9758394b838
# ╠═88912bea-7ca2-11eb-1cde-db111d594c20
# ╠═846e21da-7ca2-11eb-05a3-51e22ed04147
# ╟─846e8170-7ca2-11eb-351e-eb45c629f6a6
# ╠═84730f04-7ca2-11eb-3b13-3fd82a9e2109
# ╠═84763e9c-7ca2-11eb-356c-c391966cdc98
# ╠═847711e6-7ca2-11eb-326a-15f6bfc05347


text:
using DitherPunk

A = [1 2; 3 4]

@test upscale(A, 3) == [
    1 1 1 2 2 2
    1 1 1 2 2 2
    1 1 1 2 2 2
    3 3 3 4 4 4
    3 3 3 4 4 4
    3 3 3 4 4 4
]


text:
module AssemblyView

export AssemblyOperand, AssemblyRegister, AssemblyMemoryOperand,
    AssemblyImmediate, AssemblyOffset, AssemblyStatement, AssemblyComment,
    AssemblyLabel, AssemblyInstruction,
    parsed_asm, asm_offsets, view_asm

using InteractiveUtils: _dump_function

################################################################# REGISTER NAMES

const X86_REGISTER_NAMES = Dict{String,Symbol}(
    "ah" => :A, "al" => :A, "ax" => :A, "eax" => :A, "rax" => :A,
    "ch" => :C, "cl" => :C, "cx" => :C, "ecx" => :C, "rcx" => :C,
    "dh" => :D, "dl" => :D, "dx" => :D, "edx" => :D, "rdx" => :D,
    "bh" => :B, "bl" => :B, "bx" => :B, "ebx" => :B, "rbx" => :B,
    "spl" => :SP, "sp" => :SP, "esp" => :SP, "rsp" => :SP,
    "bpl" => :BP, "bp" => :BP, "ebp" => :BP, "rbp" => :BP,
    "sil" => :SI, "si" => :SI, "esi" => :SI, "rsi" => :SI,
    "dil" => :DI, "di" => :DI, "edi" => :DI, "rdi" => :DI,
    "r8b" => :R8, "r8w" => :R8, "r8d" => :R8, "r8" => :R8,
    "r9b" => :R9, "r9w" => :R9, "r9d" => :R9, "r9" => :R9,
    "r10b" => :R10, "r10w" => :R10, "r10d" => :R10, "r10" => :R10,
    "r11b" => :R11, "r11w" => :R11, "r11d" => :R11, "r11" => :R11,
    "r12b" => :R12, "r12w" => :R12, "r12d" => :R12, "r12" => :R12,
    "r13b" => :R13, "r13w" => :R13, "r13d" => :R13, "r13" => :R13,
    "r14b" => :R14, "r14w" => :R14, "r14d" => :R14, "r14" => :R14,
    "r15b" => :R15, "r15w" => :R15, "r15d" => :R15, "r15" => :R15,
    "xmm0" => :SSE0, "ymm0" => :SSE0, "zmm0" => :SSE0,
    "xmm1" => :SSE1, "ymm1" => :SSE1, "zmm1" => :SSE1,
    "xmm2" => :SSE2, "ymm2" => :SSE2, "zmm2" => :SSE2,
    "xmm3" => :SSE3, "ymm3" => :SSE3, "zmm3" => :SSE3,
    "xmm4" => :SSE4, "ymm4" => :SSE4, "zmm4" => :SSE4,
    "xmm5" => :SSE5, "ymm5" => :SSE5, "zmm5" => :SSE5,
    "xmm6" => :SSE6, "ymm6" => :SSE6, "zmm6" => :SSE6,
    "xmm7" => :SSE7, "ymm7" => :SSE7, "zmm7" => :SSE7,
    "xmm8" => :SSE8, "ymm8" => :SSE8, "zmm8" => :SSE8,
    "xmm9" => :SSE9, "ymm9" => :SSE9, "zmm9" => :SSE9,
    "xmm10" => :SSE10, "ymm10" => :SSE10, "zmm10" => :SSE10,
    "xmm11" => :SSE11, "ymm11" => :SSE11, "zmm11" => :SSE11,
    "xmm12" => :SSE12, "ymm12" => :SSE12, "zmm12" => :SSE12,
    "xmm13" => :SSE13, "ymm13" => :SSE13, "zmm13" => :SSE13,
    "xmm14" => :SSE14, "ymm14" => :SSE14, "zmm14" => :SSE14,
    "xmm15" => :SSE15, "ymm15" => :SSE15, "zmm15" => :SSE15,
    "xmm16" => :SSE16, "ymm16" => :SSE16, "zmm16" => :SSE16,
    "xmm17" => :SSE17, "ymm17" => :SSE17, "zmm17" => :SSE17,
    "xmm18" => :SSE18, "ymm18" => :SSE18, "zmm18" => :SSE18,
    "xmm19" => :SSE19, "ymm19" => :SSE19, "zmm19" => :SSE19,
    "xmm20" => :SSE20, "ymm20" => :SSE20, "zmm20" => :SSE20,
    "xmm21" => :SSE21, "ymm21" => :SSE21, "zmm21" => :SSE21,
    "xmm22" => :SSE22, "ymm22" => :SSE22, "zmm22" => :SSE22,
    "xmm23" => :SSE23, "ymm23" => :SSE23, "zmm23" => :SSE23,
    "xmm24" => :SSE24, "ymm24" => :SSE24, "zmm24" => :SSE24,
    "xmm25" => :SSE25, "ymm25" => :SSE25, "zmm25" => :SSE25,
    "xmm26" => :SSE26, "ymm26" => :SSE26, "zmm26" => :SSE26,
    "xmm27" => :SSE27, "ymm27" => :SSE27, "zmm27" => :SSE27,
    "xmm28" => :SSE28, "ymm28" => :SSE28, "zmm28" => :SSE28,
    "xmm29" => :SSE29, "ymm29" => :SSE29, "zmm29" => :SSE29,
    "xmm30" => :SSE30, "ymm30" => :SSE30, "zmm30" => :SSE30,
    "xmm31" => :SSE31, "ymm31" => :SSE31, "zmm31" => :SSE31,
)

############################################################## ASSEMBLY OPERANDS

abstract type AssemblyOperand end

struct AssemblyRegister <: AssemblyOperand
    name::String
end

Base.print(io::IO, reg::AssemblyRegister) = print(io, reg.name)

struct AssemblyMemoryOperand <: AssemblyOperand
    type::String
    address::String
end

Base.print(io::IO, mem::AssemblyMemoryOperand) =
    print(io, "*($(mem.address))")

struct AssemblyImmediate <: AssemblyOperand
    value::String
end

Base.print(io::IO, imm::AssemblyImmediate) = print(io, imm.value)

struct AssemblyOffset <: AssemblyOperand
    name::String
end

Base.print(io::IO, off::AssemblyOffset) = print(io, "offset ", off.name)

############################################################ ASSEMBLY STATEMENTS

abstract type AssemblyStatement end

struct AssemblyComment <: AssemblyStatement
    contents::String
end

Base.print(io::IO, comment::AssemblyComment) = print(io, comment.contents)

struct AssemblyLabel <: AssemblyStatement
    name::String
end

Base.print(io::IO, label::AssemblyLabel) = print(io, label.name, ':')

struct AssemblyInstruction <: AssemblyStatement
    opcode::String
    operands::Vector{AssemblyOperand}
    comment::String
end

AssemblyInstruction(instruction::AbstractString) =
    AssemblyInstruction(String(instruction), AssemblyOperand[], "")

AssemblyInstruction(instruction::AbstractString,
                    operands::Vector{T}) where {T <: AssemblyOperand} =
    AssemblyInstruction(String(instruction), operands, "")

const PRINT_HANDLERS = Dict{String,Function}()

function Base.print(io::IO, instr::AssemblyInstruction)
    if haskey(PRINT_HANDLERS, instr.opcode)
        PRINT_HANDLERS[instr.opcode](io, instr)
    else
        print(io, rpad('{' * instr.opcode * '}', 16))
        for (i, op) in enumerate(instr.operands)
            (i > 1) && print(io, ", ")
            print(io, op)
        end
        if !isempty(instr.comment)
            print(io, " // ", instr.comment)
        end
    end
    return nothing
end

############################################################### ASSEMBLY PARSING

function parse_assembly_operand(op::AbstractString)
    op = strip(op)
    if any(op == reg_name for reg_name in keys(X86_REGISTER_NAMES))
        return AssemblyRegister(op)
    end
    mem_op_match = match(r"([a-z]+) ptr \[(.*)\]", op)
    if !isnothing(mem_op_match)
        return AssemblyMemoryOperand(mem_op_match[1], mem_op_match[2])
    end
    offset_match = match(r"offset (.*)", op)
    if !isnothing(offset_match)
        return AssemblyOffset(offset_match[1])
    end
    # TODO: How to parse immediates?
    return AssemblyImmediate(op)
end

function parse_assembly_statement(stmt::AbstractString)

    # Return nothing for non-statements.
    if (stmt == "\t.text") || isempty(strip(stmt))
        return nothing
    end

    # Assembly comments begin with a semicolon.
    if startswith(strip(stmt), ';')
        return AssemblyComment(stmt)
    end

    # Tabs separate parts of an assembly statement.
    tokens = split(stmt, '\t')

    # If the statement contains no tabs, then it is a label.
    if length(tokens) == 1
        @assert endswith(tokens[1], ':')
        label_name = tokens[1][1:end-1]
        @assert !isempty(label_name)
        return AssemblyLabel(label_name)

    # If the statement begins with a tab, then it is an instruction.
    else
        @assert length(tokens) > 1
        @assert isempty(tokens[1])

        # If there is no second tab, then the instruction takes no operands.
        if length(tokens) == 2
            return AssemblyInstruction(tokens[2])

        # If there is a second tab, then the instruction operands follow it.
        else
            @assert length(tokens) == 3

            # Some instructions are output with a #-delimited comment.
            arg_tokens = split(tokens[3], " # ")
            if length(arg_tokens) == 1
                return AssemblyInstruction(tokens[2],
                    parse_assembly_operand.(split(arg_tokens[1], ',')))
            else
                @assert length(arg_tokens) == 2
                return AssemblyInstruction(tokens[2],
                    parse_assembly_operand.(split(arg_tokens[1], ',')),
                    strip(arg_tokens[2]))
            end
        end
    end

    error("AssemblyView.jl INTERNAL ERROR: Unable to parse" *
          " assembly statement <<<$stmt>>>.")
end

function parsed_asm(@nospecialize(func), @nospecialize(types...);
                    keep_comments::Bool=false)::Vector{AssemblyStatement}

    # Support either a tuple of argument types or varargs.
    if (length(types) == 1) && !(types[1] isa Type)
        types = types[1]
    end

    if (VERSION.major > 1) || (VERSION.major == 1 && VERSION.minor >= 7)

        # Call internal Julia API to generate x86 assembly code.
        code = _dump_function(func, types,
            true,    # Generate native code (as opposed to LLVM IR).
            false,   # Don't generate wrapper code.
            true,    # (strip_ir_metadata) Ignored when dumping native code.
            true,    # (dump_module) Ignored when dumping native code.
            :intel,  # I prefer Intel assembly syntax.
            true,    # (optimize) Ignored when dumping native code.
            :source, # TODO: What does debuginfo=:source mean?
            false    # Don't display binary machine code.
        )

    else

        # Call internal Julia API to generate x86 assembly code.
        code = _dump_function(func, types,
            true,   # Generate native code (as opposed to LLVM IR).
            false,  # Don't generate wrapper code.
            true,   # (strip_ir_metadata) Ignored when dumping native code.
            true,   # (dump_module) Ignored when dumping native code.
            :intel, # I prefer Intel assembly syntax.
            true,   # (optimize) Ignored when dumping native code.
            :source # TODO: What does debuginfo=:source mean?
        )

    end

    # Parse each line of code, discarding comments if requested.
    result = AssemblyStatement[]
    for stmt in split(code, '\n')
        parsed_stmt = parse_assembly_statement(stmt)
        if !isnothing(parsed_stmt) && (keep_comments ||
                                       !(parsed_stmt isa AssemblyComment))
            push!(result, parsed_stmt)
        end
    end

    return result
end

function asm_offsets(
        @nospecialize(func), @nospecialize(types...))::Vector{String}
    result = String[]
    stmts = parsed_asm(func, types...)
    for stmt in stmts
        if stmt isa AssemblyInstruction
            for op in stmt.operands
                if op isa AssemblyOffset
                    push!(result, op.name)
                end
            end
        end
    end
    return result
end

################################################################ PRINT UTILITIES

function assert_num_operands(instr::AssemblyInstruction, n::Int)::Nothing
    if length(instr.operands) != n
        throw(AssertionError("AssemblyView.jl INTERNAL ERROR: Encountered" *
                             " $(instr.opcode) instruction with wrong" *
                             " number of operands (expected $n; found" *
                             " $(length(instr.operands)))."))
    end
end

################################################################ IGNORED OPCODES

const X86_IGNORED_OPCODES = [
    "nop",
    "vzeroupper",
    "ud2",
]

for opcode in X86_IGNORED_OPCODES
    PRINT_HANDLERS[opcode] = (io::IO, instr::AssemblyInstruction) -> print(io)
end

################################################################ CONTROL OPCODES

const X86_CONTROL_OPCODES = [
    "push",
    "pop",
    "call",
    "ret",
    "jmp",
]

function verbatim_print_hander(io::IO, instr::AssemblyInstruction)::Nothing
    assert_num_operands(instr, 1)
    print(io, instr.opcode, ' ', instr.operands[1])
end

PRINT_HANDLERS["push"] =
PRINT_HANDLERS["pop"] =
PRINT_HANDLERS["call"] = verbatim_print_hander

PRINT_HANDLERS["ret"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 0)
    print(io, "return")
end

PRINT_HANDLERS["jmp"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 1)
    label = instr.operands[1]
    @assert label isa Union{AssemblyImmediate,AssemblyMemoryOperand}
    print(io, "goto $label")
end

####################################################### CONDITIONAL JUMP OPCODES

const X86_COMPARISON_OPCODES = [
    "cmp",
    "vucomisd",
]

PRINT_HANDLERS["cmp"] =
PRINT_HANDLERS["vucomisd"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 2)
    a, b = instr.operands
    print(io, "$a <=> $b")
end

const X86_CONDITIONAL_JUMP_OPCODES = [
    "ja", "jnbe", "jae", "jnb", "jb", "jnae", "jbe", "jna",
    "jg", "jnle", "jge", "jnl", "jl", "jnge", "jle", "jng",
    "je", "jne",
]

function make_conditional_jump_handler(op::String)::Function
    return (io::IO, instr::AssemblyInstruction) -> begin
        if length(instr.operands) == 1
            label = instr.operands[1]
            print(io, "if ($op) goto $label")
        elseif length(instr.operands) == 3
            a, b, label = instr.operands
            print(io, "if ($a $op $b) goto $label")
        else
            @assert false
        end
    end
end

PRINT_HANDLERS["ja"] =
PRINT_HANDLERS["jnbe"] =
PRINT_HANDLERS["jg"] =
PRINT_HANDLERS["jnle"] = make_conditional_jump_handler(">")

PRINT_HANDLERS["jae"] =
PRINT_HANDLERS["jnb"] =
PRINT_HANDLERS["jge"] =
PRINT_HANDLERS["jnl"] = make_conditional_jump_handler(">=")

PRINT_HANDLERS["jb"] =
PRINT_HANDLERS["jnae"] =
PRINT_HANDLERS["jl"] =
PRINT_HANDLERS["jnge"] = make_conditional_jump_handler("<")

PRINT_HANDLERS["jbe"] =
PRINT_HANDLERS["jna"] =
PRINT_HANDLERS["jle"] =
PRINT_HANDLERS["jng"] = make_conditional_jump_handler("<=")

PRINT_HANDLERS["je"] = make_conditional_jump_handler("==")

PRINT_HANDLERS["jne"] = make_conditional_jump_handler("!=")

#################################################################### MOV OPCODES

const X86_MOV_OPCODES = [
    "mov",
    "movzx",
    "movabs",
    "vmovss",
    "vmovsd",
    "vmovups",
    "vmovupd",
    "vmovdqu",
    "vmovaps",
    "vmovapd",
    "vmovdqa",
    "vextractf128",
]

PRINT_HANDLERS["mov"] =
PRINT_HANDLERS["movzx"] =
PRINT_HANDLERS["movabs"] =
PRINT_HANDLERS["vmovss"] =
PRINT_HANDLERS["vmovsd"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 2)
    dst, src = instr.operands
    print(io, "$dst = $src")
end

PRINT_HANDLERS["vmovups"] =
PRINT_HANDLERS["vmovupd"] =
PRINT_HANDLERS["vmovdqu"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 2)
    dst, src = instr.operands
    print(io, "$dst .= $src")
end

PRINT_HANDLERS["vmovaps"] =
PRINT_HANDLERS["vmovapd"] =
PRINT_HANDLERS["vmovdqa"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 2)
    dst, src = instr.operands
    print(io, "$dst .= $src [aligned]")
end

PRINT_HANDLERS["vextractf128"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 3)
    dst, src, imm = instr.operands
    @assert imm isa AssemblyImmediate
    if imm.value == "0"
        print(io, "$dst = $src[0:127]")
    elseif imm.value == "1"
        print(io, "$dst = $src[128:255]")
    else
        @assert false
    end
end

############################################################# ARITHMETIC OPCODES

const X86_ARITHMETIC_OPCODES = [
    "inc",
    "dec",
    "add",
    "sub",
    "and",
    "andn",
    "xor",
    "vxorps",
    "vxorpd",
    "shl",
    "sar",
    "lea",
    "vcvtsi2sd",
    "vaddss",
    "vaddsd",
    "vaddps",
    "vaddpd",
    "vsubss",
    "vsubsd",
    "vsubps",
    "vsubpd",
    "vmulss",
    "vmulsd",
    "vmulps",
    "vmulpd",
    "vdivss",
    "vdivsd",
    "vdivps",
    "vdivpd",
    "vpermilpd",
    "vfmadd213sd",
    "vfmadd213pd",
    "vfmadd231sd",
    "vfmadd231pd",
    "vfmsub213sd",
    "vfmsub213pd",
    "vfmsub132sd",
    "vfmsub132pd",
    "vfnmsub213sd",
    "vfnmsub213pd",
]

PRINT_HANDLERS["inc"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 1)
    dst = instr.operands[1]
    print(io, "$dst++")
end

PRINT_HANDLERS["dec"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 1)
    dst = instr.operands[1]
    print(io, "$dst--")
end

PRINT_HANDLERS["add"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 2)
    dst, src = instr.operands
    print(io, "$dst += $src")
end

PRINT_HANDLERS["sub"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 2)
    dst, src = instr.operands
    print(io, "$dst -= $src")
end

PRINT_HANDLERS["and"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 2)
    dst, src = instr.operands
    print(io, "$dst &= $src")
end

PRINT_HANDLERS["andn"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 3)
    dst, a, b = instr.operands
    print(io, "$dst = ~$a & $b")
end

PRINT_HANDLERS["xor"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 2)
    dst, src = instr.operands
    if dst == src
        print(io, "$dst = 0")
    else
        print(io, "$dst ^= src")
    end
end

PRINT_HANDLERS["vxorps"] =
PRINT_HANDLERS["vxorpd"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 3)
    dst, a, b = instr.operands
    if a == b
        print(io, "$dst .= 0")
    else
        print(io, "$dst .= $a .^ $b")
    end
end

PRINT_HANDLERS["shl"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 2)
    dst, src = instr.operands
    print(io, "$dst <<= $src")
end

PRINT_HANDLERS["sar"] = (io::IO, instr::AssemblyInstruction) -> begin
    if length(instr.operands) == 1
        dst = instr.operands[1]
        print(io, "$dst >>= 1")
    elseif length(instr.operands) == 2
        dst, src = instr.operands
        print(io, "$dst >>= $src")
    else
        @assert false
    end
end

PRINT_HANDLERS["lea"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 2)
    dst, src = instr.operands
    @assert dst isa AssemblyRegister
    @assert src isa AssemblyImmediate
    @assert startswith(src.value, '[')
    @assert endswith(src.value, ']')
    print(io, "$dst = $(src.value[2:end-1])")
end

PRINT_HANDLERS["vcvtsi2sd"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 3)
    dst, a, b = instr.operands
    print(io, "$dst = (double) $b[0], $a[1]")
end

PRINT_HANDLERS["vaddss"] =
PRINT_HANDLERS["vaddsd"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 3)
    dst, a, b = instr.operands
    print(io, "$dst = $a + $b")
end

PRINT_HANDLERS["vaddps"] =
PRINT_HANDLERS["vaddpd"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 3)
    dst, a, b = instr.operands
    print(io, "$dst .= $a .+ $b")
end

PRINT_HANDLERS["vsubss"] =
PRINT_HANDLERS["vsubsd"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 3)
    dst, a, b = instr.operands
    print(io, "$dst = $a - $b")
end

PRINT_HANDLERS["vsubps"] =
PRINT_HANDLERS["vsubpd"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 3)
    dst, a, b = instr.operands
    print(io, "$dst .= $a .- $b")
end

PRINT_HANDLERS["vmulss"] =
PRINT_HANDLERS["vmulsd"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 3)
    dst, a, b = instr.operands
    print(io, "$dst = $a * $b")
end

PRINT_HANDLERS["vmulps"] =
PRINT_HANDLERS["vmulpd"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 3)
    dst, a, b = instr.operands
    print(io, "$dst .= $a .* $b")
end

PRINT_HANDLERS["vdivss"] =
PRINT_HANDLERS["vdivsd"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 3)
    dst, a, b = instr.operands
    print(io, "$dst = $a / $b")
end

PRINT_HANDLERS["vdivps"] =
PRINT_HANDLERS["vdivpd"] = (io::IO, instr::AssemblyInstruction) -> begin
    assert_num_operands(instr, 3)
    dst, a, b = instr.operands
    print(io, "$dst .= $a ./ $b")
end

function comment_print_handler(io::IO, instr::AssemblyInstruction)
    @assert !isempty(instr.comment)
    parts = split(instr.comment, " = ")
    @assert length(parts) == 2
    lhs, rhs = parts
    print(io, "$lhs = $rhs")
end

PRINT_HANDLERS["vpermilpd"] =
PRINT_HANDLERS["vfmadd213sd"] =
PRINT_HANDLERS["vfmadd213pd"] =
PRINT_HANDLERS["vfmadd231sd"] =
PRINT_HANDLERS["vfmadd231pd"] =
PRINT_HANDLERS["vfmsub213sd"] =
PRINT_HANDLERS["vfmsub213pd"] =
PRINT_HANDLERS["vfmsub132sd"] =
PRINT_HANDLERS["vfmsub132pd"] =
PRINT_HANDLERS["vfnmsub213sd"] =
PRINT_HANDLERS["vfnmsub213pd"] = comment_print_handler

############################################################### COVERAGE TESTING

@assert isempty(symdiff(
    Set{String}(keys(PRINT_HANDLERS)),
    union(
        Set{String}(X86_IGNORED_OPCODES),
        Set{String}(X86_CONTROL_OPCODES),
        Set{String}(X86_COMPARISON_OPCODES),
        Set{String}(X86_CONDITIONAL_JUMP_OPCODES),
        Set{String}(X86_MOV_OPCODES),
        Set{String}(X86_ARITHMETIC_OPCODES),
    )
))

######################################################### ASSEMBLY PREPROCESSING

is_opcode(stmt::AssemblyStatement, instr::AbstractString)::Bool =
    (stmt isa AssemblyInstruction) && (stmt.opcode == instr)

function remove_nops(
        stmts::Vector{AssemblyStatement})::Vector{AssemblyStatement}
    result = AssemblyStatement[]
    for stmt in stmts
        if all(!is_opcode(stmt, nop) for nop in X86_IGNORED_OPCODES)
            push!(result, stmt)
        end
    end
    return result
end

function remove_prologue_epilogue(
        stmts::Vector{AssemblyStatement})::Vector{AssemblyStatement}

    # Check whether the first two instructions are:
    #     push rbp
    #     mov rbp, rsp
    # Exit early if this is not the case.
    if ((length(stmts) < 2)
        || !is_opcode(stmts[1], "push")
        || (length(stmts[1].operands) != 1)
        || (stmts[1].operands[1] != AssemblyRegister("rbp"))
        || !is_opcode(stmts[2], "mov")
        || (stmts[2].operands != [AssemblyRegister("rbp"),
                                  AssemblyRegister("rsp")]))
        return copy(stmts)
    end

    # At this point, we know we have at least two prologue instructions.
    # The rest of the prologue consists of "push register" instructions...
    prologue_len = 3
    saved_regs = [AssemblyRegister("rbp")]
    while ((prologue_len <= length(stmts))
           && is_opcode(stmts[prologue_len], "push")
           && (length(stmts[prologue_len].operands) == 1)
           && (stmts[prologue_len].operands[1] isa AssemblyRegister))
        push!(saved_regs, stmts[prologue_len].operands[1])
        prologue_len += 1
    end
    reverse!(saved_regs)

    # ...possibly followed by "add/sub rsp, immediate".
    if ((prologue_len <= length(stmts))
        && (is_opcode(stmts[prologue_len], "add")
            || is_opcode(stmts[prologue_len], "sub"))
        && (length(stmts[prologue_len].operands) == 2)
        && (stmts[prologue_len].operands[1] isa AssemblyRegister)
        && (stmts[prologue_len].operands[1].name == "rsp")
        && (stmts[prologue_len].operands[2] isa AssemblyImmediate))
    else
        prologue_len -= 1
    end

    # Delete all prologue instructions.
    deletion_indices = collect(1 : prologue_len)

    # To identify the epilogue, we search for all occurences of "ret".
    ret_indices = [i for i = 1 : length(stmts) if is_opcode(stmts[i], "ret")]

    # We scan backwards from each occurrence, looking for a sequence of "pop"
    # instructions that matches the "push" sequence from the prologue.
    for i in ret_indices
        epilogue = stmts[i-length(saved_regs) : i-1]
        for (stmt, reg) in zip(epilogue, saved_regs)
            if (!is_opcode(stmt, "pop")
                || (length(stmt.operands) != 1)
                || (stmt.operands[1] != reg))
                return copy(stmts)
            end
        end

        # As before, the epilogue might optionally be prefixed with an
        # "add/sub rsp, immediate" instruction. Delete this if it exists.
        j = i - length(saved_regs) - 1
        if ((j > 0)
            && (is_opcode(stmts[j], "add")
                || is_opcode(stmts[j], "sub"))
            && (length(stmts[j].operands) == 2)
            && (stmts[j].operands[1] isa AssemblyRegister)
            && (stmts[j].operands[1].name == "rsp")
            && (stmts[j].operands[2] isa AssemblyImmediate))
            append!(deletion_indices, j : i-1)
        elseif ((j > 0)
                && (is_opcode(stmts[j], "lea"))
                && (length(stmts[j].operands) == 2)
                && (stmts[j].operands[1] isa AssemblyRegister)
                && (stmts[j].operands[1].name == "rsp")
                && (stmts[j].operands[2] isa AssemblyImmediate)
                && startswith(stmts[j].operands[2].value, "[rbp"))
            append!(deletion_indices, j : i-1)
        else
            append!(deletion_indices, j+1 : i-1)
        end
    end

    return deleteat!(copy(stmts), deletion_indices)
end

function fuse_conditional_jumps(
        stmts::Vector{AssemblyStatement})::Vector{AssemblyStatement}
    stmts = copy(stmts)
    for i = 1 : length(stmts)
        if (stmts[i] isa AssemblyInstruction
            && stmts[i].opcode in X86_CONDITIONAL_JUMP_OPCODES
            && stmts[i-1] isa AssemblyInstruction
            && stmts[i-1].opcode in X86_COMPARISON_OPCODES)
            prepend!(stmts[i].operands, stmts[i-1].operands)
            stmts[i-1] = AssemblyInstruction("nop")
        end
    end
    return remove_nops(stmts)
end

is_register_arithmetic_instruction(stmt::AssemblyStatement)::Bool = (
    stmt isa AssemblyInstruction
    && (stmt.opcode in AssemblyView.X86_ARITHMETIC_OPCODES
        || stmt.opcode in AssemblyView.X86_MOV_OPCODES)
    && all(((op isa AssemblyRegister)
            || (op isa AssemblyImmediate))
           for op in stmt.operands))

is_memory_arithmetic_instruction(stmt::AssemblyStatement)::Bool = (
    stmt isa AssemblyInstruction
    && (stmt.opcode in AssemblyView.X86_ARITHMETIC_OPCODES
        || stmt.opcode in AssemblyView.X86_MOV_OPCODES)
    && all(((op isa AssemblyRegister)
            || (op isa AssemblyImmediate)
            || (op isa AssemblyMemoryOperand))
           for op in stmt.operands))

function fuse_instructions(predicate, pseudo_opcode::String, min_length::Int,
        stmts::Vector{AssemblyStatement})::Vector{AssemblyStatement}
    stmts = copy(stmts)
    i = 1
    while i <= length(stmts)
        if predicate(stmts[i])
            j = i
            affected_operands = Set{AssemblyOperand}()
            while predicate(stmts[j])
                for op in stmts[j].operands
                    if !(op isa AssemblyImmediate)
                        push!(affected_operands, op)
                    end
                end
                j += 1
            end
            j -= 1
            if (j - i + 1) >= min_length
                stmts[i] = AssemblyInstruction(pseudo_opcode,
                    push!(collect(affected_operands),
                        AssemblyImmediate(string(length(affected_operands)))))
                for k = i+1 : j
                    stmts[k] = AssemblyInstruction("nop")
                end
            end
            i = j + 1
        else
            i += 1
        end
    end
    return remove_nops(stmts)
end

################################################################ ASSEMBLY OUTPUT

const INSTRUCTION_PREFIX = "\t"
const INSTRUCTION_SUFFIX = ";"

function view_asm(io::IO, @nospecialize(func), @nospecialize(types...))::Nothing

    parsed_stmts = parsed_asm(func, types...)
    parsed_stmts = remove_nops(parsed_stmts)
    parsed_stmts = remove_prologue_epilogue(parsed_stmts)
    parsed_stmts = fuse_conditional_jumps(parsed_stmts)

    unknown_opcodes = Dict{String,Int}()
    for stmt in parsed_stmts
        if stmt isa AssemblyInstruction
            println(io, INSTRUCTION_PREFIX, stmt, INSTRUCTION_SUFFIX)
            if !haskey(PRINT_HANDLERS, stmt.opcode)
                if haskey(unknown_opcodes, stmt.opcode)
                    unknown_opcodes[stmt.opcode] += 1
                else
                    unknown_opcodes[stmt.opcode] = 1
                end
            end
        else
            println(io, stmt)
        end
    end

    frequency_table = [(count, opcode) for (opcode, count) in unknown_opcodes]
    sort!(frequency_table, rev=true)

    if !isempty(frequency_table)
        println(io)
        println(io, "Unknown opcodes:")
        for (count, opcode) in frequency_table
            println(io, INSTRUCTION_PREFIX, "$opcode ($count occurrences)")
        end
    end

    println(io)
    return nothing
end

view_asm(@nospecialize(func), @nospecialize(types...))::Nothing =
    view_asm(stdout, func, types...)

end # module AssemblyView


text:
#
# Copyright (c) 2021 Tobias Thummerer, Lars Mikelsons, Josef Kircher
# Licensed under the MIT license. See LICENSE file in the project root for details.
#

using OrdinaryDiffEq: ODESolution
using FMIImport: FMU2Solution

"""
Plots data from a ME-FMU.

Optional `t_in_solution` controls if the first state in the solution is interpreted as t(ime).
Optional keyword argument `maxLabelLength` controls the maximum length for legend labels (too long labels are cut from front).
"""
function fmiPlot(solution::FMU2Solution; kwargs...)
    fig = Plots.plot(; xlabel="t [s]")
    fmiPlot!(fig, solution; kwargs...)
    return fig
end
function fmiPlot!(fig, solution::FMU2Solution;
    states::Union{Bool, Nothing}=nothing,
    values::Union{Bool, Nothing}=nothing,
    stateEvents::Union{Bool, Nothing}=nothing,
    timeEvents::Union{Bool, Nothing}=nothing,
    stateIndices=nothing,
    valueIndices=nothing,
    maxLabelLength=64,
    plotkwargs...)

    numStateEvents = 0
    numTimeEvents = 0
    for e in solution.events
        if e.indicator > 0
            numStateEvents += 1
        else
            numTimeEvents += 1
        end
    end

    if states === nothing
        states = (solution.states !== nothing)
    end

    if values === nothing
        values = (solution.values !== nothing)
    end

    if stateEvents === nothing
        stateEvents = false
        for e in solution.events
            if e.indicator > 0
                stateEvents = true
                break
            end
        end

        if numStateEvents > 100
            @info "fmiPlot(...): Number of state events ($(numStateEvents)) exceeding 100, disabling automatic plotting of state events (can be forced with keyword `stateEvents=true`)."
            stateEvents = false
        end
    end

    if timeEvents === nothing
        timeEvents = false
        for e in solution.events
            if e.indicator == 0
                timeEvents = true
                break
            end
        end

        if numTimeEvents > 100
            @info "fmiPlot(...): Number of time events ($(numTimeEvents)) exceeding 100, disabling automatic plotting of time events (can be forced with keyword `timeEvents=true`)."
            timeEvents = false
        end
    end

    if stateIndices === nothing
        stateIndices = 1:length(solution.fmu.modelDescription.stateValueReferences)
    end

    if valueIndices === nothing
        if solution.values !== nothing
            valueIndices = 1:length(solution.values.saveval[1])
        end
    end

    plot_min = Inf
    plot_max = -Inf

    # plot states
    if states
        t = solution.states.t
        numValues = length(solution.states.u[1])

        for v in 1:numValues
            if v ∈ stateIndices
                vr = solution.fmu.modelDescription.stateValueReferences[v]
                vrNames = fmi2ValueReferenceToString(solution.fmu, vr)
                vrName = vrNames[1]

                vals = collect(data[v] for data in solution.states.u)

                plot_min = min(plot_min, vals...)
                plot_max = max(plot_max, vals...)

                # prevent legend labels from getting too long
                label = "$vrName ($vr)"
                labelLength = length(label)
                if labelLength > maxLabelLength
                    label = "..." * label[labelLength-maxLabelLength:end]
                end

                Plots.plot!(fig, t, vals; label=label, plotkwargs...)
            end
        end
    end

    # plot recorded values
    if values
        t = solution.values.t
        numValues = length(solution.values.saveval[1])

        for v in 1:numValues
            if v ∈ valueIndices
                vr = solution.valueReferences[v]
                vrNames = fmi2ValueReferenceToString(solution.fmu, vr)
                vrName = vrNames[1]

                vals = collect(data[v] for data in solution.values.saveval)

                plot_min = min(plot_min, vals...)
                plot_max = max(plot_max, vals...)

                # prevent legend labels from getting too long
                label = "$vrName ($vr)"
                labelLength = length(label)
                if labelLength > maxLabelLength
                    label = "..." * label[labelLength-maxLabelLength:end]
                end

                Plots.plot!(fig, t, vals; label=label, plotkwargs...)
            end
        end
    end

    if stateEvents
        first = true
        for e in solution.events
            if e.indicator > 0
                Plots.plot!(fig, [e.t, e.t], [plot_min, plot_max]; label=(first ? "State event(s)" : nothing), style=:dash, color=:blue)
                first = false
            end
        end
    end

    if timeEvents
        first = true
        for e in solution.events
            if e.indicator == 0
                Plots.plot!(fig, [e.t, e.t], [plot_min, plot_max]; label=(first ? "Time event(s)" : nothing), style=:dash, color=:red)
                first = false
            end
        end
    end

    return fig
end

"""
Extended the original plot-command by plotting FMUs.
"""
function Plots.plot(fmu::FMU2, args...; kwargs...)
    fmiPlot(fmu, args...; kwargs...)
end
function Plots.plot!(fig, fmu::FMU2, args...; kwargs...)
    fmiPlot!(fig, fmu, args...; kwargs...)
end


text:
@testset "Events/Attributes" begin
    evaluate(s) = collect(E.events(L.State(IOBuffer(s))))

    E = Events
    L = Events.Lexical

    AS = E.AttributeSpecification
    ES = E.ElementStart
    DC = E.DataContent
    ME = E.MarkupError

    @testset "Events/Attributes (Positive)" begin
        @test (evaluate("") == [ ES("a", [ AS("b", [ DC("c", L.Location("a buffer", -1)) ],
                                                         L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("abc", [ DC("def", L.Location("a buffer", -1)) ],
                                                             L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("é", [ DC("è", L.Location("a buffer", -1)) ],
                                                         L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("aé", [ DC("aè", L.Location("a buffer", -1)) ],
                                                           L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("éb", [ DC("èb", L.Location("a buffer", -1)) ],
                                                           L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("aéb", [ DC("aèb", L.Location("a buffer", -1)) ],
                                                             L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])

        @test (evaluate("") == [ ES("a", [ AS("b", [ DC("c", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ],
                                             L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("abc", [ DC("def", L.Location("a buffer", -1)) ],
                                                           L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("é", [ DC("è", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ],
                                             L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("aé", [ DC("aè", L.Location("a buffer", -1)) ],
                                                         L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("éb", [ DC("èb", L.Location("a buffer", -1)) ],
                                                         L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("aéb", [ DC("aèb", L.Location("a buffer", -1)) ],
                                                           L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])

        @test (evaluate("") == [ ES("a", [ AS("b", [ DC("c", L.Location("a buffer", -1)) ],
                                                           L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("abc", [ DC("def", L.Location("a buffer", -1)) ],
                                                               L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("é", [ DC("è", L.Location("a buffer", -1)) ],
                                                           L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("aé", [ DC("aè", L.Location("a buffer", -1)) ],
                                                             L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("éb", [ DC("èb", L.Location("a buffer", -1)) ],
                                                             L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("") == [ ES("a", [ AS("aéb", [ DC("aèb", L.Location("a buffer", -1)) ],
                                                               L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])

        @test (evaluate("")
               == [ ES("a", [ AS("b", [ DC("c", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)),
                              AS("d", [ DC("e", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ],
                       L.Location("a buffer", -1)) ])

        @test (evaluate("")
               == [ ES("a", [ AS("b", [ DC("c", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)),
                              AS("d", [ DC("e", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)),
                              AS("f", [ DC("g", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ],
                       L.Location("a buffer", -1)) ])

        # Make sure any white space can be used as a separator.
        #
        @test (evaluate("")
               == [ ES("a", [ AS("b", [ DC("c", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)),
                              AS("d", [ DC("e", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)),
                              AS("f", [ DC("g", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ],
                       L.Location("a buffer", -1)) ])

        # Make sure trailing white space is consumed.
        #
        @test (evaluate("")
               == [ ES("a", [ AS("b", [ DC("c", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ],
                       L.Location("a buffer", -1)) ])
        @test (evaluate("")
               == [ ES("a", [ AS("b", [ DC("c", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ],
                       L.Location("a buffer", -1)) ])
        @test (evaluate("")
               == [ ES("a", [ AS("b", [ DC("c", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ],
                       L.Location("a buffer", -1)) ])

        # Check that attribute values can contain allowed special characters. Well, try a few representative cases.
        #
        @test (evaluate("")
               == [ ES("a", [ AS("b", [ DC(" ", true, L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ],
                       L.Location("a buffer", -1)) ])
        @test (evaluate("")
               == [ ES("a", [ AS("b", [ E.EntityReferenceGeneral("amp", L.Location("a buffer", -1)) ],
                                 L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ])
        @test (evaluate("")
               == [ ES("a", [ AS("b", [ E.CharacterReference("20", L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ],
                       L.Location("a buffer", -1)) ])

        # Make sure that the "other" quotes can be used inside the attribute value.
        #
        @test (evaluate("")
               == [ ES("a", [ AS("b", [ DC("'", false, L.Location("a buffer", -1)),
                                        DC("Hello", false, L.Location("a buffer", -1)),
                                        DC("'", false, L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ],
                       L.Location("a buffer", -1)) ])
        @test (evaluate("")
               == [ ES("a", [ AS("b", [ DC("\"", false, L.Location("a buffer", -1)),
                                        DC("Hello", false, L.Location("a buffer", -1)),
                                        DC("\"", false, L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ],
                       L.Location("a buffer", -1)) ])
    end

    @testset "Events/Attributes (Negative ... no VI)" begin
        @test (evaluate("' to end an element open tag.", L.Location("a buffer", -1)) ])

        events = evaluate("")
        @test length(events) == 7
        @test events == [ ME("ERROR: Expecting '=' after an attribute name.", L.Location("a buffer", -1)),
                          ES(true, false, "a", [ AS("b", [ ], L.Location("a buffer", -1)) ], L.Location("a buffer", -1)),
                          ME("ERROR: Expecting '>' to end an element open tag.", L.Location("a buffer", -1)),
                          DC("\"", false, L.Location("a buffer", -1)),
                          DC("c", false, L.Location("a buffer", -1)),
                          DC("\"", false, L.Location("a buffer", -1)),
                          DC(">", false, L.Location("a buffer", -1)) ]

        events = evaluate("")
        @test length(events) == 7
        @test events == [ ME("ERROR: Expecting '=' after an attribute name.", L.Location("a buffer", -1)),
                          ES(true, false, "a",
                             [ AS("b", [ DC("c", false, L.Location("a buffer", -1)) ], L.Location("a buffer", -1)),
                               AS("d", [ ], L.Location("a buffer", -1)) ], L.Location("a buffer", -1)),
                          ME("ERROR: Expecting '>' to end an element open tag.", L.Location("a buffer", -1)),
                          DC("\"", false, L.Location("a buffer", -1)),
                          DC("e", false, L.Location("a buffer", -1)),
                          DC("\"", false, L.Location("a buffer", -1)),
                          DC(">", false, L.Location("a buffer", -1)) ]
    end

    @testset "Events/Attributes (Negative ... no quoted value)" begin
        events = evaluate("' to end an element open tag.", L.Location("a buffer", -1)) ]

        events = evaluate("' to end an element open tag.", L.Location("a buffer", -1)) ]

        events = evaluate("' to end an element open tag.", L.Location("a buffer", -1)),
                          DC("c", false, L.Location("a buffer", -1)) ]

        events = evaluate("' to end an element open tag.", L.Location("a buffer", -1)),
                          DC("c", false, L.Location("a buffer", -1)) ]
    end

    @testset "Events/Attributes (Negative ... incomplete quoted value)" begin
        events = evaluate("' to end an element open tag.", L.Location("a buffer", -1)) ]

        events = evaluate("' to end an element open tag.", L.Location("a buffer", -1)) ]
    end

    @testset "Events/Attributes (Negative ... unescaped special characters)" begin
        events = evaluate("<a b=\"<\">")
        @test length(events) == 2
        @test events == [ ME("ERROR: A '<' must be escaped inside an attribute value.", L.Location("a buffer", -1)),
                          ES(false, false, "a", [ AS("b", [ ], L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ]

        events = evaluate("<a b=\"<<<\">")
        @test length(events) == 4
        @test events == [ ME("ERROR: A '<' must be escaped inside an attribute value.", L.Location("a buffer", -1)),
                          ME("ERROR: A '<' must be escaped inside an attribute value.", L.Location("a buffer", -1)),
                          ME("ERROR: A '<' must be escaped inside an attribute value.", L.Location("a buffer", -1)),
                          ES(false, false, "a", [ AS("b", [ ], L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ]

        events = evaluate("<a b=\"Hello, <World!\">")
        @test length(events) == 2
        @test events == [ ME("ERROR: A '<' must be escaped inside an attribute value.", L.Location("a buffer", -1)),
                          ES(false, false, "a", [ AS("b", [ DC("Hello", false, L.Location("a buffer", -1)),
                                                            DC(",", false, L.Location("a buffer", -1)),
                                                            DC(" ", true, L.Location("a buffer", -1)),
                                                            DC("World!", false, L.Location("a buffer", -1)) ],
                                                     L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ]

        events = evaluate("")
        @test length(events) == 2
        @test events == [ ME("ERROR: A '&' must be escaped inside an attribute value.", L.Location("a buffer", -1)),
                          ES(false, false, "a", [ AS("b", [ DC("&", false, L.Location("a buffer", -1)) ],
                                                     L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ]

        events = evaluate("")
        @test length(events) == 4
        @test events == [ ME("ERROR: A '&' must be escaped inside an attribute value.", L.Location("a buffer", -1)),
                          ME("ERROR: A '&' must be escaped inside an attribute value.", L.Location("a buffer", -1)),
                          ME("ERROR: A '&' must be escaped inside an attribute value.", L.Location("a buffer", -1)),
                          ES(false, false, "a", [ AS("b", [ DC("&", false, L.Location("a buffer", -1)),
                                                            DC("&", false, L.Location("a buffer", -1)),
                                                            DC("&", false, L.Location("a buffer", -1)) ],
                                                     L.Location("a buffer", -1)) ],
                             L.Location("a buffer", -1)) ]

        events = evaluate("")
        @test length(events) == 2
        @test events == [ ME("ERROR: A '&' must be escaped inside an attribute value.", L.Location("a buffer", -1)),
                          ES(false, false, "a", [ AS("b", [ DC("Jack", false, L.Location("a buffer", -1)),
                                                            DC(" ", true, L.Location("a buffer", -1)),
                                                            DC("&", false, L.Location("a buffer", -1)),
                                                            DC(" ", true, L.Location("a buffer", -1)),
                                                            DC("Jill", false, L.Location("a buffer", -1)) ],
                                                     L.Location("a buffer", -1)) ], L.Location("a buffer", -1)) ]
    end
end

max_stars_repo_name: jlegare/NativeXML

max_issues_repo_name: jlegare/NativeXML

max_forks_repo_name: jlegare/NativeXML


text:
module NaivePDE
using NaiveBoundary

# package code goes here
include("oneDopeartor.jl")
include("TimeIntergral/TimeIntergral.jl")
include("hypobolicPDE/schemes.jl")

export eval

end # module


text:
using Oceananigans.Models: ShallowWaterModel
using Oceananigans.Grids: Periodic, Bounded

function time_stepping_shallow_water_model_works(arch, topo, coriolis, advection; timestepper=:RungeKutta3)
    grid = RegularCartesianGrid(size=(1, 1, 1), extent=(2π, 2π, 2π), topology=topo)
    model = ShallowWaterModel(grid=grid, gravitational_acceleration=1, architecture=arch, coriolis=coriolis,
                              advection=advection, timestepper=:RungeKutta3)
    set!(model, h=1)

    simulation = Simulation(model, Δt=1.0, stop_iteration=1)
    run!(simulation)

    return model.clock.iteration == 1
end

function shallow_water_model_tracers_and_forcings_work(arch)
    grid = RegularCartesianGrid(size=(1, 1, 1), extent=(2π, 2π, 2π))
    model = ShallowWaterModel(grid=grid, gravitational_acceleration=1, architecture=arch, tracers=(:c, :d))
    set!(model, h=1)

    @test model.tracers.c isa Field
    @test model.tracers.d isa Field

    @test haskey(model.forcing, :uh)
    @test haskey(model.forcing, :vh)
    @test haskey(model.forcing, :h)
    @test haskey(model.forcing, :c)
    @test haskey(model.forcing, :d)

    simulation = Simulation(model, Δt=1.0, stop_iteration=1)
    run!(simulation)

    @test model.clock.iteration == 1

    return nothing
end

@testset "Shallow Water Models" begin
    @info "Testing shallow water models..."

    @testset "Model constructor errors" begin
        grid = RegularCartesianGrid(size=(1, 1, 1), extent=(1, 1, 1))
        @test_throws TypeError ShallowWaterModel(architecture=CPU, grid=grid, gravitational_acceleration=1)
        @test_throws TypeError ShallowWaterModel(architecture=GPU, grid=grid, gravitational_acceleration=1)
    end

    topos = (
             (Periodic, Periodic,  Bounded),
             (Periodic,  Bounded,  Bounded),
             (Bounded,   Bounded,  Bounded),
            )

    for topo in topos
        @testset "$topo model construction" begin
            @info "  Testing $topo model construction..."
            for arch in archs, FT in float_types
		        arch isa GPU && topo == (Bounded, Bounded, Bounded) && continue

                grid = RegularCartesianGrid(FT, topology=topo, size=(1, 1, 1), extent=(1, 2, 3))
                model = ShallowWaterModel(grid=grid, gravitational_acceleration=1, architecture=arch)

                # Just testing that the model was constructed with no errors/crashes.
                @test model isa ShallowWaterModel

                # Test that the grid didn't get mangled (sort of)
                @test size(grid) === size(model.grid)

                too_big_grid = RegularCartesianGrid(FT, topology=topo, size=(1, 1, 2), extent=(1, 2, 3))

                @test_throws ArgumentError ShallowWaterModel(grid=too_big_grid, gravitational_acceleration=1, architecture=arch)
            end
        end
    end

    @testset "Setting ShallowWaterModel fields" begin
        @info "  Testing setting shallow water model fields..."
        for arch in archs, FT in float_types
            N = (4, 4, 1)
            L = (2π, 3π, 5π)

            grid = RegularCartesianGrid(FT, size=N, extent=L)
            model = ShallowWaterModel(grid=grid, gravitational_acceleration=1, architecture=arch)

            x, y, z = nodes((Face, Cell, Cell), model.grid, reshape=true)

            uh₀(x, y, z) = x * y^2
            uh_answer = @. x * y^2

            h₀ = rand(size(grid)...)
            h_answer = deepcopy(h₀)

            set!(model, uh=uh₀, h=h₀)

            uh, vh, h = model.solution

            @test all(interior(uh) .≈ uh_answer)
            @test all(interior(h) .≈ h_answer)
        end
    end

    for arch in archs
        for topo in topos
            @testset "Time-stepping ShallowWaterModels [$arch, $topo]" begin
                @info "  Testing time-stepping ShallowWaterModels [$arch, $topo]..."
                @test time_stepping_shallow_water_model_works(arch, topo, nothing, nothing)
            end
        end

        for coriolis in (nothing, FPlane(f=1), BetaPlane(f₀=1, β=0.1))
            @testset "Time-stepping ShallowWaterModels [$arch, $(typeof(coriolis))]" begin
                @info "  Testing time-stepping ShallowWaterModels [$arch, $(typeof(coriolis))]..."
                @test time_stepping_shallow_water_model_works(arch, topos[1], coriolis, nothing)
            end
        end

        for advection in (nothing, CenteredSecondOrder(), WENO5())
            @testset "Time-stepping ShallowWaterModels [$arch, $(typeof(advection))]" begin
                @info "  Testing time-stepping ShallowWaterModels [$arch, $(typeof(advection))]..."
                @test time_stepping_shallow_water_model_works(arch, topos[1], nothing, advection)
            end
        end

        for timestepper in (:RungeKutta3, :QuasiAdamsBashforth2)
            @testset "Time-stepping ShallowWaterModels [$arch, $timestepper]" begin
                @info "  Testing time-stepping ShallowWaterModels [$arch, $timestepper]..."
                @test time_stepping_shallow_water_model_works(arch, topos[1], nothing, nothing, timestepper=timestepper)
            end
        end

        @testset "ShallowWaterModel with tracers and forcings [$arch]" begin
            @info "  Testing ShallowWaterModel with tracers and forcings [$arch]..."
            shallow_water_model_tracers_and_forcings_work(arch)
        end
    end
end


text:
print("Hello World!")

function greet(name)
  print("Greetings, " * name)
end

greet("TechnoDot")
greet("himybro")
greet("HawkEye")
greet("Test_User")

max_stars_repo_name: TechDudie/julia

max_issues_repo_name: TechDudie/julia

max_forks_repo_name: TechDudie/julia


text:
module kernels
  using CUDAnative, CuArrays
  export sync, dispatch
  macro sync(type, expr)
    ex = quote
      if $type == CuArray
        CuArrays.@sync begin
        $expr
        end
      end
      if $type == Array
        $expr
      end
    end
    return esc(ex)
  end

  macro dispatch(type, threads, blocks, expr)
    cuda = Meta.parse("kernels.cuda_$expr")
    cpu = Meta.parse("kernels.cpu_$expr")
    ex = nothing
    ex = quote
      if $type == CuArray
          @cuda $threads $blocks $cuda
      end
      if $type == Array
          $cpu
      end
    end
    return esc(ex)
  end
  function cuda_term1(viewToR, cuVm, colptrVm, nzvalVm, colptrYftR, nzvalYftR,
                              cuVa, colptrVa, nzvalVa, colptrYftI, nzvalYftI,
                              sizeFromLines)
      index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
      stride = blockDim().x * gridDim().x
      for b in index:stride:size(viewToR,1)
        for (i,c) in enumerate(colptrVm[b]:colptrVm[b+1]-1)
          @inbounds viewToR[b] += cuVm[b] * nzvalVm[c] *
                                  (  nzvalYftR[c] * CUDAnative.cos(cuVa[b] - nzvalVa[c])
                                   + nzvalYftI[c] * CUDAnative.sin(cuVa[b] - nzvalVa[c])
                                  )
        end
      end
      return nothing
  end

  function cuda_term2(viewFromR, cuVm, colptrVm, nzvalVm, colptrYtfR, nzvalYtfR,
                                  cuVa, colptrVa, nzvalVa, colptrYtfI, nzvalYtfI,
                                  sizeToLines)
      index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
      stride = blockDim().x * gridDim().x
      for b in index:stride:size(viewFromR,1)
        for (i,c) in enumerate(colptrVm[b]:colptrVm[b+1]-1)
          @inbounds viewFromR[b] += cuVm[b] * nzvalVm[c] *
                                    (  nzvalYtfR[c] * CUDAnative.cos(cuVa[b] - nzvalVa[c])
                                     + nzvalYtfI[c] * CUDAnative.sin(cuVa[b] - nzvalVa[c])
                                    )
        end
      end
      return nothing
  end

  function cuda_term3(viewToI, cuVm, colptrVm, nzvalVm, colptrYftI, nzvalYftI,
                              cuVa, colptrVa, nzvalVa, colptrYftR, nzvalYftR,
                              sizeFromLines)
      index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
      stride = blockDim().x * gridDim().x
      for b in index:stride:size(viewToI,1)
        for (i,c) in enumerate(colptrVm[b]:colptrVm[b+1]-1)
          @inbounds viewToI[b] += cuVm[b] * nzvalVm[c] *
                                  (  ( -nzvalYftI[c]) * CUDAnative.cos(cuVa[b] - nzvalVa[c])
                                      + nzvalYftR[c] * CUDAnative.sin(cuVa[b] - nzvalVa[c])
                                  )
        end
      end
      return nothing
  end

  function cuda_term4(viewFromI, cuVm, colptrVm, nzvalVm, colptrYtfI, nzvalYtfI,
                                  cuVa, colptrVa, nzvalVa, colptrYtfR, nzvalYtfR,
                                  sizeToLines)
      index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
      stride = blockDim().x * gridDim().x
      for b in index:stride:size(viewFromI,1)
        for (i,c) in enumerate(colptrVm[b]:colptrVm[b+1]-1)
          @inbounds viewFromI[b] += cuVm[b] * nzvalVm[c] *
                                    (  ( -nzvalYtfI[c]) * CUDAnative.cos(cuVa[b] - nzvalVa[c])
                                        + nzvalYtfR[c] * CUDAnative.sin(cuVa[b] - nzvalVa[c])
                                    )
        end
      end
      return nothing
  end
  function cuda_sumPg(sumPg, colptr, nzval)
      index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
      stride = blockDim().x * gridDim().x
      for b in index:stride:size(sumPg,1)
        for c in colptr[b]:colptr[b+1]-1
            @inbounds sumPg[b] += nzval[c]
        end
      end
      return nothing
  end
  function cuda_sumQg(sumQg, colptr, nzval)
      index = (blockIdx().x - 1) * blockDim().x + threadIdx().x
      stride = blockDim().x * gridDim().x
      for b in index:stride:size(sumQg,1)
        for c in colptr[b]:colptr[b+1]-1
            @inbounds sumQg[b] += nzval[c]
        end
      end
      return nothing
  end
  function cpu_term1(viewToR, cuVm, colptrVm, nzvalVm, colptrYftR, nzvalYftR,
                              cuVa, colptrVa, nzvalVa, colptrYftI, nzvalYftI,
                              sizeFromLines)
      Threads.@threads for b in 1:size(viewToR,1)
        for (i,c) in enumerate(colptrVm[b]:colptrVm[b+1]-1)
          viewToR[b] += cuVm[b] * nzvalVm[c] *
                                  (  nzvalYftR[c] * cos(cuVa[b] - nzvalVa[c])
                                   + nzvalYftI[c] * sin(cuVa[b] - nzvalVa[c])
                                  )
        end
      end
      return nothing
  end

  function cpu_term2(viewFromR, cuVm, colptrVm, nzvalVm, colptrYtfR, nzvalYtfR,
                                  cuVa, colptrVa, nzvalVa, colptrYtfI, nzvalYtfI,
                                  sizeToLines)
      for b in 1:size(viewFromR,1)
        for (i,c) in enumerate(colptrVm[b]:colptrVm[b+1]-1)
          viewFromR[b] += cuVm[b] * nzvalVm[c] *
                                    (  nzvalYtfR[c] * cos(cuVa[b] - nzvalVa[c])
                                     + nzvalYtfI[c] * sin(cuVa[b] - nzvalVa[c])
                                    )
        end
      end
      return nothing
  end

  function cpu_term3(viewToI, cuVm, colptrVm, nzvalVm, colptrYftI, nzvalYftI,
                              cuVa, colptrVa, nzvalVa, colptrYftR, nzvalYftR,
                              sizeFromLines)
      for b in 1:size(viewToI,1)
        for (i,c) in enumerate(colptrVm[b]:colptrVm[b+1]-1)
          viewToI[b] += cuVm[b] * nzvalVm[c] *
                                  (  ( -nzvalYftI[c]) * cos(cuVa[b] - nzvalVa[c])
                                      + nzvalYftR[c] * sin(cuVa[b] - nzvalVa[c])
                                  )
        end
      end
      return nothing
  end

  function cpu_term4(viewFromI, cuVm, colptrVm, nzvalVm, colptrYtfI, nzvalYtfI,
                                  cuVa, colptrVa, nzvalVa, colptrYtfR, nzvalYtfR,
                                  sizeToLines)
      for b in 1:size(viewFromI,1)
        for (i,c) in enumerate(colptrVm[b]:colptrVm[b+1]-1)
          viewFromI[b] += cuVm[b] * nzvalVm[c] *
                                    (  ( -nzvalYtfI[c]) * cos(cuVa[b] - nzvalVa[c])
                                        + nzvalYtfR[c] * sin(cuVa[b] - nzvalVa[c])
                                    )
        end
      end
      return nothing
  end
  function cpu_sumPg(sumPg, colptr, nzval)
      for b in 1:size(sumPg,1)
        for c in colptr[b]:colptr[b+1]-1
            @inbounds sumPg[b] += nzval[c]
        end
      end
      return nothing
  end
  function cpu_sumQg(sumQg, colptr, nzval)
      for b in 1:size(sumQg,1)
        for c in colptr[b]:colptr[b+1]-1
            sumQg[b] += nzval[c]
        end
      end
      return nothing
  end
end


text:
tmpOriginData = Dict(
    # ---------- Constants
    :Smax => 10,  # maximum age
    :Sret => 1,  # retirement age
    :alpha => 1.5, # leisure preference than consumption
    :gamma => 0.5, # the inter-temporal elasticity of substitution
    :k1 => 0.0, # capital when born
    # ---------- Vectors
    :Survival => fill(0.001, 10),  # survival probabilities between two years
    :q => fill(0.15, 10),  # the ratio of health expenditure on consumption
    :r => fill(0.08, 10),  # interest rate
    :w => fill(3.25, 1),  # wage level
    :z => fill(0.85, 1),  # the collection rate of PAYG pension
    :θ => fill(0.08, 1),  # contribution: agent -> PAYG pension
    :η => fill(0.20, 1),  # contribution: firm  -> PAYG pension
    :ϕ => fill(0.02, 1),  # contribution: agent -> UEBMI
    :ζ => fill(0.06, 1),  # contribution: firm  -> UEBMI
    :cpB => fill(0.30, 10),  # copayment rate of UEBMI (inhospital)
    :p => fill(1.10, 10),  # the ratio of outpatient expenditure on inpatient expenditure
    :Λ => fill(0.95, 10-1),  # the benefits of PAYG pension
    :𝕡 => fill(0.10, 10-1),  # the amount of the transfer payment from this year's firm contribution to UEBMI to those have retired in this year
    :𝕒 => fill(0.30, 1),  # the rate of the money transferred from this year's firm contribution to those working men's individual account of UEBMI
    # ------------ Constant in this paper but converted to vectors
    :σ => fill(0.24, 1),  # wage taxation
    :μ => fill(0.10, 10),  # consumption taxation
    :δ => fill(1/0.99 - 1, 10), # the discounting rate of utility
)
# 正常working+retire
# tmpOriginData = House.SampleOrigindata;
tmpConst, tmpDict = House.lev0Abbr( tmpOriginData )
# tmpConst = House.SampleConst; tmpDict = copy(House.SampleLev0Abbr);
House.lev1Abbr!(tmpDict,tmpConst)
House.lev2Abbr!(tmpDict,tmpConst)
tmpcs, tmpls = House.getcls( 0.1, tmpDict, tmpConst )
tmpks = House.getks( tmpcs, tmpls, tmpDict, tmpConst )
House.G( 0.1, tmpDict, tmpConst )
@time tmpRes = House.HHSolve( tmpOriginData , ReturnData = true )
House.ExtractAPhi!( tmpRes, tmpOriginData , a1 = 0.0 )

# 只有1期working
tmpOriginData1 = House.SampleOrigindata1;
tmpConst1, tmpDict1 = House.lev0Abbr( tmpOriginData1 )
# tmpConst1 = House.SampleConst1; tmpDict1 = copy(House.SampleLev0Abbr1);
House.lev1Abbr!(tmpDict1,tmpConst1)
House.lev2Abbr!(tmpDict1,tmpConst1)
tmpcs1, tmpls1 = House.getcls( 0.1, tmpDict1, tmpConst1 )
tmpks1 = House.getks( tmpcs1, tmpls1, tmpDict1, tmpConst1 )
House.G( 0.1, tmpDict1, tmpConst1 )
@time tmpRes1 = House.HHSolve( tmpOriginData1 , ReturnData = true )
House.ExtractAPhi!( tmpRes1, tmpOriginData1 , a1 = 0.0 )

# 只有退休期
tmpOriginData_Retired = House.SampleOrigindata_Retired;
tmpConst_Retired, tmpDict_Retired = House.lev0Abbr_Retired( tmpOriginData_Retired )
# tmpConst_Retired = House.SampleConst_Retired; tmpDict_Retired = copy(House.SampleLev0Abbr_Retired);
House.lev1Abbr_Retired!(tmpDict_Retired,tmpConst_Retired)
House.lev2Abbr_Retired!(tmpDict_Retired,tmpConst_Retired)
tmpcs_Retired = House.getcls_Retired( 0.1, tmpDict_Retired, tmpConst_Retired )
tmpks_Retired = House.getks_Retired( tmpcs_Retired, tmpDict_Retired, tmpConst_Retired )
House.G_Retired( 0.1, tmpDict_Retired, tmpConst_Retired )
@time tmpRes_Retired = House.HHSolve_Retired( tmpOriginData_Retired , ReturnData = true )
House.ExtractAPhi_Retired!( tmpRes_Retired, tmpOriginData_Retired , a1 = 0.0 )

# 只有1期退休期
tmpOriginData_Retired1 = House.SampleOrigindata_Retired1;
tmpConst_Retired1, tmpDict_Retired1 = House.lev0Abbr_Retired( tmpOriginData_Retired1 )
# tmpConst_Retired1 = House.SampleConst_Retired1; tmpDict_Retired1 = copy(House.SampleLev0Abbr_Retired1);
House.lev1Abbr_Retired!(tmpDict_Retired1,tmpConst_Retired1)
House.lev2Abbr_Retired!(tmpDict_Retired1,tmpConst_Retired1)
tmpcs_Retired1 = House.getcls_Retired( 0.1, tmpDict_Retired1, tmpConst_Retired1 )
tmpks_Retired1 = House.getks_Retired( tmpcs_Retired1, tmpDict_Retired1, tmpConst_Retired1 )
House.G_Retired( 0.1, tmpDict_Retired1, tmpConst_Retired1 )
@time tmpRes_Retired1 = House.HHSolve_Retired( tmpOriginData_Retired1 , ReturnData = true )
House.ExtractAPhi_Retired!( tmpRes_Retired1, tmpOriginData_Retired1 , a1 = 0.0 )

#


text: using Pkg;Pkg.add(["IJulia", "PyPlot", "Plots", "SparseIR", "Revise", "FFTW", "Roots", "OMEinsum", "GR"])

max_stars_repo_name: SpM-lab/sparse-ir-doc

max_issues_repo_name: SpM-lab/sparse-ir-doc

max_forks_repo_name: SpM-lab/sparse-ir-doc


text:
#=
Stack
=#

"""
Space complexity: O(n)
"""
struct SimpleStack{T}
    l::SimpleLinkedList{T}

    SimpleStack{T}() where T = new{T}(SimpleLinkedList{T}())
end

#=
Core methods
=#

"""
Time complexity: O(1)
"""
function Base.push!(ss::SimpleStack{T}, item::T) where T
    push!(ss.l, item)
    return ss
end

"""
Time complexity: O(1)
"""
Base.pop!(ss::SimpleStack) = pop!(ss.l)


text:
function simulate_waveforms( mcevents::TypedTables.Table, sim::Simulation{T},
                             output_dir::AbstractString,
                             output_base_name::AbstractString = "generated_waveforms";
                             chunk_n_physics_events::Int = 1000,
                             Δt::RealQuantity = 4u"ns",
                             max_nsteps::Int = 1000,
                             diffusion::Bool = false,
                             self_repulsion::Bool = false,
                             number_of_carriers::Int = 1,
                             number_of_shells::Int = 1,
                             verbose = false) where {T <: SSDFloat}
    n_total_physics_events = length(mcevents)
    Δtime = T(to_internal_units(Δt))
    n_contacts = length(sim.detector.contacts)
    @info "Detector has $(n_contacts) contact(s)"
    contacts = sim.detector.contacts
    if !ispath(output_dir) mkpath(output_dir) end
    nfmt(i::Int) = format(i, zeropadding = true, width = length(digits(n_total_physics_events)))
    evt_ranges = chunked_ranges(n_total_physics_events, chunk_n_physics_events)
    @info "-> $(length(flatview(mcevents.edep))) energy depositions to simulate."

    for evtrange in evt_ranges
        ofn = joinpath(output_dir, "$(output_base_name)_evts_$(nfmt(first(evtrange)))-$(nfmt(last(evtrange))).h5")
        @info "Now simulating $(evtrange) and storing it in\n\t \"$ofn\""
        mcevents_sub = simulate_waveforms(mcevents[evtrange], sim; Δt, max_nsteps, diffusion, self_repulsion, number_of_carriers, number_of_shells, verbose)

        HDF5.h5open(ofn, "w") do output
            LegendHDF5IO.writedata(output, "generated_waveforms", mcevents_sub)
        end
    end
end


text:
__precompile__(true)

module FixedEffectModels

##############################################################################
##
## Dependencies
##
##############################################################################
using Compat
import Base.BLAS: axpy!
import Base: A_mul_B!, Ac_mul_B!, size, sumabs2, copy!, getindex, length, fill!, norm, scale!, eltype, length, view, start, next, done
import Distributions: TDist, ccdf, FDist, Chisq, AliasTable, Categorical
import DataArrays: RefArray, PooledDataArray, PooledDataVector, DataArray, DataVector, compact, NAtype
import DataFrames: @~, DataFrame, AbstractDataFrame, ModelMatrix, ModelFrame, Terms, coefnames, Formula, complete_cases, names!, pool
import StatsBase: coef, nobs, coeftable, vcov, predict, residuals, var, RegressionModel, model_response, stderr, confint, fit, CoefTable, df_residual
##############################################################################
##
## Exported methods and types
##
##############################################################################

export group,
reg,
partial_out,
residualize!,
getfe!,
decompose!,
allvars,

Ones,
FixedEffect,
FixedEffectProblem,

AbstractRegressionResult,
title,
top,
RegressionResult,
RegressionResultIV,
RegressionResultFE,
RegressionResultFEIV,

AbstractVcovMethod,
AbstractVcovMethodData,
vcov!,
shat!,
VcovMethodData,
VcovData,
VcovSimple,
VcovWhite,
VcovCluster

##############################################################################
##
## Load files
##
##############################################################################
include("utils/group.jl")
include("utils/formula.jl")
include("utils/lsmr.jl")
include("utils/basecol.jl")
include("utils/Ones.jl")
include("utils/combinations.jl")

include("RegressionResult.jl")

include("fixedeffect/FixedEffect.jl")
include("fixedeffect/FixedEffectProblem.jl")
include("fixedeffect/FixedEffectProblem_LSMR.jl")
include("fixedeffect/FixedEffectProblem_Factorization.jl")

include("vcov/types.jl")
include("vcov/vcovsimple.jl")
include("vcov/vcovwhite.jl")
include("vcov/vcovcluster.jl")
include("vcov/ranktest.jl")

include("reg.jl")
include("partial_out.jl")

end  # module FixedEffectModels


text:
export iauPvtob
"""
Position and velocity of a terrestrial observing station.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
    elong   double       longitude (radians, east +ve, Note 1)
    phi     double       latitude (geodetic, radians, Note 1)
    hm      double       height above ref. ellipsoid (geodetic, m)
    xp,yp   double       coordinates of the pole (radians, Note 2)
    sp      double       the TIO locator s' (radians, Note 2)
    theta   double       Earth rotation angle (radians, Note 3)

Returned:
    pv      double[2][3] position/velocity vector (m, m/s, CIRS)

Notes:

    1. The terrestrial coordinates are with respect to the WGS84
    reference ellipsoid.

    2. xp and yp are the coordinates (in radians) of the Celestial
    Intermediate Pole with respect to the International Terrestrial
    Reference System (see IERS Conventions), measured along the
    meridians 0 and 90 deg west respectively.  sp is the TIO locator
    s', in radians, which positions the Terrestrial Intermediate
    Origin on the equator.  For many applications, xp, yp and
    (especially) sp can be set to zero.

    3. If theta is Greenwich apparent sidereal time instead of Earth
    rotation angle, the result is with respect to the true equator
    and equinox of date, i.e. with the x-axis at the equinox rather
    than the celestial intermediate origin.

    4. The velocity units are meters per UT1 second, not per SI second.
    This is unlikely to have any practical consequences in the modern
    era.

    5. No validation is performed on the arguments.  Error cases that
    could lead to arithmetic exceptions are trapped by the iauGd2gc
    function, and the result set to zeros.

References:

    McCarthy, D. D., Petit, G. (eds.), IERS Conventions (2003),
    IERS Technical Note No. 32, BKG (2004)

    Urban, S. & Seidelmann, P. K. (eds), Explanatory Supplement to
    the Astronomical Almanac, 3rd ed., University Science Books
    (2013), Section 7.4.3.3.

Called:
    iauGd2gc     geodetic to geocentric transformation
    iauPom00     polar motion matrix
    iauTrxp      product of transpose of r-matrix and p-vector

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPvtob(double elong, double phi, double hm,
#               double xp, double yp, double sp, double theta,
#               double pv[2][3])

function iauPvtob(elong::Real, phi::Real, hm::Real,
                  xp::Real, yp::Real, sp::Real, theta::Real)

    pv = zeros(Float64, 3, 2)

    status = ccall((:iauPvtob, libsofa_c), Cvoid,
         (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble,
          Ptr{Float64}),
         convert(Float64, elong),
         convert(Float64, phi),
         convert(Float64, hm),
         convert(Float64, xp),
         convert(Float64, yp),
         convert(Float64, sp),
         convert(Float64, theta), pv)

    return SMatrix{2,3}(pv')
end


text:
COBREXA.@_make_logging_tag TEST "testing stuff"

log_TEST()
@testset "Logging on" begin
    @test _TEST_log_enabled
    @test_logs (:warn, "qweasdzxc") @_TEST_log @warn "qweasdzxc"
end

log_TEST(false)
@testset "Logging off" begin
    @test !_TEST_log_enabled
    @test_logs @_TEST_log @warn "all okay!"
end


text:
using LinearAlgebra
using Plots
using Revise

includet("src/IterativeMethods.jl")
includet("src/InternalArray.jl")

IArray = InternalArray

t0 = 0
x0 = [0, 100]
v0 = [300, 0]
m = 1
c = 0.1
g = [0, -9.8]

# Solving without friction
x_nofriction = eulermethod(t0, x0) do t, x
	v0 + t*g
end
plot_nofriction = plot(IArray(x_nofriction, 1), IArray(x_nofriction, 2), label = "without friction")

# Solving with friction
x_friction = eulermethod(t0, [x0; v0]) do t, x
	v = [x[3], x[4]]
	return [v; g-c/m*norm(v)*v]
end
plot_friction = plot(IArray(x_friction, 1), IArray(x_friction, 2), label = "with friction")

plot(plot_nofriction, plot_friction, plot_title = "comparing with and without friction")

# Solving for collision with floor
t, x = targetting_eulermethod(x -> x[2], t0, [x0; v0]; h = 5e-3) do t, x
	v = [x[3], x[4]]
	return [v; g-c/m*norm(v)*v]
end
p = plot(IArray(x, 1), IArray(x, 2), label = "landing")

# Needed to stop showing due to error
using IJulia

x1 = IArray(x, 1)
x2 = IArray(x, 2)
v1 = IArray(x, 3)
v2 = IArray(x, 4)

ppp = 35 # points per frame
lims(v) = (minimum(v), maximum(v))
p1 = plot(1, xlim = lims(x1), ylim = lims(x2), label = "position")
p2 = plot(1, xlim = lims(v1), ylim = lims(v2), label = "speed")
anim = Animation()
for i = ppp:ppp:length(x)
	for j = i-ppp+1:i
    	push!(p1, x1[j], x2[j])
    	push!(p2, v1[j], v2[j])
    end
	plot(p1, p2)
	frame(anim)
end
gif(anim, "media/dropped-object-trayectory-to-collision.gif", fps = 24)

max_stars_repo_name: useredsa/numeric-differential-equations

max_issues_repo_name: useredsa/numeric-differential-equations

max_forks_repo_name: useredsa/numeric-differential-equations