text:
using LinearAlgebra, LoopVectorization, Libdl

const REGISTER_SIZE = LoopVectorization.VectorizationBase.register_size()

# const LOOPVECBENCHDIR = joinpath(pkgdir(LoopVectorization), "benchmark")
include(joinpath(LOOPVECBENCHDIR, "looptests.jl"))

const LIBCTEST = joinpath(LOOPVECBENCHDIR, "libctests.so")
const LIBFTEST = joinpath(LOOPVECBENCHDIR, "libftests.so")
const LIBICTEST = joinpath(LOOPVECBENCHDIR, "libictests.so")
const LIBIFTEST = joinpath(LOOPVECBENCHDIR, "libiftests.so")
const LIBEIGENTEST = joinpath(LOOPVECBENCHDIR, "libetest.so")
const LIBIEIGENTEST = joinpath(LOOPVECBENCHDIR, "libietest.so")

# requires Clang with polly to build
cfile = joinpath(LOOPVECBENCHDIR, "looptests.c")
if !isfile(LIBCTEST) || mtime(cfile) > mtime(LIBCTEST)
  if (Sys.ARCH === :aarch64) && Sys.isapple() # assume no `-march=native` support
    run(
      `clang -Ofast -mprefer-vector-width=$(8REGISTER_SIZE) -lm -shared -fPIC $cfile -o $LIBCTEST`,
    )
  else
    run(
      `clang -Ofast -march=native -mprefer-vector-width=$(8REGISTER_SIZE) -lm -shared -fPIC $cfile -o $LIBCTEST`,
    )
  end

  # run(`/usr/local/bin/clang -Ofast -march=native -mprefer-vector-width=$(8REGISTER_SIZE) -lm -mllvm -polly -mllvm -polly-vectorizer=stripmine -shared -fPIC $cfile -o $LIBCTEST`)
end
ffile = joinpath(LOOPVECBENCHDIR, "looptests.f90")
if !isfile(LIBFTEST) || mtime(ffile) > mtime(LIBFTEST)
  # --param max-unroll-times defaults to ≥8, which is generally excessive
  if (Sys.ARCH === :x86_64)
    run(
      `gfortran -Ofast -march=native -funroll-loops -mprefer-vector-width=$(8REGISTER_SIZE) -fvariable-expansion-in-unroller --param max-variable-expansions-in-unroller=4 -shared -fPIC $ffile -o $LIBFTEST`,
    )
  else
    run(
      `gfortran -Ofast -march=native -funroll-loops -fvariable-expansion-in-unroller --param max-variable-expansions-in-unroller=4 -shared -fPIC $ffile -o $LIBFTEST`,
    )
  end
  # run(`gfortran -Ofast -march=native -funroll-loops -floop-nest-optimize -mprefer-vector-width=$(8REGISTER_SIZE) -shared -fPIC $ffile -o $LIBFTEST`)
end

const INTEL_BENCH = try
  if !isfile(LIBIFTEST) || mtime(ffile) > mtime(LIBIFTEST)
    run(
      `ifort -fast -qopt-zmm-usage=high -qoverride-limits -shared -fPIC $ffile -o $LIBIFTEST`,
    )
  end
  if !isfile(LIBICTEST) || mtime(cfile) > mtime(LIBICTEST)
    run(
      `icc -fast -qopt-zmm-usage=high -fargument-noalias-global -qoverride-limits -shared -fPIC $cfile -o $LIBICTEST`,
    )
  end
  true
catch ProcessFailedException
  false
end

# g++ -> ICE, so Clang and ICPC
eigenfile = joinpath(LOOPVECBENCHDIR, "looptestseigen.cpp")
if !isfile(LIBEIGENTEST) || mtime(eigenfile) > mtime(LIBEIGENTEST)
  # Clang seems to have trouble finding includes
  if Bool(LoopVectorization.VectorizationBase.has_feature(Val(:x86_64_avx512f)))
    run(
      `g++ -O3 -march=native -mprefer-vector-width=$(8REGISTER_SIZE) -DEIGEN_VECTORIZE_AVX512 -I/usr/include/eigen3 -shared -fPIC $eigenfile -o $LIBEIGENTEST`,
    )
  elseif (Sys.ARCH === :aarch64) && Sys.isapple() # assume homebrew
    run(
      `g++-10 -O3 -march=native -I/opt/homebrew/Cellar/eigen/3.3.9/include/eigen3 -shared -fPIC $eigenfile -o $LIBEIGENTEST`,
    )
  else
    run(
      `g++ -O3 -march=native -mprefer-vector-width=$(8REGISTER_SIZE) -I/usr/include/eigen3 -shared -fPIC $eigenfile -o $LIBEIGENTEST`,
    )
  end
end
if !isfile(LIBIEIGENTEST) || mtime(eigenfile) > mtime(LIBIEIGENTEST)
  # run(`/usr/bin/clang++ -Ofast -march=native -mprefer-vector-width=$(8REGISTER_SIZE) -DEIGEN_VECTORIZE_AVX512 -I/usr/include/c++/9 -I/usr/include/c++/9/x86_64-generic-linux -I/usr/include/eigen3 -shared -fPIC $eigenfile -o $LIBEIGENTEST`)
  if Bool(LoopVectorization.VectorizationBase.has_feature(Val(:x86_64_avx512f)))
    run(
      `clang++ -Ofast -march=native -mprefer-vector-width=$(8REGISTER_SIZE) -DEIGEN_VECTORIZE_AVX512 -I/usr/include/eigen3 -shared -fPIC $eigenfile -o $LIBIEIGENTEST`,
    )
  elseif (Sys.ARCH === :aarch64) && Sys.isapple() # assume homebrew and no `-march=native`
    run(
      `clang++ -Ofast -I/opt/homebrew/Cellar/eigen/3.3.9/include/eigen3 -shared -fPIC $eigenfile -o $LIBIEIGENTEST`,
    )
  else
    run(
      `clang++ -Ofast -march=native -mprefer-vector-width=$(8REGISTER_SIZE) -I/usr/include/eigen3 -shared -fPIC $eigenfile -o $LIBIEIGENTEST`,
    )
  end
  # run(`icpc -fast -qopt-zmm-usage=high -fargument-noalias-global -qoverride-limits -I/usr/include/eigen3 -shared -fPIC $eigenfile -o $LIBIEIGENTEST`)
end

# MKL_ROOT = "/home/chriselrod/intel"
# directcalljitfile = joinpath(LOOPVECBENCHDIR, "directcalljit.f90")
# if !isfile(LIBDIRECTCALLJIT) || mtime(directcalljitfile) > mtime(LIBDIRECTCALLJIT)
#     run(`ifort -fast -DMKL_DIRECT_CALL_SEQ_JIT -fpp -qopt-zmm-usage=high -Wl,--start-group $(joinpath(MKL_ROOT,"mkl/lib/intel64/libmkl_intel_lp64.a")) $(joinpath(MKL_ROOT,"mkl/lib/intel64/libmkl_sequential.a")) $(joinpath(MKL_ROOT,"mkl/lib/intel64/libmkl_core.a")) -Wl,--end-group -I$(joinpath(MKL_ROOT, "mkl/include")) -I$(joinpath(MKL_ROOT, "compilers_and_libraries_2020.1.217/linux/mkl/include/intel64/lp64")) -shared -fPIC $directcalljitfile -o $LIBDIRECTCALLJIT`)
#     # run(`gfortran -Ofast -march=native -DMKL_DIRECT_CALL_SEQ_JIT -cpp -mprefer-vector-width=$(8REGISTER_SIZE) -Wl,--start-group $(joinpath(MKL_ROOT,"mkl/lib/intel64/libmkl_intel_lp64.a")) $(joinpath(MKL_ROOT,"mkl/lib/intel64/libmkl_sequential.a")) $(joinpath(MKL_ROOT,"mkl/lib/intel64/libmkl_core.a")) -Wl,--end-group -I$(joinpath(MKL_ROOT, "mkl/include")) -I$(joinpath(MKL_ROOT, "compilers_and_libraries_2020.1.217/linux/mkl/include/intel64/lp64")) -shared -fPIC $directcalljitfile -o $LIBDIRECTCALLJIT`)

#     # run(`gfortran -Ofast -march=native -DMKL_DIRECT_CALL_SEQ_JIT -cpp -mprefer-vector-width=$(8REGISTER_SIZE) -shared -fPIC $directcalljitfile -o $LIBDIRECTCALLJIT`)
# end

randa(::Type{T}, dim...) where {T} = rand(T, dim...)
randa(::Type{T}, dim...) where {T<:Signed} = rand(T(-100):T(200), dim...)

using MKL_jll, OpenBLAS_jll

const MKL_BENCH = isdefined(MKL_jll, :libmkl_rt)

if MKL_BENCH
  const libMKL = Libdl.dlopen(MKL_jll.libmkl_rt)
  const DGEMM_MKL = Libdl.dlsym(libMKL, :dgemm)
  const SGEMM_MKL = Libdl.dlsym(libMKL, :sgemm)
  const DGEMV_MKL = Libdl.dlsym(libMKL, :dgemv)
  const MKL_SET_NUM_THREADS = Libdl.dlsym(libMKL, :MKL_Set_Num_Threads)
  true
else
  false
end

const libOpenBLAS = Libdl.dlopen(OpenBLAS_jll.libopenblas)
const DGEMM_OpenBLAS = Libdl.dlsym(libOpenBLAS, :dgemm_64_)
const SGEMM_OpenBLAS = Libdl.dlsym(libOpenBLAS, :sgemm_64_)
const DGEMV_OpenBLAS = Libdl.dlsym(libOpenBLAS, :dgemv_64_)
const OPENBLAS_SET_NUM_THREADS = Libdl.dlsym(libOpenBLAS, :openblas_set_num_threads64_)

istransposed(x) = 'N'
istransposed(x::Adjoint{<:Real}) = 'T'
istransposed(x::Adjoint) = 'C'
istransposed(x::Transpose) = 'T'
for (lib, f) ∈ [(:GEMM_MKL, :gemmmkl!), (:GEMM_OpenBLAS, :gemmopenblas!)]
  for (T, prefix) ∈ [(Float32, :S), (Float64, :D)]
    fm = Symbol(prefix, lib)
    @eval begin
      function $f(C::AbstractMatrix{$T}, A::AbstractMatrix{$T}, B::AbstractMatrix{$T})
        transA = istransposed(A)
        transB = istransposed(B)
        M, N = size(C)
        K = size(B, 1)
        pA = parent(A)
        pB = parent(B)
        ldA = stride(pA, 2)
        ldB = stride(pB, 2)
        ldC = stride(C, 2)
        α = one($T)
        β = zero($T)
        ccall(
          $fm,
          Cvoid,
          (
            Ref{UInt8},
            Ref{UInt8},
            Ref{Int64},
            Ref{Int64},
            Ref{Int64},
            Ref{$T},
            Ref{$T},
            Ref{Int64},
            Ref{$T},
            Ref{Int64},
            Ref{$T},
            Ref{$T},
            Ref{Int64},
          ),
          transA,
          transB,
          M,
          N,
          K,
          α,
          pA,
          ldA,
          pB,
          ldB,
          β,
          C,
          ldC,
        )
      end
    end
  end
end
if MKL_BENCH
  mkl_set_num_threads(N::Integer) = ccall(MKL_SET_NUM_THREADS, Cvoid, (Int32,), N % Int32)
  mkl_set_num_threads(1)
end
openblas_set_num_threads(N::Integer) = ccall(OPENBLAS_SET_NUM_THREADS, Cvoid, (Int64,), N)
openblas_set_num_threads(1)

function dgemvmkl!(
  y::AbstractVector{Float64},
  A::AbstractMatrix{Float64},
  x::AbstractVector{Float64},
  α = 1.0,
  β = 0.0,
)
  transA = istransposed(A)
  pA = parent(A)
  M, N = size(pA)
  ldA = stride(pA, 2)
  incx = LinearAlgebra.stride1(x)
  incy = LinearAlgebra.stride1(y)
  ccall(
    DGEMV_MKL,
    Cvoid,
    (
      Ref{UInt8},
      Ref{Int64},
      Ref{Int64},
      Ref{Float64},
      Ref{Float64},
      Ref{Int64},
      Ref{Float64},
      Ref{Int64},
      Ref{Float64},
      Ref{Float64},
      Ref{Int64},
    ),
    transA,
    M,
    N,
    α,
    pA,
    ldA,
    x,
    incx,
    β,
    y,
    incy,
  )
end
function dgemvopenblas!(
  y::AbstractVector{Float64},
  A::AbstractMatrix{Float64},
  x::AbstractVector{Float64},
)
  transA = istransposed(A)
  pA = parent(A)
  M, N = size(pA)
  ldA = stride(pA, 2)
  incx = LinearAlgebra.stride1(x)
  incy = LinearAlgebra.stride1(y)
  α = 1.0
  β = 0.0
  ccall(
    DGEMV_OpenBLAS,
    Cvoid,
    (
      Ref{UInt8},
      Ref{Int64},
      Ref{Int64},
      Ref{Float64},
      Ref{Float64},
      Ref{Int64},
      Ref{Float64},
      Ref{Int64},
      Ref{Float64},
      Ref{Float64},
      Ref{Int64},
    ),
    transA,
    M,
    N,
    α,
    pA,
    ldA,
    x,
    incx,
    β,
    y,
    incy,
  )
end

for (prefix, Eshared) ∈ ((Symbol(""), LIBEIGENTEST), (:i, LIBIEIGENTEST))
  @eval function $(Symbol(prefix, :egemm!))(C, A, B)
    M, N = size(C)
    K = size(B, 1)
    ccall(
      (:AmulB, $Eshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
      C,
      A,
      B,
      M,
      K,
      N,
    )
  end
  let (p, s) = (:e, Eshared)
    @eval function $(Symbol(prefix, p, :gemm!))(C, A::Adjoint, B)
      M, N = size(C)
      K = size(B, 1)
      ccall(
        (:AtmulB, $s),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
        C,
        parent(A),
        B,
        M,
        K,
        N,
      )
    end
    @eval function $(Symbol(prefix, p, :gemm!))(C, A, B::Adjoint)
      M, N = size(C)
      K = size(B, 1)
      ccall(
        (:AmulBt, $s),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
        C,
        A,
        parent(B),
        M,
        K,
        N,
      )
    end
    @eval function $(Symbol(prefix, p, :gemm!))(C, A::Adjoint, B::Adjoint)
      M, N = size(C)
      K = size(B, 1)
      ccall(
        (:AtmulBt, $s),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
        C,
        parent(A),
        parent(B),
        M,
        K,
        N,
      )
    end
    @eval function $(Symbol(prefix, p, :dot))(a, b)
      N = length(a)
      ccall((:dot, $s), Float64, (Ptr{Float64}, Ptr{Float64}, Clong), a, b, N)
    end
    @eval function $(Symbol(prefix, p, :gemv!))(y, A, x)
      M, K = size(A)
      ccall(
        ($(QuoteNode(p === :c ? :gemv : :Amulvb)), $s),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        y,
        A,
        x,
        M,
        K,
      )
    end
    @eval function $(Symbol(prefix, p, :selfdot))(a)
      N = length(a)
      ccall((:selfdot, $s), Float64, (Ptr{Float64}, Clong), a, N)
    end
    @eval function $(Symbol(prefix, p, :dot3))(x, A, y)
      M, N = size(A)
      ccall(
        (:dot3, $s),
        Float64,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        x,
        A,
        y,
        M,
        N,
      )
    end
    @eval function $(Symbol(prefix, p, :gemv!))(y, A::Adjoint, x)
      M, K = size(A)
      ccall(
        (:Atmulvb, $s),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        y,
        parent(A),
        x,
        M,
        K,
      )
    end
    @eval function $(Symbol(prefix, p, :aplusBc!))(D, a, B, c)
      M, K = size(B)
      ccall(
        (:aplusBc, $s),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        D,
        a,
        B,
        c,
        M,
        K,
      )
    end
    @eval function $(Symbol(prefix, p, :OLSlp))(y, X, β)
      N, P = size(X)
      ccall(
        (:OLSlp, $s),
        Float64,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        y,
        X,
        β,
        N,
        P,
      )
    end
    @eval function $(Symbol(prefix, p, :AplusAt!))(B, A)
      N = size(B, 1)
      ccall((:AplusAt, $s), Cvoid, (Ptr{Float64}, Ptr{Float64}, Clong), B, A, N)
    end
    @eval function $(Symbol(prefix, p, :logdettriangle))(
      T::Union{LowerTriangular,UpperTriangular},
    )
      N = size(T, 1)
      Tp = parent(T)
      ccall((:logdettriangle, $s), Float64, (Ptr{Float64}, Clong), Tp, N)
    end
  end
end
funcs_to_define = if INTEL_BENCH
  ((Symbol(""), LIBCTEST, LIBFTEST), (:i, LIBICTEST, LIBIFTEST))
else
  ((Symbol(""), LIBCTEST, LIBFTEST),)
end
for (prefix, Cshared, Fshared) ∈ funcs_to_define
  for order ∈ (:kmn, :knm, :mkn, :mnk, :nkm, :nmk)
    gemm = Symbol(:gemm_, order)
    @eval function $(Symbol(prefix, :c, gemm, :!))(C, A, B)
      M, N = size(C)
      K = size(B, 1)
      ccall(
        ($(QuoteNode(gemm)), $Cshared),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
        C,
        A,
        B,
        M,
        K,
        N,
      )
    end
    @eval function $(Symbol(prefix, :f, gemm, :!))(C, A, B)
      M, N = size(C)
      K = size(B, 1)
      ccall(
        ($(QuoteNode(gemm)), $Fshared),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
        C,
        A,
        B,
        Ref(M),
        Ref(K),
        Ref(N),
      )
    end
  end
  @eval $(Symbol(prefix, :cgemm!))(C, A, B) = $(Symbol(prefix, :cgemm_nkm!))(C, A, B)
  @eval $(Symbol(prefix, :fgemm!))(C, A, B) = $(Symbol(prefix, :fgemm_nkm!))(C, A, B)
  @eval function $(Symbol(prefix, :fgemm_builtin!))(C, A, B)
    M, N = size(C)
    K = size(B, 1)
    ccall(
      (:gemmbuiltin, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
      C,
      A,
      B,
      Ref(M),
      Ref(K),
      Ref(N),
    )
  end
  let (p, s) = (:c, Cshared)
    @eval function $(Symbol(prefix, p, :gemm!))(C, A::Adjoint, B)
      M, N = size(C)
      K = size(B, 1)
      ccall(
        (:AtmulB, $s),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
        C,
        parent(A),
        B,
        M,
        K,
        N,
      )
    end
  end
  @eval function $(Symbol(prefix, :fgemm!))(C, A::Adjoint, B)
    M, N = size(C)
    K = size(B, 1)
    ccall(
      (:AtmulB, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
      C,
      parent(A),
      B,
      Ref(M),
      Ref(K),
      Ref(N),
    )
  end
  @eval function $(Symbol(prefix, :fgemm_builtin!))(C, A::Adjoint, B)
    M, N = size(C)
    K = size(B, 1)
    ccall(
      (:AtmulBbuiltin, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
      C,
      parent(A),
      B,
      Ref(M),
      Ref(K),
      Ref(N),
    )
  end
  let (p, s) = (:c, Cshared)# (:e,Eshared)]
    @eval function $(Symbol(prefix, p, :gemm!))(C, A, B::Adjoint)
      M, N = size(C)
      K = size(B, 1)
      ccall(
        (:AmulBt, $s),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
        C,
        A,
        parent(B),
        M,
        K,
        N,
      )
    end
  end
  @eval function $(Symbol(prefix, :fgemm!))(C, A, B::Adjoint)
    M, N = size(C)
    K = size(B, 1)
    ccall(
      (:AmulBt, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
      C,
      A,
      parent(B),
      Ref(M),
      Ref(K),
      Ref(N),
    )
  end
  @eval function $(Symbol(prefix, :fgemm_builtin!))(C, A, B::Adjoint)
    M, N = size(C)
    K = size(B, 1)
    ccall(
      (:AmulBtbuiltin, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
      C,
      A,
      parent(B),
      Ref(M),
      Ref(K),
      Ref(N),
    )
  end
  @eval function $(Symbol(prefix, :fgemm!))(C, A::Adjoint, B::Adjoint)
    M, N = size(C)
    K = size(B, 1)
    ccall(
      (:AtmulBt, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
      C,
      parent(A),
      parent(B),
      Ref(M),
      Ref(K),
      Ref(N),
    )
  end
  @eval function $(Symbol(prefix, :fgemm_builtin!))(C, A::Adjoint, B::Adjoint)
    M, N = size(C)
    K = size(B, 1)
    ccall(
      (:AtmulBtbuiltin, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
      C,
      parent(A),
      parent(B),
      Ref(M),
      Ref(K),
      Ref(N),
    )
  end
  @eval function $(Symbol(prefix, :fdot))(a, b)
    N = length(a)
    ccall((:dot, $Fshared), Float64, (Ptr{Float64}, Ptr{Float64}, Ref{Clong}), a, b, Ref(N))
  end

  @eval function $(Symbol(prefix, :fselfdot))(a)
    N = length(a)
    ccall((:selfdot, $Fshared), Float64, (Ptr{Float64}, Ref{Clong}), a, Ref(N))
  end
  @eval function $(Symbol(prefix, :fdot3))(x, A, y)
    M, N = size(A)
    ccall(
      (:dot3, $Fshared),
      Float64,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
      x,
      A,
      y,
      Ref(M),
      Ref(N),
    )
  end
  @eval function $(Symbol(prefix, :fgemv!))(y, A, x)
    M, K = size(A)
    ccall(
      (:gemv, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
      y,
      A,
      x,
      Ref(M),
      Ref(K),
    )
  end
  @eval function $(Symbol(prefix, :fgemv_builtin!))(y, A, x)
    M, K = size(A)
    ccall(
      (:gemvbuiltin, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
      y,
      A,
      x,
      Ref(M),
      Ref(K),
    )
  end
  @eval function $(Symbol(prefix, :fgemv!))(y, A::Adjoint, x)
    M, K = size(A)
    ccall(
      (:Atmulvb, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
      y,
      parent(A),
      x,
      Ref(M),
      Ref(K),
    )
  end
  @eval function $(Symbol(prefix, :fgemv_builtin!))(y, A::Adjoint, x)
    M, K = size(A)
    ccall(
      (:Atmulvbbuiltin, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
      y,
      parent(A),
      x,
      Ref(M),
      Ref(K),
    )
  end
  let (p, s) = (:c, Cshared)# (:e,Eshared)]
    @eval function $(Symbol(prefix, p, :gemm!))(C, A::Adjoint, B::Adjoint)
      M, N = size(C)
      K = size(B, 1)
      ccall(
        (:AtmulBt, $s),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
        C,
        parent(A),
        parent(B),
        M,
        K,
        N,
      )
    end
    @eval function $(Symbol(prefix, p, :dot))(a, b)
      N = length(a)
      ccall((:dot, $s), Float64, (Ptr{Float64}, Ptr{Float64}, Clong), a, b, N)
    end
    @eval function $(Symbol(prefix, p, :gemv!))(y, A, x)
      M, K = size(A)
      ccall(
        ($(QuoteNode(p === :c ? :gemv : :Amulvb)), $s),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        y,
        A,
        x,
        M,
        K,
      )
    end
    @eval function $(Symbol(prefix, p, :selfdot))(a)
      N = length(a)
      ccall((:selfdot, $s), Float64, (Ptr{Float64}, Clong), a, N)
    end
    @eval function $(Symbol(prefix, p, :dot3))(x, A, y)
      M, N = size(A)
      ccall(
        (:dot3, $s),
        Float64,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        x,
        A,
        y,
        M,
        N,
      )
    end
    @eval function $(Symbol(prefix, p, :gemv!))(y, A::Adjoint, x)
      M, K = size(A)
      ccall(
        (:Atmulvb, $s),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        y,
        parent(A),
        x,
        M,
        K,
      )
    end
    @eval function $(Symbol(prefix, p, :aplusBc!))(D, a, B, c)
      M, K = size(B)
      ccall(
        (:aplusBc, $s),
        Cvoid,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        D,
        a,
        B,
        c,
        M,
        K,
      )
    end
    @eval function $(Symbol(prefix, p, :OLSlp))(y, X, β)
      N, P = size(X)
      ccall(
        (:OLSlp, $s),
        Float64,
        (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
        y,
        X,
        β,
        N,
        P,
      )
    end
    @eval function $(Symbol(prefix, p, :AplusAt!))(B, A)
      N = size(B, 1)
      ccall((:AplusAt, $s), Cvoid, (Ptr{Float64}, Ptr{Float64}, Clong), B, A, N)
    end
    @eval function $(Symbol(prefix, p, :logdettriangle))(
      T::Union{LowerTriangular,UpperTriangular},
    )
      N = size(T, 1)
      Tp = parent(T)
      ccall((:logdettriangle, $s), Float64, (Ptr{Float64}, Clong), Tp, N)
    end
  end
  @eval function $(Symbol(prefix, :faplusBc!))(D, a, B, c)
    M, K = size(B)
    ccall(
      (:aplusBc, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
      D,
      a,
      B,
      c,
      Ref(M),
      Ref(K),
    )
  end
  @eval function $(Symbol(prefix, :fOLSlp))(y, X, β)
    N, P = size(X)
    ccall(
      (:OLSlp, $Fshared),
      Float64,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
      y,
      X,
      β,
      Ref(N),
      Ref(P),
    )
  end
  @eval function $(Symbol(prefix, :cvexp!))(b, a)
    N = length(b)
    ccall((:vexp, $Cshared), Cvoid, (Ptr{Float64}, Ptr{Float64}, Clong), b, a, N)
  end
  @eval function $(Symbol(prefix, :fvexp!))(b, a)
    N = length(b)
    ccall((:vexp, $Fshared), Cvoid, (Ptr{Float64}, Ptr{Float64}, Ref{Clong}), b, a, Ref(N))
  end
  @eval function $(Symbol(prefix, :fvexpsum))(a)
    N = length(a)
    s = Ref{Float64}()
    ccall((:svexp, $Fshared), Cvoid, (Ref{Float64}, Ptr{Float64}, Ref{Clong}), s, a, Ref(N))
    s[]
  end
  @eval function $(Symbol(prefix, :fAplusAt!))(B, A)
    N = size(B, 1)
    ccall(
      (:AplusAt, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ref{Clong}),
      B,
      A,
      Ref(N),
    )
  end
  @eval function $(Symbol(prefix, :fAplusAt_builtin!))(B, A)
    N = size(B, 1)
    ccall(
      (:AplusAtbuiltin, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ref{Clong}),
      B,
      A,
      Ref(N),
    )
  end

  @eval function $(Symbol(prefix, :crandomaccess))(P, basis, coefs)
    A, C = size(P)
    ccall(
      (:randomaccess, $Cshared),
      Float64,
      (Ptr{Float64}, Ptr{Clong}, Ptr{Float64}, Clong, Clong),
      P,
      basis,
      coefs,
      A,
      C,
    )
  end
  @eval function $(Symbol(prefix, :frandomaccess))(P, basis, coefs)
    A, C = size(P)
    ccall(
      (:randomaccess, $Fshared),
      Float64,
      (Ptr{Float64}, Ptr{Clong}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
      P,
      basis,
      coefs,
      Ref(A),
      Ref(C),
    )
  end
  @eval function $(Symbol(prefix, :flogdettriangle))(
    T::Union{LowerTriangular,UpperTriangular},
  )
    N = size(T, 1)
    Tp = parent(T)
    ccall((:logdettriangle, $Fshared), Float64, (Ptr{Float64}, Ref{Clong}), Tp, Ref(N))
  end
  @eval function $(Symbol(prefix, :cfilter2d!))(
    B::OffsetArray,
    A::AbstractArray,
    K::OffsetArray,
  )
    Ma, Na = size(A)
    offset = first(B.offsets)
    ccall(
      (:filter2d, $Cshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong, Clong),
      parent(B),
      A,
      parent(K),
      Ma,
      Na,
      offset,
    )
  end
  @eval function $(Symbol(prefix, :ffilter2d!))(
    B::OffsetArray,
    A::AbstractArray,
    K::OffsetArray,
  )
    Ma, Na = size(A)
    offset = first(B.offsets)
    ccall(
      (:filter2d, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}, Ref{Clong}),
      parent(B),
      A,
      parent(K),
      Ref(Ma),
      Ref(Na),
      Ref(offset),
    )
  end
  @eval function $(Symbol(prefix, :cfilter2d!))(
    B::OffsetArray,
    A::AbstractArray,
    K::SizedOffsetMatrix{Float64,-1,1,-1,1},
  )
    Ma, Na = size(A)
    ccall(
      (:filter2d3x3, $Cshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
      parent(B),
      A,
      K,
      Ma,
      Na,
    )
  end
  @eval function $(Symbol(prefix, :ffilter2d!))(
    B::OffsetArray,
    A::AbstractArray,
    K::SizedOffsetMatrix{Float64,-1,1,-1,1},
  )
    Ma, Na = size(A)
    ccall(
      (:filter2d3x3, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
      parent(B),
      A,
      K,
      Ref(Ma),
      Ref(Na),
    )
  end
  @eval function $(Symbol(prefix, :cfilter2dunrolled!))(
    B::OffsetArray,
    A::AbstractArray,
    K::SizedOffsetMatrix{Float64,-1,1,-1,1},
  )
    Ma, Na = size(A)
    ccall(
      (:filter2d3x3unrolled, $Cshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Clong, Clong),
      parent(B),
      A,
      K,
      Ma,
      Na,
    )
  end
  @eval function $(Symbol(prefix, :ffilter2dunrolled!))(
    B::OffsetArray,
    A::AbstractArray,
    K::SizedOffsetMatrix{Float64,-1,1,-1,1},
  )
    Ma, Na = size(A)
    ccall(
      (:filter2d3x3unrolled, $Fshared),
      Cvoid,
      (Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ref{Clong}, Ref{Clong}),
      parent(B),
      A,
      K,
      Ref(Ma),
      Ref(Na),
    )
  end

end


text:
"""
Draw a network with timedepth 1 and the greedy-paths chosen between 3 end user pairs.
"""

timedepth = 3
grid_size = 10

# Index offset for asynchronus nodes
off = grid_size^2 * timedepth
# Choose asynchronus endusers
userpairs = [(1 + off, 100 + off), (50 + off, 81 + off), (87 + off, 22 + off)]

net = GridNetwork(grid_size, grid_size)
T = QuNet.TemporalGraph(net, timedepth, memory_prob=1.0)
QuNet.add_async_nodes!(T, userpairs)
T_copy = deepcopy(T)
user_paths, dum1, dum2 = QuNet.greedy_multi_path!(T_copy, QuNet.purify, userpairs)
QuNet.plot_network(T.graph["Z"], user_paths, T.locs_x, T.locs_y)

ext: jl

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
using Printf

function romanencode(n::Integer)
    if n < 1 || n > 4999 throw(DomainError()) end

    DR = [["I", "X", "C", "M"] ["V", "L", "D", "MMM"]]
    rnum = ""
    for (omag, d) in enumerate(digits(n))
        if d == 0
            omr = ""
        elseif d <  4
            omr = DR[omag, 1] ^ d
        elseif d == 4
            omr = DR[omag, 1] * DR[omag, 2]
        elseif d == 5
            omr = DR[omag, 2]
        elseif d <  9
            omr = DR[omag, 2] * DR[omag, 1] ^ (d - 5)
        else
            omr = DR[omag, 1] * DR[omag + 1, 1]
        end
        rnum = omr * rnum
    end
    return rnum
end

testcases = [1990, 2008, 1668]
append!(testcases, rand(1:4999, 12))
testcases = unique(testcases)

println("Test romanencode, arabic => roman:")
for n in testcases
    @printf("%-4i => %s\n", n, romanencode(n))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module ZeroRing

using LinearAlgebra
using Random

export zeroelem, ZeroElem

struct ZeroElem <: Number end

const zeroelem = ZeroElem()

Base.:(==)(::ZeroElem, ::ZeroElem) = true
Base.:(<)(::ZeroElem, ::ZeroElem) = false
Base.hash(::ZeroElem, h::UInt) = hash(UInt(0x56212a80), h)

Base.rand(rng::AbstractRNG, ::Random.SamplerType{ZeroElem}) = zeroelem

Base.:+(::ZeroElem) = zeroelem
Base.:-(::ZeroElem) = zeroelem

Base.:+(::ZeroElem, ::ZeroElem) = zeroelem
Base.:-(::ZeroElem, ::ZeroElem) = zeroelem

Base.inv(::ZeroElem) = zeroelem

Base.:*(::ZeroElem, ::ZeroElem) = zeroelem
Base.:/(::ZeroElem, ::ZeroElem) = zeroelem
Base.:\(::ZeroElem, ::ZeroElem) = zeroelem

Base.:^(::ZeroElem, ::ZeroElem) = zeroelem

Base.copy(::ZeroElem) = zeroelem

LinearAlgebra.adjoint(::ZeroElem) = zeroelem
LinearAlgebra.transpose(::ZeroElem) = zeroelem

end


text:
module Aref_Lagrangian_Descriptor
export Aref_Lagrangian_Descriptor_Function
export final_T

using Roots
using DifferentialEquations
using Printf
include("leap_frog_definitions.jl")

barrier=5
condition_escape(u, t, integrator) =  maximum([abs(u[1]),abs(u[2]),abs(u[3]),abs(u[4])]) > barrier
affect_stop!(integrator) = terminate!(integrator)
callback_escape = DiscreteCallback(condition_escape, affect_stop!)
function condition_hits_PSS_q1(u, t, integrator) # Event when event_f(u,t) == 0
    u[1]
end
cb=ContinuousCallback(condition_hits_PSS_q1, affect_stop!, rootfind=true)

function Aref_Lagrangian_Descriptor_Function(mesh, H,  t_end)
    Q=mesh[2]
    P=mesh[1]
     #for a given Q,P,H with X=0
    Y=Yfind_Aref(Q,P,H)
    if ~isempty(Y)
        u0=zeros(5)
        u0[1]=0 #X
        u0[2]=P #P
        u0[3]=Q #Q
        u0[4]=Y #Y
        u0[5]=0 # counter variable

        #constructor for ODE

        prob_b = ODEProblem(Eq_of_M_Lagrangian_Descriptors_Backwards,u0,(0., t_end))
        prob_f = ODEProblem(Eq_of_M_Lagrangian_Descriptors,u0,(0., t_end))

        sol_f=solve(prob_f, Tsit5(),maxiters=1e20,reltol=1e-6,abstol=1e-9,save_idxs = [5],save_every_step=false, save_end=true, dense=false)#,abstol=1e-9)
        # sol_b=solve(prob_b, Tsit5(),maxiters=1e20,reltol=1e-6,abstol=1e-9,save_idxs = [5],save_every_step=false, save_end=true, dense=false)#,abstol=1e-9)

        # sol_f=solve(prob_f,  Vern9(),maxiters=1e20,reltol=1e-14,abstol=1e-14,save_idxs = [5],save_every_step=false, save_end=true, dense=false)#,abstol=1e-9)
        # sol_b=solve(prob_b,  Vern9(),maxiters=1e20,reltol=1e-14,abstol=1e-14,save_idxs = [5],save_every_step=false, save_end=true, dense=false)#,abstol
        # sol_b=solve(prob_b, RK4(),maxiters=1e20,reltol=1e-6,abstol=1e-8)

        # sol_f=solve(prob_f, RK4(),maxiters=1e20,reltol=1e-8,abstol=1e-10)
        # sol_b=solve(prob_b, RK4(),maxiters=1e20,reltol=1e-8,abstol=1e-10)

        # sol_b=solve(prob_b, RK4(),maxiters=1e20)
        # sol_f=solve(prob_f, RK4(),maxiters=1e20)

        # sol_b=solve(prob_b, RK4(),maxiters=1e20,callback=callback_escape)
        # sol_f=solve(prob_f, RK4(),maxiters=1e20,callback=callback_escape)

        # sol_f=solve(prob_f, Vern9(),maxiters=1e20,reltol=1e-13,abstol=1e-15)
        # sol_b=solve(prob_b, Vern9(),maxiters=1e20,reltol=1e-13,abstol=1e-15)

        #find final value to check if energy was conserved
        # uf=zeros(5)
        # uf[1]=sol[1,end] #X
        # uf[2]=sol[2,end] #P
        # uf[3]=sol[3,end] #Q
        # uf[4]=sol[4,end] #Y
        # dH=abs(H_test(uf)-H)

        LD=sol_f[1,end]#+sol_b[1,end]
        # if maximum([abs(sol_f[1,end]),abs(sol_f[2,end]),abs(sol_f[3,)nd]),abs(sol_f[4,end]),abs(sol_b[1,end]),abs(sol_b[2,end]),abs(sol_b[3,end]),abs(sol_b[4,end])])<barrier
        return LD
    else
        return NaN
    end
end

function final_T(H)
    Q=1e-8
    P=1e-8
    Y=Yfind_Aref(Q,P,H)
    if ~isempty(Y)
        u0=zeros(5)
        u0[1]=0 #X
        u0[2]=P #P
        u0[3]=Q #Q
        u0[4]=Y #Y
        u0[5]=0 # counter variable

        #constructor for ODE

        prob = ODEProblem(Eq_of_M_Lagrangian_Descriptors_Backwards,u0,(0., 20))
        sol=solve(prob, Vern9(),maxiters=1e20,reltol=1e-13,abstol=1e-15,callback=cb)

        return sol.t[end]
    end
end

end


text:
# ============================================================================ #
# Triangle
# ============================================================================ #
abstract type AbstractTriangle end

struct Triangle{PointType <: AbstractPoint} <: AbstractTriangle
    p1::PointType
    p2::PointType
    p3::PointType
end

Triangle(x::A, y::A) where {A<:AbstractArray} = Triangle( Point(x[1], y[1]), Point(x[2], y[2]),  Point(x[3], y[3]) )
Triangle(p1::A, p2::A, p3::A) where {A<:AbstractArray} = Triangle( Point(p1[1], p1[2]), Point(p2[1], p2[2]),  Point(p3[1], p3[2]) )

mutable struct MutableTriangle{T} <: AbstractTriangle
    p1::MutablePoint{T}
    p2::MutablePoint{T}
    p3::MutablePoint{T}
end

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

Base.:(==)(t1::AbstractTriangle, t2::AbstractTriangle) = t1.p1 == t2.p1 && t1.p2 == t2.p2  && t1.p3 == t2.p3

function Base.getindex(triangle::AbstractTriangle, index::Int)
    if index == 1
        return triangle.p1
    elseif index == 2
        return triangle.p2
    elseif index == 3
        return triangle.p3
    else
        return NaN
    end
end

function shift(triangle::T, v) where {T <: AbstractTriangle}
    return T( triangle.p1 + v, triangle.p2 + v, triangle.p3 + v )
end
function scale(line::T, factor) where {T <: AbstractTriangle}
    mid = mid_point(triangle)

    vm1 = triangle.p1 - mid
    vm2 = triangle.p2 - mid
    vm3 = triangle.p3 - mid

    return T( vm1, vm2, vm3 )
end

function circumference(triangle::AbstractTriangle)
    v12 = triangle.p2 - triangle.p1
    v23 = triangle.p3 - triangle.p2
    v31 = triangle.p1 - triangle.p3
    return norm(v12) + norm(v23) + norm(v31)
end

function area(triangle::AbstractTriangle)
    v12 = triangle.p2 - triangle.p1
    v13 = triangle.p3 - triangle.p1
    return cross2(v12, v13)/2
end

function centroid(triangle::AbstractTriangle)
    return (triangle.p1 + triangle.p2 + triangle.p3)./3
end
function mid_point(triangle)
    return centroid(triangle)
end

function is_in_triangle(triangle::AbstractTriangle, point::AbstractPoint)::Bool
    line_12 = Line(triangle.p1, triangle.p2)
    line_23 = Line(triangle.p2, triangle.p3)
    line_31 = Line(triangle.p3, triangle.p1)

    # b12 = is_left_on(line_12, point)
    b12 = is_left_on_inlined(line_12, point)
    if !b12
        return false
    end

    # b23 = is_left_on(line_23, point)
    b23 = is_left_on_inlined(line_23, point)
    if !b23
        return false
    end

    # b31 = is_left_on(line_31, point)
    b31 = is_left_on_inlined(line_31, point)
    return b31
end

function is_in_triangle(p1, p2, p3, query_point)::Bool

    # b12 = is_left(p1, p2, query_point)
    b12 = is_left_on_inlined(p1, p2, query_point)
    if !b12
        return false
    end

    # b23 = is_left(p2, p3, query_point)
    b23 = is_left_on_inlined(p2, p3, query_point)
    if !b23
        return false
    end

    # b31 = is_left(p3, p1, query_point)
    b31 = is_left_on_inlined(p3, p1, query_point)
    return b31
end

function side_lengths(triangle::AbstractTriangle)
    v12 = triangle.p2 - triangle.p1
    v23 = triangle.p3 - triangle.p2
    v31 = triangle.p1 - triangle.p3
    a = norm(v12)
    b = norm(v23)
    c = norm(v31)

    return (a, b, c)
end

function in_circle(triangle::AbstractTriangle)

    (a, b, c) = side_lengths(triangle)
    s = 0.5 * (a + b + c)
    r = sqrt( (s-a) * (s-b) * (s-c) / s )

    A = triangle.p3
    B = triangle.p1
    C = triangle.p2

    center = (a * A + b * B + c * C) / (a + b + c)

    return Circle(center, r)
end

function out_circle(triangle::AbstractTriangle)
    (a, b, c) = side_lengths(triangle)
    Area = heron_formula(a, b, c)
    R = a * b * c / (4 * Area)
    asq = a * a
    bsq = b * b
    csq = c * c

    w_a = asq * (bsq + csq - asq)
    w_b = bsq * (asq + csq - bsq)
    w_c = csq * (asq + bsq - csq)

    A = triangle.p3
    B = triangle.p1
    C = triangle.p2

    center = (w_a * A + w_b * B + w_c * C ) / (16 * Area * Area)

    return Circle(center, R)
end

function out_circle_radius(triangle::AbstractTriangle)
    (a, b, c) = side_lengths(triangle)
    A = heron_formula(a, b, c)
    R = a * b * c / (4 * A)
end

function heron_formula(a::T, b::T, c::T) where {T <: Real}
    s = 0.5 * (a + b + c)
    A = sqrt( s * (s-a) * (s-b) * (s-c) )
    return A
end
heron_formula(a, b, c) = heron_formula( promote(a, b, c)... )

function discretize_for_plot(t::AbstractTriangle)
    x = [t.p1.x, t.p2.x, t.p3.x, t.p1.x]
    y = [t.p1.y, t.p2.y, t.p3.y, t.p1.y]
    return (x, y)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using MarkovTransitionMatrices
using Test
using Plots
gr()

nP = 31
nC = 37
n = nP * nC

σp = [0x1.269c263357e05p-4, 0x1.9bec08eb57abfp-5]
σc = [0x1.2c9a3a5b8b5bap-5, 0x1.09d5d411d8a33p-7]

Πk = [0x1.d7e13e2f38db1p-1 0x1.40f60e8638b53p-4;
      0x1.098244df0fddfp-4 0x1.decfb7641df78p-1]
all(sum(Πk, dims=2) .≈ 1.0) || throw(error("each row of π must sum to 1"))

K = size(Πk,2)
extrema_p = [0.8776572, 2.485073]
extrema_c = [0.2437302, 1.529529]

# state-spaces
pspace = range(extrema_p[1]-log(1.5), stop=extrema_p[2] + log(1.5), length=nP)
cspace = range(extrema_c[1]-log(1.5), stop=extrema_c[2] + log(1.5), length=nC)

# make deviations
zrandwalk(x::Real, st::Real, σ::Real) = (x - st) / σ

# allocate giant transition matrix
P = Matrix{eltype(pspace)}(n*K,n*K)

# match moments for logp
Pp1, JN, Λ, L_p1, approxErr = discreteNormalApprox(pspace, pspace, (x::Real,st::Real) -> zrandwalk(x,st,σp[1]), 15)
Pp2, JN, Λ, L_p2, approxErr = discreteNormalApprox(pspace, pspace, (x::Real,st::Real) -> zrandwalk(x,st,σp[2]), 15)
plot(pspace, [L_p1,L_p2], yticks = 1:15, labels=["Hi vol", "Lo vol"], xlabel="Log p", ylabel="Moments matched")

# match moments for logc
Pc1, JN, Λ, L_c1, approxErr = discreteNormalApprox(cspace, cspace, (x::Real,st::Real) -> zrandwalk(x,st,σc[1]), 13)
Pc2, JN, Λ, L_c2, approxErr = discreteNormalApprox(cspace, cspace, (x::Real,st::Real) -> zrandwalk(x,st,σc[2]), 5)
plot(cspace, [L_c1,L_c2], yticks = 1:15, labels=["Hi vol", "Lo vol"], xlabel="Log c", ylabel="Moments matched")

# sparse versions
minp = 1e-6
sPp1 = MarkovTransitionMatrices.sparsify!(Pp1, minp)
sPp2 = MarkovTransitionMatrices.sparsify!(Pp2, minp)
sPc1 = MarkovTransitionMatrices.sparsify!(Pc1, minp)
sPc2 = MarkovTransitionMatrices.sparsify!(Pc2, minp)

# do the big matrix
P12 = (kron(sPc1, sPp1), kron(sPc2, sPp2))
for r2 in 1:K
  for r1 in 1:K
    P[(r1-1)*n+1:r1*n,  (r2-1)*n+1:r2*n] .= Πk[r1, r2] .* P12[r2]
  end
end

sP = sparse(P)


text:
function ljpotential(xij, qi, qj, ti, tj, mu, eta, kappa)

A = mu[1];
B = mu[2];

r2 = xij[1,:].*xij[1,:] .+ xij[2,:].*xij[2,:] .+ xij[3,:].*xij[3,:];
#r1 = sqrt.(r2);
r4 = r2.*r2;
r6 = r2.*r4;
r12 = r6.*r6;
eij = (A./r12 - B./r6);

r8 = r6.*r2;
r14 = r12.*r2;
fij = xij*1.0;
for i = 1:3
    fij[i,:] = xij[i,:].*((6*B)./r8 - (12*A)./r14);
end
# r7 = r6.*r2;
# r13 = r12.*r1;
# dr = xij./r1;
# fij = dr.*((6*B)./r7 - (12*A)./r13);

return eij, fij

end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# test with the network in HGTinconsistency
# Claudia April 2015

include("../src/types.jl")
include("../src/functions.jl")

truenetwork = "((((1,2),((3,4))#H1),(#H1,5)),6);"
treefile = "1.ms"

d = readInputData(treefile, "HGTtableCF.txt"); #will get list of all quartets (allQuartets.txt), and obsCF (HGTtableCF.txt)
length(d.quartet)

# for expCF
net = readTopologyUpdate(truenetwork);
printEdges(net)
extractQuartet!(net,d) #extract qnet from that list and calculate expCF
df2 = writeExpCF(d.quartet)
writetable("HGT_truenet_expCF.csv",df2)

# compare table with expCF and obsCF: very different!

# ----------------- optTopLevel with expCF and starting tree 1_astral.out, no branches updated --------
include("../src/types.jl")
include("../src/functions.jl")

df2 = readtable("HGT_truenet_expCF.csv")
d2 = readTableCF(df2); #expCF

truenetwork = "((((1,2),((3,4))#H1),(#H1,5)),6);"
net = readTopologyUpdate(truenetwork);
printEdges(net)
@time optBL!(net,d2) #loglik~1e.-15 in 0.17secs

currT0 = readTopologyUpdate("1_astral.out");
printEdges(currT0)
srand(1234) #found right network in 135secs, wrong BL, right gamma (debug3hgtGood.txt)
srand(4568) # bug with movedownlevel, will leave for later (debug4hgt.txt)
srand(11233) #found right network in 119.34secs (debug4hgtBad.txt, by mistake)
currT = deepcopy(currT0);
addHybridizationUpdate!(currT); #add hybrid at random (different test would be to start with the tree)
printEdges(currT)

@time optTopLevel!(currT,d2,1)
@time optTop!(currT,d2,1)

# ----------------- optTopLevel with expCF and starting tree 1_astral.out, branches updated --------
include("../src/types.jl")
include("../src/functions.jl")

df2 = readtable("HGT_truenet_expCF.csv")
d2 = readTableCF(df2); #expCF

truenetwork = "((((1,2),((3,4))#H1),(#H1,5)),6);"
net = readTopologyUpdate(truenetwork);
printEdges(net)
@time optBL!(net,d2) #loglik~1e.-15 in 0.17secs

currT0 = readTopologyUpdate("1_astral.out");
x = updateBL!(currT0,d2)
printEdges(currT0)
srand(1234) #right network in 20secs(debug12hgt)
srand(4568) #movedownlevel: debug13hgtBad
srand(11233) #very close to right network in 135secs (debug8hgtgood.txt)
currT = deepcopy(currT0);
addHybridizationUpdate!(currT); #add hybrid at random (different test would be to start with the tree)
printEdges(currT)

@time optTopLevel!(currT,d2,1);

# ----------------- optTopLevel with obsCF and starting tree 1_astral.out, no branches updated --------
include("../src/types.jl")
include("../src/functions.jl")

df = readtable("HGTtableCF.txt") #from 1.ms
d = readTableCF(df); #obsCF

truenetwork = "((((1,2),((3,4))#H1),(#H1,5)),6);"
net = readTopologyUpdate(truenetwork);
printEdges(net)
@time optBL!(net,d,true) #loglik ~144.74 for true net

currT0 = readTopologyUpdate("1_astral.out");
printEdges(currT0)
srand(1234) #local max found debug5hgtBad
srand(4568) #local max2 found debug6hgtBad
srand(11233) #local max found debug7hgtBad
currT = deepcopy(currT0);
addHybridizationUpdate!(currT); #add hybrid at random (different test would be to start with the tree)
printEdges(currT)

@time optTopLevel!(currT,d,1);

# ----------------- optTopLevel with obsCF and starting tree 1_astral.out, branches updated --------
include("../src/types.jl")
include("../src/functions.jl")

df = readtable("HGTtableCF.txt") #from 1.ms
d = readTableCF(df); #obsCF

truenetwork = "((((1,2),((3,4))#H1),(#H1,5)),6);"
net = readTopologyUpdate(truenetwork);
printEdges(net)
@time optBL!(net,d,true) #loglik ~144.74 for true net

currT0 = readTopologyUpdate("1_astral.out");
x = updateBL!(currT0,d)
printEdges(currT0)
srand(1234) #local max found debug9hgt.txt
srand(4568) #local max2 found debug9hgt.txt
srand(11233) #local max3 found debug9hgt.txt
currT = deepcopy(currT0);
addHybridizationUpdate!(currT); #add hybrid at random (different test would be to start with the tree)
printEdges(currT)

@time optTopLevel!(currT,d,1);

# -----------------------
# bug in srand(4568)
net = readTopologyUpdate("((4,#H7:9.99670403892172::0.43454301575229803):1.5467254857425556,((6,(5)#H7:2.512064322645178::0.565456984247702):9.221085796210835,(2,1):0.38003642076628485):0.0,3);");
srand(4568)
flag = proposedTop!(:nni,net,true,0,100, zeros(Int,18), zeros(Int,6))
df = readtable("HGTtableCF.txt") #from 1.ms
d = readTableCF(df); #obsCF

qnet=QuartetNetwork(net);
extractQuartet!(net,d)
printEdges(qnet)
printNodes(qnet)
deleteLeaf!(qnet,qnet.node[1]) #leaf1
qnet0=deepcopy(qnet);
deleteLeaf!(qnet,qnet.node[2]) #leaf4


text:
const pyjliobasetype = pynew()
const pyjlbinaryiotype = pynew()
const pyjltextiotype = pynew()

pyjlio_close(io::IO) = (close(io); Py(nothing))
pyjl_handle_error_type(::typeof(pyjlio_close), io, exc) = exc isa MethodError && exc.f === close ? pybuiltins.ValueError : PyNULL

pyjlio_closed(io::IO) = Py(!isopen(io))
pyjl_handle_error_type(::typeof(pyjlio_closed), io, exc) = exc isa MethodError && exc.f === isopen ? pybuiltins.ValueError : PyNULL

pyjlio_fileno(io::IO) = Py(fd(io))
pyjl_handle_error_type(::typeof(pyjlio_fileno), io, exc) = exc isa MethodError && exc.f === fd ? pybuiltins.ValueError : PyNULL

pyjlio_flush(io::IO) = (flush(io); Py(nothing))
pyjl_handle_error_type(::typeof(pyjlio_flush), io, exc) = exc isa MethodError && exc.f === flush ? pybuiltins.ValueError : PyNULL

pyjlio_isatty(io::IO) = Py(io isa Base.TTY)

pyjlio_readable(io::IO) = Py(isreadable(io))
pyjl_handle_error_type(::typeof(pyjlio_readable), io, exc) = exc isa MethodError && exc.f === isreadable ? pybuiltins.ValueError : PyNULL

function pyjlio_seek(io::IO, offset_::Py, whence_::Py)
    offset = pyconvertarg(Int, offset_, "offset")
    whence = pyconvertarg(Int, whence_, "whence")
    if whence == 0
        pos = offset
    elseif whence == 1
        pos = position(io) + offset
    elseif whence == 2
        seekend(io)
        pos = position(io) + offset
    else
        errset(pybuiltins.ValueError, "Argument 'whence' must be 0, 1 or 2")
        return pynew()
    end
    seek(io, pos)
    Py(position(io))
end
pyjl_handle_error_type(::typeof(pyjlio_seek), io, exc) = exc isa MethodError && (exc.f === position || exc.f === seek || exc.f === seekend) ? pybuiltins.ValueError : PyNULL

pyjlio_tell(io::IO) = Py(position(io))
pyjl_handle_error_type(::typeof(pyjlio_tell), io, exc) = exc isa MethodError && exc.f === position ? pybuiltins.ValueError : PyNULL

function pyjlio_truncate(io::IO, size_::Py)
    size = pyconvertarg(Union{Int,Nothing}, size_, "size")
    if size === nothing
        size = position(io)
    end
    truncate(io, size)
    Py(size)
end
pyjl_handle_error_type(::typeof(pyjlio_truncate), io, exc) = exc isa MethodError && (exc.f === position || exc.f === truncate) ? pybuiltins.ValueError : PyNULL

pyjlio_writable(io::IO) = Py(iswritable(io))
pyjl_handle_error_type(::typeof(pyjlio_writable), io, exc) = exc isa MethodError && exc.f === iswritable ? pybuiltins.ValueError : PyNULL

function pyjlbinaryio_read(io::IO, size_::Py)
    size = pyconvertarg(Union{Int,Nothing}, size_, "size")
    if size === nothing || size < 0
        buf = read(io)
    else
        buf = read(io, size)
    end
    pybytes(buf)
end
pyjl_handle_error_type(::typeof(pyjlbinaryio_read), io, exc) = exc isa MethodError && exc.f === read ? pybuiltins.ValueError : PyNULL

function pyjlbinaryio_readline(io::IO, size_::Py)
    size = pyconvertarg(Union{Int,Nothing}, size_, "size")
    if size === nothing
        size = -1
    end
    buf = UInt8[]
    while !eof(io) && (size < 0 || length(buf) < size)
        c = read(io, UInt8)
        push!(buf, c)
        c == 0x0A && break
    end
    pybytes(buf)
end
pyjl_handle_error_type(::typeof(pyjlbinaryio_readline), io, exc) = exc isa MethodError && exc.f === read ? pybuiltins.ValueError : PyNULL

function pyjlbinaryio_readinto(io::IO, b::Py)
    m = pybuiltins.memoryview(b)
    c = m.c_contiguous
    if !pytruth(c)
        pydel!(c)
        errset(pybuiltins.ValueError, "input buffer is not contiguous")
        return pynew()
    end
    pydel!(c)
    buf = unsafe_load(C.PyMemoryView_GET_BUFFER(getptr(m)))
    if buf.readonly != 0
        pydel!(m)
        errset(pybuiltins.ValueError, "output buffer is read-only")
        return pynew()
    end
    data = unsafe_wrap(Array, Ptr{UInt8}(buf.buf), buf.len)
    nb = readbytes!(io, data)
    pydel!(m)
    return Py(nb)
end
pyjl_handle_error_type(::typeof(pyjlbinaryio_readinto), io, exc) = exc isa MethodError && exc.f === readbytes! ? pybuiltins.ValueError : PyNULL

function pyjlbinaryio_write(io::IO, b::Py)
    m = pybuiltins.memoryview(b)
    c = m.c_contiguous
    if !pytruth(c)
        pydel!(c)
        errset(pybuiltins.ValueError, "input buffer is not contiguous")
        return pynew()
    end
    pydel!(c)
    buf = unsafe_load(C.PyMemoryView_GET_BUFFER(getptr(m)))
    data = unsafe_wrap(Array, Ptr{UInt8}(buf.buf), buf.len)
    write(io, data)
    pydel!(m)
    return Py(buf.len)
end
pyjl_handle_error_type(::typeof(pyjlbinaryio_write), io, exc) = exc isa MethodError && exc.f === write ? pybuiltins.ValueError : PyNULL

function pyjltextio_read(io::IO, size_::Py)
    size = pyconvertarg(Union{Int,Nothing}, size_, "size")
    if size === nothing
        size = -1
    end
    buf = IOBuffer()
    total = 0
    while !eof(io) && (size < 0 || total < size)
        c = read(io, Char)
        # translate "\n", "\r" and "\r\n" to "\n"
        if c == '\r'
            !eof(io) && peek(io) == 0x0A && read(io, UInt8)
            write(buf, '\n')
        else
            write(buf, c)
        end
        total += 1
    end
    Py(String(take!(buf)))
end
pyjl_handle_error_type(::typeof(pyjltextio_read), io, exc) = exc isa MethodError && exc.f === read ? pybuiltins.ValueError : PyNULL

function pyjltextio_readline(io::IO, size_::Py)
    size = pyconvertarg(Union{Int,Nothing}, size_, "size")
    if size === nothing
        size = -1
    end
    buf = IOBuffer()
    total = 0
    while !eof(io) && (size < 0 || total < size)
        c = read(io, Char)
        # translate "\n", "\r" and "\r\n" to "\n"
        if c == '\n'
            write(buf, c)
            total += 1
            break
        elseif c == '\r'
            !eof(io) && peek(io) == 0x0A && read(io, UInt8)
            write(buf, '\n')
            total += 1
            break
        else
            write(buf, c)
            total += 1
        end
    end
    Py(String(take!(buf)))
end
pyjl_handle_error_type(::typeof(pyjltextio_readline), io, exc) = exc isa MethodError && exc.f === read ? pybuiltins.ValueError : PyNULL

function pyjltextio_write(io::IO, s_::Py)
    if pyisstr(s_)
        s = pystr_asstring(s_)
        # get the line separator
        linesep_ = pyosmodule.linesep
        linesep = pystr_asstring(linesep_)
        pydel!(linesep_)
        # write the string
        # translating '\n' to os.linesep
        i = firstindex(s)
        iend = lastindex(s)
        while i ≤ iend
            j = findnext('\n', s, i)
            if j === nothing
                write(io, SubString(s, i))
                break
            else
                write(io, SubString(s, i, prevind(s, j)))
                write(io, linesep)
                i = nextind(s, j)
            end
        end
        # return number of characters written (not number of bytes)
        # TODO: is this the number of source characters, or the number of output characters?
        Py(length(s))
    else
        errset(pybuiltins.TypeError, "Argument 's' must be a 'str', got a '$(pytype(s_).__name__)'")
        pynew()
    end
end
pyjl_handle_error_type(::typeof(pyjltextio_write), io, exc) = exc isa MethodError && exc.f === write ? pybuiltins.ValueError : PyNULL

function init_jlwrap_io()
    jl = pyjuliacallmodule
    filename = "$(@__FILE__):$(1+@__LINE__)"
    pybuiltins.exec(pybuiltins.compile("""
    class IOValueBase(AnyValue):
        __slots__ = ()
        __module__ = "juliacall"
        def close(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlio_close)))
        @property
        def closed(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlio_closed)))
        def fileno(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlio_fileno)))
        def flush(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlio_flush)))
        def isatty(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlio_isatty)))
        def readable(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlio_readable)))
        def readlines(self, hint=-1):
            lines = []
            total = 0
            while hint < 0 or total < hint:
                line = self.readline()
                if line:
                    lines.append(line)
                    total += len(line)
                else:
                    break
            return lines
        def seek(self, offset, whence=0):
            return self._jl_callmethod($(pyjl_methodnum(pyjlio_seek)), offset, whence)
        def seekable(self):
            return True
        def tell(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlio_tell)))
        def truncate(self, size=None):
            return self._jl_callmethod($(pyjl_methodnum(pyjlio_truncate)), size)
        def writable(self):
            return self._jl_callmethod($(pyjl_methodnum(pyjlio_writable)))
        def writelines(self, lines):
            for line in lines:
                self.write(line)
        def __enter__(self):
            return self
        def __exit__(self, t, v, b):
            self.close()
        def __iter__(self):
            return self
        def __next__(self):
            line = self.readline()
            if line:
                return line
            else:
                raise StopIteration
    class BinaryIOValue(IOValueBase):
        __slots__ = ()
        __module__ = "juliacall"
        def detach(self):
            raise ValueError("Cannot detach '{}'.".format(type(self)))
        def read(self, size=-1):
            return self._jl_callmethod($(pyjl_methodnum(pyjlbinaryio_read)), size)
        def read1(self, size=-1):
            return self.read(size)
        def readline(self, size=-1):
            return self._jl_callmethod($(pyjl_methodnum(pyjlbinaryio_readline)), size)
        def readinto(self, b):
            return self._jl_callmethod($(pyjl_methodnum(pyjlbinaryio_readinto)), b)
        def readinto1(self, b):
            return self.readinto(b)
        def write(self, b):
            return self._jl_callmethod($(pyjl_methodnum(pyjlbinaryio_write)), b)
    class TextIOValue(IOValueBase):
        __slots__ = ()
        __module__ = "juliacall"
        @property
        def encoding(self):
            return "UTF-8"
        @property
        def errors(self):
            return "strict"
        def detach(self):
            raise ValueError("Cannot detach '{}'.".format(type(self)))
        def read(self, size=-1):
            return self._jl_callmethod($(pyjl_methodnum(pyjltextio_read)), size)
        def readline(self, size=-1):
            return self._jl_callmethod($(pyjl_methodnum(pyjltextio_readline)), size)
        def write(self, s):
            return self._jl_callmethod($(pyjl_methodnum(pyjltextio_write)), s)
    import io
    io.IOBase.register(IOValueBase)
    io.BufferedIOBase.register(BinaryIOValue)
    io.TextIOBase.register(TextIOValue)
    del io
    """, filename, "exec"), jl.__dict__)
    pycopy!(pyjliobasetype, jl.IOValueBase)
    pycopy!(pyjlbinaryiotype, jl.BinaryIOValue)
    pycopy!(pyjltextiotype, jl.TextIOValue)
end

pyiobase(v::IO) = pyjl(pyjliobasetype, v)

"""
    pybinaryio(io::IO)

Wrap `io` as a Python binary IO object.

This is the default behaviour of `Py(io)`.
"""
pybinaryio(v::IO) = pyjl(pyjlbinaryiotype, v)
export pybinaryio

"""
    pytextio(io::IO)

Wrap `io` as a Python text IO object.
"""
pytextio(v::IO) = pyjl(pyjltextiotype, v)
export pytextio

pyjl(v::IO) = pybinaryio(v)


text:
# check_generail_ledger.jl

using AppliGeneralLedger

using DataFrames

using Query

df = DataFrame(AppliGeneralLedger.read_from_file("test_ledger.txt"))

df2 = df |> @filter(_.accountid == 1300) |> DataFrame

balance_1300 = sum(df2.debit - df2.credit)

println("Balance accounts receivable is $balance_1300. Should be 1210.0.")

df2 = df |> @filter(_.accountid == 8000) |> DataFrame

balance_8000 = sum(df2.credit - df2.debit)

println("Sales is $balance_8000. Should be 4000.0")

#stm = `rm invoicing.sqlite ledger.sqlite log_master.txt`
#stm = `rm invoicing.sqlite ledger.txt journal.txt`

#run(stm)
# end


text:
using NNFEM

viscoelasticity_idx = []

# function load_crack_domain(;option::String = "elasticity", xmax::Float64 = 8.0, ymax::Float64 = 4.0, c1=(4.0, 0.0), c2=(4.0,0.5))
#     nodes, elements = meshread("crack_vertical.msh")
function load_crack_domain(;option::String = "elasticity", xmax::Float64 = 8.0, ymax::Float64 = 4.0, c1=(3.5, 0.0), c2=(4.0,0.5))
    nodes, elements = meshread("crack_wider.msh")
    slope = (c2[2]-c1[2])/(c2[1]-c1[1])
    # slope = 0
    scale_factor = 1000
    nodes *= scale_factor

    # visualize_mesh(nodes, elements)
    ids = Int64[]
    for i = 1:size(nodes, 1)
        x, y = nodes[i,:]
        if  abs(y-slope*(x-c1[1]*scale_factor))<0.05*scale_factor && (c1[1]*0.99*scale_factor < x <= c2[1]*1.01*scale_factor)
        # if  abs(c1[2]*scale_factor <= y <= c2[2]*scale_factor) && abs(x-c1[1]*scale_factor)<0.05*scale_factor
            push!(ids, i)
            # @info x, y
            # plot([x], [y], ".g")
        end
    end

    new_nodes = []
    adic = Dict{Int64, Int64}()
    id1 = Set{Int64}([]) # crack lower surface
    id2 = Set{Int64}([]) # crack upper surface
    id3 = Set{Int64}([]) # free traction
    id4 = Set{Int64}([]) # other fixed boundaries

    for i = 1:size(elements, 1)
        e = elements[i,:]
        if !any([j in ids for j in e])
            continue
        end

        @info i
        x, y = mean(nodes[e,:], dims=1)[:]
        if y-slope*(x-c1[1]*scale_factor)<0
            # plot([x], [y], ".r")
        else
            for j = 1:4
                for e[j] in ids
                    push!(id1, e[j])
                end
            end
            continue
        end

        ## add nodes from the uppper side of crack
        for j = 1:4
            if e[j] in ids
                if haskey(adic,e[j])
                    elements[i,j] = adic[e[j]]
                else
                    push!(new_nodes, nodes[e[j],:])
                    elements[i,j] = size(new_nodes,1) + size(nodes,1)
                    adic[e[j]] = size(new_nodes,1) + size(nodes,1)
                end
                push!(id2, adic[e[j]])
            end
        end

    end

    # savefig("marker.png")

    new_nodes = hcat(new_nodes...)
    nodes = [nodes;new_nodes']

    id1 = collect(id1)
    id2 = collect(id2)

    for j = 1:size(elements, 1)
        e = elements[j,:]
        for i in e
            x, y = nodes[i, :]
            if x<0.001*scale_factor || x>xmax*scale_factor-0.001*scale_factor || y>=ymax*scale_factor-0.001*scale_factor
                push!(id4, i)
            end
            if y<0.001*scale_factor && x>=0.001*scale_factor && x<=xmax*scale_factor-0.001*scale_factor
                push!(id3, i)
            end
        end

    end

    id3 = collect(id3)
    id4 = collect(id4)

    elems = []
    # prop = Dict("name"=> "PlaneStrain", "rho"=> 2700, "E"=> 3.4e10, "nu"=> 0.2339)
    if option=="viscoelasticity"
        prop = Dict("name"=> "ViscoelasticityMaxwell", "rho"=> 2000, "E"=> 1e10, "nu"=> 0.35, "eta"=>1e20)
    elseif option=="elasticity"
        prop = Dict("name"=> "PlaneStrain", "rho"=> 2000, "E"=> 1e10, "nu"=> 0.35)
    end
    global viscoelasticity_idx = []
    g_idx = 0
    for j = 1:size(elements,1)
        elnodes = elements[j,:]
        nodes_ = nodes[elnodes, :]
        ngp = 3
        if option=="mixed"
            if mean(nodes_[:,2])<c2[2]*scale_factor
                prop = Dict("name"=> "ViscoelasticityMaxwell", "rho"=> 2000, "E"=> 1e10, "nu"=> 0.35, "eta"=>1e20)
            else
                prop = Dict("name"=> "ViscoelasticityMaxwell", "rho"=> 2000, "E"=> 1e10, "nu"=> 0.35, "eta"=>1e10)
                # prop = Dict("name"=> "PlaneStrainViscoelasticityProny", "rho"=> 2000, "E"=> 1e10, "nu"=> 0.35, "eta"=>1e10, "c"=>0.99, "tau"=>1e-3)
                # prop = Dict("name"=> "PlaneStrainViscoelasticityProny", "rho"=> 2000, "E"=> 1e10, "nu"=> 0.35, "eta"=>1e20, "c"=>0.99, "tau"=>10.0)
                # prop = Dict("name"=> "ViscoelasticityMaxwell", "rho"=> 2000, "E"=> 1e10, "nu"=> 0.35, "eta"=>1e20)
                # prop = Dict("name"=> "PlaneStrain", "rho"=> 2000, "E"=> 1e8, "nu"=> 0.35)

            end
        end
        push!(elems, SmallStrainContinuum(nodes_, elnodes, prop, ngp))

        if option=="mixed"
            if mean(nodes_[:,2])<c2[2]*scale_factor
            else
                for k = 1:length(elems[end].weights)
                    push!(viscoelasticity_idx, g_idx + k)
                end
            end
        end
        g_idx += length(elems[end].weights)
    end

    # fixed on the bottom, left, and right
    EBC = zeros(Int64, size(nodes, 1), 2)
    FBC = zeros(Int64, size(nodes, 1), 2)
    g = zeros(size(nodes, 1), 2)
    f = zeros(size(nodes, 1), 2)

    EBC[id1,:] .= -1
    EBC[id2,:] .= -1
    EBC[id4,:] .= -1
    g[id1,:] .= [5 5] # m
    g[id2,:] .= [-5 -5] # m
    # g[id1,:] .= [0. 0.5] # m
    # g[id2,:] .= [0. -0.5] # m

    ndims = 2
    domain = Domain(nodes, elems, ndims, EBC, g, FBC, f)
end

function make_patch()
end

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
filename = replace(basename(@__FILE__), ".jl" => "")
filepath = joinpath(@__DIR__, "data/", filename, "input")
input = readlines(filepath)

# part 1

measurements = parse.(Int, input)
single_differences = diff(measurements)
part_1_answer = sum(single_differences .> 0)
println(part_1_answer)

# part 2

n = length(measurements)
window_sums = measurements[1:n-2] + measurements[2:n-1] + measurements[3:n]
triple_differences = diff(window_sums)
part_2_answer = sum(triple_differences .> 0)
println(part_2_answer)


text:
using DependencyTrees, Test

const DT = DependencyTrees

using DependencyTrees: TreebankReader
using DependencyTrees: projective_only, deptype
using DependencyTrees: deprel, form, id, hashead, head, root, isroot
using DependencyTrees: MultiWordTokenError, EmptyTokenError
using DependencyTrees: stack, buffer, stacklength, bufferlength
using DependencyTrees: leftmostdep, rightmostdep
using DependencyTrees: leftdeps, rightdeps
using DependencyTrees: root, noval, token, tokens, xys

using DependencyTrees: typed, untyped
using DependencyTrees: initconfig, transition_space
using DependencyTrees: leftarc, rightarc, noarc, shift, isfinal
using DependencyTrees: LeftArc, RightArc, NoArc, Reduce, Shift
using DependencyTrees: gold_transitions

function showstr(op)
    buf = IOBuffer()
    show(buf, op)
    return String(take!(buf))
end

function test_treebank(filename)
    Treebank{CoNLLU}(joinpath(@__DIR__, "data", filename))
end

function test_sentence(filename)
    first(test_treebank(filename))
end

@testset "DependencyTrees" begin
    include("test_tokens.jl")
    include("test_graphs.jl")
    include("test_slp3_ch13.jl")
    include("test_nivre08.jl")
    include("test_kubleretal09.jl")
    include("test_treebanks.jl")
    include("test_conllu.jl")
    include("test_dynamic.jl")
    include("test_arc_hybrid.jl")
    include("test_arc_swift.jl")
    include("test_evaluation.jl")
end


text:
using BenchmarkTools
using GeoPhyInv
using Test

medium = Medium(:elastic_homo3D);

names(medium)

medium.ref

medium.bounds

mgrid = [range(0.0, stop = 10.0, step = 0.1), range(0.0, stop = 30.0, step = 0.2)];

mgrid = fill(range(-10, stop = 10.0, step = 0.1), 3)

medium = Medium(mgrid, [:vp, :rho, :vs]);

vpb = [2100.0, 2200.0];
vsb = [1500, 1700];
rhob = [2100.0, 2300.0];
update!(medium, [:vp, :vs, :rho], [vpb, vsb, rhob]);

fill!(medium);

medium[:Zp];

medium[:vp] .= 3000.0;
medium[:vs] .= 2000.0;

println(medium)

update!(medium, [:vp, :rho], randn_perc = 1.0);


text:
module PolynomialEquations

using Polynomials
using LinearAlgebra
using JuMP
using Clp
using Printf

include("transformations.jl")
include("matrices.jl")
include("analysis.jl")
include("linear.jl")
include("spectral.jl")

export scale
export shift
export conjreciprocal
export cconj
export dconj
export ltbtmatrix
export sylvestermatrix
export iscoprime
export isschurstable
export ishurwitzstable
export isconjsymmetric
export axb
export axby0
export axbyc
export axaxbb
export axbycd
export axycminl1
export spectralfactor

end


text:
using BTrDB
using Test

@testset "BTrDB.jl" begin
    # Write your own tests here.

    @testset "list collections" begin
        @test 1 == 1
    end
end


text:
using ConstrainedDynamics
using ConstrainedDynamics: vrotate
using Rotations
using StaticArrays
using LinearAlgebra

Δt = 0.01
length1 = 1.0
width, depth = 1.0, 1.0
box1 = Box(width, depth, length1, length1, color = RGBA(1., 1., 0.))

for i=1:10
    # Links
    origin = Origin{Float64}()
    link1 = Body(box1)
    link1.m = 1.0
    link1.J = diagm(ones(3))

    p1 = rand(3)
    p2 = rand(3)
    qoff = rand(UnitQuaternion)

    # Constraints
    joint1 = EqualityConstraint(Fixed(origin, link1; p1=p1, p2=p2, qoffset = qoff))

    links = [link1]
    constraints = [joint1]
    shapes = [box1]

    mech = Mechanism(origin, links, constraints, g = 0., shapes = shapes)

    setPosition!(mech,joint1,SA_F64[]; iter=false)
    setVelocity!(mech,joint1,SA_F64[])
    # setForce!(mech,joint1,SA_F64[]) # Does not do anything

    storage = simulate!(mech, 10., record = true)

    truex = [(p1 - vrotate(p2,qoff)) for i=0:999]
    trueq = [qoff for i=0:999]

    @test isapprox(norm(minimalCoordinates(mech, joint1) - SA_F64[]), 0.0; atol = 1e-8)
    @test isapprox(sum(norm.(storage.x[1].-truex))/1000, 0.0; atol = 1e-8)
    @test isapprox(sum(norm.(storage.q[1].-trueq))/1000, 0.0; atol = 1e-8)
end


text:
using HttpServer
using WebSockets

const SDEF = "Tell me more ...";
const SBYE = "Goodbye, please let's talk again soon";

d = Dict();

d["hello"] = "Hi, my name is Eliza";
d["no"]    = "Tell me more about the detail.";
d["yes"]   = "Tell me more about the detail.";
d["you"]   = "Let's not talk about me.";
d["think"] = "Why do you think that?";
d["hate"]  = "So you hate something, tell me more.";
d["what"]  = "Why do you ask ?";
d["want"]  = "I'm here help you with what you want";
d["need"]  = "We all need many thing, is this is important?";
d["why"]   = "Remember therapy is good for you.";
d["know"]  = "How do you know that?";
d["can't"] = "Don't be negative -- be positive!";
d["never"] = "Don't be negative -- be positive!";
d["unhappy"] = "Why are you unhappy?";
d["like"]  = "Why do you like that?";
d["help"]  = "I'm here to help you.";

looping = true;

wsh = WebSocketHandler() do req,clt
  while looping
    s0 = lowercase(read(client))
    s1 = SDEF
    for k in keys(d)
      s = lowercase(k)
      if s == "bye"
        looping = false
        s1 = SBYE
      else if contains(s0, s)
        s1 = d[k]
        break
      end
    end
    write(clt, "ELIZA> $s1")
  end
end

max_stars_repo_path: Module 3/Chapter09/eliza.jl

max_issues_repo_path: Module 3/Chapter09/eliza.jl

max_forks_repo_path: Module 3/Chapter09/eliza.jl


text:
const TRAVERSAL_VECTOR_PRESIZE = 1000

#####################
#
# Preorder vector
#
#####################

function _prevec!(v::Vector{T}, p::T, q::T) where T <: AbstractNode
    push!(v, q)
    for link in q.links
        link.to == p && continue
        _prevec!(v, q, link.to)
    end

    return v
end

"""
    preorder_vector(tree)
    preorder_vector(p [, q])

Get a vector with the nodes of `tree` in preorder.

The traversal is initiated either from the default origin node of the tree, or from an arbitrary node `p`. Subtrees can be traversed by giving the nodes `p` and `q` to define the stem.

See also: [`postorder`](@ref), [`postorder_vector`](@ref), [`preorder`](@ref)
"""
function preorder_vector(p::AbstractNode, q::AbstractNode)
    p == q ||
        areneighbours(p, q) ||
        throw(InvalidTopology("the nodes must be neighbours."))
    v = Vector{typeof(p)}()
    sizehint!(v, TRAVERSAL_VECTOR_PRESIZE)
    return _prevec!(v, p, q)
end
preorder_vector(p::AbstractNode) = preorder_vector(p, p)
preorder_vector(t::AbstractTree) = preorder_vector(t.origin)

#####################
#
# Preorder iterator
#
#####################

function _gather_preorder(list::Vector{Tuple{T,T}}, p::T, q::T) where T <: AbstractNode
    push!(list, (p, q))
    for link in q.links
        link.to == p && continue
        _gather_preorder(list, q, link.to)
    end

    return list
end

"""
    preorder(tree)
    preorder(p [, q])

Get a vector of parent-child tuples with the nodes of `tree` in preorder.

The traversal is initiated either from the default origin node of the tree, or from an arbitrary node `p`. Subtrees can be traversed by giving the nodes `p` and `q` to define the stem.

See also: [`postorder`](@ref), [`postorder_vector`](@ref), [`preorder_vector`](@ref)
"""
function preorder(p::T, q::T) where T <: AbstractNode
    p == q ||
        areneighbours(p, q) ||
        throw(InvalidTopology("the nodes must be neighbours."))

    return _gather_preorder(Vector{Tuple{T,T}}(), p, q)
end
preorder(p::AbstractNode) = preorder(p, p)
preorder(t::AbstractTree) = preorder(t.origin, t.origin)

#####################
#
# Postorder vector
#
#####################

function _postvec!(v::Vector{T}, p::T, q::T) where T <: AbstractNode
    for link in Iterators.reverse(q.links)
        link.to == p && continue
        _postvec!(v, q, link.to)
    end
    push!(v, q)

    return v
end

"""
    postorder_vector(tree)
    postorder_vector(p [, q])

Get a vector with the nodes of `tree` in postorder.

The traversal is initiated either from the default origin node of the tree, or from an arbitrary node `p`. Subtrees can be traversed by giving the nodes `p` and `q` to define the stem.

See also: [`postorder`](@ref), [`preorder`](@ref), [`preorder_vector`](@ref)
"""
function postorder_vector(p::AbstractNode, q::AbstractNode)
    p == q ||
        areneighbours(p, q) ||
        throw(InvalidTopology("the nodes must be neighbours."))
    v = Vector{typeof(p)}()
    sizehint!(v, TRAVERSAL_VECTOR_PRESIZE)
    return _postvec!(v, p, q)
end
postorder_vector(p::AbstractNode) = postorder_vector(p, p)
postorder_vector(t::AbstractTree) = postorder_vector(t.origin)

#####################
#
# Postorder iterator
#
#####################

function _gather_postorder(list::Vector{Tuple{T,T}}, p::T, q::T) where T <: AbstractNode
    for link in Iterators.reverse(q.links)
        link.to == p && continue
        _gather_postorder(list, q, link.to)
    end
    push!(list, (p, q))

    return list
end

"""
    postorder(tree)
    postorder(p [, q])

Get a vector of parent-child tuples with the nodes of `tree` in postorder.

The traversal is initiated either from the default origin node of the tree, or from an arbitrary node `p`. Subtrees can be traversed by giving the nodes `p` and `q` to define the stem.

See also: [`postorder_vector`](@ref), [`preorder`](@ref), [`preorder_vector`](@ref)
"""
function postorder(p::T, q::T) where T <: AbstractNode
    p == q ||
        areneighbours(p, q) ||
        throw(InvalidTopology("the nodes must be neighbours."))

    return _gather_postorder(Vector{Tuple{T,T}}(), p, q)
end
postorder(p::AbstractNode) = postorder(p, p)
postorder(t::AbstractTree) = postorder(t.origin, t.origin)


text:
parr = [
    ["Q", "u[1]", 12, 12, 12, 12],
    ["Q", "u[2]", 12, 12, 12, 12],
    ["Q", "η", 12, 12, 11, 12],
    ["Q", "θ", 12, 12, 12, 0],
    ["s_aux", "y", 12, 12, 12, 12],
    ["s_aux", "w", 12, 12, 12, 12],
    ["s_aux", "pkin", 12, 12, 12, 12],
    ["s_aux", "wz0", 12, 12, 10, 12],
    ["s_aux", "uᵈ[1]", 12, 12, 12, 12],
    ["s_aux", "uᵈ[2]", 12, 12, 12, 12],
    ["s_aux", "ΔGᵘ[1]", 12, 12, 12, 12],
    ["s_aux", "ΔGᵘ[2]", 12, 12, 12, 12],
]

#! format: off
imex = [
 [     "Q",   "u[1]", -3.79026877309730850230e-02,  3.77520487392776424307e-02, -1.40391503429245286812e-06,  5.33572816364154614566e-03 ],
 [     "Q",   "u[2]", -7.13128520899746973227e-03,  6.54348134460207026680e-03, -5.44654209456798120995e-06,  9.83670484102493409076e-04 ],
 [     "Q",      "η", -5.75156897888350494147e-03,  5.06823909948787842961e-03, -1.61345029906403004787e-05,  1.58599716799621352596e-03 ],
 [     "Q",      "θ",  1.99999999999999928946e+01,  2.00000000000002948752e+01,  2.00000000000000426326e+01,  6.30374859517183302003e-14 ],
 [ "s_aux",      "y",  0.00000000000000000000e+00,  1.00000000000000011642e+06,  5.00000000000000000000e+05,  2.92779390974978974555e+05 ],
 [ "s_aux",      "w", -1.93125254703910637023e-05,  1.89413928194037212366e-05,  6.96287201311761918867e-08,  1.82962329428446588302e-06 ],
 [ "s_aux",   "pkin", -1.60000000000000408562e+00,  0.00000000000000000000e+00, -8.00000000000002042810e-01,  4.68447025559967145103e-01 ],
 [ "s_aux",    "wz0", -1.78460606500689411704e-06,  1.42096532423882541587e-06, -3.15575167929213483841e-09,  7.00511144160156445126e-07 ],
 [ "s_aux",  "uᵈ[1]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux",  "uᵈ[2]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux", "ΔGᵘ[1]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux", "ΔGᵘ[2]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
]

explicit_cpu = [
 [     "Q",   "u[1]", -3.74270752639261211625e-02,  3.72763215301363109999e-02, -1.40392694287316997219e-06,  5.29521849931491837837e-03 ],
 [     "Q",   "u[2]", -7.13776376792300999707e-03,  6.54949226335132476257e-03, -5.45004642311143043000e-06,  9.84283148803805074678e-04 ],
 [     "Q",      "η", -5.75146380759523553894e-03,  5.06819905742867966164e-03, -1.61399463692112299070e-05,  1.58594255118538803723e-03 ],
 [     "Q",      "θ",  1.99999999999996518341e+01,  2.00000000000013891110e+01,  2.00000000000004192202e+01,  2.61373866244148208458e-13 ],
 [ "s_aux",      "y",  0.00000000000000000000e+00,  1.00000000000000011642e+06,  5.00000000000000000000e+05,  2.92779390974978974555e+05 ],
 [ "s_aux",      "w", -1.91903846873268650749e-05,  1.88201368003043060118e-05,  6.88331165208082022114e-08,  1.81657738735220659856e-06 ],
 [ "s_aux",   "pkin", -1.60000000000003339551e+00,  0.00000000000000000000e+00, -8.00000000000017919000e-01,  4.68447025559975749331e-01 ],
 [ "s_aux",    "wz0", -1.79359066322475932405e-06,  1.42978124248321576704e-06, -3.01455602004428728418e-09,  6.97246485535118604352e-07 ],
 [ "s_aux",  "uᵈ[1]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux",  "uᵈ[2]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux", "ΔGᵘ[1]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux", "ΔGᵘ[2]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
]

explicit_gpu = deepcopy(explicit_cpu)
explicit_gpu[3][5] = -1.61399463694379806270e-05 # CPU value: -1.61399463692112299070e-05 -> Δ =  2.267507199799762e-16
explicit_gpu[6][5] =  6.88331165196929377526e-08 # CPU value:  6.88331165208082022114e-08 -> Δ = -1.1152644588480957e-18

long = [
 [     "Q",   "u[1]", -1.26047572708648997208e-01,  8.73949752192989676169e-02,  4.77293205335180186562e-05,  7.01620954203814161526e-03 ],
 [     "Q",   "u[2]", -7.74039642996040555545e-02,  1.17079049312627345159e-01, -1.81789027475552530900e-04,  1.05959299008928087282e-02 ],
 [     "Q",      "η", -7.88344390091704622092e-02,  3.98295576513588017731e-02, -9.41579560194082235179e-05,  2.80408708312616314351e-02 ],
 [     "Q",      "θ",  1.99999999999542694695e+01,  2.00000000000501145792e+01,  2.00000000000023661073e+01,  1.20202839114602287074e-12 ],
 [ "s_aux",      "y",  0.00000000000000000000e+00,  4.00000000000000046566e+06,  2.00000000000000000000e+06,  1.15573129229947482236e+06 ],
 [ "s_aux",      "w", -2.39010107630067514077e-05,  6.21800814928817288142e-05,  2.57411032722755278831e-07,  4.26276984128505141513e-06 ],
 [ "s_aux",   "pkin", -1.60000000000220565788e+00,  0.00000000000000000000e+00, -8.00000000000101851860e-01,  4.61946285916540189120e-01 ],
 [ "s_aux",    "wz0", -2.51713525692209742640e-06,  9.22709404863107399641e-07, -8.84750036743458214697e-10,  5.52781925119928323386e-07 ],
 [ "s_aux",  "uᵈ[1]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux",  "uᵈ[2]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux", "ΔGᵘ[1]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
 [ "s_aux", "ΔGᵘ[2]",  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00,  0.00000000000000000000e+00 ],
]

#! format: on

refVals = (
    explicit_cpu = (explicit_cpu, parr),
    explicit_gpu = (explicit_gpu, parr),
    imex = (imex, parr),
    long = (long, parr),
)


text:
using DataFrames, Gadfly, Dates

set_default_plot_size(6inch, 3inch)

# Geom.smooth already has tests for these types:
# Floats: smooth_lm.jl
# Ints: subplot_layers.jl

# This adds tests for Date and DateTime objects

t1 = Date("2001-01-15"):Dates.Month(1):Date("2016-12-31")
t2 = DateTime("2001-01-15"):Dates.Month(1):DateTime("2016-12-31")
t = Float64.(Dates.value.(t1))
n = length(t)
D = DataFrame(t1=t1, t2=t2,
              cycle = 5*sin.(t*2π/365.25)+randn(n),
              trend = 0.1*[1.0:n;].+2*randn(n))
Dl = stack(D, setdiff(names(D), [:t1,:t2]))

plot(D,
    x=:t1, y=:trend, Geom.point,
    Geom.smooth(method=:lm),
    Theme(point_size=1.8pt, key_position=:none)
    )

plot(Dl,
    x=:t2, y=:value, color=:variable, Geom.point,
    Geom.smooth(smoothing=0.05),
    Theme(point_size=1.8pt, key_position=:none)
    )


text:
# test module for bayesian/variational_inference

module TestInference
  using Test

  # target modules
  include(joinpath(split(@__FILE__, "test")[1], "src/bayesian_inference/gauss_gamma.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/variational_inference/categorical_distribution.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/variational_inference/dirichlet_distribution.jl"))
  include(joinpath(split(@__FILE__, "test")[1], "src/variational_inference/variational_inference.jl"))

  function main()
    @testset "BayesianInference" begin
      @testset "GaussGamma" begin
        @test_nowarn GaussGamma.main(test=true)
      end
    end
    @testset "VariationalInference" begin
      @testset "CategoricalDistribution" begin
        @test_nowarn CategoricalDistribution.main()
      end
      @testset "DirichletDistribution" begin
        @test_nowarn DirichletDistribution.main()
      end
      @testset "VariationalInference" begin
        @test_nowarn VariationalInference.main(test=true)
      end
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
"""
```
measurement(m::SmetsWouters{T}, TTT::Matrix{T}, RRR::Matrix{T},
            CCC::Vector{T}) where {T<:AbstractFloat}
```

Assign measurement equation

```
y_t = ZZ*s_t + DD + u_t
```

where

```
Var(ϵ_t) = QQ
Var(u_t) = EE
Cov(ϵ_t, u_t) = 0
```
"""
function measurement(m::SmetsWouters{T},
                     TTT::Matrix{T},
                     RRR::Matrix{T},
                     CCC::Vector{T}; reg::Int = 1) where {T<:AbstractFloat}
    endo      = m.endogenous_states
    endo_addl = m.endogenous_states_augmented
    exo       = m.exogenous_shocks
    obs       = m.observables

    _n_observables = n_observables(m)
    _n_states = n_states_augmented(m)
    _n_shocks_exogenous = n_shocks_exogenous(m)

    ZZ = zeros(_n_observables, _n_states)
    DD = zeros(_n_observables)
    EE = zeros(_n_observables, _n_observables)
    QQ = zeros(_n_shocks_exogenous, _n_shocks_exogenous)

    for para in m.parameters
        if !isempty(para.regimes)
            if (haskey(get_settings(m), :model2para_regime) ? haskey(get_setting(m, :model2para_regime), para.key) : false)
                ModelConstructors.toggle_regime!(para, reg, get_setting(m, :model2para_regime)[para.key])
            else
                ModelConstructors.toggle_regime!(para, reg)
            end
        end
    end

    ## Output growth - Quarterly!
    ZZ[obs[:obs_gdp], endo[:y_t]]       = 1.0
    ZZ[obs[:obs_gdp], endo_addl[:y_t1]] = -1.0
    ZZ[obs[:obs_gdp], endo[:z_t]]       = 1.0
    DD[obs[:obs_gdp]]                   = 100*(exp(m[:zstar])-1)

    ## Hours growth
    ZZ[obs[:obs_hours], endo[:L_t]] = 1.0
    DD[obs[:obs_hours]]             = m[:Lmean]

    ## Labor Share/real wage growth
    ZZ[obs[:obs_wages], endo[:w_t]]       = 1.0
    ZZ[obs[:obs_wages], endo_addl[:w_t1]] = -1.0
    ZZ[obs[:obs_wages], endo[:z_t]]       = 1.0
    DD[obs[:obs_wages]]                   = 100*(exp(m[:zstar])-1)

    ## Inflation (GDP Deflator)
    ZZ[obs[:obs_gdpdeflator], endo[:π_t]]  = 1.0
    DD[obs[:obs_gdpdeflator]]              = 100*(m[:π_star]-1)

    ## Nominal interest rate
    ZZ[obs[:obs_nominalrate], endo[:R_t]]       = 1.0
    DD[obs[:obs_nominalrate]]                   = m[:Rstarn]

    ## Consumption Growth
    ZZ[obs[:obs_consumption], endo[:c_t]]       = 1.0
    ZZ[obs[:obs_consumption], endo_addl[:c_t1]] = -1.0
    ZZ[obs[:obs_consumption], endo[:z_t]]       = 1.0
    DD[obs[:obs_consumption]]                   = 100*(exp(m[:zstar])-1)

    ## Investment Growth
    ZZ[obs[:obs_investment], endo[:i_t]]       = 1.0
    ZZ[obs[:obs_investment], endo_addl[:i_t1]] = -1.0
    ZZ[obs[:obs_investment], endo[:z_t]]       = 1.0
    DD[obs[:obs_investment]]                   = 100*(exp(m[:zstar])-1)

    #Measurement error
    EE[obs[:obs_gdp],1] = m[:e_y]^2
    EE[obs[:obs_hours],2] = m[:e_L]^2
    EE[obs[:obs_wages],3] = m[:e_w]^2
    EE[obs[:obs_gdpdeflator],4] = m[:e_π]^2
    EE[obs[:obs_nominalrate],5] = m[:e_R]^2
    EE[obs[:obs_consumption],6] = m[:e_c]^2
    EE[obs[:obs_investment],7] = m[:e_i]^2

    #Variance of innovations
    QQ[exo[:g_sh], exo[:g_sh]]           = m[:σ_g]^2
    QQ[exo[:b_sh], exo[:b_sh]]           = m[:σ_b]^2
    QQ[exo[:μ_sh], exo[:μ_sh]]           = m[:σ_μ]^2
    QQ[exo[:z_sh], exo[:z_sh]]           = m[:σ_z]^2
    QQ[exo[:λ_f_sh], exo[:λ_f_sh]]       = m[:σ_λ_f]^2
    QQ[exo[:λ_w_sh], exo[:λ_w_sh]]       = m[:σ_λ_w]^2
    QQ[exo[:rm_sh], exo[:rm_sh]]         = m[:σ_rm]^2

    # These lines set the standard deviations for the anticipated
    # shocks to be equal to the standard deviation for the
    # unanticipated policy shock
    for i = 1:n_anticipated_shocks(m)
        ZZ[obs[Symbol("obs_nominalrate$i")], :] = ZZ[obs[:obs_nominalrate], :]' * (TTT^i)
        DD[obs[Symbol("obs_nominalrate$i")]]    = m[:Rstarn]
        if subspec(m) == "ss1"
            QQ[exo[Symbol("rm_shl$i")], exo[Symbol("rm_shl$i")]] = m[Symbol("σ_rm")]^2 / n_anticipated_shocks(m)
        else
            QQ[exo[Symbol("rm_shl$i")], exo[Symbol("rm_shl$i")]] = m[Symbol("σ_rm")]^2 / 16
        end
    end

    # Adjustment to DD because measurement equation assumes CCC is the zero vector
    if any(CCC .!= 0)
        DD += ZZ*((UniformScaling(1) - TTT)\CCC)
    end

    for para in m.parameters
        if !isempty(para.regimes)
            ModelConstructors.toggle_regime!(para, 1)
        end
    end

    return Measurement(ZZ, DD, QQ, EE)
end


text:
@enum Orientation begin
    AlongX
    AlongY
    AlongZ
end

collect_grid(grid) = grid

function collect_grid(grid::NTuple{N, T}) where {N, T <: AbstractRange}
    map(idxs->SVector{N}(idxs...), Iterators.product(grid...))
end

function collect_grid(grid::NTuple{N, T}) where {N, T <: AbstractVector}
    map(idxs->SVector{N}(idxs...), zip(grid...))
end

to_cylindrical(::Val{AlongZ}) = CylindricalFromCartesian()

to_cylindrical(f::AbstractPICDataStructure{3}, o::Orientation) = to_cylindrical(scalarness(f), f, o)

function to_cylindrical(::ScalarQuantity, f, orientation)
    # TODO: Add coordinate system to fields
    # assume cartesian for now
    dense_grid = collect_grid(f.grid)
    t = to_cylindrical(Val(orientation))
    transformed_grid = t.(dense_grid)

    parameterless_type(f)(f.data, transformed_grid)
end


text:
@testset "oversample" begin
    x = rand(2, 5)
    ya = ["a", "c", "c", "a", "b"]
    y2 = ["c", "c", "c", "a", "b"]

    o = oversample((x, ya), fraction=1, shuffle=false)
    @test o == oversample((x, ya), ya, shuffle=false)
    ox, oy = getobs(o)
    @test ox isa Matrix
    @test oy isa Vector
    @test size(ox) == (2, 6)
    @test size(oy) == (6,)
    @test ox[:,1:5] == x
    @test ox[:,6] == ox[:,5]
    @test oy[1:5] == ya
    @test oy[6] == ya[5]

    o = oversample((x, ya), y2, shuffle=false)
    ox, oy = getobs(o)
    @test ox isa Matrix
    @test oy isa Vector
    @test size(ox) == (2, 9)
    @test size(oy) == (9,)
    @test ox[:,1:5] == x
    @test ox[:,6] == ox[:,7] == x[:,5]
    @test ox[:,8] == ox[:,9] == x[:,4]
    @test oy[1:5] == ya
    @test oy[6] == oy[7] == ya[5]
    @test oy[8] == oy[9] == ya[4]
end

@testset "undersample" begin
    x = rand(2, 5)
    ya = ["a", "c", "c", "a", "b"]
    y2 = ["c", "c", "c", "a", "b"]

    o = undersample((x, ya), shuffle=false)
    ox, oy = getobs(o)
    @test ox isa Matrix
    @test oy isa Vector
    @test size(ox) == (2, 3)
    @test size(oy) == (3,)
    @test ox[:,3] == x[:,5]

    o = undersample((x, ya), y2, shuffle=false)
    ox, oy = getobs(o)
    @test ox isa Matrix
    @test oy isa Vector
    @test size(ox) == (2, 3)
    @test size(oy) == (3,)
    @test ox[:,2:3] == x[:,4:5]
end


text:
using RoME
using Test
using Manifolds: ProductRepr

##

@testset "Basic PartialPriorYawPose2 test" begin
##

fg = initfg()

# fg.solverParams.graphinit=false
addVariable!(fg, :x, Pose2)
addVariable!(fg, :l, Point2)

addFactor!(fg, [:x], PartialPriorYawPose2(MvNormal([0.,0, 0], [0.01, 0.01, 0.01])))
addFactor!(fg, [:l], PriorPoint2(MvNormal([1.,1], [0.01, 0.01])))

addFactor!(fg, [:x; :l], Pose2Point2(MvNormal([1.,1], [0.1,0.1])))

initAll!(fg)

tree = solveTree!(fg)

M = getManifold(Pose2)
me_ = mean(M, getVal(fg, :x))
@test isapprox(M.manifold[1], me_.parts[1], [0,0], atol=0.2)
@test isapprox(M.manifold[2], me_.parts[2], [1 0; 0 1], atol=0.05)
@test isapprox(mean(getVal(fg, :l)), [1,1], atol = 0.05)

##
end