text:
# ------------------------------------
# advection term by AUSM
# ------------------------------------
function AUSM(E_adv_hat, F_adv_hat, Qbase, Qcon, cellxmax, cellymax, vecAx, vecAy, specific_heat_ratio, volume, nval, Minf, ad_scheme)
    g         = specific_heat_ratio
    temp_vecX = zeros(nval)   # vecAx used for the pressure term
    temp_vecY = zeros(nval)   # vecAx used for the pressure term
    Lpsi = zeros(nval)
    Rpsi = zeros(nval)

    mdot = 0.0
    ph = 0.0

    tvAxx = 0.0
    tvAxy = 0.0
    vAxx = 0.0
    vAxy = 0.0
    tvAyx = 0.0
    tvAyy = 0.0
    vAyx = 0.0
    vAyy = 0.0
    for j in 2:cellymax -1
        for i in 2:cellxmax+1 -1
            volume_av = 0.5*(volume[i-1,j] + volume[i,j])

            # i-1 cell
            # 規格化
            tvAxx = 0.5*(vecAx[i-1,j,1]+vecAx[i,j,1])
            tvAxy = 0.5*(vecAx[i-1,j,2]+vecAx[i,j,2])
            vAxx = tvAxx / (tvAxx^2+tvAxy^2)^0.5
            vAxy = tvAxy / (tvAxx^2+tvAxy^2)^0.5

            rhoL = Qbase[i-1,j,1]
            UL = Qbase[i-1,j,2]*vAxx + Qbase[i-1,j,3]*vAxy
            pL = Qbase[i-1,j,4]

            # i cell
            # 規格化
            tvAxx = 0.5*(vecAx[i,j,1]+vecAx[i+1,j,1])
            tvAxy = 0.5*(vecAx[i,j,2]+vecAx[i+1,j,2])
            vAxx = tvAxx / (tvAxx^2+tvAxy^2)^0.5
            vAxy = tvAxy / (tvAxx^2+tvAxy^2)^0.5

            rhoR = Qbase[i,j,1]
            UR = Qbase[i,j,2]*vAxx + Qbase[i,j,3]*vAxy
            pR = Qbase[i,j,4]

            #=
            if i==100 && j== 200
                mdot, ph = AUSM_plus_half(rhoL, rhoR, UL, UR, pL, pR, g,i,j)
                println("AUSM+")
                println(mdot)
                println(ph)

                mdot, ph = AUSM_plusup_half(rhoL, rhoR, UL, UR, pL, pR, Minf, g,i,j)
                println("AUSM+up")
                println(mdot)
                println(ph)

                velocity = 0.5 * ((Qbase[i-1,j,2]^2 + Qbase[i-1,j,3]^2)^0.5 + (Qbase[i,j,2]^2 + Qbase[i,j,3]^2)^0.5)
                mdot, ph = SLAU_half(rhoL, rhoR, UL, UR, pL, pR, velocity, specific_heat_ratio)
                println("SLAU")
                println(mdot)
                println(ph)
            end
            =#

            # scheme
            if  ad_scheme == 1
                mdot, ph = AUSM_plus_half(rhoL, rhoR, UL, UR, pL, pR, g,i,j)
            elseif ad_scheme == 2
                mdot, ph = AUSM_plusup_half(rhoL, rhoR, UL, UR, pL, pR, Minf, g, i,j)
            elseif ad_scheme == 4
                velocity = (0.5 * (Qbase[i-1,j,2]^2 + Qbase[i-1,j,3]^2 + Qbase[i,j,2]^2 + Qbase[i,j,3]^2))^0.5
                mdot, ph = SLAU_half(rhoL, rhoR, UL, UR, pL, pR, velocity, specific_heat_ratio)
            end

            # flux half
            sqAx = (vecAx[i,j,1]^2 + vecAx[i,j,2]^2)^0.5
            temp_vecX[2] = vecAx[i,j,1] / sqAx
            temp_vecX[3] = vecAx[i,j,2] / sqAx

            for l in 1:nval
                Lpsi[l] = Qcon[i-1,j,l] / Qcon[i-1,j,1]
                Rpsi[l] = Qcon[i,j,l]   / Qcon[i,j,1]
            end
            Lpsi[4] = (Qcon[i-1,j,4] + Qbase[i-1,j,4]) / Qcon[i-1,j,1]
            Rpsi[4] = (Qcon[i,j,4]   + Qbase[i,j,4])   / Qcon[i,j,1]

            if mdot > 0
                for l in 1:nval
                    E_adv_hat[i,j,l] = (mdot * Lpsi[l] + ph * temp_vecX[l]) * sqAx
                end
            else
                for l in 1:nval
                    E_adv_hat[i,j,l] = (mdot * Rpsi[l] + ph * temp_vecX[l]) * sqAx
                end
            end
        end
    end

    for j in 2:cellymax+1 -1
        for i in 2:cellxmax -1
            volume_av = 0.5*(volume[i,j-1] + volume[i,j])

            # j-1 cell
            # 規格化
            tvAyx = 0.5*(vecAy[i,j-1,1]+vecAy[i,j,1])
            tvAyy = 0.5*(vecAy[i,j-1,2]+vecAy[i,j,2])
            vAyx = tvAyx / (tvAyx^2+tvAyy^2)^0.5
            vAyy = tvAyy / (tvAyx^2+tvAyy^2)^0.5

            rhoL = Qbase[i,j-1,1]
            VL = Qbase[i,j-1,2]*vAyx + Qbase[i,j-1,3]*vAyy
            pL = Qbase[i,j-1,4]

            # j cell
            tvAyx = 0.5*(vecAy[i,j,1]+vecAy[i,j+1,1])
            tvAyy = 0.5*(vecAy[i,j,2]+vecAy[i,j+1,2])
            vAyx = tvAyx / (tvAyx^2+tvAyy^2)^0.5
            vAyy = tvAyy / (tvAyx^2+tvAyy^2)^0.5

            rhoR = Qbase[i,j,1]
            VR = Qbase[i,j,2]*vAyx + Qbase[i,j,3]*vAyy
            pR = Qbase[i,j,4]

            # scheme
            if  ad_scheme == 1
                mdot, ph = AUSM_plus_half(rhoL, rhoR, VL, VR, pL, pR, g,i,j)
            elseif ad_scheme == 2
                mdot, ph = AUSM_plusup_half(rhoL, rhoR, VL, VR, pL, pR, Minf, g, i,j)
            elseif ad_scheme == 4
                velocity = (0.5 * (Qbase[i,j-1,2]^2 + Qbase[i,j-1,3]^2 + Qbase[i,j,2]^2 + Qbase[i,j,3]^2))^0.5
                mdot, ph = SLAU_half(rhoL, rhoR, VL, VR, pL, pR, velocity, specific_heat_ratio)
            end
            #=
            if i==170 && j== 2
                mdot, ph = AUSM_plus_half(rhoL, rhoR, VL,VR, pL, pR, g)
                println("AUSM+")
                println(mdot)
                println(ph)

                mdot, ph = AUSM_plusup_half(rhoL, rhoR, VL,VR, pL, pR, Minf, g)
                println("AUSM+up")
                println(mdot)
                println(ph)

                velocity = (0.5 * (Qbase[i,j-1,2]^2 + Qbase[i,j-1,3]^2 + Qbase[i,j,2]^2 + Qbase[i,j,3]^2))^0.5
                mdot, ph = SLAU_half(rhoL, rhoR, VL, VR, pL, pR, velocity, specific_heat_ratio)
                println("SLAU")
                println(mdot)
                println(ph)
            end
            =#

            # flux half
            sqAy = (vecAy[i,j,1]^2 + vecAy[i,j,2]^2)^0.5
            temp_vecY[2] = vecAy[i,j,1] / sqAy
            temp_vecY[3] = vecAy[i,j,2] / sqAy

            for l in 1:nval
                Lpsi[l] = Qcon[i,j-1,l] / Qcon[i,j-1,1]
                Rpsi[l] = Qcon[i,j,l]   / Qcon[i,j,1]
            end
            Lpsi[4] = (Qcon[i,j-1,4] + Qbase[i,j-1,4]) / Qcon[i,j-1,1]
            Rpsi[4] = (Qcon[i,j,4]   + Qbase[i,j,4])   / Qcon[i,j,1]

            if mdot > 0
                for l in 1:nval
                    F_adv_hat[i,j,l] = (mdot * Lpsi[l] + ph * temp_vecY[l]) * sqAy
                end
            else
                for l in 1:nval
                    F_adv_hat[i,j,l] = (mdot * Rpsi[l] + ph * temp_vecY[l]) * sqAy
                end
            end
        end
    end

    return E_adv_hat, F_adv_hat
end

function AUSM_plus_half(rhoL, rhoR, UL, UR, pL, pR, g,i,j)
    # param
    beta  = 1/8
    alpha = 3/16

    # L, R
    aL = (g * pL / rhoL)^0.5
    aR = (g * pR / rhoR)^0.5

    # half
    ah   = 0.5*( aL + aR )
    rhoh = 0.5*( rhoL + rhoR )

    ML = UL/ah
    MR = UR/ah

    M_p4 = 0
    M_m4 = 0
    p_p5 = 0
    p_m5 = 0
    if abs(ML) >= 1
        M_p4 = 0.5*(ML + abs(ML))
        p_p5 = 0.5*(ML + abs(ML)) / ML
    else
        M_p4 = 0.25*(ML+1)^2 + beta*(ML^2-1)^2
        p_p5 = 0.25*(ML+1)^2 * (2-ML) + alpha*ML*(ML^2-1)^2
    end

    if abs(MR) >= 1
        M_m4 = 0.5*(MR - abs(MR))
        p_m5 = 0.5*(MR - abs(MR)) / MR
    else
        M_m4 = -0.25*(MR-1)^2 - beta*(MR^2-1)^2
        p_m5 = 0.25*(MR-1)^2 * (2+MR) - alpha*MR*(MR^2-1)^2
    end

    # Mh = M_p4 + M_m4 + Mp/fa
    Mh = M_p4 + M_m4
#=
    if i==100 && j==200
        println(" MMM_AUSM+ ")
        println(M_p4)
        println(M_m4)
    end
  =#
    # mdot half
    mdot = ah * Mh
    if Mh > 0
        mdot = mdot * rhoL
    else
        mdot = mdot * rhoR
    end

    # p half
    ph = p_p5*pL + p_m5*pR

    return mdot, ph
end

function AUSM_plusup_half(rhoL, rhoR, UL, UR, pL, pR, Minf, g, i,j)
    # param
    beta  = 1/8
    kp    = 0.25
    sigma = 1.0
    ku    = 0.75

    # L, R
    aL = (g * pL / rhoL)^0.5
    aR = (g * pR / rhoR)^0.5

    # half
    ah   = 0.5*( aL + aR )
    rhoh = 0.5*( rhoL + rhoR )

    Mbar = (( UL^2 + UR^2 ) / ( 2 * ah^2 ))^0.5
    Mo   = (min(1,max( Mbar^2, Minf^2 )))^0.5
    fa   = Mo * (2-Mo)

    alpha = 3/16 * ( -4 + 5*fa^2 )

    ML = UL/ah
    MR = UR/ah

    M_p4 = 0
    M_m4 = 0
    p_p5 = 0
    p_m5 = 0
    if abs(ML) >= 1
        M_p4 = 0.5*(ML + abs(ML))
        p_p5 = 0.5*(ML + abs(ML)) / ML
    else
        Mtp  = 0.25*(ML + 1)^2
        Mtm  = -0.25*(ML - 1)^2
        M_p4 = Mtp * (1 - 16*beta*Mtm)
        p_p5 = Mtp * ((2-ML) - 16*alpha*ML*Mtm)
    end

    if abs(MR) >= 1
        M_m4 = 0.5*(MR - abs(MR))
        p_m5 = 0.5*(MR - abs(MR)) / MR
    else
        Mtp  = 0.25*(MR + 1)^2
        Mtm  = -0.25*(MR - 1)^2
        M_m4 = Mtm * (1 + 16*beta*Mtp)
        p_m5 = Mtm * ((-2-MR) + 16*alpha*MR*Mtp)
    end

    # M half
    Mp = -kp * max( 1 - sigma*Mbar^2, 0.0) * (pR - pL)/(rhoh*ah^2)

    Mh = M_p4 + M_m4 + Mp/fa

    #=
    if i==100 && j==200
        println(" MMM_AUSM+up ")
        println(M_p4)
        println(M_m4)
        println(Mp/fa)
    end
    =#

    # mdot half
    mdot = ah * Mh
    if Mh > 0
        mdot = mdot * rhoL
    else
        mdot = mdot * rhoR
    end

    # p half
    pu = -ku * p_p5 * p_m5 * (rhoL + rhoR) * ah *(UR - UL)

    ph = p_p5*pL + p_m5*pR + fa * pu
    return mdot, ph
end

function set_Minf(bdcon, specific_heat_ratio, Rd, nval)
    Minf = 0.0
    check = 0
    for i in 1:4
        if Int(bdcon[i][1]) == 0 || Int(bdcon[i][1]) == 5
            rho = bdcon[i][2]
            u = bdcon[i][3]
            v = bdcon[i][4]
            p = bdcon[i][5]

            a = (specific_heat_ratio * p / rho)^0.5
            Minf = (u^2+v^2)^0.5 / a

            println(" Minf = " * string(Minf))

            check = 1

        elseif Int(bdcon[i][1]) == 6
            rho = bdcon[i][2]
            u = bdcon[i][3]
            v = bdcon[i][4]
            T = bdcon[i][nval+2]

            p = (rho*Rd) * T

            a = (specific_heat_ratio * p / rho)^0.5
            Minf = (u^2+v^2)^0.5 / a

            println(" Minf = " * string(Minf))

            check = 1
        end
    end

    if check == 0
        println(" Minf error ")
        println("  ")
        println(" if you don't use inlet condition, ")
        println(" don't use AUSM+up ")
        println("  ")
        throw(UndefVarError(:x))
    end

    return Minf
end

function SLAU_half(rhoL, rhoR, UL, UR, pL, pR, velocity, specific_heat_ratio)
    mdot =1.0
    ph =1.0

    # L, R
    aL = (specific_heat_ratio * pL / rhoL)^0.5
    aR = (specific_heat_ratio * pR / rhoR)^0.5

    # half
    ah   = 0.5*( aL + aR )

    ML = UL/ah
    MR = UR/ah

    temp1 = -max( min( ML, 0.0 ), -1.0)
    temp2 =  min( max( MR, 0.0 ), 1.0)
    g = temp1 * temp2

    barUL = (1-g) * (rhoL*abs(UL) + rhoR*abs(UR))/(rhoL + rhoR) + g*abs(UL)
    barUR = (1-g) * (rhoL*abs(UL) + rhoR*abs(UR))/(rhoL + rhoR) + g*abs(UR)

    hatM = min(1.0, 1/ah * velocity)
    chi  = (1-hatM)^2

    beta_a = 0
    beta_b = 0
    if abs(MR) >= 1
        beta_a = 0.5 * (1+sign(-MR))
    else
        beta_a = 0.25 * (2+MR) * (MR-1)^2
    end
    if abs(ML) >= 1
        beta_b = 0.5 * (1+sign(ML))
    else
        beta_b = 0.25 * (2-ML) * (ML+1)^2
    end

    # mdot half
    mdot = 0.5 * (rhoL*(UL+abs(barUL)) + rhoR*(UR-abs(barUR)) - chi/ah*(pR-pL))

    # p half
    ph = 0.5*(pR+pL) + 0.5*(beta_b-beta_a)*(pL-pR) + (1-chi) * (beta_b+beta_a-1) * 0.5*(pR+pL)

    # SLAU2
    #ph = 0.5*(pR+pL) + 0.5*(beta_b-beta_a) + ((UR^2+UL^2)/2)^0.5 * (beta_b+beta_a-1) * 0.5*(rhoR+rhoL)*ah

    return mdot, ph
end


text:
# ============================================================================ #
# Circle
# ============================================================================ #

abstract type AbstractCircle end

struct Circle{T <: Real, PointType <: AbstractPoint} <: AbstractCircle
    center::PointType
    r::T
end

Circle(r::T) where {T} = Circle(Point(0, 0), r)
Circle(a::AbstractArray, r::T) where {T} = Circle(Point(a[1], a[2]), r)
Circle(x0::T, y0::T, r::S) where {T, S} = Circle(Point(x0, y0), r)

# ---------------------------------------------------------------------------- #
# Methods
# ---------------------------------------------------------------------------- #

Base.:(==)(c1::Circle, c2::Circle) = c1.center == c2.center && c1.r == c2.r

function shift(circle::AbstractCircle, v)
    return Circle( circle.x + v[1], circle.y + v[2], circle.r )
end

function shift(circle::T, v) where {T <: AbstractCircle}
    return T( circle.x + v[1], circle.y + v[2], circle.r )
end
function scale(c::T, factor) where {T <: AbstractCircle}
    return T( circle.x, circle.y, circle.r * factor )
end
function rotate(c::T, rot_mat) where {T <: AbstractCircle}
    return c
end

diameter(circle::Circle) = circle.r * 2
circumference(circle::Circle) = 2 * circle.r * π
area(circle::Circle) = circle.r^2 * π
is_in_circle(c::Circle, p::AbstractPoint) = distance( c.center, p ) <= c.r

function discretize_for_plot(c::Circle, n::Int = 100)
    phi = LinRange(0, 2*pi, n)
    x = c.center.x .+ c.r .* cos.(phi)
    y = c.center.y .+ c.r .* sin.(phi)
    return (x, y)
end

function mid_point(circle::T) where {T <: AbstractCircle}
    return circle.center
end

ext: jl

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
module GlobTest
include("../src/Glob.jl")
using .Glob
using Base.Test
using Compat

if VERSION < v"0.3-"
    macro test_throws(a,b)
        :( Base.Test.@test_throws($b) )
    end
end

macro test_types(arr, types)
    quote
        @test length($arr) == length($types)
        for i in 1:length($arr)
            @test isa($arr[i], $types[i])
        end
    end
end

# various unit tests, in no particular order

@test !ismatch(fn"abc*", "ABXABXAB")
@test ismatch(fn"AB*AB*AB", "ABXABXAB")
@test ismatch(fn"AB*AB*AB", "ABXABAB")
@test !ismatch(fn"AB*AB*AB", "AABAB")
@test ismatch(fn"AB*AB*AB", "ABABAB")
@test ismatch(fn"AB*AB*B", "ABABAB")
@test ismatch(fn"AB*AB*", "ABABAB")
@test ismatch(fn"AB*AB*??", "ABABAB")
@test !ismatch(fn"AB*AB*???", "ABABAB")
@test ismatch(fn"AB*AB*??***", "ABABAB")
@test ismatch(fn"AB*AB*??***", "ABABABC")
@test ismatch(fn"AB*AB*??***", "ABABABCDEFG")
@test ismatch(fn"?AB*AB*??***", ".ABABABCDEFG")
@test !ismatch(fn"?AB*AB*??***"p, ".ABABABCDEFG")
@test_throws ErrorException Glob.FilenameMatch("?AB*AB*??***","z")
@test ismatch(fn"[abc]", "a")
@test !ismatch(fn"[abc]", "A")
@test ismatch(fn"[abc]"i, "A")
@test ismatch(fn"[abc]", "b")
@test ismatch(fn"[abc]", "c")
@test !ismatch(fn"[abc]", "x")
@test !ismatch(fn"[a-d]", "x")
@test ismatch(fn"[a-d]", "a")
@test ismatch(fn"[a-d]", "d")
@test !ismatch(fn"[a--]", "e")
@test ismatch(fn"[--e]", "e")
@test !ismatch(fn"[f-e]", "e")
@test !ismatch(fn"[g-e]", "e")
@test !ismatch(fn"[g-f]", "e")
@test ismatch(fn"[a-f]", "e")
@test ismatch(fn"[a-e]", "e")
@test ismatch(fn"[f-ee]", "e")
@test !ismatch(fn"[A-Z]", "e")
@test ismatch(fn"[A-z]", "e")
@test ismatch(fn"[\]", "\\")
@test_throws ErrorException ismatch(fn"[[:a:]]", "e")
@test_throws ErrorException !ismatch(fn"[\[:a:]]", "e")
@test !ismatch(fn"[\[:a:]]"x, "e")
@test_throws ErrorException ismatch(fn"[\\[:a:]]", "e")
@test ismatch(fn"[\[:a-e:]]"x, "e]")
@test ismatch(fn"[\[:a-e:]"x, "e")
@test !ismatch(fn"[\[:a-e:xxxx"x, "e")
@test ismatch(fn"[\[:a-e:xxxx-]"x, "e")
@test ismatch(fn"[a-]", "a")
@test ismatch(fn"[a-]", "-")
@test !ismatch(fn"[a-]", "b")
@test ismatch(fn"[!a-]", "b")
@test !ismatch(fn"[!a-]", "a")
@test ismatch(fn"[!a]", "!")
@test !ismatch(fn"[!!]", "!")
@test !ismatch(fn"[!a]", "a")
@test ismatch(fn"[!a]", "b")
@test !ismatch(fn"[][]", "")
@test ismatch(fn"[]", "[]")
@test ismatch(fn"[][]", "[")
@test ismatch(fn"[][]", "]")
@test !ismatch(fn"[][]", "x")
@test !ismatch(fn"[]-[]", "x")
@test !ismatch(fn"[]-[]", "-")
@test !ismatch(fn"[]-[]", "\\")
@test ismatch(fn"[\[-\]]*"x, "][")
@test ismatch(fn"[\]]*"x, "][")
@test ismatch(fn"[\[-\]]*"x, "][")
@test ismatch(fn"[[-\]]*"x, "][")
@test ismatch(fn"base/[\[-\]]*"dpx,"base/][x")
@test ismatch(fn"[\[-\]]"x, "\\")
@test ismatch(fn"[[-\]]"x, "\\")
@test ismatch(fn"[---]", "-")
@test !ismatch(fn"[!---]", "-")
@test ismatch(fn"[!---]", "0")
@test !ismatch(fn"[---a-zA-Z]", "0")
@test !ismatch(fn"[---a-zA-Z:]", "0")
@test !ismatch(fn"[---!]", "0")
@test ismatch(fn"[---!]", "!")
@test !ismatch(fn"[---!]", "a")
@test !ismatch(fn"[---!]", "\0")
@test ismatch(fn"ab/c/d"dp, "ab/c/d")
@test !ismatch(fn"ab/c/d"dp, "ab/c?d")
@test !ismatch(fn"ab/./d"dp, "ab/?/d")
@test !ismatch(fn"ab*d"dp, "aba/d")
@test !ismatch(fn"ab*d"dp, "ab/d")
@test ismatch(fn"ab*d", "ab/d")
@test ismatch(fn"ab*d", "aba/d")
@test ismatch(fn"[a-z]"i, "B")
@test !ismatch(fn"[a-z]"i, "_")
@test ismatch(fn"[A-z]"i, "_")
@test !ismatch(fn"[a-Z]"i, "_")
@test !ismatch(fn"#[a-Z]%"i, "#a%")
@test ismatch(fn"#[α-ω]%"i, "#Γ%")
@test !ismatch(fn"#[α-ω]%", "#Γ%")
@test ismatch(fn"#[α-ω]%", "#γ%")
@test !ismatch(fn"a?b"d, "a/b")
@test ismatch(fn"a?b", "a/b")
@test !ismatch(fn"?b"p, ".b")
@test ismatch(fn"?b", ".b")
@test ismatch(fn"?/?b", "./.b")
@test !ismatch(fn"?/?b"p, "./.b")
@test ismatch(fn"./?b"p, "./.b")
@test !ismatch(fn"./?b"pd, "./.b")
@test ismatch(fn"./.b"pd, "./.b")
@test !ismatch(fn"?/.b"pd, "./.b")
@test ismatch(fn"""./.b"""pd, "./.b")
@test !ismatch(fn"""?/.b"""pd, "./.b")
@test ismatch(fn"_[[:blank:][.a.]-c]_", "_b_")
@test !ismatch(fn"_[[:blank:][.a.]-c]_", "_-_")
@test ismatch(fn"_[[:blank:][.a.]-c]_", "_ _")
@test ismatch(fn"_[[:alnum:]]_", "_a_")
@test !ismatch(fn"_[[:alnum:]]_", "_[_")
@test !ismatch(fn"_[[:alnum:]]_", "_]_")
@test !ismatch(fn"_[[:alnum:]]_", "_:_")
@test ismatch(fn"_[[:alpha:]]_", "_z_")
@test !ismatch(fn"_[[:alpha:]]_", "_[_")
@test ismatch(fn"_[[:cntrl:]]_", "_\0_")
@test ismatch(fn"_[[:cntrl:]]_", "_\b_")
@test !ismatch(fn"_[[:cntrl:]]_", "_:_")
@test !ismatch(fn"_[[:cntrl:]]_", "_ _")
@test ismatch(fn"_[[:digit:]]_", "_0_")
@test !ismatch(fn"_[[:digit:]]_", "_:_")
@test ismatch(fn"_[[:graph:]]_", "_._")
@test !ismatch(fn"_[[:graph:]]_", "_ _")
@test ismatch(fn"_[[:lower:]]_", "_a_")
@test ismatch(fn"_[[:lower:]]_"i, "_A_")
@test !ismatch(fn"_[[:lower:]]_", "_A_")
@test !ismatch(fn"_[[:lower:]]_", "_:_")
@test ismatch(fn"_[[:print:]]_", "_a_")
@test !ismatch(fn"_[[:print:]]_", "_\7_")
@test ismatch(fn"_[[:punct:]]_", "_:_")
@test !ismatch(fn"_[[:punct:]]_", "_p_")
@test ismatch(fn"_[[:space:]]_", "_\f_")
@test !ismatch(fn"_[[:space:]]_", "_:_")
@test !ismatch(fn"_[[:space:]]_", "_\r\n_")
@test ismatch(fn"_[[:upper:]]_", "_A_")
@test ismatch(fn"_[[:upper:]]_"i, "_a_")
@test !ismatch(fn"_[[:upper:]]_", "_a_")
@test !ismatch(fn"_[[:upper:]]_", "_:_")
@test ismatch(fn"_[[:xdigit:]]_", "_a_")
@test !ismatch(fn"_[[:xdigit:]]_", "_:_")
@test ismatch(fn"_[[.a.]-[.z.]]_", "_c_")
@test !ismatch(fn"_[[.a.]-[.z.]]_", "_-_")
@test !ismatch(fn"_[[.a.]-[.z.]]_", "_]_")
@test !ismatch(fn"_[[.a.]-[.z.]]_", "_[_")
@test ismatch(fn"_[[=a=]]_", "_a_")
@test !ismatch(fn"_[[=a=]]_", "_=_")
@test !ismatch(fn"_[[=a=]]_", "_á_")
@test ismatch(fn"[[=a=]-z]", "-")
@test_throws ErrorException ismatch(fn"[a-[=z=]]", "e")

@test !ismatch(fn"\?", "\\?")
@test ismatch(fn"\?", "?")
@test ismatch(fn"\?"e, "\\!")
@test !ismatch(fn"\?"e, "?")

@test_types glob"ab/?/d".pattern (AbstractString, Glob.FilenameMatch, AbstractString)
@test_types glob"""ab/*/d""".pattern (AbstractString, Glob.FilenameMatch, AbstractString)
@test length(glob"ab/[/d".pattern) == 3
@test length(glob"ab/[/]d".pattern) == 3
@test_types glob"ab/[/]d".pattern (AbstractString, AbstractString, AbstractString)
@test_types glob"ab/[/d".pattern (AbstractString, AbstractString, AbstractString)
@test_types glob"ab/[]/d".pattern (AbstractString, AbstractString, AbstractString)
@test_types glob"ab/[]]/d".pattern (AbstractString, Glob.FilenameMatch, AbstractString)

@test glob("*") == filter(x->!startswith(x,'.'), readdir()) == readdir(glob"*")
@test glob(".*") == filter(x->startswith(x,'.'), readdir()) == readdir(glob".*")
@test isempty(Glob.glob("[.]*"))
@test glob([r".*"]) == readdir()
@test glob([".", r".*"]) == map(x->joinpath(".",x), readdir())
@test all([!startswith(x,'.') for x in Glob.glob("*.*")])

function test_string(x1)
    x2 = string(eval(parse(x1)))
    x1 == x2 ? nothing : error(string(
        "string test failed:",
        "\noriginal: ", x1,
        "\n\nstringify: ", x2))
end

if VERSION > v"0.4-dev"
    test_string("""Glob.GlobMatch(Any["base",r"h\.+"])""")
else
    test_string("""Glob.GlobMatch({"base",r"h\.+"})""")
end

test_string("""glob"base/*/a/[b]\"""")
test_string("""fn"base/*/a/[b]\"ipedx""")
test_string("""fn"base/*/a/[b]\"""")

@test_throws ErrorException Glob.GlobMatch("")
@test_throws ErrorException Glob.GlobMatch("/a/b/c")

end


text:
using Dates
using GaussQuadrature
using JLD2
using LinearAlgebra
using Logging
using MPI
using Printf
using StaticArrays
using Test

using ClimateMachine
using ClimateMachine.Atmos: NoReferenceState
using ClimateMachine.MPIStateArrays
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Geometry
using ClimateMachine.Mesh.Topologies
using ClimateMachine.Mesh.Grids
using ClimateMachine.DGMethods
using ClimateMachine.DGMethods.NumericalFluxes
using ClimateMachine.BalanceLaws
using ClimateMachine.ODESolvers
using ClimateMachine.SystemSolvers
using ClimateMachine.Orientations
using ClimateMachine.VTK

# to be removed
using ClimateMachine.GenericCallbacks: EveryXWallTimeSeconds, EveryXSimulationSteps
using ClimateMachine.Thermodynamics: soundspeed_air
using ClimateMachine.VariableTemplates: flattenednames

# to be removed: needed for updating ref state
import ClimateMachine.ODESolvers: update_backward_Euler_solver!
import ClimateMachine.DGMethods: update_auxiliary_state!

# to be removed
using CLIMAParameters#: AbstractEarthParameterSet
struct PlanetParameterSet <: AbstractEarthParameterSet end
get_planet_parameter(p::Symbol) = getproperty(CLIMAParameters.Planet, p)(PlanetParameterSet())

# utils
include("../utilities/operations.jl")
include("../utilities/sphere_utils.jl")
# grid
include("../grid/domains.jl")
include("../grid/grids.jl")
# numerics
include("../models/models.jl")
include("../numerics/filters.jl")
# backends
#include("../balance_law_interface.jl")
include("../esdg_balance_law_interface.jl")
include("../boundary_conditions/boundary_conditions.jl")
include("../numerics/numerical_volume_fluxes.jl")
include("../numerics/numerical_interface_fluxes.jl")
include("../numerics/timestepper_abstractions.jl")
include("../simulations.jl")
include("../utilities/callbacks.jl")
# physics
include("../grid/orientations.jl")
include("../physics/physics.jl")
include("../physics/advection.jl")
include("../physics/thermodynamics.jl")
include("../physics/coriolis.jl")
include("../physics/pressure_force.jl")
include("../physics/gravity.jl")
include("../physics/diffusion.jl")
include("../physics/microphysics.jl")
include("../physics/temperature_profiles.jl")
# coupler
include("../coupler/put_get_calls.jl")
include("../coupler/surface_fluxes.jl")

ClimateMachine.init()

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
const DEFAULT_MODELDIRS = Dict(DEFAULT_MODEL_NAME => "$(homedir())/.config/JustSayIt/models/default/",
                               TYPE_MODEL_NAME    => "$(homedir())/.config/JustSayIt/models/type/")
const DEFAULT_NOISES    = Dict(DEFAULT_MODEL_NAME => NOISES_ENGLISH,
                               TYPE_MODEL_NAME    => NOISES_ENGLISH)
const DEFAULT_COMMANDS  = Dict("help"     => Help.help,
                               "type"     => Keyboard.type,
                               "ma"       => Mouse.click_left,
                               "select"   => Mouse.press_left,
                               "okay"     => Mouse.release_left,
                               "middle"   => Mouse.click_middle,
                               "right"    => Mouse.click_right,
                               "double"   => Mouse.click_double,
                               "triple"   => Mouse.click_triple,
                               "email"    => Email.email,
                               "internet" => Internet.internet)

"""
    just_say_it()
    just_say_it(; <keyword arguments>)

Start offline, low latency, highly accurate speech to command translation.

# Keyword arguments
- `modeldirs::Dict{String, String}=DEFAULT_MODELDIRS`: the directories where the speech recognition models to be used are located.
- `noises::Dict{String, <:AbstractArray{String}}=DEFAULT_NOISES`: for each model, an array of strings with noises (tokens that are to be ignored in the speech as, e.g., "huh").
- `commands::Dict{String, Function}=DEFAULT_COMMANDS`: the commands to be recognized with their mapping to a function.
- `subset::NTuple{N,String}=nothing`: a subset of the `commands` to be recognised and executed (instead of the complete `commands` list).

# Default `modeldirs`
```
$(pretty_dict_string(DEFAULT_MODELDIRS))
```

# Default `noises`
```
$(pretty_dict_string(DEFAULT_NOISES))
```

# Default `commands`
```
$(pretty_dict_string(DEFAULT_COMMANDS))
```

# Examples
```
# Listen to all commands with exception of the mouse click commands.
using JustSayIt
just_say_it(; subset=("help", "type", "email", "internet"))
```

```
# Listen only to the mouse click commands.
using JustSayIt
just_say_it(; subset=("ma", "select", "okay", "middle", "right", "double", "triple"))
```

```
# Define custom modeldirs and commands
using JustSayIt
modeldirs = Dict(DEFAULT_MODEL_NAME => "$(homedir())/.config/JustSayIt/models/default/",
                 TYPE_MODEL_NAME    => "$(homedir())/.config/JustSayIt/models/type/")
commands = Dict("help"     => Help.help,
                "type"     => Keyboard.type,
                "double"   => Mouse.click_double,
                "triple"   => Mouse.click_triple,
                "email"    => Email.email,
                "internet" => Internet.internet)
just_say_it(; modeldirs=modeldirs, commands=commands)
```

"""
function just_say_it(; modeldirs::Dict{String,String}=DEFAULT_MODELDIRS, noises::Dict{String,<:AbstractArray{String}}=DEFAULT_NOISES, commands::Dict{String, Function}=DEFAULT_COMMANDS, subset::Union{Nothing, NTuple{N,String}}=nothing) where N
    if (!isnothing(subset) && !issubset(subset, keys(commands))) error("obtained command name subset ($(subset)) is not a subset of the command names ($(keys(commands))).") end
    if !isnothing(subset) commands = filter(x -> x[1] in subset, commands) end

    # Initializations
    @info "Initializing JustSayIt..."
    init_jsi(modeldirs, noises, commands)
    start_recording()

    # Interprete commands
    @info "Listening for commands..."
    valid_cmd_names = command_names()
    cmd_name = ""
    is_sleeping = false
    try
        while true
            if is_sleeping
                if cmd_name == COMMAND_NAME_AWAKE
                    if is_confirmed() is_sleeping = false end
                    if (is_sleeping) @info("I heard \"awake\". If you want to awake me, say \"awake JustSayIt\".")
                    else             @info("... awake again. Listening for commands...")
                    end
                end
            else
                if cmd_name in valid_cmd_names
                    cmd = command(string(cmd_name))
                    if (t0_latency() > 0.0) @debug "Latency of command `$cmd_name`: $(round(toc(t0_latency()),sigdigits=2))." end
                    try
                        @info "Starting command: $cmd (latency: $(round(Int,toc(t0_latency())*1000)) ms)"
                        cmd()
                    catch e
                        if isa(e, InsecureRecognitionException)
                            @info("Command `$cmd_name` aborted: insecure command argument recognition.")
                        else
                            rethrow(e)
                        end
                    end
                elseif cmd_name == COMMAND_NAME_SLEEP
                    if is_confirmed() is_sleeping = true end
                    if (!is_sleeping) @info("I heard \"sleep\". If you want me to sleep, say \"sleep JustSayIt\".")
                    else              @info("Going to sleep... (To awake me, just say \"awake JustSayIt\".)")
                    end
                elseif cmd_name != ""
                    @info "Invalid command: $cmd_name." # NOTE: this might better go to stderr or @debug later.
                end
            end
            try
                force_reset_previous(recognizer(COMMAND_RECOGNIZER_ID))
                cmd_name = next_token(recognizer(COMMAND_RECOGNIZER_ID), _noises(DEFAULT_MODEL_NAME); use_partial_recognitions = true)
            catch e
                if isa(e, InsecureRecognitionException)
                    @info(e.msg)
                    cmd_name = ""
                else
                    rethrow(e)
                end
            end
        end
    finally
        stop_recording()
        @info "Stopped listening for commands."
    end

    finalize_jsi()
end

@voiceargs words=>(valid_input=["just say it"], use_max_accuracy=true, vararg_timeout=2.0, vararg_max=3) function is_confirmed(words::String...)
    return (join(words, " ") == "just say it")
end


text:
import IRTools
import IRTools: IR, @dynamo, self, insertafter!
import UUIDs: UUID, uuid1

function module_functions(modl)
    res = Vector{Function}()
    for s in Base.names(modl; all=true)
        isdefined(modl, s) || continue
        fn = getfield(modl, s)
        # && match(r"^[a-z#]+$", string(s)) != nothing
        if fn isa Function && !in(string(fn)[1], "#@")
            push!(res, fn)
        end
    end
    return res
end

# include("primitives.jl")

const BASE_PRIMITIVE_FUNCTIONS = vcat(
    module_functions(Base),
    module_functions(Core),
    module_functions(Core.Intrinsics),
    [Broadcast.materialize, Broadcast.broadcasted, Colon(), (:),
     Base.not_int,
     # our own special functions
     __new__]);

const PRIMITIVES = FunctionResolver{Bool}(
    collect(Tuple{typeof(f), Vararg} => true for f in BASE_PRIMITIVE_FUNCTIONS)
)

"""
    is_primitive(sig)

The default implementation of `is_primitive` argument in [`trace()`](@ref).
Returns `true` if the method with the provided signature is defined
in one of the Julia's built-in modules, e.g. `Base`, `Core`, `Broadcast`, etc.
"""
function is_primitive(sig)
    return sig in PRIMITIVES
end

################################################################################
#                                Frame                                         #
################################################################################

"""Frame of a call stack"""
mutable struct Frame
    # IR var ID to Tape var. Note that IR ID refers to the _original_ IR,
    # not the transformed one. In fact, any unique name of SSA instructions
    # would fit, using IR IDs is just convenient.
    # Note that a constant value can be passed instead of a Tape var
    ir2tape::Dict{Int, Any}
    # result ID - tape ID corresponding to latest return value
    # from this call frame
    result::V
    # function and arguments - used purely for debugging
    fargs::Tuple
end

function Base.show(io::IO, fr::Frame)
    map_str = join(["$k=>$v" for (k, v) in fr.ir2tape], ",")
    print(io, "Frame($map_str, $(fr.result))")
end

################################################################################
#                              Tracer Options                                  #
################################################################################

const TRACING_OPTIONS = Ref(Dict())

"""
    should_trace_loops!(val=false)

Turn on/off loop tracing. Without parameters, resets the flag to the default value
"""
should_trace_loops!(val::Bool=false) = (TRACING_OPTIONS[][:trace_loops] = val)

"""
    should_trace_loops()

Check the current value of the loop tracing option.
"""
should_trace_loops() = get(TRACING_OPTIONS[], :trace_loops, false)

"""
    should_assert_branches!(val=false)

Turn on/off loop tracing. Without parameters, resets the flag to the default value
"""
should_assert_branches!(val::Bool=false) = (TRACING_OPTIONS[][:assert_branches] = val)

"""
    should_assert_branches()

Check the current value of the loop tracing option.
"""
should_assert_branches() = get(TRACING_OPTIONS[], :assert_branches, false)

"""
Tracer options. Configured globally via the following methods:

- should_trace_loops!()
- should_assert_branches!()
"""
struct TracerOptions
    trace_loops::Bool
    assert_branches::Bool
end

TracerOptions() = TracerOptions(should_trace_loops(),
                                should_assert_branches())

################################################################################
#                        IRTracer (defs and utils)                             #
################################################################################

mutable struct IRTracer
    is_primitive::Function
    tape::Tape
    frames::Vector{Frame}
    options::TracerOptions
    last_block_jmp_target::Union{V, Nothing}
end

function IRTracer(; ctx=Dict(), is_primitive=is_chainrules_primitive)
    tape = Tape(ctx)
    return IRTracer(is_primitive, tape, [], TracerOptions(), nothing)
end

Base.show(io::IO, t::IRTracer) = print(io, "IRTracer($(length(t.tape)))")

promote_const_value(x::QuoteNode) = x.value
promote_const_value(x::GlobalRef) = getproperty(x.mod, x.name)
promote_const_value(x) = x

function get_tape_vars(t::IRTracer, farg_irvars::Union{Vector, Tuple})
    result = Vector{Any}(undef, length(farg_irvars))
    for (i, arg) in enumerate(farg_irvars)
        if arg isa IRTools.Variable
            # cat be either tape var or constant value if a function
            # wall called with a constant
            x = t.frames[end].ir2tape[arg.id]
            result[i] = x isa V ? x : promote_const_value(x)
        else
            result[i] = promote_const_value(arg)
        end
    end
    return result
end

"""Push a new call frame to tracer, setting function params accordingly"""
function push_frame!(t::IRTracer, fargs, farg_irvars...)
    tape_vars = get_tape_vars(t, farg_irvars)
    frame = Frame(
        Dict(i => v for (i, v) in enumerate(tape_vars)),
        V(0),
        fargs,
    )
    push!(t.frames, frame)
end

"""Pop call frame from tracer"""
function pop_frame!(t::IRTracer, res_sid::Int)
    frame = pop!(t.frames)
    # create mapping from the current SSA ID to the last instruction on the tape
    t.frames[end].ir2tape[res_sid] =
        (frame.result.id == 0 ? bound(t.tape, V(length(t.tape))) : frame.result)
end

"""Set target branch parameters to variables corresponding to SSA args"""
function set_branch_params!(t::IRTracer, ssa_args, target_params)
    tape_vars = get_tape_vars(t, ssa_args)
    ir2tape = t.frames[end].ir2tape
    for (v, p) in zip(tape_vars, target_params)
        ir2tape[p] = v
    end
end

"""Set return variable for the current frame"""
function set_return!(t::IRTracer, arg_sid_ref)
    tape_var = get_tape_vars(t, [arg_sid_ref[]])[1]
    if tape_var isa V
        t.frames[end].result = tape_var
    else
        # if the function returns a constant (e.g. nothing), record it
        res = push!(t.tape, Constant(tape_var))
        t.frames[end].result = res
    end
end

function check_and_set_branch!(t::IRTracer, condition, block)
    condition = get_tape_vars(t, [condition[]])[1]
    t.last_block_jmp_target = push!(t.tape, mkcall(_check_and_set_branch!, condition, block, t.last_block_jmp_target))
end
function _check_and_set_branch!(condition, block, last_block_jmp_target)
    if last_block_jmp_target == nothing && (condition == nothing || condition == false)
        return block
    else
        return last_block_jmp_target
    end
end

function check_block(t::IRTracer, block_id)
    t.last_block_jmp_target = push!(t.tape, mkcall(_check_block, block_id, t.last_block_jmp_target))
end
function _check_block(block_id, next_block)
    @assert next_block !== nothing ? next_block == block_id : true
    return nothing
end

################################################################################
#                                  Loop utils                                  #
################################################################################

function is_loop(block::IRTools.Block)
    for br in IRTools.branches(block)
        # if a branch refers to an earlier block and is not return
        # then it must be a loop
        if br.block <= block.id && br.block != 0
            return true
        end
    end
    return false
end

"""Get SSA IDs of a block's inputs"""
function block_input_ir_ids(block::IRTools.Block)
    result = []
    for (stmt_id, (block_id, arg_id)) in enumerate(block.ir.defs)
        if block_id == block.id && arg_id < 0
            push!(result, stmt_id)
        end
    end
    return result
end

"""
Get SSA IDs of arguments which are not part of this block
(e.g. come from outside of a loop)
"""
function block_outsider_ir_ids(block::IRTools.Block)
    result = []
    min_id = minimum(block_input_ir_ids(block))
    for (_, stmt) in block
        ex = stmt.expr
        @assert Meta.isexpr(ex, :call)
        for arg in ex.args[2:end]
            if arg isa IRTools.Variable && arg.id < min_id
                push!(result, arg.id)
            end
        end
    end
    return result
end

function loop_exit_branch(block::IRTools.Block)
    branches = IRTools.branches(block)
    target_branch_id = findfirst([br.block > block.id for br in branches])
    return target_branch_id !== nothing ? branches[target_branch_id] : nothing
end

function loop_continue_branch(block::IRTools.Block)
    branches = IRTools.branches(block)
    target_branch_id = findfirst([br.block <= block.id for br in branches])
    return return target_branch_id !== nothing ? branches[target_branch_id] : nothing
end

function loop_condition_ir_id(block::IRTools.Block)
    br = loop_exit_branch(block)
    return br !== nothing ? br.condition.id : nothing
end

"""Get SSA IDs of exit arguments of a loop block"""
function loop_exit_ir_ids(block::IRTools.Block)
    br = loop_exit_branch(block)
    return br !== nothing ? [arg.id for arg in br.args] : nothing
end

function loop_continue_ir_ids(block::IRTools.Block)
    br = loop_continue_branch(block)
    return br !== nothing ? [arg.id for arg in br.args] : nothing
end

"""Pseudo op to designate loop end. Removed after Loop op is created"""
mutable struct _LoopEnd <: AbstractOp
    id::Int
    tape
end
_LoopEnd() = _LoopEnd(0, nothing)

const LOOP_EXIT_TAPE_IDS = "loop_exit_tape_ids"
const LOOP_COND_ID = "loop_cond_id"
const LOOP_CONTINUE_TAPE_IDS = "loop_continue_tape_ids"

is_loop_traced(t::IRTracer, loop_id::UUID) = haskey(t.tape.meta, "loop_already_traced_$loop_id")
loop_traced!(t::IRTracer, loop_id::UUID) = (t.tape.meta["loop_already_traced_$loop_id"] = true)

"""
Trigger loop start operations.

Arguments:
----------

 * t :: IRTracer
    Current tracer
 * loop_id :: Int
    Unique ID of a loop being entered
 * loop_input_ir_ids :: Vector{Int}
    IR IDs of variables which will be used as loop inputs. Includes
    loop block inputs and any outside IDs

This function is added to the very beginning of the loop block(s).
During the first iteration we initialize a subtape which will be used
later to create the Loop operation on the parent tape. Since all iterations
of the loop produce identical operations, we only need to trace it once.
However, it turns out to be easier to record all iterations (seprated by
a special _LoopEnd op) and then prune unused iterations.

Another important detail is that Loop's subtape is a completely valid
and independent tape with its own call frame and inputs which include
all explicit and implicit inputs to the loop's block in the original IR.
"""
function enter_loop!(t::IRTracer, loop_id, loop_input_ir_ids::Vector)
    t.options.trace_loops || return
    # skip if it's not the first iteration
    is_loop_traced(t, loop_id) && return
    # create subtape, with the current tape as parent
    C = typeof(t.tape.c)
    subtape = Tape(C())
    subtape.parent = t.tape
    # create a new frame and push to the list
    frame = Frame(Dict(), V(0), ())
    push!(t.frames, frame)
    # record inputs to the subtape & populate the new frame's ir2tape
    for ir_id in loop_input_ir_ids
        parent_tape_var = t.frames[end - 1].ir2tape[ir_id]
        val = subtape.parent[parent_tape_var].val
        tape_var = push!(subtape, Input(val))
        t.frames[end].ir2tape[ir_id] = tape_var
    end
    # replace IRTracer.tape with subtape
    t.tape = subtape
end

"""
Set flags designating the end of the first run of the loop.
"""
function stop_loop_tracing!(t::IRTracer,
                            loop_id::Any,
                            input_ir_ids::Vector,
                            cond_ir_id::Any,
                            cont_ir_ids::Vector,
                            exit_ir_ids::Vector,
                            exit_target_ir_ids::Vector)
    t.options.trace_loops || return
    if !is_loop_traced(t, loop_id)
        # record exit tape IDs as of first iteration
        # we will use them later
        t.tape.meta[LOOP_EXIT_TAPE_IDS] =
            [t.frames[end].ir2tape[ir_id] for ir_id in exit_ir_ids]
        t.tape.meta[LOOP_COND_ID] = t.frames[end].ir2tape[cond_ir_id]
        t.tape.meta[LOOP_CONTINUE_TAPE_IDS] =
            [t.frames[end].ir2tape[ir_id] for ir_id in cont_ir_ids]
        # set flag to stop creatnig new subtapes
        loop_traced!(t, loop_id)
    end
    # record a special op to designate the end of the loop code
    # tracer will continue to record ops, but later we truncate
    # the tape to get only ops before _LoopEnd
    push!(t.tape, _LoopEnd())
end

"""
Trigget loop end operations.

This function is added just before the end of the loop block.
Since we record all iterations of the loop, we must remember tape IDs
of continuation condition and exit variables during the first run.
"""
function exit_loop!(t::IRTracer,
                    input_ir_ids::Vector,
                    cond_ir_id::Any,
                    cont_ir_ids::Vector,
                    exit_ir_ids::Vector,
                    exit_target_ir_ids::Vector)
    t.options.trace_loops || return
    # loop subtape already contains a variable designating condition
    # of loop continuation; if this condition is false,
    # we are ready to exit the loop and record Loop operation
    cond_op = t.tape[t.frames[end].ir2tape[cond_ir_id]]
    if !cond_op.val
        # swap active tape back
        pop!(t.frames)
        parent_ir2tape = t.frames[end].ir2tape
        subtape = t.tape
        t.tape = t.tape.parent
        # remove repeating blocks
        first_loop_end = findfirst(op -> isa(op, _LoopEnd), subtape.ops)
        subtape.ops = subtape.ops[1:first_loop_end-1]
        # record output tuple
        exit_tape_vars = subtape.meta[LOOP_EXIT_TAPE_IDS]
        exit_val = tuple([subtape[id].val for id in exit_tape_vars]...)
        # exit_var = push!(subtape, mkcall(tuple, exit_tape_vars...))
        # exit_val = subtape[exit_var].val
        subtape.result = bound(t.tape, V(length(t.tape)))   # dummy, not used in practice
        # record the loop operation
        parent_input_vars = [parent_ir2tape[ir_id] for ir_id in input_ir_ids]
        condition = subtape.meta[LOOP_COND_ID]
        cont_vars = subtape.meta[LOOP_CONTINUE_TAPE_IDS]
        loop_id = push!(t.tape, Loop(0, parent_input_vars,
                          condition, cont_vars,
                          exit_tape_vars, subtape, exit_val))
        # destructure loop return values to separate vars on the main tape
        # and map branch arguments to these vars
        for i=1:length(exit_val)
            res_id = push!(t.tape, mkcall(getfield, loop_id, i))
            parent_ir2tape[exit_target_ir_ids[i]] = res_id
        end
    end
end

################################################################################
#                        IRTracer (body) + irtrace()                           #
################################################################################

function rewrite_special_cases!(ir::IR)
    for (v, st) in ir
        if Meta.isexpr(st.expr, :new)
            ir[v] = Expr(:call, __new__, st.expr.args...)
        end
    end
end

"""
Record function call onto a tape or recurse into it.

Params:
-------
* t::IRTracer - current tracer
* res_id::Int - IR ID of the operation
* farg_irvars - IR variables of the operation
* fargs - values of the operation
"""
function record_or_recurse!(t::IRTracer, res_sid::Int, farg_irvars, fargs...)
    fn, args = fargs[1], fargs[2:end]
    if t.is_primitive(Tuple{map(typeof, fargs)...})
        tape_vars = get_tape_vars(t, farg_irvars)
        # record corresponding op to the tape
        res = push!(t.tape, mkcall(tape_vars...))
        # update mapping from SSA var to tape var
        # note that in functions with loops this mapping may change over time
        t.frames[end].ir2tape[res_sid] = res
        val = t.tape[res].val
    else
        push_frame!(t, fargs, farg_irvars...)
        val = t(fn, args...)
        pop_frame!(t, res_sid)
    end
    return val
end

function record_const!(t::IRTracer, res_sid, val)
    val = val isa QuoteNode ? val.value : val
    res = push!(t.tape, Constant(val))
    t.frames[end].ir2tape[res_sid] = res
    return val
end

"""
Get SSA IDs of the branch target parameters.
For example, given code like this:

  2: (%9, %10, %11)
    ...
    br 3 (%14, %15) unless %18
    br 2 (%16, %14, %15)
  3: (%19, %20)
    ...

This function will return:

  branch_target_params(ir, <br 3>) ==> [19, 20]

"""
function branch_target_params(ir:: IR, branch::IRTools.Branch)
    return [v.id for v in ir.blocks[branch.block].args]
end

function trace_branches!(ir::IR)
    # if a block ends with a branch, we map its parameters to tape IDs
    # which currently correspond to argument SSA IDs
    for block in IRTools.blocks(ir)
        if should_assert_branches()
            pushfirst!(block, Expr(:call, check_block, self, block.id))
        end
        for branch in IRTools.branches(block)
            if should_assert_branches()
                push!(block, Expr(:call, check_and_set_branch!, self, Ref(branch.condition), branch.block))
            end
            if IRTools.isreturn(branch)
                ret_v = branch.args[1]
                if should_assert_branches()
                    # Returns are equivalent to branches to 0 in IRTools
                    push!(block, Expr(:call, check_block, self, 0))
                end
                push!(block, Expr(:call, set_return!, self, Ref(ret_v)))
            else
                ssa_args = branch.args
                target_params = [v.id for v in ir.blocks[branch.block].args]
                push!(block, Expr(:call, set_branch_params!, self, ssa_args, target_params))
            end
        end
    end
end

function loop_start_end_blocks(ir::IR, block::IRTools.Block)
    start_block = IRTools.block(ir, loop_continue_branch(block).block)
    exit_branch = loop_exit_branch(block)
    # check if this block contains the exit branch
    if exit_branch !== nothing && exit_branch.block > block.id
        end_block = block
    else
        exit_branch = loop_exit_branch(start_block)
        if exit_branch !== nothing && exit_branch.block > block.id
            end_block = start_block
        else
            error("Cannot find end block of a loop")
        end
    end
    return start_block, end_block
end

function trace_loops!(ir::IR)
    for block in IRTools.blocks(ir)
        if is_loop(block)
            loop_id = uuid1()   # unique ID of this loop
            start_block, end_block = loop_start_end_blocks(ir, block)
            # loop start - the first block of the loop
            loop_input_ir_ids = vcat(
                block_input_ir_ids(start_block),
                # block_outsider_ir_ids(start_block),
            )
            pushfirst!(start_block, Expr(:call, enter_loop!, self, loop_id, loop_input_ir_ids))
            # loop tracing border - at this point all operations of the loop
            # have been executed at least once, even if continuation condition
            # is in another block
            push!(block, Expr(:call, stop_loop_tracing!, self,
                              loop_id,
                              loop_input_ir_ids,
                              loop_condition_ir_id(end_block),
                              loop_continue_ir_ids(block),
                              loop_exit_ir_ids(end_block),
                              branch_target_params(ir, loop_exit_branch(end_block))))
            # loop end - continuation condition is checked here
            push!(end_block, Expr(:call, exit_loop!, self,
                              loop_input_ir_ids,
                              loop_condition_ir_id(end_block),
                              loop_continue_ir_ids(block),
                              loop_exit_ir_ids(end_block),
                              branch_target_params(ir, loop_exit_branch(end_block))))
        end
    end
end

function expand!(ir::IR)
    # TODO delete this function
    # Variables orders may be switched in IRTools.expand!
    # depending on details of the julia version/environment.
    # This bandaid forces a uniform order
    # See also https://github.com/dfdx/Ghost.jl/pull/7
    worklist = IRTools.blocks(ir)
    spats = Dict(b => Dict() for b in IRTools.blocks(ir))
    while !isempty(worklist)
      b = pop!(worklist)
      b.id == 1 && continue
      defs = IRTools.definitions(b)
      uses = IRTools.usages(b)
      uses = sort(collect(uses), by=v->v.id, rev=true)
      for v in setdiff(uses, defs)
        haskey(spats[b], v) && continue
        spats[b][v] = IRTools.argument!(b, v)
        for c in IRTools.predecessors(b)
          c in worklist || push!(worklist, c)
        end
      end
      ir.blocks[b.id] = IRTools.prewalk(x -> get(spats[b], x, x), ir.blocks[b.id])
    end
    return ir
  end

@dynamo function (t::IRTracer)(fargs...)
    ir = IR(fargs...)
    ir === nothing && return   # intrinsic functions
    expand!(ir)
    rewrite_special_cases!(ir)
    for (v, st) in ir
        ex = st.expr
        # note the difference:
        # * `ex.args` is an array and thus will be passed to a function as is,
        # including definitions of IRTools.Variable;
        # * `ex.args...` is top-level to this expression and thus all Variable's
        # will be replaced with actual values during runtime
        if Meta.isexpr(ex, :call)
            ir[v] = IRTools.xcall(record_or_recurse!, self, v.id, ex.args, ex.args...)
        else
            # e.g. GlobalRef
            # note: using insertafter!() due to
            # https://github.com/FluxML/IRTools.jl/issues/78
            # ir[v] = Expr(:call, record_const!, self, v.id, v)
            insertafter!(ir, v, IRTools.xcall(record_const!, self, v.id, v))
        end
    end
    trace_branches!(ir)
    trace_loops!(ir)
    return ir
end

"""
    trace(f, args...; is_primitive, primitives)

Trace function call, produce call value and a Tape.

`trace` records to the tape primitive methods and recursively dives into
non-primitives. There are 2 ways to tell `trace` that a particular method
is a primitive:

* provide `is_primitive(sig) -> Bool` function, where `sig` is
    is a method signature, e.g. `Tuple{map(typeof, (f, args...))...}`
* provide an iterable `primitives`; in this case `trace` matches
    all methods of this function
"""
function trace(f, args...; is_primitive=is_primitive, primitives=nothing, ctx=Dict())
    if primitives !== nothing
        sigs = FunctionResolver{Bool}([Tuple{typeof(f), Vararg} => true for f in primitives])
        is_primitive = sig -> sig in sigs
    end
    t = IRTracer(; ctx=ctx, is_primitive=is_primitive)
    arg_vars = inputs!(t.tape, f, args...)
    frame = Frame(
        Dict(i => a for (i, a) in enumerate(arg_vars)), V(0), (f, args...)
    )
    push!(t.frames, frame)
    val = t(f, args...)
    t.tape.result = t.frames[1].result
    tape = t.tape
    return val, tape
end


text:
function digitalroot(n::Integer, bs::Integer=10)
    if n < 0 || bs < 2 throw(DomainError()) end
    ds, pers = n, 0
    while bs ≤ ds
        ds = sum(digits(ds, bs))
        pers += 1
    end
    return pers, ds
end

for i in [627615, 39390, 588225, 393900588225, big(2) ^ 100]
    pers, ds = digitalroot(i)
    println(i, " has persistence ", pers, " and digital root ", ds)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    CylinderNeighborhood{N}(radius, height)

A `N`-dimensional cylinder neighborhood with `radius` and `height`.

### Notes

The `height` parameter is only half of the actual cylinder height.
"""
struct CylinderNeighborhood{T,N} <: AbstractNeighborhood{T,N}
  radius::T
  height::T
end

function CylinderNeighborhood{N}(radius::T, height::T) where {N,T}
  @assert radius > 0 "cylinder radius must be positive"
  @assert height > 0 "cylinder height must be positive"
  CylinderNeighborhood{T,N}(radius, height)
end

"""
    radius(cylinder)

Return the radius of the `cylinder`.
"""
radius(cylinder::CylinderNeighborhood) = cylinder.radius

"""
    height(cylinder)

Return (half of) the height of the `cylinder`.
"""
height(cylinder::CylinderNeighborhood) = cylinder.height

function isneighbor(cylinder::CylinderNeighborhood, xₒ::AbstractVector, x::AbstractVector)
  l = length(xₒ)
  @inbounds abs(x[l] - xₒ[l]) ≤ cylinder.height && norm(view(x,1:l-1) - view(xₒ,1:l-1)) ≤ cylinder.radius
end

# ------------
# IO methods
# ------------
function Base.show(io::IO, cylinder::CylinderNeighborhood{T,N}) where {N,T}
  r = cylinder.radius
  h = cylinder.height
  print(io, "$(N)D CylinderNeighborhood($r, $h)")
end

function Base.show(io::IO, ::MIME"text/plain",
                   cylinder::CylinderNeighborhood{T,N}) where {N,T}
  println(io, "$(N)D CylinderNeighborhood")
  println(io, "  radius: ", cylinder.radius)
  print(  io, "  height: ", cylinder.height)
end


text:
# small scrip to compute the solution of Lipman Schinwer equation
# We test the types introduced in FastConvolution.jl
# we test that the application is fast and that the construction
# is performed fast.

using PyPlot
using IterativeSolvers
using SpecialFunctions

@everywhere include("../src/FastConvolution.jl")
@everywhere include("../src/Preconditioner.jl")

FFTW.set_num_threads(4);
BLAS.set_num_threads(4)

#Defining Omega
h = 0.001
k = 1/h

# size of box
a  = 1
x = -a/2:h:a/2
y = -a/2:h:a/2
(n,m) = length(x), length(y)
N = n*m
X = repmat(x, 1, m)[:]
Y = repmat(y', n,1)[:]
# we solve \triangle u + k^2(1 + nu(x))u = 0

# We use the modified quadrature in Ruan and Rohklin
(ppw,D) = referenceValsTrapRule();
D0 = D[1];

spline(y) = (y.<0) +  (y.>=0).*(y.<1).*(2*y.^3 - 3*y.^2 + 1)
window(a1,b1,b2,a2,y) = (y.>=a1).*( (y.<b1).*spline(1/(abs(b1-a1)).*(-y+b1)) + (y.>=b1).*(y.<b2) + (y.>=b2).*spline(1/(abs(b2-a2)).*(y-b2)))

# Defining the smooth perturbation of the slowness
nu(x,y) = 0.05*sin(2*pi*x).*window(-0.48,-0.4,0.4,0.48,y).*window(-0.48,-0.05,0.05,0.48,x);

figure(1); clf();
imshow(real(reshape( (1-nu(X,Y)),n,m)), extent=[x[1], x[end], y[1], y[end]]);cb =  colorbar();

# Sampling the Green's function for the Toeplitz form
Ge = buildGConv(x,y,h,n,m,D0,k);
GFFT = fft(Ge);

fastconv = FastM(GFFT,nu(X,Y),2*n-1,2*m-1,n, m, k);

# assembling the sparsifiying preconditioner
As = buildSparseA(k,X,Y,D0, n ,m);

# assembling As*( I + k^2G*nu)
Mapproxsp = As + k^2*(buildSparseAG(k,X,Y,D0, n ,m)*spdiagm(nu(X,Y)));

# defining the preconditioner
precond = SparsifyingPreconditioner(Mapproxsp, As)

# building the RHS from the incident field
u_plane = exp(k*im*sqrt(1-nu(0*Y,X)).*X);

u_inc = exp(k*im*X);

rhs = -(fastconv*u_inc - u_inc);

# allocating the solution
u = zeros(Complex128,N);

# solving the system using GMRES
@time info =  gmres!(u, fastconv, rhs, Pl = precond)
# println(info[2].residuals[:])

# plotting the solution
figure(2)
clf()
imshow(real(reshape(u+u_inc,n,m)))

figure(3); clf()
imshow(real(reshape((u+u_inc)./u_inc, n,m))); colorbar()
title("\$u(x)/e^{ikx}\$ ")

intx = cumsum(sqrt(1-nu(x,0*x)))*h;
u_gem = exp(k*im*intx);

UUgem = repmat(u_gem, 1,n)

figure(4); clf()
imshow(real(reshape((u+u_inc)./UUgem[:], n,m))); colorbar()
title("\$u(x)/e^{ik \\phi (x)}\$ where \$ \\phi \$ comes from \$ \\int_{0}^{x} \\sqrt{1 + n(y)} dy \$ ")

refFact = 10;
xRef = -a/2:h/refFact:a/2

intxRef = cumsum(sqrt(1-nu(xRef,0*xRef)))*h/refFact;
u_gem = exp(k*im*intxRef)[1:refFact:end];

UUgem = repmat(u_gem, 1,n)

figure(14); clf()
imshow(real(reshape(-(u+u_inc)./UUgem[:], n,m))); colorbar()
title("\$u(x)/e^{ik \\phi (x)}\$ where \$ \\phi \$ comes from \$ \\int_{0}^{x} \\sqrt{1 + n(y)} dy \$ ")

UU = reshape(u+u_inc,n,m)
ucomp = UU[:,round(Integer,floor(end/2))]

# there is a two that had to disappear... I need to find the bug... :/
phi_1 = -1im*rhs./u_inc

figure(5); clf()
imshow(real(reshape(phi_1, n,m))); colorbar()
title("phase function \$phi_1(x)\$ ")

u_rytov = u_inc.*exp(1im*phi_1)

figure(6); clf()
imshow(real(reshape(u_rytov, n,m))); colorbar()
title("\$u_{rytov}\$" )

figure(7); clf()
imshow(real(reshape((u+u_inc)./ u_rytov, n,m))); colorbar()
title("\$u(x)/u_{rytov}(x)\$ ")

#Defining the smooth perturbation of the slowness
nu1(x,y) = nu(x,y).*window(-0.28,-0.10,0.10,0.28,y);

figure(24); clf();
imshow(real(reshape( (1-nu1(X,Y)),n,m)), extent=[x[1], x[end], y[1], y[end]]);cb =  colorbar();

fastconvWindow = FastM(GFFT,nu1(X,Y),2*n-1,2*m-1,n, m, k);

rhsWindow = -(fastconvWindow*u_inc - u_inc);

# there is a two that had to disappear... I need to find the bug... :/
phi_Window = -1im*rhsWindow./u_inc

figure(15); clf()
imshow(real(reshape(phi_1, n,m))); colorbar()
title("phase function phi_1(x) ")

u_rytovWindow = u_inc.*exp(1im*phi_Window)

figure(16); clf()
imshow(real(reshape(u_rytovWindow, n,m))); colorbar()
title("u_rytov ")

figure(17); clf()
imshow(real(reshape((u+u_inc)./ u_rytovWindow, n,m))); colorbar()
title("u/u_rytov ")

UU = reshape((u+u_inc)./u_rytov,n,m)
ucomp = UU[:,round(Integer,floor(end/2))]

figure(20); clf();
plot(real(ucomp[:]))

# u_inc = exp(k*im*X);

# rhs1 = -(fastconv1*u_inc - u_inc);

# # allocating the solution
# u1 = zeros(Complex128,N);

# # solving the system using GMRES
# @time info =  gmres!(u1, fastconv1, rhs1, precond1)
# println(info[2].residuals[:])

# # plotting the solution
# figure(5)
# clf()
# imshow(real(reshape(u1+u_inc,n,m)))

# figure(6); clf()
# imshow(real(reshape((u1+u_inc)./u_inc, n,m))); colorbar()

max_stars_repo_name: Forgotten/Lippmann-Schingwer

max_issues_repo_name: Forgotten/Lippmann-Schingwer

max_forks_repo_name: Forgotten/Lippmann-Schingwer


text:
using ArgParse
include("../specfem_gll.jl/tasks/numerical_operation_on_bin_file.jl")
include("../specfem_gll.jl/src/utils/readfiles.jl")

function parse_command_line()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "--target_basedir"
            help = "the target gll directory"
        "--reference_basedir"
            help = "the reference gll directory"
        "--mesh_basedir"
            help = "the mesh directory"
        "--output_basedir"
            help = "the output gll directory for perturbation"
        "--tags"
            help = "the tags, eg: alpha,beta"
        "--nproc"
            help = "number of processors the gll directory correspons to"
    end
    return parse_args(s)
    end

function main()
    # parse args
    parsed_args = parse_command_line()
    target_basedir = parsed_args["target_basedir"]
    mesh_basedir = parsed_args["mesh_basedir"]
    reference_basedir = parsed_args["reference_basedir"]
    output_basedir = parsed_args["output_basedir"]
    tags = parsed_args["tags"]
    nproc = parse(Int64, parsed_args["nproc"])
    # get nspec
    mesh_info = sem_mesh_read(mesh_basedir, 0)
    nspec = mesh_info.nspec
    # run generate_perturbation
    generate_difference(target_basedir, reference_basedir, output_basedir, tags, nproc, nspec)
end

main()

max_stars_repo_name: ziyixi/seisflow

max_issues_repo_name: ziyixi/seisflow

max_forks_repo_name: ziyixi/seisflow


text:
include("ycreator_initiate.jl")
include("ycreator_create.jl")

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
function _checkfunc(x, xs)
    x in [:(<), :(<=), :(>), :(>=)] && error("Only Boolean combinations of equalities allowed, got $x")
    :($x($(xs...)))
end

function _checksym(x, lc)
    issymbol(x) && x != lc && return :($x($lc))
    x
end

function preprocess_invariant(expr::Expr, lc::Symbol)
    atom_walk(expr) do ex
        @match ex begin
            x_ && y_ => :($x & $y)
            x_ || y_ => :($x | $y)
            x_ != 0  => :(~($x == 0))
            0 != x_  => :(~($x == 0))
            x_ != y_ => :(~($x-$y == 0))
            x_ == 0  => :($x == 0)
            0 == x_  => :($x == 0)
            x_ == y_ => :($x-$y == 0)
            x_(xs__) => _checkfunc(x, xs)
            _        => _checksym(ex, lc)
        end
    end
end

function check_loop_counter(expr::Expr)
    lc = nothing
    function_walk(expr) do _, args
        length(args) != 1 && error("More than one argument, got $args")
        symbol_walk(args[1]) do s
            !isnothing(lc) && lc != s && error("More than one loop counter found, got ($lc, $s)")
            lc = s
        end
    end
    lc
end

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

struct Invariant
    x::Expr
    lc::Symbol
    function Invariant(x::Expr)
        lc = check_loop_counter(x)
        unique_lc = gensym_unhashed(isnothing(lc) ? :n : lc)
        x = replace(x, lc, unique_lc)
        new(preprocess_invariant(x, unique_lc), unique_lc)
    end
end

macro invariant(x)
    Invariant(x)
end

function program_variables(i::Invariant)
    ls = Symbol[]
    postwalk(i.x) do x
        if @capture(x, f_(_)) && issymbol(f)
            push!(ls, f)
        end
        x
    end
    unique(ls)
end

"Return true if the invariant is a relation among a single loop iteration."
function isloopinvariant(i::Invariant)
    flag = true
    function_walk(i) do _, args
        # assume there is just one argument (should be checked on construction)
        if args[1] != i.lc && args[1] != 0
            flag = false
        end
    end
    flag
end

function loopinvariant(i::Invariant)
    @assert isloopinvariant(i)
    function_walk(i) do f, args
        # assume there is just one argument (should be checked on construction)
        args[1] == 0 ? init_var(f) : f
    end
end

constraint_walk(f, i::Invariant) = constraint_walk(f, i.x)
function_walk(f, i::Invariant) = function_walk(f, i.x)

Base.show(io::IO, i::Invariant) = print(io, string(i.x))


text:
module ResizableArraysTests

using Test
using ResizableArrays
using ResizableArrays: checkdimension, checkdimensions, _same_elements
using Base: unsafe_convert, elsize

# FIXME: used @generated
slice(A::AbstractArray{<:Any,2}, I) = A[:,I]
slice(A::AbstractArray{<:Any,3}, I) = A[:,:,I]
slice(A::AbstractArray{<:Any,4}, I) = A[:,:,:,I]
slice(A::AbstractArray{<:Any,5}, I) = A[:,:,:,:,I]

sum_v1(iter::AbstractArray) = (s = zero(eltype(iter));
                               for x in iter; s += x; end;
                               return s)
sum_v2(iter::AbstractArray) = (s = zero(eltype(iter));
                               @inbounds for x in iter; s += x; end;
                               return s)
sum_v3(iter::AbstractArray) = (s = zero(eltype(iter));
                               @inbounds @simd for x in iter; s += x; end;
                               return s)

@testset "Basic methods" begin
    @testset "Utilities" begin
        @test checkdimension(Bool, π) == false
        @test checkdimensions(Bool, ()) == true
        @test checkdimensions(Bool, (1,)) == true
        @test checkdimensions(Bool, (1,2,0)) == true
        @test checkdimensions(Bool, (1,-2,0)) == false
        @test_throws ErrorException checkdimensions((1,-2,0))
        @test isgrowable(π) == false
        @test isgrowable((1,2,3)) == false
        @test isgrowable([1,2,3]) == true
        for T in (Int16, Float32, Tuple{Float64,Float64})
            @test elsize(ResizableArray{T,3,Vector{T}}) == sizeof(T)
        end

        # Make sure all variants of _same_elements are tested.
        A = randn(3,4)
        indexstyles = (IndexLinear(), IndexCartesian())
        for indexstyle1 in indexstyles, indexstyle2 in indexstyles
            @test _same_elements(indexstyle1, A, indexstyle2, A, length(A))
        end

    end
    @testset "Dimensions: $dims" for dims in ((), (3,), (2,3), (2,3,4))
        altdims = map(UInt, dims) # used later
        N = length(dims)
        if N > 0
            A = rand(dims...)
            T = eltype(A)
        else
            T = Float64
            A = Array{T}(undef, dims)
            A[1] = rand(dims...)
        end
        B = ResizableArray{T}(undef, size(A))
        @test isgrowable(B) == (N > 0)
        @test IndexStyle(typeof(B)) == IndexLinear()
        @test eltype(B) == eltype(A)
        @test elsize(B) == elsize(A)
        @test sizeof(B) == sizeof(A)
        @test ndims(B) == ndims(A) == N
        @test size(B) == size(A)
        @test all(d -> size(B,d) == size(A,d), 1:(N+2))
        @test axes(B) == axes(A)
        @test Base.axes1(B) == axes(B,1)
        @test length(B) == length(A) == prod(dims)
        @test maxlength(B) == length(B)
        @test all(d -> axes(B,d) == axes(A,d), 1:(N+2))
        @test unsafe_convert(Ptr{T}, B) == unsafe_convert(Ptr{T}, parent(B))
        @test pointer(B) == pointer(parent(B))
        @test pointer(B, 2) == pointer(parent(B), 2)
        copyto!(B, A)
        @test all(i -> A[i] == B[i], 1:length(A))
        @test all(i -> A[i] == B[i], CartesianIndices(A))
        @test all(i -> A[i] == parent(B)[i], 1:length(A))
        @test A == B
        for i in eachindex(B)
            B[i] = rand()
        end
        copyto!(A, B)
        @test A == B
        if N > 0
            # Extend array B.
            tmpdims = collect(size(B))
            tmpdims[end] += 1
            resize!(B, tmpdims...)
            for i in length(A)+1:length(B); B[i] = 0; end
            @test maxlength(B) == length(B) == prod(tmpdims)
            @test A != B
            @test B != A
            @test all(i -> B[i] == A[i], 1:length(A))
            C = view(B, axes(A)...)
            @test C != B && B != C
            @test C == A && A == C
            # Shrink array B.
            oldmaxlen = maxlength(B)
            resize!(B, dims)
            @test B == A
            @test C == B && B == C
            @test maxlength(B) == oldmaxlen
            # Use copy to make a fresh resizable copy
            C = copy(ResizableArray, B)
            @test C == B
            @test pointer(C) != pointer(B)
            @test maxlength(C) == length(C)
            C = copy(ResizableArray{T}, B)
            @test C == B
            @test pointer(C) != pointer(B)
            @test maxlength(C) == length(C)
            C = copy(ResizableArray{T,N}, B)
            @test C == B
            @test pointer(C) != pointer(B)
            @test maxlength(C) == length(C)
            shrink!(B)
            @test B == A
            @test maxlength(B) == length(B)
        end

        # Check errors.
        @test_throws BoundsError B[0]
        @test_throws BoundsError B[length(B) + 1]
        @test_throws ErrorException resize!(B, (dims..., 5))
        @test_throws DimensionMismatch ResizableArray{T,N+1}(A)
        @test_throws DimensionMismatch copy(ResizableArray{T,N+1}, A)
        @test_throws ErrorException ResizableArray{T,N,Vector{Char}}(A)
        @test_throws ErrorException copy(ResizableArray{T,N,Vector{Char}}, A)

        # Make a copy of A using a resizable array.
        C = copyto!(similar(ResizableArray{T}, axes(A)), A)
        @test C == A

        # Check equality for a different list of dimensions.
        C = rand(7)
        @test (B == C) == false
        @test (C == B) == false
        @test (B == ResizableArray(C)) == false
        @test (ResizableArray(C) == B) == false

        # Check various constructors and custom buffer
        # (do not splat dimensions if N=0).
        buf = Vector{T}(undef, length(A))
        for arg in (undef, buf), sz in (dims, altdims)
            if isa(arg, Vector)
                # No parameters.
                C = copyto!(ResizableArray(arg, sz), A)
                @test eltype(C) == eltype(A) && C == A
                if N > 0
                    C = copyto!(ResizableArray(arg, sz...), A)
                    @test eltype(C) == eltype(A) && C == A
                end
            end
            # Parameter {T}.
            C = copyto!(ResizableArray{T}(arg, sz), A)
            @test eltype(C) == eltype(A) && C == A
            if N > 0
                C = copyto!(ResizableArray{T}(arg, sz...), A)
                @test eltype(C) == eltype(A) && C == A
            end
            # Parameters {T,N}.
            C = copyto!(ResizableArray{T,N}(arg, sz), A)
            if N > 0
                @test eltype(C) == eltype(A) && C == A
                C = copyto!(ResizableArray{T,N}(arg, sz...), A)
            end
        end

        # Use constructor to convert array.
        C = ResizableArray{T,N}()
        resize!(C, dims)
        copyto!(C, A)
        @test eltype(C) == eltype(A) && C == A

        # Use constructor to convert ordinary array.
        C = ResizableArray(A)
        @test eltype(C) == eltype(A) && C == A
        C = ResizableArray{T}(A)
        @test eltype(C) == eltype(A) && C == A
        C = ResizableArray{T,N}(A)
        @test eltype(C) == eltype(A) && C == A

        # Use convert to convert ordinary array.
        C = convert(ResizableArray, A)
        @test eltype(C) == eltype(A) && C == A
        C = convert(ResizableArray{T}, A)
        @test eltype(C) == eltype(A) && C == A
        C = convert(ResizableArray{T,N}, A)
        @test eltype(C) == eltype(A) && C == A

        # Use convert to convert resizable array, result should be identical.
        for pass in 1:4
            C = (pass == 1 ? convert(ResizableArray,             B) :
                 pass == 2 ? convert(ResizableArray{T},          B) :
                 pass == 3 ? convert(ResizableArray{T,N},        B) :
                 pass == 4 ? convert(ResizableArray{T,N,Vector}, B) : nothing)
            @test eltype(C) == eltype(B)
            @test C == B
            @test pointer(C) == pointer(B)
        end
        # Use convert to convert resizable array, result should be different.
        B = ResizableArray{Int16,N}(undef, dims)
        for i in 1:length(B); B[i] = i; end
        for pass in 1:3
            C = (pass == 1 ? convert(ResizableArray{Int32},   B)        :
                 pass == 2 ? convert(ResizableArray{Int32,N}, B)        :
                 pass == 3 ? convert(ResizableArray{Int32,N,Vector}, B) :
                 nothing)
            @test eltype(C) == Int32
            @test elsize(C) == sizeof(eltype(C))
            @test sizeof(C) == elsize(C)*length(C)
            @test C == B
            @test pointer(C) != pointer(B)
        end
    end
end

@testset "Queue methods" begin
    T = Float32
    @testset "Dimensions: $dims" for dims in ((3,), (2,3), (2,3,4))
        N = length(dims)
        m = 5
        altdims = map(UInt16, dims) # used later
        extdims = (dims..., m)
        extaltdims = (altdims..., Int16(m))
        A = rand(T, dims..., m)
        B = rand(T, dims)
        C = rand(T, dims)
        R = ResizableArray(A)
        append!(R, B)
        @test slice(R, 1:m) == A
        @test slice(R, m+1) == B
        prepend!(R, C)
        @test slice(R, 1) == C
        @test slice(R, 2:m+1) == A
        @test slice(R, m+2) == B
        @test length(extdims) == length(extaltdims) == N+1
        for hint in (prod(extdims), extdims, extaltdims)
            R = sizehint!(ResizableArray{T}(undef, dims..., 0), hint...)
            for k in 1:m
                if isodd(k)
                    append!(R, B)
                else
                    prepend!(R, C)
                end
            end
            @test slice(R, 1) == C
            @test slice(R, m) == B
        end
    end
end

@testset "Iterations" begin
    T = Float64
    @testset "Dimensions: $dims" for dims in ((3,), (2,3), (2,3,4))
        N = length(dims)
        A = rand(T, dims)
        B = ResizableArray(A)
        val = sum(A)
        @test val ≈ sum(B)
        @test val ≈ sum_v1(B)
        @test val ≈ sum_v2(B)
        @test val ≈ sum_v3(B)
    end
end

end # module


text:
# Autogenerated wrapper script for atomsk_jll for aarch64-linux-gnu
export atomsk

using LAPACK_jll
JLLWrappers.@generate_wrapper_header("atomsk")
JLLWrappers.@declare_executable_product(atomsk)
function __init__()
    JLLWrappers.@generate_init_header(LAPACK_jll)
    JLLWrappers.@init_executable_product(
        atomsk,
        "bin/atomsk",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#=
    Following the description in
        _Simulation of two-dimensional quantum systems on an infinite lattice revisited:
        Corner Transfer matrix for tensor contraction_
    by Orus and Vidal

    Implemented here is the directional variant of the CTMRG as outlined
    on page 2 and in Fig2.
    DOI: 10.1103/PhysRevB.80.094403
=#

struct CTMIterable{T,
        TA <: AbstractTensor{T,4},
        TC <: AbstractTensor{T,2},
        TT <: AbstractTensor{T,3}} <: AbstractCTMIterable
    A::TA
    χ::Int
    Csinit::Union{NTuple{4,TC},Nothing}
    Tsinit::Union{NTuple{4,TT},Nothing}
end

function ctmiterable(A::DTensor{T,4}, χ::Int, Csinit = nothing, Tsinit = nothing) where T
    CTMIterable{T,DTensor{T,4},DTensor{T,2},DTensor{T,3}}(A, χ, Csinit, Tsinit)
end

function ctmiterable(A::DASTensor{T,4,SYM,CHS,SS,CH}, χ::Int,
        Csinit = nothing, Tsinit = nothing) where {T,N,SYM,CHS,SS,CH}
    TC = DASTensor{T,2,SYM,CHS,SS,CH}
    TT = DASTensor{T,3,SYM,CHS,SS,CH}
    TA = DASTensor{T,4,SYM,CHS,SS,CH}
    CTMIterable{T,TA,TC,TT}(A, χ, Csinit, Tsinit)
end

struct CTMState{S,
        TA <: AbstractTensor,
        TC <: AbstractTensor,
        TT <: AbstractTensor} <: AbstractCTMState
    Cs::NTuple{4,TC}
    Ts::NTuple{4,TT}
    oldsvdvals::Vector{S}
    diffs::Vector{S}
    n_it::Ref{Int}
end

function iterate(iter::CTMIterable{S,TA,TC,TT}) where {S,TA,TC,TT}
    @unpack A, χ, Tsinit, Csinit = iter
    Cs = Csinit == nothing ? ntuple(i -> initializeC(A,χ), 4) : deepcopy(Csinit)
    Ts = Tsinit == nothing ? ntuple(i -> initializeT(A,χ,i), 4) : deepcopy(Tsinit)

    l = ifelse(TA <: DASTensor, 2χ, χ)
    oldsvdvals = zeros(real(S),l)
    state = CTMState{real(S),TA,TC,TT}(Cs, Ts, oldsvdvals, [], Ref(0))
    return state, state
end

function iterate(iter::CTMIterable, state::CTMState{S,TA,TC,TT}) where {S,TA,TC,TT}
    @unpack A, χ = iter
    @unpack Cs, Ts, oldsvdvals, diffs, n_it = state
    C1, C2, C3, C4 = Cs
    T1, T2, T3, T4 = Ts
    #=
        [C1] - [T1] - [C2]
          |      |      |
        [T4] - [ A] - [T2]
          |      |      |
        [C4] - [T3] - [C3]
    =#

    #leftmove
    leftmove!((C1,T4,C4),(T1,A,T3),χ)

    # upmove
    upmove!((C1,T1,C2),(T4,A,T2),χ)

    #rightmove
    rightmove!((C3,T2,C2),(T3,A,T1),χ)

    #downmove
    downmove!((C4,T3,C3),(T4,A,T2),χ)

    @tensor C[1,2] := C1[1,-1] * C2[-1,-2] * C3[-2,-3] * C4[-3,2]
    vals  = diag(tensorsvd(C)[2])
    normalize!(vals,1)

    #compare
    push!(diffs, sum(abs, oldsvdvals - vals))
    oldsvdvals[:] = vals
    n_it[] += 1

    return state, state
end

function leftmove!((C1,T4,C4), (T1,A,T3), χ)
    #=
    leftmove:
        [C1]--[T1]--
          |    |
        [T4]--[ A]--
          |    |
        [C4]--[T3]--
    =#

    @tensor begin
        C1p[1,2,3]     := C1[1,-1]   * T1[-1,2,3]
        T4p[1,2,3,4,5] := T4[1,-1,5] * A[2,3,4,-1]
        C4p[1,2,3]     := C4[-1,3]   * T3[1,2,-1]
    end

    #renormalize
    @tensor begin
        C1pC1p[1,2,3,4] := C1p[1,2,-1] * C1p'[4,3,-1]
        C4pC4p[1,2,3,4] := C4p[-1,2,1] * C4p'[-1,3,4]
        CCpCC[1,2,4,3]  := C1pC1p[1,2,3,4] + C4pC4p[1,2,3,4]
    end

    CCpCCmat, reshaper = fuselegs(CCpCC,((1,2),(3,4)))
    E, = tensoreig(CCpCCmat,truncfun = svdtrunc_maxχ(χ))
    U = splitlegs(E, ((1,1,1),(1,1,2),2), reshaper...)

    @tensor begin
        C1[1,2]   = U'[-1,-2,1] * C1p[-1,-2,2]
        T4[1,2,3] = U'[-1,-2,1] * T4p[-1,-2,2,-4,-3] * U[-3,-4,3]
        C4[1,2]   = C4p[1,-1,-2] * U[-2,-1,2]
    end

    mval = maximum(diag(tensorsvd(C1)[2]))
    apply!(C1, x -> x .= x ./ mval)
    mval = maximum(diag(tensorsvd(C4)[2]))
    apply!(C4, x -> x .= x ./ mval)
    apply!(T4, x -> x .= x ./ sqrt(norm(T4)))
end

function rightmove!((C3,T2,C2),(T3,A,T1),χ)
    #=
    rightmove:
      --[T1]---[C2]
         |      |
      --[ A]---[T2]
         |      |
      --[T3]---[C3]
    =#

    @tensor begin
        C2p[1,2,3]     := T1[1,2,-1]  * C2[-1,3]
        T2p[1,2,3,4,5] := A[4,-1,2,3] * T2[1,-1,5]
        C3p[1,2,3]     := C3[1,-1]    * T3[-1,2,3]
    end

    #renormalize
    @tensor begin
        C2pC2p[1,2,3,4] := C2p[-1,2,1] * C2p'[-1,3,4]
        C3pC3p[1,2,3,4] := C3p[1,2,-1] * C3p'[4,3,-1]
        CCpCC[1,2,4,3]  := C2pC2p[1,2,3,4] + C3pC3p[1,2,3,4]
    end

    CCpCCmat, reshaper = fuselegs(CCpCC,((1,2),(3,4)))
    E, = tensoreig(CCpCCmat,truncfun = svdtrunc_maxχ(χ))
    U = splitlegs(E, ((1,1,1),(1,1,2),2), reshaper...)

    @tensor begin
        C2[1,2]   = C2p[1,-1,-2] * U'[-2,-1,2]
        T2[1,2,3] = U[-1,-2,1] * T2p[-1,-2,2,-4,-3] * U'[-3,-4,3]
        C3[1,2]   = C3p[-1,-2,2] * U[-1,-2,1]
    end

    mval = maximum(diag(tensorsvd(C2)[2]))
    apply!(C2, x -> x .= x ./ mval)
    mval = maximum(diag(tensorsvd(C3)[2]))
    apply!(C3, x -> x .= x ./ mval)
    apply!(T2, x -> x .= x ./ sqrt(norm(T2)))
end

function upmove!((C1,T1,C2),(T4,A,T2),χ)
    #=
    upmove:
        [C1]--[T1]--[C2]
          |     |     |
        [T4]--[ A]--[T2]
          |     |     |
    =#

    @tensor begin
        C1p[1,2,3]     := T4[1,2,-1]  * C1[-1,3]
        T1p[1,2,3,4,5] := A[3,4,-1,2] * T1[1,-1,5]
        C2p[1,2,3]     := C2[1,-1]    * T2[-1,2,3]
    end

    #renormalize
    @tensor begin
        C1pC1p[1,2,3,4] := C1p[-1,2,1] * C1p'[-1,3,4]
        C2pC2p[1,2,3,4] := C2p[1,2,-1] * C2p'[4,3,-1]
        CCpCC[1,2,4,3]  := C1pC1p[1,2,3,4] + C2pC2p[1,2,3,4]
    end

    CCpCCmat, reshaper = fuselegs(CCpCC,((1,2),(3,4)))
    E, = tensoreig(CCpCCmat,truncfun = svdtrunc_maxχ(χ))
    U = splitlegs(E, ((1,1,1),(1,1,2),2), reshaper...)

    @tensor begin
        C1[1,2]   = C1p[1,-1,-2] * U'[-2,-1,2]
        T1[1,2,3] = U'[-1,-2,1] * T1p[-1,-2,2,-4,-3] * U[-3,-4,3]
        C2[1,2]   = C2p[-1,-2,2] * U[-1,-2,1]
    end

    mval = maximum(diag(tensorsvd(C1)[2]))
    apply!(C1, x -> x .= x ./ mval)
    mval = maximum(diag(tensorsvd(C2)[2]))
    apply!(C2, x -> x .= x ./ mval)
    apply!(T1, x -> x .= x ./ sqrt(norm(T1)))
end

function downmove!((C4,T3,C3),(T4,A,T2),χ)
    #=
    downmove:
          |     |     |
        [T4]--[ A]--[T2]
          |     |     |
        [C4]--[T3]--[C3]
    =#

    @tensor begin
        C4p[1,2,3]     := C4[1,-1]  * T4[-1,2,3]
        T3p[1,2,3,4,5] := T3[1,-1,5] * A[-1,2,3,4]
        C3p[1,2,3]     := T2[1,2,-1] * C3[-1,3]
    end

    #renormalize
    @tensor begin
        C3pC3p[1,2,3,4] := C3p[-1,2,1] * C3p'[-1,3,4]
        C4pC4p[1,2,3,4] := C4p[1,2,-1] * C4p'[4,3,-1]
        CCpCC[1,2,4,3]  := C3pC3p[1,2,3,4] + C4pC4p[1,2,3,4]
    end

    CCpCCmat, reshaper = fuselegs(CCpCC,((1,2),(3,4)))
    E, = tensoreig(CCpCCmat,truncfun = svdtrunc_maxχ(χ))
    U = splitlegs(E, ((1,1,1),(1,1,2),2), reshaper...)

    @tensor begin
        C3[1,2]   = C3p[1,-1,-2] * U[-2,-1,2]
        T3[1,2,3] = U'[-1,-2,1] * T3p[-1,-2,2,-3,-4] * U[-4,-3,3]
        C4[1,2]   = C4p[-1,-2,2] * U'[-1,-2,1]
    end

    mval = maximum(diag(tensorsvd(C4)[2]))
    apply!(C4, x -> x .= x ./ mval)
    mval = maximum(diag(tensorsvd(C3)[2]))
    apply!(C3, x -> x .= x ./ mval)
    apply!(T3, x -> x .= x ./ sqrt(norm(T3)))
end


text:
using LiveServer
servedocs()


text:
using MDPTools, NestedFixedPoint
# Set RNG seed
srand(124)

# State space
n = 175
o_max = 450
Xval = linspace(0, o_max, n)
p = [0.0937; 0.4475; 0.4459; 0.0127; 0.0002];

F1 = zeros(n, n)
offset = -1
for i = 1:n, j = 1:length(p)
    i+offset+j > n && continue
    F1[i,i+offset+j] = i+offset+j == n ? sum(p[j:end]) : p[j]
end

# We can handle sparse matrices
sparse_F = [sparse(F1), sparse(F1)[ones(Int64, n), :]]
F = [F1, F1[ones(Int64, n), :]]

S = State(Xval, F)
sparse_S = State(Xval, sparse_F)
bdisc = 0.9999
Z1 = [zeros(n) -0.001*Xval]
Z2 = [-ones(n) zeros(n)]

U = LinearUtility([Z1, Z2], bdisc, [11.;2.5])
solve(U,S)
T = 120
N = 50

D = simulate(U, S, 1, T, N)

U.Θ*=0.
rust_nfxp = fit_nfxp(U, S, D)


text:
"""
Core functionality

## Basic usage

The primary export of the package is `Mat{T,N}`, which is the Julia type for
`cv::Mat`. `Mat{T,N}` is designed to be a subtype of `AbstractArray{T,N}`. It
has element type (`T`) and dimension (`N`) as type parameters, whereas
`cv::Mat` doesn't. Note that matrix construction interface is different between
Julia and C++. `cv::Mat(3,3,CV_8U)` in C++ can be translated in Julia as
`Mat{UInt8}(3,3)`.

### Create uninitialized matrix:

```julia
using CVCore

julia> A = Mat{Float64}(3,3)
3×3 CVCore.Mat{Float64,2}:
   3.39519e-313  4.94066e-324   2.122e-314
 NaN             1.72723e-77   -2.32036e77
   1.97626e-323  0.0            0.0
```

### Create matrix filled with Scalar (`cv::Scalar`)

```julia
julia> A = Mat{Float64}(3,3,Scalar(1))
3×3 CVCore.Mat{Float64,2}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
```

### Matirx operations

```julia
julia> A * A
CVCore.MatExpr{Float64,2}
3×3 CVCore.Mat{Float64,2}:
 3.0  3.0  3.0
 3.0  3.0  3.0
 3.0  3.0  3.0

julia> A + A
CVCore.MatExpr{Float64,2}
3×3 CVCore.Mat{Float64,2}:
 2.0  2.0  2.0
 2.0  2.0  2.0
 2.0  2.0  2.0

julia> A - A
CVCore.MatExpr{Float64,2}
3×3 CVCore.Mat{Float64,2}:
 0.0  0.0  0.0
 0.0  0.0  0.0
 0.0  0.0  0.0
```

### Create multi-channel matrix

```julia
julia> A = Mat{Float64}(3,3,3,Scalar(1,2,3))
3×3×3 CVCore.Mat{Float64,3}:
[:, :, 1] =
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

[:, :, 2] =
 2.0  2.0  2.0
 2.0  2.0  2.0
 2.0  2.0  2.0

[:, :, 3] =
 3.0  3.0  3.0
 3.0  3.0  3.0
 3.0  3.0  3.0
```

### Conversion between `Mat{T,N}` and `Array{T,N}`

```julia
julia> B = Array(A)
3×3 Array{Float64,2}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> C = Mat(B)
3×3 CVCore.Mat{Float64,2}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0
```

Note that `Mat(B)` where `B` is an array shares the underlying data.

"""
module CVCore

export AbstractCvMat, MatExpr, Mat, UMat, depth, channels, flags, dims, rows,
    cols, clone, total, isContinuous, elemSize, Scalar, TermCriteria,
    convertScaleAbs!, convertScaleAbs, addWeighted!, addWeighted

#=
Naming convention:
  1. Cxx types should have prefix cv: e.g. cv::Mat -> cvMat
  2. Julia types should not have prefix cv: e.g. cv::Mat -> Mat
=#

using LibOpenCV
using Cxx

"""
Returns a refrence of the interal structure
"""
function handle end

include("macros.jl")
include("const.jl")

import Base: convert, eltype, size

### Cxx types ###

cvScalar_{T} =  cxxt"cv::Scalar_<$T>"

const cvScalar = cxxt"cv::Scalar"
cvScalar(v) = icxx"return cv::Scalar($v);"

const AbstractCvScalar = Union{cvScalar, cvScalar_}

cvPoint_{T} =  cxxt"cv::Point_<$T>"
cvPoint_{T}(x, y) where {T} = icxx"cv::Point_<$T>($x, $y);"
eltype(p::cvPoint_{T}) where {T} = T

cvPoint3_{T} =  cxxt"cv::Point3_<$T>"
cvPoint3_{T}(x, y, z) where {T} = icxx"cv::Point3_<$T>($x, $y, $z);"
eltype(p::cvPoint3_{T}) where {T} = T

const cvPoint = cxxt"cv::Point"
cvPoint(x, y) = icxx"cv::Point($x, $y);"

const AbstractCvPoint = Union{cvPoint, cvPoint_}

cvSize_{T} =  cxxt"cv::Size_<$T>"
cvSize_{T}(x, y) where {T} = icxx"cv::Size_<$T>($x, $y);"
eltype(s::cvSize_{T}) where {T} = T

const cvSize = cxxt"cv::Size"
cvSize(x, y) = icxx"cv::Size($x, $y);"

const AbstractCvSize = Union{cvSize, cvSize_}

height(s::AbstractCvSize) = Int(icxx"$s.height;")
width(s::AbstractCvSize) = Int(icxx"$s.width;")
area(s::AbstractCvSize) = Int(icxx"$s.area();")

"""Determine julia type from the depth of cv::Mat
"""
function jltype(depth)
    if depth == CV_8U
        return UInt8
    elseif depth == CV_8S
        return Int8
    elseif depth == CV_16U
        return UInt16
    elseif depth == CV_16S
        return Int16
    elseif depth == CV_32S
        return Int32
    elseif depth == CV_32F
        return Float32
    elseif depth == CV_64F
        return Float64
    else
        error("This shouldn't happen")
    end
end

"""Determine cv::Mat depth from Julia type
"""
function cvdepth(T)
    if T == UInt8
        return CV_8U
    elseif T == Int8
        return CV_8S
    elseif T == UInt16
        return CV_16U
    elseif T == Int16
        return CV_16S
    elseif T == Int32
        return CV_32S
    elseif T == Float32
        return CV_32F
    elseif T == Float64
        return CV_64F
    else
        error("$T: not supported in cv::Mat")
    end
end

mat_depth(flags) = flags & CV_MAT_DEPTH_MASK
mat_channel(flags) = (flags & CV_MAT_CN_MASK) >> CV_CN_SHIFT + 1
maketype(depth, cn) = mat_depth(depth) + ((cn-1) << CV_CN_SHIFT)

### Scalar ###

# TODO: need to be subtype of cv::Vec
const Scalar = cvScalar_
Scalar{T}(s1=01,s2=0,s3=0,s4=0) where {T} =
    icxx"return cv::Scalar_<$T>($s1,$s2,$s3,$s4);"
Scalar(s1=0,s2=0,s3=0,s4=0) = Scalar{Float64}(s1,s2,s3,s4)
eltype(s::Scalar{T}) where {T} = T

### cv::TermCriteria ###

const TERM_CRITERIA_COUNT = icxx"cv::TermCriteria::COUNT;"
const TERM_CRITERIA_MAX_ITER = icxx"cv::TermCriteria::MAX_ITER;"
const TERM_CRITERIA_EPS = icxx"cv::TermCriteria::EPS;"

const TermCriteria = cxxt"cv::TermCriteria"
TermCriteria(typ, maxCount, epsilon) =
    icxx"cv::TermCriteria($typ, $maxCount, $epsilon);"

include("mat.jl")

@gencxxf(convertScaleAbs!(src::AbstractCvMat, dst::AbstractCvMat,
    alpha=1,beta=0), "cv::convertScaleAbs")
function convertScaleAbs(src::AbstractCvMat, alpha=1, beta=0)
    dst = similar_empty(src)
    convertScaleAbs!(src, dst, alpha, beta)
    dst
end

@gencxxf(addWeighted!(src1::AbstractCvMat, alpha, src2::AbstractCvMat, beta,
    gamma, dst::AbstractCvMat, dtype=-1), "cv::addWeighted")
function addWeighted(src1::AbstractCvMat, alpha, src2::AbstractCvMat, beta,
    gamma, dtype=-1)
    dst = similar_empty(src)
    addWeighted!(src1, alpha, src2, beta, gamma, dst, dtype)
    dst
end

end # module


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    struct ForwardDiffAD

*Experimental feature, not part of stable public API.*

Constructors:

* ```ForwardDiffAD()```
"""
struct ForwardDiffAD <: DifferentiationAlgorithm end
export ForwardDiffAD

function unshaped_gradient!(grad_f_x::AbstractVector{<:Real}, ::Type{R}, f::Function, x::AbstractVector{<:Real}, diffalg::ForwardDiffAD) where {R<:Real}
    result = DiffResults.MutableDiffResult(zero(R), (grad_f_x,))

    # chunk = ForwardDiff.Chunk(x)
    # config = ForwardDiff.GradientConfig(f, x, chunk)
    # ForwardDiff.gradient!(result, f, x, config)
    ForwardDiff.gradient!(result, f, x)

    @assert DiffResults.gradient(result) === grad_f_x
    convert(R, DiffResults.value(result))::R
end


text:
push!(LOAD_PATH,"../src/")

using Documenter
using CalculatedABC

makedocs(
    sitename = "CalculatedABC.jl",
    authors = "Tobias Frilling",
    format = Documenter.HTML(),
    modules = [CalculatedABC],
    pages = [
        "Home" => "index.md",
        "ABC Analysis" => "abcanalysis.md",
        "Plotting" => "plotting.md",
        "Gini Coefficient" => "gini.md"
    ]

)

deploydocs(
    repo = "github.com/ckafi/CalculatedABC.jl"
)


text:
using Kinetic, Plots
using KitBase.JLD2
using Base.Threads: @threads
using KitBase.ProgressMeter: @showprogress
cd(@__DIR__)

X = Float32.([[1e-4, 0.0, 0.0, 1e-4]; zeros(4); 1.0])
Y = Int32(0)

begin
    set = Setup(space = "1d1f3v", boundary = ["fix", "fix"], collision = "fsm")
    ps = PSpace1D(-0.5, 0.5, 250, 1)
    vs = VSpace3D(-6.0, 6.0, 28, -6.0, 6.0, 28, -6.0, 6.0, 28)
    gas = begin
        Kn = 5e-3
        Gas(Kn = Kn, K = 0.0, fsm = fsm_kernel(vs, ref_vhs_vis(Kn, 1.0, 0.5)))
    end
    fw = function(x)
        prim = zeros(5)
        if x <= 0
            prim .= [1.0, 0.0, 1.0, 0.0, 1.0]
        else
            prim .= [1.0, 0.0, -1.0, 0.0, 2.0]
        end

        return prim_conserve(prim, ks.gas.γ)
    end
    ib = IB1F(fw, vs, gas)
    ks = SolverSet(set, ps, vs, gas, ib)
    ctr, face = init_fvm(ks, ks.ps)
end

τ0 = vhs_collision_time(ctr[1].prim, ks.gas.μᵣ, ks.gas.ω)
tmax = 50 * τ0
t = 0.0
dt = timestep(ks, ctr, t)
nt = Int(tmax ÷ dt)
res = zero(ctr[1].w)

@showprogress for iter = 1:nt
    evolve!(ks, ctr, face, dt)
    update!(ks, ctr, face, dt, res)

    if iter%5 == 0
        for i = 1:ks.ps.nx
            if abs(ctr[i].prim[2]) > 1e-6
                sw = (ctr[i+1].w .- ctr[i-1].w) / ks.ps.dx[i] / 2
                x, y = regime_data(ks, ctr[i].w, sw, zero(sw), zero(sw), ctr[i].f)
                rx = [x[1:3]; x[5:8]; x[10:11]]
                global X = hcat(X, rx)
                global Y = hcat(Y, y)
            end
        end
    end

    global t += dt
end
@save "data/dataset.jld2" X Y

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
function test(signalVector)

    memory = length(signalVector)

    # we use binary 0,1 under the covers - so lets get rid of the -1 and replace w/0 -
    binarySignal = replace(signalVector, -1 => 0)
    iv = range(memory, stop=1, step=-1) |> collect

    # convert the binarySignal to an int -
    strategy_index = 0
    for index in iv
        value = binarySignal[index]
        strategy_index += value * 2^(memory - index)
    end

    return (strategy_index)
end

signalVector = [-1,-1, 1]
r = test(signalVector)


text:
# Quick test to make sure that import works as intended
# Solve problem defined by imported file using ConvexJL and SCS
# by Michael Garstka
# University of Oxford, Control Group
# January 2018
# updated in April 2019

# We work with a semidefinite programming problem that has been written
# in the following standard form:

# (P)    min c1*x1+c2*x2+...+cm*xm
#        st  F1*x1+F2*x2+...+Fm*xn-F0=X
#                                  X >= 0

using FileIO, COSMO, SparseArrays, LinearAlgebra, Test

# import data from .jld2 file and create variables
data = load("./sdplib/maxG11.jld2");
F = data["F"]
c = data["c"]
m = data["m"]
n = data["n"]
obj_true = data["optVal"]

# Describe primal problem using JuMP and solve with SCS Solver
model = JuMP.Model(with_optimizer(COSMO.Optimizer, verbose = true, decompose = true));
@variable(model, x[1:m]);
@objective(model, Min, c' * x);
@constraint(model, con1,  Symmetric(-Matrix(F[1]) + sum(Matrix(F[k + 1]) .* x[k] for k in 1:m))  in JuMP.PSDCone());
JuMP.optimize!(model);

# solve with native interface

d = div(n * (n + 1), 2)
A_man = zeros(d, m)
# now manually create A and B
for (i, Fi) in enumerate(F[2:end])
  t = zeros(d)
  COSMO.extract_upper_triangle!(Fi, t, sqrt(2))
  A_man[:, i] = t
end

b_man = zeros(d)
COSMO.extract_upper_triangle!(-F[1], b_man, sqrt(2))

cs1 = COSMO.Constraint(A_man, b_man, COSMO.PsdConeTriangle)
model_direct = COSMO.Model();
assemble!(model_direct, spzeros(m, m), c, cs1, settings = COSMO.Settings(verbose = true, decompose = false));
res = COSMO.optimize!(model_direct);


text:
@testset "LogDet" begin
    bridged_mock = MOIB.LogDet{Float64}(mock)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [0, 1, 0, 1, 1, 0, 1, 0, 0, 1])
    MOIT.logdett1vtest(bridged_mock, config)
    MOIT.logdett1ftest(bridged_mock, config)
    # Dual is not yet implemented for LogDet bridge
    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.LogDetConeTriangle}()))
    test_delete_bridge(bridged_mock, ci, 5, ((MOI.VectorAffineFunction{Float64}, MOI.ExponentialCone, 0), (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle, 0)))
end

@testset "RootDet" begin
    bridged_mock = MOIB.RootDet{Float64}(mock)
    mock.optimize! = (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(mock, [1, 1, 0, 1, 1, 0, 1])
    MOIT.rootdett1vtest(bridged_mock, config)
    MOIT.rootdett1ftest(bridged_mock, config)
    # Dual is not yet implemented for RootDet bridge
    ci = first(MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorAffineFunction{Float64}, MOI.RootDetConeTriangle}()))
    test_delete_bridge(bridged_mock, ci, 4, ((MOI.VectorAffineFunction{Float64}, MOI.RotatedSecondOrderCone, 0),
                                            (MOI.VectorAffineFunction{Float64}, MOI.GeometricMeanCone, 0),
                                            (MOI.VectorAffineFunction{Float64}, MOI.PositiveSemidefiniteConeTriangle, 0)))
end


text:
"""
Returns: (n, perm) where n is the number of columns in
the beginning of `cols`, `perm` is one possible permutation of those
first `n` columns.
"""
function best_perm_estimate(perms, cols)
    bestperm = nothing
    bestmatch = 0
    for p in perms
        matched_cols = 0
        l = min(length(cols), length(p.columns))
        for i in 1:l
            cols[i] != p.columns[i] && break
            matched_cols += 1
        end
        if matched_cols == length(cols)
            return (matched_cols, p.perm)
        elseif bestmatch < matched_cols
            bestmatch = matched_cols
            bestperm = p.perm
        end
    end
    return (bestmatch, bestperm)
end

function sortpermby(t, by; cache=false, return_keys=false)
    by = lowerselection(t, by)
    if !isa(by, Tuple)
        by = (by,)
    end
    canonorder = map((x, y)->y isa Pair ? 0 : x, colindex(t, by), by)
    canonorder_vec = Int[canonorder...]
    perms = permcache(t)
    matched_cols, partial_perm = best_perm_estimate(perms, canonorder_vec)

    if matched_cols == length(by)
        # first n index columns
        return return_keys ? (partial_perm, compact_mem(rows(t, by))) : partial_perm
    end

    byrows = compact_mem(rows(t, by))
    bycols = columns(byrows)
    perm = if matched_cols > 0
        nxtcol = bycols[matched_cols+1]
        p = convert(Array{UInt32}, partial_perm)
        refine_perm!(p, bycols, matched_cols,
                     rows(byrows, Tuple(1:matched_cols)),
                     nxtcol, 1, length(t))
        p
    else
        sortperm(byrows)
    end

    if cache
        cacheperm!(t, Perm(canonorder_vec, perm))
    end

    return return_keys ? (perm, byrows) : perm
end

function sortpermby(t, by::AbstractArray; cache=true)
    sortperm_fast(by)
end