text:
# RUN: julia -e "import Brutus; Brutus.lit(:emit_optimized)" --startup-file=no %s 2>&1 | FileCheck %s

# Intentionally bad code
function matmul!(C, A, B)
    for i in 1:size(A, 1)
        for k in 1:size(A, 2)
            for j in 1:size(B, 2)
                C[i, j] += A[i, k] * B[k, j]
            end
        end
    end
   return nothing
end

emit(matmul!, Matrix{Float64}, Matrix{Float64}, Matrix{Float64})

# CHECK: Core.MethodMatch(Tuple{typeof(Main.Main.matmul!), Matrix{Float64}, Matrix{Float64}, Matrix{Float64}}, svec(), matmul!(C, A, B) in Main.Main at /home/mccoy/Dev/brutus/test/Codegen/canonicalize/matmul.jl:4, true)after translating to MLIR in JLIR dialect:module  {
# CHECK-NEXT:   func nested @"Tuple{typeof(Main.matmul!), Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}}"(%arg0: !jlir<"typeof(Main.matmul!)">, %arg1: !jlir<"Array{Float64, 2}">, %arg2: !jlir<"Array{Float64, 2}">, %arg3: !jlir<"Array{Float64, 2}">) -> !jlir.Nothing attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir.Core.arraysize} : () -> !jlir<"typeof(Core.arraysize)">
# CHECK-NEXT:     %1 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %2 = "jlir.call"(%0, %arg2, %1) : (!jlir<"typeof(Core.arraysize)">, !jlir<"Array{Float64, 2}">, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %3 = "jlir.constant"() {value = #jlir<"#<intrinsic #29 sle_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %4 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %5 = "jlir.call"(%3, %4, %2) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     %6 = "jlir.constant"() {value = #jlir.ifelse} : () -> !jlir<"typeof(ifelse)">
# CHECK-NEXT:     %7 = "jlir.constant"() {value = #jlir<"0">} : () -> !jlir.Int64
# CHECK-NEXT:     %8 = "jlir.call"(%6, %5, %2, %7) : (!jlir<"typeof(ifelse)">, !jlir.Bool, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %9 = "jlir.constant"() {value = #jlir<"#<intrinsic #27 slt_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %10 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %11 = "jlir.call"(%9, %8, %10) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%11)[^bb3, ^bb2] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb2:  // pred: ^bb1
# CHECK-NEXT:     %12 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %13 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %14 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%12, %13, %14)[^bb4] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb3:  // pred: ^bb1
# CHECK-NEXT:     %15 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     %16 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %17 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%15, %16, %17)[^bb4] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb4(%18: !jlir.Bool, %19: !jlir.Int64, %20: !jlir.Int64):  // 2 preds: ^bb2, ^bb3
# CHECK-NEXT:     %21 = "jlir.constant"() {value = #jlir<"#<intrinsic #44 not_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %22 = "jlir.call"(%21, %18) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%22, %19, %20)[^bb28, ^bb5] {operand_segment_sizes = dense<[1, 0, 2]> : vector<3xi32>} : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb5(%23: !jlir.Int64, %24: !jlir.Int64):  // 2 preds: ^bb4, ^bb27
# CHECK-NEXT:     %25 = "jlir.constant"() {value = #jlir.Core.arraysize} : () -> !jlir<"typeof(Core.arraysize)">
# CHECK-NEXT:     %26 = "jlir.constant"() {value = #jlir<"2">} : () -> !jlir.Int64
# CHECK-NEXT:     %27 = "jlir.call"(%25, %arg2, %26) : (!jlir<"typeof(Core.arraysize)">, !jlir<"Array{Float64, 2}">, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %28 = "jlir.constant"() {value = #jlir<"#<intrinsic #29 sle_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %29 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %30 = "jlir.call"(%28, %29, %27) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     %31 = "jlir.constant"() {value = #jlir.ifelse} : () -> !jlir<"typeof(ifelse)">
# CHECK-NEXT:     %32 = "jlir.constant"() {value = #jlir<"0">} : () -> !jlir.Int64
# CHECK-NEXT:     %33 = "jlir.call"(%31, %30, %27, %32) : (!jlir<"typeof(ifelse)">, !jlir.Bool, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %34 = "jlir.constant"() {value = #jlir<"#<intrinsic #27 slt_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %35 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %36 = "jlir.call"(%34, %33, %35) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%36)[^bb7, ^bb6] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb6:  // pred: ^bb5
# CHECK-NEXT:     %37 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %38 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %39 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%37, %38, %39)[^bb8] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb7:  // pred: ^bb5
# CHECK-NEXT:     %40 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     %41 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %42 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%40, %41, %42)[^bb8] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb8(%43: !jlir.Bool, %44: !jlir.Int64, %45: !jlir.Int64):  // 2 preds: ^bb6, ^bb7
# CHECK-NEXT:     %46 = "jlir.constant"() {value = #jlir<"#<intrinsic #44 not_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %47 = "jlir.call"(%46, %43) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%47, %44, %45)[^bb23, ^bb9] {operand_segment_sizes = dense<[1, 0, 2]> : vector<3xi32>} : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb9(%48: !jlir.Int64, %49: !jlir.Int64):  // 2 preds: ^bb8, ^bb22
# CHECK-NEXT:     %50 = "jlir.constant"() {value = #jlir.Core.arraysize} : () -> !jlir<"typeof(Core.arraysize)">
# CHECK-NEXT:     %51 = "jlir.constant"() {value = #jlir<"2">} : () -> !jlir.Int64
# CHECK-NEXT:     %52 = "jlir.call"(%50, %arg3, %51) : (!jlir<"typeof(Core.arraysize)">, !jlir<"Array{Float64, 2}">, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %53 = "jlir.constant"() {value = #jlir<"#<intrinsic #29 sle_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %54 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %55 = "jlir.call"(%53, %54, %52) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     %56 = "jlir.constant"() {value = #jlir.ifelse} : () -> !jlir<"typeof(ifelse)">
# CHECK-NEXT:     %57 = "jlir.constant"() {value = #jlir<"0">} : () -> !jlir.Int64
# CHECK-NEXT:     %58 = "jlir.call"(%56, %55, %52, %57) : (!jlir<"typeof(ifelse)">, !jlir.Bool, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %59 = "jlir.constant"() {value = #jlir<"#<intrinsic #27 slt_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %60 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %61 = "jlir.call"(%59, %58, %60) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%61)[^bb11, ^bb10] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb10:  // pred: ^bb9
# CHECK-NEXT:     %62 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %63 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %64 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%62, %63, %64)[^bb12] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb11:  // pred: ^bb9
# CHECK-NEXT:     %65 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     %66 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %67 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%65, %66, %67)[^bb12] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb12(%68: !jlir.Bool, %69: !jlir.Int64, %70: !jlir.Int64):  // 2 preds: ^bb10, ^bb11
# CHECK-NEXT:     %71 = "jlir.constant"() {value = #jlir<"#<intrinsic #44 not_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %72 = "jlir.call"(%71, %68) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%72, %69, %70)[^bb18, ^bb13] {operand_segment_sizes = dense<[1, 0, 2]> : vector<3xi32>} : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb13(%73: !jlir.Int64, %74: !jlir.Int64):  // 2 preds: ^bb12, ^bb17
# CHECK-NEXT:     %75 = "jlir.constant"() {value = #jlir.Core.arrayref} : () -> !jlir<"typeof(Core.arrayref)">
# CHECK-NEXT:     %76 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %77 = "jlir.call"(%75, %76, %arg1, %23, %73) : (!jlir<"typeof(Core.arrayref)">, !jlir.Bool, !jlir<"Array{Float64, 2}">, !jlir.Int64, !jlir.Int64) -> !jlir.Float64
# CHECK-NEXT:     %78 = "jlir.constant"() {value = #jlir.Core.arrayref} : () -> !jlir<"typeof(Core.arrayref)">
# CHECK-NEXT:     %79 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %80 = "jlir.call"(%78, %79, %arg2, %23, %48) : (!jlir<"typeof(Core.arrayref)">, !jlir.Bool, !jlir<"Array{Float64, 2}">, !jlir.Int64, !jlir.Int64) -> !jlir.Float64
# CHECK-NEXT:     %81 = "jlir.constant"() {value = #jlir.Core.arrayref} : () -> !jlir<"typeof(Core.arrayref)">
# CHECK-NEXT:     %82 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %83 = "jlir.call"(%81, %82, %arg3, %48, %73) : (!jlir<"typeof(Core.arrayref)">, !jlir.Bool, !jlir<"Array{Float64, 2}">, !jlir.Int64, !jlir.Int64) -> !jlir.Float64
# CHECK-NEXT:     %84 = "jlir.constant"() {value = #jlir<"#<intrinsic #14 mul_float>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %85 = "jlir.call"(%84, %80, %83) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Float64, !jlir.Float64) -> !jlir.Float64
# CHECK-NEXT:     %86 = "jlir.constant"() {value = #jlir<"#<intrinsic #12 add_float>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %87 = "jlir.call"(%86, %77, %85) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Float64, !jlir.Float64) -> !jlir.Float64
# CHECK-NEXT:     %88 = "jlir.constant"() {value = #jlir.Core.arrayset} : () -> !jlir<"typeof(Core.arrayset)">
# CHECK-NEXT:     %89 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %90 = "jlir.call"(%88, %89, %arg1, %87, %23, %73) : (!jlir<"typeof(Core.arrayset)">, !jlir.Bool, !jlir<"Array{Float64, 2}">, !jlir.Float64, !jlir.Int64, !jlir.Int64) -> !jlir<"Array{Float64, 2}">
# CHECK-NEXT:     %91 = "jlir.constant"() {value = #jlir<"===">} : () -> !jlir<"typeof(===)">
# CHECK-NEXT:     %92 = "jlir.call"(%91, %74, %58) : (!jlir<"typeof(===)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%92)[^bb15, ^bb14] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb14:  // pred: ^bb13
# CHECK-NEXT:     %93 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %94 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %95 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%93, %94, %95)[^bb16] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb15:  // pred: ^bb13
# CHECK-NEXT:     %96 = "jlir.constant"() {value = #jlir<"#<intrinsic #2 add_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %97 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %98 = "jlir.call"(%96, %74, %97) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %99 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%98, %98, %99)[^bb16] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb16(%100: !jlir.Int64, %101: !jlir.Int64, %102: !jlir.Bool):  // 2 preds: ^bb14, ^bb15
# CHECK-NEXT:     %103 = "jlir.constant"() {value = #jlir<"#<intrinsic #44 not_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %104 = "jlir.call"(%103, %102) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%104)[^bb18, ^bb17] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb17:  // pred: ^bb16
# CHECK-NEXT:     "jlir.goto"(%100, %101)[^bb13] : (!jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb18:  // 2 preds: ^bb12, ^bb16
# CHECK-NEXT:     %105 = "jlir.constant"() {value = #jlir<"===">} : () -> !jlir<"typeof(===)">
# CHECK-NEXT:     %106 = "jlir.call"(%105, %49, %33) : (!jlir<"typeof(===)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%106)[^bb20, ^bb19] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb19:  // pred: ^bb18
# CHECK-NEXT:     %107 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %108 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %109 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%107, %108, %109)[^bb21] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb20:  // pred: ^bb18
# CHECK-NEXT:     %110 = "jlir.constant"() {value = #jlir<"#<intrinsic #2 add_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %111 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %112 = "jlir.call"(%110, %49, %111) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %113 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%112, %112, %113)[^bb21] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb21(%114: !jlir.Int64, %115: !jlir.Int64, %116: !jlir.Bool):  // 2 preds: ^bb19, ^bb20
# CHECK-NEXT:     %117 = "jlir.constant"() {value = #jlir<"#<intrinsic #44 not_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %118 = "jlir.call"(%117, %116) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%118)[^bb23, ^bb22] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb22:  // pred: ^bb21
# CHECK-NEXT:     "jlir.goto"(%114, %115)[^bb9] : (!jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb23:  // 2 preds: ^bb8, ^bb21
# CHECK-NEXT:     %119 = "jlir.constant"() {value = #jlir<"===">} : () -> !jlir<"typeof(===)">
# CHECK-NEXT:     %120 = "jlir.call"(%119, %24, %8) : (!jlir<"typeof(===)">, !jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%120)[^bb25, ^bb24] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb24:  // pred: ^bb23
# CHECK-NEXT:     %121 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %122 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %123 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%121, %122, %123)[^bb26] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb25:  // pred: ^bb23
# CHECK-NEXT:     %124 = "jlir.constant"() {value = #jlir<"#<intrinsic #2 add_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %125 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %126 = "jlir.call"(%124, %24, %125) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %127 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%126, %126, %127)[^bb26] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb26(%128: !jlir.Int64, %129: !jlir.Int64, %130: !jlir.Bool):  // 2 preds: ^bb24, ^bb25
# CHECK-NEXT:     %131 = "jlir.constant"() {value = #jlir<"#<intrinsic #44 not_int>">} : () -> !jlir<"typeof(Core.IntrinsicFunction)">
# CHECK-NEXT:     %132 = "jlir.call"(%131, %130) : (!jlir<"typeof(Core.IntrinsicFunction)">, !jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%132)[^bb28, ^bb27] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb27:  // pred: ^bb26
# CHECK-NEXT:     "jlir.goto"(%128, %129)[^bb5] : (!jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb28:  // 2 preds: ^bb4, ^bb26
# CHECK-NEXT:     %133 = "jlir.constant"() {value = #jlir.nothing} : () -> !jlir.Nothing
# CHECK-NEXT:     "jlir.return"(%133) : (!jlir.Nothing) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

# CHECK:   func nested @"Tuple{typeof(Main.matmul!), Array{Float64, 2}, Array{Float64, 2}, Array{Float64, 2}}"(%arg0: !jlir<"typeof(Main.matmul!)">, %arg1: !jlir<"Array{Float64, 2}">, %arg2: !jlir<"Array{Float64, 2}">, %arg3: !jlir<"Array{Float64, 2}">) -> !jlir.Nothing attributes {llvm.emit_c_interface} {
# CHECK-NEXT:     "jlir.goto"()[^bb1] : () -> ()
# CHECK-NEXT:   ^bb1:  // pred: ^bb0
# CHECK-NEXT:     %0 = "jlir.constant"() {value = #jlir<"1">} : () -> !jlir.Int64
# CHECK-NEXT:     %1 = "jlir.arraysize"(%arg2, %0) : (!jlir<"Array{Float64, 2}">, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %2 = "jlir.sle_int"(%0, %1) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     %3 = "jlir.constant"() {value = #jlir<"0">} : () -> !jlir.Int64
# CHECK-NEXT:     %4 = "jlir.ifelse"(%2, %1, %3) : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %5 = "jlir.slt_int"(%4, %0) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%5)[^bb3, ^bb2] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb2:  // pred: ^bb1
# CHECK-NEXT:     %6 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %7 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %8 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%6, %7, %8)[^bb4] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb3:  // pred: ^bb1
# CHECK-NEXT:     %9 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%9, %0, %0)[^bb4] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb4(%10: !jlir.Bool, %11: !jlir.Int64, %12: !jlir.Int64):  // 2 preds: ^bb2, ^bb3
# CHECK-NEXT:     %13 = "jlir.not_int"(%10) : (!jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%13, %11, %12)[^bb28, ^bb5] {operand_segment_sizes = dense<[1, 0, 2]> : vector<3xi32>} : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb5(%14: !jlir.Int64, %15: !jlir.Int64):  // 2 preds: ^bb4, ^bb27
# CHECK-NEXT:     %16 = "jlir.constant"() {value = #jlir<"2">} : () -> !jlir.Int64
# CHECK-NEXT:     %17 = "jlir.arraysize"(%arg2, %16) : (!jlir<"Array{Float64, 2}">, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %18 = "jlir.sle_int"(%0, %17) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     %19 = "jlir.ifelse"(%18, %17, %3) : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %20 = "jlir.slt_int"(%19, %0) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%20)[^bb7, ^bb6] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb6:  // pred: ^bb5
# CHECK-NEXT:     %21 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %22 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %23 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%21, %22, %23)[^bb8] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb7:  // pred: ^bb5
# CHECK-NEXT:     %24 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%24, %0, %0)[^bb8] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb8(%25: !jlir.Bool, %26: !jlir.Int64, %27: !jlir.Int64):  // 2 preds: ^bb6, ^bb7
# CHECK-NEXT:     %28 = "jlir.not_int"(%25) : (!jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%28, %26, %27)[^bb23, ^bb9] {operand_segment_sizes = dense<[1, 0, 2]> : vector<3xi32>} : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb9(%29: !jlir.Int64, %30: !jlir.Int64):  // 2 preds: ^bb8, ^bb22
# CHECK-NEXT:     %31 = "jlir.arraysize"(%arg3, %16) : (!jlir<"Array{Float64, 2}">, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %32 = "jlir.sle_int"(%0, %31) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     %33 = "jlir.ifelse"(%32, %31, %3) : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %34 = "jlir.slt_int"(%33, %0) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%34)[^bb11, ^bb10] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb10:  // pred: ^bb9
# CHECK-NEXT:     %35 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %36 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %37 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%35, %36, %37)[^bb12] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb11:  // pred: ^bb9
# CHECK-NEXT:     %38 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%38, %0, %0)[^bb12] : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb12(%39: !jlir.Bool, %40: !jlir.Int64, %41: !jlir.Int64):  // 2 preds: ^bb10, ^bb11
# CHECK-NEXT:     %42 = "jlir.not_int"(%39) : (!jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%42, %40, %41)[^bb18, ^bb13] {operand_segment_sizes = dense<[1, 0, 2]> : vector<3xi32>} : (!jlir.Bool, !jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb13(%43: !jlir.Int64, %44: !jlir.Int64):  // 2 preds: ^bb12, ^bb17
# CHECK-NEXT:     %45 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     %46 = "jlir.arrayref"(%45, %arg1, %14, %43) : (!jlir.Bool, !jlir<"Array{Float64, 2}">, !jlir.Int64, !jlir.Int64) -> !jlir.Float64
# CHECK-NEXT:     %47 = "jlir.arrayref"(%45, %arg2, %14, %29) : (!jlir.Bool, !jlir<"Array{Float64, 2}">, !jlir.Int64, !jlir.Int64) -> !jlir.Float64
# CHECK-NEXT:     %48 = "jlir.arrayref"(%45, %arg3, %29, %43) : (!jlir.Bool, !jlir<"Array{Float64, 2}">, !jlir.Int64, !jlir.Int64) -> !jlir.Float64
# CHECK-NEXT:     %49 = "jlir.mul_float"(%47, %48) : (!jlir.Float64, !jlir.Float64) -> !jlir.Float64
# CHECK-NEXT:     %50 = "jlir.add_float"(%46, %49) : (!jlir.Float64, !jlir.Float64) -> !jlir.Float64
# CHECK-NEXT:     %51 = "jlir.arrayset"(%45, %arg1, %50, %14, %43) : (!jlir.Bool, !jlir<"Array{Float64, 2}">, !jlir.Float64, !jlir.Int64, !jlir.Int64) -> !jlir<"Array{Float64, 2}">
# CHECK-NEXT:     %52 = "jlir.==="(%44, %33) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%52)[^bb15, ^bb14] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb14:  // pred: ^bb13
# CHECK-NEXT:     %53 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %54 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     "jlir.goto"(%53, %54, %45)[^bb16] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb15:  // pred: ^bb13
# CHECK-NEXT:     %55 = "jlir.add_int"(%44, %0) : (!jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %56 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%55, %55, %56)[^bb16] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb16(%57: !jlir.Int64, %58: !jlir.Int64, %59: !jlir.Bool):  // 2 preds: ^bb14, ^bb15
# CHECK-NEXT:     %60 = "jlir.not_int"(%59) : (!jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%60)[^bb18, ^bb17] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb17:  // pred: ^bb16
# CHECK-NEXT:     "jlir.goto"(%57, %58)[^bb13] : (!jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb18:  // 2 preds: ^bb12, ^bb16
# CHECK-NEXT:     %61 = "jlir.==="(%30, %19) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%61)[^bb20, ^bb19] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb19:  // pred: ^bb18
# CHECK-NEXT:     %62 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %63 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %64 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%62, %63, %64)[^bb21] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb20:  // pred: ^bb18
# CHECK-NEXT:     %65 = "jlir.add_int"(%30, %0) : (!jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %66 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%65, %65, %66)[^bb21] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb21(%67: !jlir.Int64, %68: !jlir.Int64, %69: !jlir.Bool):  // 2 preds: ^bb19, ^bb20
# CHECK-NEXT:     %70 = "jlir.not_int"(%69) : (!jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%70)[^bb23, ^bb22] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb22:  // pred: ^bb21
# CHECK-NEXT:     "jlir.goto"(%67, %68)[^bb9] : (!jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb23:  // 2 preds: ^bb8, ^bb21
# CHECK-NEXT:     %71 = "jlir.==="(%15, %4) : (!jlir.Int64, !jlir.Int64) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%71)[^bb25, ^bb24] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb24:  // pred: ^bb23
# CHECK-NEXT:     %72 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %73 = "jlir.undef"() : () -> !jlir.Int64
# CHECK-NEXT:     %74 = "jlir.constant"() {value = #jlir.true} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%72, %73, %74)[^bb26] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb25:  // pred: ^bb23
# CHECK-NEXT:     %75 = "jlir.add_int"(%15, %0) : (!jlir.Int64, !jlir.Int64) -> !jlir.Int64
# CHECK-NEXT:     %76 = "jlir.constant"() {value = #jlir.false} : () -> !jlir.Bool
# CHECK-NEXT:     "jlir.goto"(%75, %75, %76)[^bb26] : (!jlir.Int64, !jlir.Int64, !jlir.Bool) -> ()
# CHECK-NEXT:   ^bb26(%77: !jlir.Int64, %78: !jlir.Int64, %79: !jlir.Bool):  // 2 preds: ^bb24, ^bb25
# CHECK-NEXT:     %80 = "jlir.not_int"(%79) : (!jlir.Bool) -> !jlir.Bool
# CHECK-NEXT:     "jlir.gotoifnot"(%80)[^bb28, ^bb27] {operand_segment_sizes = dense<[1, 0, 0]> : vector<3xi32>} : (!jlir.Bool) -> ()
# CHECK-NEXT:   ^bb27:  // pred: ^bb26
# CHECK-NEXT:     "jlir.goto"(%77, %78)[^bb5] : (!jlir.Int64, !jlir.Int64) -> ()
# CHECK-NEXT:   ^bb28:  // 2 preds: ^bb4, ^bb26
# CHECK-NEXT:     %81 = "jlir.constant"() {value = #jlir.nothing} : () -> !jlir.Nothing
# CHECK-NEXT:     "jlir.return"(%81) : (!jlir.Nothing) -> ()
# CHECK-NEXT:   }
# CHECK-NEXT: }

max_stars_repo_name: JuliaLabs/brutus

max_issues_repo_name: JuliaLabs/brutus

max_forks_repo_name: JuliaLabs/brutus


text:
function areequal(theory::Vector{<:AbstractRule}, exprs...; params=SaturationParams())
    g = EGraph(exprs[1])
    areequal(g, theory, exprs...; params=params)
end

function areequal(g::EGraph, t::Vector{<:AbstractRule}, exprs...; params=SaturationParams())
    @log "Checking equality for " exprs
    if length(exprs) == 1; return true end
    # rebuild!(G)

    @log "starting saturation"

    n = length(exprs)
    ids = Vector{EClassId}(undef, n)
    nodes = Vector{AbstractENode}(undef, n)
    for i ∈ 1:n
        ec, node = addexpr!(g, exprs[i])
        ids[i] = ec.id
        nodes[i] = node
    end

    goal = EqualityGoal(collect(exprs), ids)

    # alleq = () -> (all(x -> in_same_set(G.uf, ids[1], x), ids[2:end]))

    params.goal = goal
    # params.stopwhen = alleq

    report = saturate!(g, t, params)

    # display(g.classes); println()
    if !(report.reason isa Saturated) && !reached(g, goal)
        return missing # failed to prove
    end
    return reached(g, goal)
end

import ..gettheory

macro areequal(theory, exprs...)
    t = gettheory(theory, __module__)
    areequal(t, exprs...)
end

macro areequalg(G, theory, exprs...)
    t = gettheory(theory, __module__)
    areequal(getfield(__module__, G), t, exprs...)
end

ext: jl


text:
# ---
# title: 1541. Minimum Insertions to Balance a Parentheses String
# id: problem1541
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Stack
# link: <https://leetcode.com/problems/minimum-insertions-to-balance-a-parentheses-string/description/>
# hidden: true
# ---
#
# Given a parentheses string `s` containing only the characters `'('` and `')'`.
# A parentheses string is **balanced** if:
#
#   * Any left parenthesis `'('` must have a corresponding two consecutive right parenthesis `'))'`.
#   * Left parenthesis `'('` must go before the corresponding two consecutive right parenthesis `'))'`.
#
# In other words, we treat `'('` as openning parenthesis and `'))'` as closing
# parenthesis.
#
# For example, `"())"`, `"())(())))"` and `"(())())))"` are balanced, `")()"`,
# `"()))"` and `"(()))"` are not balanced.
#
# You can insert the characters `'('` and `')'` at any position of the string to
# balance it if needed.
#
# Return _the minimum number of insertions_ needed to make `s` balanced.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "(()))"
#     Output: 1
#     Explanation: The second '(' has two matching '))', but the first '(' has only ')' matching. We need to to add one more ')' at the end of the string to be "(())))" which is balanced.
#
#
# **Example 2:**
#
#
#
#     Input: s = "())"
#     Output: 0
#     Explanation: The string is already balanced.
#
#
# **Example 3:**
#
#
#
#     Input: s = "))())("
#     Output: 3
#     Explanation: Add '(' to match the first '))', Add '))' to match the last '('.
#
#
# **Example 4:**
#
#
#
#     Input: s = "(((((("
#     Output: 12
#     Explanation: Add 12 ')' to balance the string.
#
#
# **Example 5:**
#
#
#
#     Input: s = ")))))))"
#     Output: 5
#     Explanation: Add 4 '(' at the beginning of the string and one ')' at the end. The string becomes "(((())))))))".
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 10^5`
#   * `s` consists of `'('` and `')'` only.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
"""

```
edit!(atoms::Vector{Atom})
```

Opens a temporary PDB file in which the fields of the vector of atoms can be edited.

"""
function edit!(atoms::AbstractVector{Atom})
    tmp_file_name = tempname()
    writePDB(atoms, tmp_file_name)
    InteractiveUtils.edit(tmp_file_name)
    read_again = readPDB(tmp_file_name)
    resize!(atoms, length(read_again))
    @. atoms = read_again
    rm(tmp_file_name)
end

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
# unary.jl: Unary Array->Array operations.
# The following list comes from the NVIDIA math docs with some extras.
# http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#mathematical-functions-appendix
# The entry format is (cudaname, julianame, kernelcode)
# With single name entries cudaname=julianame and kernelcode=name(xi).
# I commented out functions if I don't know the Julia equivalent.
unary_ops = [
("abs2", "abs2", "(xi*xi)"),
("abs", "abs", "(xi<0?-xi:xi)"),
("htanh", "htanh", "(xi<-1?-1 : (xi > 1 ? 1 :xi))"),
"acos",
"acosh",
"asin",
"asinh",
"atan",
"atanh",
"cbrt",
"ceil",
"cos",
"cosh",
"cospi",
# "cyl_bessel_i0",
# "cyl_bessel_i1",
"erf",
"erfc",
"erfcinv",
"erfcx",
"erfinv",
"exp",
"exp10",
"exp2",
"expm1",
"floor",
# "ilogb",
("invx", "invx", "1/xi"),
# "j0",
# "j1",
# "lgamma", # missing digamma for derivative
# "llrint",
# "llround",
"log",
"log10",
"log1p",
"log2",
# "logb",
# "lrint",
# "lround",
# "nearbyint",
("neg", "-", "-xi"),
# "normcdf",
# "normcdfinv",
# "rcbrt",
("relu", "relu", "(xi>0?xi:0)"),
("lzorelu", "lzorelu", "(xi<1e-6?1e-6:xi)"),
# "rint",
"round",
# "rsqrt",
("sigm", "sigm", "(xi>=0?1/(1+exp(-xi)):(exp(xi)/(1+exp(xi))))"),
("sign", "sign", "(xi>0?1:xi<0?-1:0)"),
"sin",
"sinh",
"sinpi",
"sqrt",
"tan",
"tanh",
# "tgamma",
"trunc",
# "y0",
# "y1",
("G", "G", "exp(-(xi*xi)/2) / 2.5066"),
("H", "H", "0.5 * erfc(xi / 1.4142)"),
("GH", "GH", "xi > 30 ? xi + 1/xi * (1 - 2/(xi*xi) * (1 - 5/(xi*xi) * (1 - 7.4/(xi*xi)))) : exp(-(xi*xi)/2) / (1.2533 * erfc(xi / 1.4142))"),
"atanh2Hm1",
]

function unary_op(f, j=f, o...)
    J=Symbol(j)
    if isdefined(Base, J); eval(Expr(:import,:Base,J)); end
    for S in (32,64)
        T = Symbol("Float$S")
        F = "$(f)_$S"
        @eval begin
            function $J(x::KnetArray{$T})
                y = similar(x)
                ccall(($F,$libknet8),Void,(Cint,Ptr{$T},Ptr{$T}),length(y),x,y)
                return y
            end
        end
    end
end

for f in unary_ops
    if !isa(f,Tuple); f=(f,); end
    unary_op(f...)
end

# Define some common operations as primitives for efficiency:
# 1. Avoid creating intermediate arrays
# 2. Avoid taking derivatives of intermediate operations

for (f,g,y,dx) in
    ((:invx, :invxback, :(one(T)/xi), :(-yi*yi*dyi)),
     (:relu, :reluback, :(max(zero(T),xi)), :(ifelse(yi>0,dyi,zero(T)))),
     #(:H, :G, :(H(xi)), :(dyi*G(xi))),
     (:tanx, :tanhback, :(tanh(xi)), :(dyi*(one(T)-yi*yi))),
     (:htanx, :htanhback, :(htanh(xi)), :(dyi*(-1.0<=yi<=1.0))),
     (:sigm, :sigmback,
      # Numerically stable implementation from
      # http://timvieira.github.io/blog/post/2014/02/11/exp-normalize-trick
      :(if xi>=0; z=exp(-xi); one(T)/(one(T)+z); else; z=exp(xi); z/(one(T)+z); end),
      :(dyi*yi*(one(T)-yi))),
     )
    @eval begin
        function $f{T<:AbstractFloat}(x::Array{T})
            y = similar(x)
            @inbounds for i=1:length(y)
                xi = x[i]
                y[i] = $y
            end
            return y
        end
        function $g{T<:AbstractFloat}(dy::Array{T},y::Array{T})
            dx = similar(dy)
            @inbounds for i=1:length(dx)
                yi = y[i]
                dyi = dy[i]
                dx[i] = $dx
            end
            return dx
        end
        $f{T<:Number}(xi::T)=$y
        $g{T<:Number}(dyi::T,yi::T)=$dx
        @primitive $f(x),dy,y $g(dy,y)
    end
end

"`invx(x) = (1./x)`" invx
"`relu(x) = max(0,x)`" relu
"`sigm(x) = (1./(1+exp(-x)))`" sigm

# To avoid conflict with AutoGrad:
@primitive tanh(x::Array),dy,y     tanhback(dy,y)
@primitive tanh(x::KnetArray),dy,y tanhback(dy,y)
@primitive tanhback(dy,y),ddx  ddx.*(1.-y.*y)  ddx.*(-2.*dy.*y)
@primitive H(x),dy dy*G(x)

"""

    logp(x,[dims])

Treat entries in `x` as as unnormalized log probabilities and return
normalized log probabilities.

`dims` is an optional argument, if not specified the normalization is
over the whole `x`, otherwise the normalization is performed over the
given dimensions.  In particular, if `x` is a matrix, `dims=1`
normalizes columns of `x` and `dims=2` normalizes rows of `x`.

"""

# Math for the cross-entropy loss: x is unnormalized input, p is
# target probabilities, q is estimated probabilities. Read left column
# down, right column (loss gradients) back up.

# x			dx = -p + qz/z = -p + exp(logq)
# xmax  = max(x,1)	-sum(db)=0
# logqz = x .- xmax	-p + qz/z
# qz    = exp(logqz)	rep(1/z)
# z     = sum(qz,1)	1/z
# logz  = log(z)	sum(p)=1
# logq  = logqz.-logz	-p
# plogq = p .* logq	-1
# loss  = -sum(plogq)	1

function logp(x,d...)
    if isa(x,Number)
        return zero(x)
    elseif isempty(x)
        return x
    else
        x = x .- maximum(x,d...)
        return (x .- log(sum(exp(x),d...)))
    end
end

function logpback(x,y,dy,d...)
    x = AutoGrad.getval(x)
    if isa(x,Number)
        return zero(x)
    elseif isempty(x)
        return x
    else
        return (dy - exp(y).*sum(dy,d...))
    end
end

# dy should be -p and y=logq so this should give us -p+q
@primitive  logp(x,d...),dy,y  logpback(x,y,dy,d...)


text:
using LoopVectorization
function Hd0(αindex::Array{Int64,1}, βindex::Array{Int64,1}, h::Array{T, 2}, V::Array{T, 4}) where T <: AbstractFloat
    """
    Σ [m|h|m] + 1/2 ΣΣ [mm|nn] - [mn|nm]
    """
    # One electron energy
    E1 = 0.0
    # Two electron energy
    E2 = 0.0
    la = length(αindex)
    lb = length(βindex)
    for n in 1:la
        N = αindex[n]
        E1 += h[N,N]
        for m in n:la
            M = αindex[m]
            E2 += V[M,M,N,N] - V[M,N,N,M]
        end
    end
    for m in 1:la
        M = αindex[m]
        for n in m:la
            N = αindex[n]
            E2 += V[N,N,M,M] - V[N,M,M,N]
        end
        for n in 1:lb
            N = βindex[n]
            E2 += 2V[M,M,N,N]
        end
    end

    for m in 1:lb
        M = βindex[m]
        E1 += h[M,M]
        for n in m:lb
            N = βindex[n]
            E2 += V[M,M,N,N] - V[M,N,N,M]
            E2 += V[N,N,M,M] - V[N,M,M,N]
        end
    end
    return E1 + 0.5*E2
end

function Hd1(αindex::Array{Int64,1}, βindex::Array{Int64,1}, D1::Determinant, D2::Determinant, h::Array{T,2}, V::Array{T, 4}, αexc::Int) where T <: AbstractFloat
    """
    differ m -> p
    [m|h|p] + Σ[mp|nn] - [mn|np]
    """

    # if m and p are α
    if αexc == 1
        m = first_αexclusive(D1, D2)
        p = first_αexclusive(D2, D1)

        # Compute phase by counting the number of occupied orbitals between m and p

        i,f = minmax(m,p)
        r = count_ones(D1.α >> i)
        l = count_ones(D1.α << (65-f))
        t = count_ones(D1.α)
        t -= r
        t -= l
        #t = l
        ph = (-1)^t

        # Compute matrix element
        E = h[m,p]
        @avx for i in 1:length(αindex)
            I = αindex[i]
            E += V[m,p,I,I] - V[m,I,I,p]
        end
        @avx for i in 1:length(βindex)
            I = βindex[i]
            E += V[m,p,I,I]
        end
        return ph*E

    else
        m = first_βexclusive(D1, D2)
        p = first_βexclusive(D2, D1)

        # Compute phase by counting the number of occupied orbitals between m and p

        i,f = minmax(m,p)
        r = count_ones(D1.β >> i)
        l = count_ones(D1.β << (65-f))
        t = count_ones(D1.β)
        t -= r
        t -= l
        #t = l
        ph = (-1)^t

        # Compute matrix element
        E = h[m,p]
        @avx for i in 1:length(αindex)
            I = αindex[i]
            E += V[m,p,I,I]
        end
        @avx for i in 1:length(βindex)
            I = βindex[i]
            E += V[m,p,I,I] - V[m,I,I,p]
        end
        return ph*E
        return E
    end
end
function Hd2(D1::Determinant, D2::Determinant, V::Array{T, 4}, αexc::Int) where T <: AbstractFloat
    """
    mn -> pq
    [mp|nq] - [mq|np]
    """

    # If α excitation is one, it means m and n have different spins
    if αexc == 1
        m = first_αexclusive(D1, D2)
        n = first_βexclusive(D1, D2)
        p = first_αexclusive(D2, D1)
        q = first_βexclusive(D2, D1)

        # For the phase factor, there is no interference between the two cases since they have difference spin
        # Move m <-> p
        #i = min(m,p)
        #f = max(m,p)-1
        i,f = minmax(m,p)
        #println("$i $f")
        #println(bitstring(D1.α))
        r = count_ones(D1.α >> i)
        l = count_ones(D1.α << (65-f))
        t = count_ones(D1.α)
        t -= r
        t -= l
        ph1 = (-1)^t

        # Move n <-> q
        i,f = minmax(n,q)
        #i = min(n,q)
        #f = max(n,q)-1
        r = count_ones(D1.β >> i)
        l = count_ones(D1.β << (65-f))
        t = count_ones(D1.β)
        t -= r
        t -= l
        ph2 = (-1)^t

        return ph1*ph2*V[m,p,n,q]

    # If α excitation is two, it means m,n,p and q are all α.
    elseif αexc == 2
        m = first_αexclusive(D1, D2)
        n = second_αexclusive(D1, D2)
        p = first_αexclusive(D2, D1)
        q = second_αexclusive(D2, D1)

        #Transform D1 -> D2. First take m into p
        #i = min(m,p)
        #f = max(m,p)-1
        i,f = minmax(m,p)
        #println("$i $f")
        #println(bitstring(D1.α))
        r = count_ones(D1.α >> i)
        l = count_ones(D1.α << (65-f)) # 65 because f used to be f = f - 1
        t = count_ones(D1.α)
        t -= r
        t -= l
        #t = l
        #iseven(t) ? ph1 = 1 : ph1 = -1
        ph1 = (-1)^t

        # Update bits
        newα = (D1.α ⊻ (1 << ((m-1)&63))) | (1 << ((p-1)&63))

        # Take n into q using the updated bit

        #i = (min(n,q))
        #f = ((max(n,q) - 1))
        i,f = minmax(n,q)

        r = count_ones(newα >> (i&63))
        l = count_ones(newα << ((65-f)&63))
        t = count_ones(newα)
        t -= r
        t -= l
        ph2 = (-1)^t
        #iseven(t) ? nothing : ph1 *= -1

        return ph1*ph2*(V[m,p,n,q] - V[m,q,n,p])

    # If α excitation is zero, it means m,n,p and q are all β.
    elseif αexc == 0
        m = first_βexclusive(D1, D2)
        n = second_βexclusive(D1, D2)
        p = first_βexclusive(D2, D1)
        q = second_βexclusive(D2, D1)

        #Transform D1 -> D2. First take m into p
        i,f = minmax(m,p)
        #println("$i $f")
        #println(bitstring(D1.α))
        r = count_ones(D1.β >> i)
        l = count_ones(D1.β << (65-f))
        t = count_ones(D1.β)
        t -= r
        t -= l
        #t = l
        ph1 = (-1)^t

        # Update bits
        newβ = (D1.β ⊻ (1 << ((m-1)&63))) | (1 << ((p-1)&63))

        # Take n into q using the updated bit

        i,f = minmax(n,q)

        r = count_ones(newβ >> (i&63))
        l = count_ones(newβ << ((65-f)&63))
        t = count_ones(newβ)
        t -= r
        t -= l
        ph2 = (-1)^t

        return ph1*ph2*(V[m,p,n,q] - V[m,q,n,p])
    end
end


text:
using Reproject: wcs_to_celestial_frame
@testset "wcs to celestial frame" begin
    wcs1 = WCSTransform(2;
                       ctype = ["RA---AIR", "DEC--AIR"],
                       )
    wcs2 = WCSTransform(2;
                       ctype = ["RA---AIR", "DEC--AIR"],
                       equinox = 1888.67
                       )
    wcs3 = WCSTransform(2;
                       ctype = ["RA---AIR", "DEC--AIR"],
                       equinox = 2000
                       )
    wcs4 = WCSTransform(2;
                       ctype = ["GLON--", "GLAT--"],
                       )
    wcs5 = WCSTransform(2;
                       ctype = ["TLON", "TLAT"],
                      )
    wcs6 = WCSTransform(2;
                       ctype = ["RA---AIR", "DEC--AIR"],
                       radesys = "UNK"
                      )

    @test wcs_to_celestial_frame(wcs1) == "ICRS"
    @test wcs_to_celestial_frame(wcs2) == "FK4"
    @test wcs_to_celestial_frame(wcs3) == "FK5"
    @test wcs_to_celestial_frame(wcs4) == "Gal"
    @test wcs_to_celestial_frame(wcs5) == "ITRS"
    @test wcs_to_celestial_frame(wcs6) == "UNK"
end


text:
 import BITEModel
const BM=BITEModel
using Plots, StatPlots, Measures, Parameters

# fl = ["figs/f2-5-34-42-h.jls"][1]
# sol, sol_prior = open(fl, "r") do io
#     deserialize(io)
# end
# gl = sol.theta0.gb.gl

function plotinv1d(sol::BM.MCMCSol; ylims1=(), ylims2=(), ylims3=(0,100), ele_only=false, kwargs...)
    @unpack mc_hs1d, mc_ivs1d = sol
    mc_hs1d = mc_hs1d[:, .!bad]
    mc_ivs1d = mc_ivs1d[:, .!bad]

    gb = sol.theta0.gb
    ivs1d_measured = gb.fields1d_original[:iv]
    xl = (gb.x[1]/1e3,gb.x[end]/1e3)
    hs1d = mean(mc_hs1d,2)
    hs1d_max = maximum(mc_hs1d,2)
    hs1d_min = minimum(mc_hs1d,2)
    hs1d_std = std(mc_hs1d,2)

    ivs1d = mean(mc_ivs1d,2)
    ivs1d_max = maximum(mc_ivs1d,2)
    ivs1d_min = minimum(mc_ivs1d,2)
    ivs1d_std = std(mc_ivs1d,2)

    p1 = plot(gb.x/1e3, gb.bands_; label="surface")#, title="Flow line")
    plot!(gb.x/1e3, gb.bands_-hs1d_max, linecolor=:red, linewidth=0,
              fillrange=gb.bands_-hs1d_min, fillcolor=:red, fillalpha=0.5, label="min/max predicted bed",
              ylabel="elevation (m)")
    plot!(gb.x/1e3, gb.bands_-hs1d, label="mean predicted bed", xlims=xl, ylims=ylims1, linecolor=:black)
    if ele_only
        return p1
    end

    p2 = plot(gb.x/1e3, hs1d_max,  linecolor=:red, linewidth=0,
              fillrange=hs1d_min, fillcolor=:red, fillalpha=0.3,
              ylabel="h (m)", label="min/max")
    plot!(gb.x/1e3, hs1d+hs1d_std,  linecolor=:blue, linewidth=0,
          fillrange=hs1d-hs1d_std, fillcolor=:blue, fillalpha=0.5,
          label="std")

    plot!(gb.x/1e3, hs1d, label="mean", xlims=xl, ylims=ylims2, linecolor=:black)

    p3 = plot(gb.x/1e3, ivs1d_max,  linecolor=:red, fillrange=ivs1d_min,
              fillcolor=:red, fillalpha=0.5, label="min/max", linewidth=0,
              ylabel="v (m/a)", xlabel="χ (km)", xlims=xl, ylims=ylims3)
    plot!(gb.x/1e3, ivs1d+ivs1d_std,  linecolor=:blue, linewidth=0,
          fillrange=ivs1d-ivs1d_std, fillcolor=:blue, fillalpha=0.5,
          label="std")
    plot!(gb.x/1e3, ivs1d, label="mean", linecolor=:black)
    im = copy(ivs1d_measured)
    im[im.<0] = NaN
    im[end]=NaN
    # scatter!(gb.xmid/1e3, im, label="v 1D measured", xlims=xl)

    plot(p1,p2,p3;layout=(3,1), size=(750,900), kwargs...)
end

function plotinv1d_v(sol; kwargs...)
    @unpack mc_hs1d, mc_ivs1d = sol
    mc_hs1d = mc_hs1d[:, .!bad]
    mc_ivs1d = mc_ivs1d[:, .!bad]
    gb = sol.theta0.gb
    ivs1d_measured = gb.fields1d_original[:iv]
    xl = (gb.x[1]/1e3,gb.x[end]/1e3)
    hs1d = mean(mc_hs1d,2)
    hs1d_max = maximum(mc_hs1d,2)
    hs1d_min = minimum(mc_hs1d,2)
    hs1d_std = std(mc_hs1d,2)

    ivs1d = mean(mc_ivs1d,2)
    ivs1d_max = maximum(mc_ivs1d,2)
    ivs1d_min = minimum(mc_ivs1d,2)
    ivs1d_std = std(mc_ivs1d,2)

    # p1 = plot(gb.x/1e3, gb.bands_; label="surface")#, title="Flow line")
    # plot!(gb.x/1e3, gb.bands_-hs1d_max, linecolor=:red, linewidth=0,
    #           fillrange=gb.bands_-hs1d_min, fillcolor=:red, fillalpha=0.5, label="min/max predicted bed",
    #           ylabel="elevation (m)")
    # plot!(gb.x/1e3, gb.bands_-hs1d, label="mean predicted bed", xlims=xl, ylims=ylims1, linecolor=:black)
    # if ele_only
    #     return p1
    # end

    # p2 = plot(gb.x/1e3, hs1d_max,  linecolor=:red, linewidth=0,
    #           fillrange=hs1d_min, fillcolor=:red, fillalpha=0.3,
    #           ylabel="h (m)", label="min/max")
    # plot!(gb.x/1e3, hs1d+hs1d_std,  linecolor=:blue, linewidth=0,
    #       fillrange=hs1d-hs1d_std, fillcolor=:blue, fillalpha=0.5,
    #       label="std")

    # plot!(gb.x/1e3, hs1d, label="mean", xlims=xl, ylims=ylims2, linecolor=:black)

    p3 = plot(gb.x/1e3, ivs1d_max,  linecolor=:red, fillrange=ivs1d_min,
              fillcolor=:red, fillalpha=0.5, label="min/max", linewidth=0,
              ylabel="v (m/a)", xlabel="χ (km)", xlims=xl)
    plot!(gb.x/1e3, ivs1d+ivs1d_std,  linecolor=:blue, linewidth=0,
          fillrange=ivs1d-ivs1d_std, fillcolor=:blue, fillalpha=0.5,
          label="std")
    plot!(gb.x/1e3, ivs1d, label="mean", linecolor=:black)
    # im = copy(ivs1d_measured)
    # im[im.<0] = NaN
    # im[end]=NaN
    # # scatter!(gb.xmid/1e3, im, label="v 1D measured", xlims=xl)

    # plot(p1,p2,p3;layout=(3,1), size=(750,900), kwargs...)
end

p1 = plotinv1d(sol, reuse=false, ele_only=true, xlabel="χ (km)", ylabel="elevation (m a.s.l.)");
annotate!(0.5,2000, text("A", 25), xlabel="χ (km)", ylabel="elevation (m a.s.l.)");

p2 = plotinv1d_v(sol, reuse=false, xlabel="χ (km)");
annotate!(0.5,30, text("B", 25), xlabel="χ (km)", ylabel="speed (ma⁻¹)");

mh = sol.vols./BM.area(gl)
p3 = histogram(mh[.!bad], xlabel="mean thickness (m)", label="", bins=30,
               xlims=(minimum(mh[.!bad])-1, maximum(mh+1)), ylabel="frequency",
               grid=false);
yticks!(Int[]);
annotate!(minimum(mh)+1,300, text("C", 25));

p4 = BM.plottheta(sol.thetas[:, .!bad], sol.theta0, reuse=false, plt=corrplot2, left_margin=0mm, nsize=5,
                  #labels=["b̃₂", "fₛₗ,₂", "T", "d", "dᵥ", "σₕₘ"],
                  labels=["b̃₂", "fₛₗ,₂", "T", "d", "dᵥ", "σₕₘ", "σᵥₘ"],
#                  xticks=([-0.4:0.2:0.4;], [-0.4:0.2:0.4;], [-0.4:0.2:0.4;], [-0.4:0.2:0.4;], [-0.4:0.2:0.4;], [-0.4:0.2:0.4;], ),
                  );

# pC = plot(right_margin=0mm, left_margin=0mm, grid=:none, framestyle=:none);
# annotate!(0,0.9, text("D", 25));

# l = @layout [grid(3,1){0.4w}  a{0.001w} b{0.5w}]
# p = plot(p1,p2, p3,pC, p4, layout=l, size=(1350,900), reuse=false)

l = @layout [grid(3,1){0.48w} a{0.02w} b{0.5w}]
kwargsj = [(:tickfontsize, 13), (:legendfontsize, 13),
           (:xguidefontsize, 15), (:yguidefontsize, 15),
           (:xtickfontrotation, pi/2)]

p = plot(p1,p2, p3, plot(framestyle=:none), p4, layout=l, size=(1350,900), reuse=true; kwargsj...)

#savefig("figs/f3.png")
#p


text:
# Before far field patterns and the pair field patterns are calculated here. These are needed for the effective wavenumber and reflection for low volume fraction.

d2D(x,m) = diffbesselj(m,x)*diffhankelh1(m,x) + (1.0 - (m/x)^2)*besselj(m,x)*hankelh1(m,x)

d3D(x,m) = x * diffsbesselj(m,x) * (x * diffshankelh1(m,x) + shankelh1(m,x)) + (x^2 + - m * (m+1)) * sbesselj(m,x)^2

function far_field_pattern(ω::T, medium::Acoustic{T,2}, species::Species{T,2}; basis_order = 2) where T<:Number

    Zs = - get_t_matrices(medium, species, ω, basis_order)
    num_density_inv = one(T)/sum(number_density.(species))

    far_field(θ::T) where T <: Number = -num_density_inv*sum(
        number_density(species[i])*Zs[i][n+basis_order+1,n+basis_order+1]*exp(im*θ*n)
    for i in eachindex(species), n = -basis_order:basis_order)

    return far_field
end

function diff_far_field_pattern(ω::T, medium::Acoustic{T,2}, species::Species{T}; tol=1e-6, basis_order = 2, verbose = false, kws...) where T<:Number

    Zs = - get_t_matrices(medium, species, ω, basis_order)
    num_density_inv = one(T) / sum(number_density.(species))

    far_field(θ::T) where T <: Number = -num_density_inv*sum(
        number_density(species[i])*Zs[i][n+basis_order+1,n+basis_order+1]*im*n*exp(im*θ*n)
    for i in eachindex(species), n=-basis_order:basis_order)

    return far_field
end

function pair_field_pattern(ω::T, medium::Acoustic{T,2}, species::Species{T}; tol::T = T(1e-6), basis_order = 2) where T<:Number

    Zs = - get_t_matrices(medium, species, ω, basis_order)
    num_density_inv = one(T)/sum(number_density.(species))

    pair_field(θ::T) where T <: Number = -T(π)*num_density_inv^(2.0)*sum(
        begin
            a12 = (outer_radius(species[i]) * species[i].exclusion_distance + outer_radius(species[j]) * species[j].exclusion_distance )
            number_density(species[i]) * number_density(species[j]) * a12^2.0 * d2D(a12*ω/medium.c,m-n) * Zs[i][n+basis_order+1,n+basis_order+1] * Zs[j][m+basis_order+1,m+basis_order+1]*exp(im*m*θ)
        end
    for i in eachindex(species), j in eachindex(species),
    n = -basis_order:basis_order, m = -basis_order:basis_order)

    return pair_field
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
# Dense algorithms relating to the evaluation of scalar/matrix phi functions
# that are used by the exponential integrators.

"""
    phi(z,k[;cache]) -> [phi_0(z),phi_1(z),...,phi_k(z)]

Compute the scalar phi functions for all orders up to k.

The phi functions are defined as

```math
\\varphi_0(z) = \\exp(z),\\quad \\varphi_{k+1}(z) = \\frac{\\varphi_k(z) - 1}{z}
```

Instead of using the recurrence relation, which is numerically unstable, a
formula given by Sidje is used (Sidje, R. B. (1998). Expokit: a software
package for computing matrix exponentials. ACM Transactions on Mathematical
Software (TOMS), 24(1), 130-156. Theorem 1).
"""
function phi(z::T, k::Integer; cache=nothing,expmethod=ExpMethodHigham2005()) where {T <: Number}
    # Construct the matrix
    if cache == nothing
        cache = fill(zero(T), k+1, k+1)
    else
        fill!(cache, zero(T))
    end
    cache[1,1] = z
    for i = 1:k
        cache[i,i+1] = one(T)
    end
    P = _exp!(cache,expmethod)
    return P[1,:]
end

"""
    phiv_dense(A,v,k[;cache]) -> [phi_0(A)v phi_1(A)v ... phi_k(A)v]

Compute the matrix-phi-vector products for small, dense `A`. `k`` >= 1.

The phi functions are defined as

```math
\\varphi_0(z) = \\exp(z),\\quad \\varphi_{k+1}(z) = \\frac{\\varphi_k(z) - 1}{z}
```

Instead of using the recurrence relation, which is numerically unstable, a
formula given by Sidje is used (Sidje, R. B. (1998). Expokit: a software
package for computing matrix exponentials. ACM Transactions on Mathematical
Software (TOMS), 24(1), 130-156. Theorem 1).
"""
function phiv_dense(A, v, k; kwargs...)
    w = Matrix{eltype(A)}(undef, length(v), k+1)
    phiv_dense!(w, A, v, k; kwargs...)
end
"""
    phiv_dense!(w,A,v,k[;cache]) -> w

Non-allocating version of `phiv_dense`.
"""
function phiv_dense!(w::AbstractMatrix{T}, A::AbstractMatrix{T},
                     v::AbstractVector{T}, k::Integer;
                     cache=nothing, expmethod=ExpMethodHigham2005()
                     ) where {T <: Number}
    @assert size(w, 1) == size(A, 1) == size(A, 2) == length(v) "Dimension mismatch"
    @assert size(w, 2) == k+1 "Dimension mismatch"
    m = length(v)
    # Construct the extended matrix
    if cache == nothing
        cache = fill(zero(T), m+k, m+k)
    else
        @assert size(cache) == (m+k, m+k) "Dimension mismatch"
        fill!(cache, zero(T))
    end
    cache[1:m, 1:m] = A
    cache[1:m, m+1] = v
    for i = m+1:m+k-1
        cache[i, i+1] = one(T)
    end
    P = _exp!(cache,expmethod)
    # Extract results
    @views mul!(w[:, 1], P[1:m, 1:m], v)
    @inbounds for i = 1:k
        @inbounds for j = 1:m
            w[j, i+1] = P[j, m+i]
        end
    end
    return w
end

"""
    phi(A,k[;cache]) -> [phi_0(A),phi_1(A),...,phi_k(A)]

Compute the matrix phi functions for all orders up to k. `k` >= 1.

The phi functions are defined as

```math
\\varphi_0(z) = \\exp(z),\\quad \\varphi_{k+1}(z) = \\frac{\\varphi_k(z) - 1}{z}
```

Calls `phiv_dense` on each of the basis vectors to obtain the answer. If A
is `Diagonal`, instead calls the scalar `phi` on each diagonal element and the
return values are also `Diagonal`s
"""
function phi(A::AbstractMatrix{T}, k; kwargs...) where {T <: Number}
    m = size(A, 1)
    if A isa Diagonal
        out = [similar(A) for i = 1:k+1]
    else
        out = [Matrix{T}(undef, m, m) for i = 1:k+1]
    end
    phi!(out, A, k; kwargs...)
end
"""
    phi!(out,A,k[;caches]) -> out

Non-allocating version of `phi` for non-diagonal matrix inputs.
"""
function phi!(out::Vector{Matrix{T}}, A::AbstractMatrix{T}, k::Integer; caches=nothing,expmethod=ExpMethodHigham2005()) where {T <: Number}
    m = size(A, 1)
    @assert length(out) == k + 1 && all(P -> size(P) == (m,m), out) "Dimension mismatch"
    if caches == nothing
        e = Vector{T}(undef, m)
        W = Matrix{T}(undef, m, k+1)
        C = Matrix{T}(undef, m+k, m+k)
    else
        e, W, C = caches
        @assert size(e) == (m,) && size(W) == (m, k+1) && size(C) == (m+k, m+k) "Dimension mismatch"
    end
    @inbounds for i = 1:m
        fill!(e, zero(T)); e[i] = one(T) # e is the ith basis vector
        phiv_dense!(W, A, e, k; cache=C, expmethod=expmethod) # W = [phi_0(A)*e phi_1(A)*e ... phi_k(A)*e]
        @inbounds for j = 1:k+1
            @inbounds for s = 1:m
                out[j][s, i] = W[s, j]
            end
        end
    end
    return out
end
function phi!(out::Vector{Diagonal{T,V}}, A::Diagonal{T,V}, k::Integer;
              caches=nothing) where {T <: Number, V <: AbstractVector{T}}
    for i = 1:size(A, 1)
        phiz = phi(A[i, i], k; cache=caches)
        for j = 1:k+1
            out[j][i, i] = phiz[j]
        end
    end
    return out
end


text:
using JLD

include("../src/NirsGadgets.jl")

data = JLD.load("./test/resting_state.jld")
x = data["resting_state"]

y = x |>
    x -> NirsGadgets.to_1ds(x,             sample_rate = 55.55) |>
    x -> NirsGadgets.rm_spk(x, :fill_mean, window_width = 5*55) |>
    x -> NirsGadgets.rm_glb(x, :pca      , remove_ratio = 0.8 ) |>
    x -> NirsGadgets.rm_glb(x, :glb_mean                      ) |>
    x -> NirsGadgets.rm_trd(x, :linear                        ) |>
    x -> NirsGadgets.fl_pas(x, :low      , w = 1              ) |>
    x -> NirsGadgets.re_smp(x, :decimate , rate = 1/55.55     )


text:
using Oceananigans.Grids: short_show
using .TurbulenceClosures: ν₀, κ₀

mutable struct Model{TS, E, A<:AbstractArchitecture, G, T, B, R, SW, U, C, Φ, F,
                     BCS, S, K, OW, DI, Θ} <: AbstractModel

           architecture :: A         # Computer `Architecture` on which `Model` is run
                   grid :: G         # Grid of physical points on which `Model` is solved
                  clock :: Clock{T}  # Tracks iteration number and simulation time of `Model`
               buoyancy :: B         # Set of parameters for buoyancy model
               coriolis :: R         # Set of parameters for the background rotation rate of `Model`
          surface_waves :: SW        # Set of parameters for surfaces waves via the Craik-Leibovich approximation
             velocities :: U         # Container for velocity fields `u`, `v`, and `w`
                tracers :: C         # Container for tracer fields
              pressures :: Φ         # Container for hydrostatic and nonhydrostatic pressure
                forcing :: F         # Container for forcing functions defined by the user
                closure :: E         # Diffusive 'turbulence closure' for all model fields
    boundary_conditions :: BCS       # Container for 3d bcs on all fields
            timestepper :: TS        # Object containing timestepper fields and parameters
         poisson_solver :: S         # Poisson Solver
          diffusivities :: K         # Container for turbulent diffusivities
         output_writers :: OW        # Objects that write data to disk
            diagnostics :: DI        # Objects that calc diagnostics on-line during simulation
             parameters :: Θ         # Container for arbitrary user-defined parameters
end

"""
    Model(;
                   grid,
           architecture = CPU(),
             float_type = Float64,
                tracers = (:T, :S),
                closure = ConstantIsotropicDiffusivity(float_type, ν=ν₀, κ=κ₀),
                  clock = Clock{float_type}(0, 0),
               buoyancy = SeawaterBuoyancy(float_type),
               coriolis = nothing,
          surface_waves = nothing,
                forcing = ModelForcing(),
    boundary_conditions = HorizontallyPeriodicSolutionBCs(),
         output_writers = OrderedDict{Symbol, AbstractOutputWriter}(),
            diagnostics = OrderedDict{Symbol, AbstractDiagnostic}(),
             parameters = nothing,
             velocities = VelocityFields(architecture, grid),
              pressures = PressureFields(architecture, grid),
          diffusivities = TurbulentDiffusivities(architecture, grid, tracernames(tracers), closure),
            timestepper = :AdamsBashforth,
         poisson_solver = PoissonSolver(architecture, PoissonBCs(boundary_conditions), grid)
    )

Construct an `Oceananigans.jl` model on `grid`.

Keyword arguments
=================
- `grid`: (required) The resolution and discrete geometry on which `model` is solved.
- `architecture`: `CPU()` or `GPU()`. The computer architecture used to time-step `model`.
- `float_type`: `Float32` or `Float64`. The floating point type used for `model` data.
- `closure`: The turbulence closure for `model`. See `TurbulenceClosures`.
- `buoyancy`: Buoyancy model parameters.
- `coriolis`: Parameters for the background rotation rate of the model.
- `forcing`: User-defined forcing functions that contribute to solution tendencies.
- `boundary_conditions`: User-defined boundary conditions for model fields. Can be either`SolutionBoundaryConditions`
  or `ModelBoundaryConditions`. See `BoundaryConditions`, `HorizontallyPeriodicSolutionBCs`, and `ChannelSolutionBCs`.
- `parameters`: User-defined parameters for use in user-defined forcing functions and boundary condition functions.
"""
function Model(;
                   grid,
           architecture = CPU(),
             float_type = Float64,
                tracers = (:T, :S),
                closure = ConstantIsotropicDiffusivity(float_type, ν=ν₀, κ=κ₀),
                  clock = Clock{float_type}(0, 0),
               buoyancy = SeawaterBuoyancy(float_type),
               coriolis = nothing,
          surface_waves = nothing,
                forcing = ModelForcing(),
    boundary_conditions = HorizontallyPeriodicSolutionBCs(),
         output_writers = OrderedDict{Symbol, AbstractOutputWriter}(),
            diagnostics = OrderedDict{Symbol, AbstractDiagnostic}(),
             parameters = nothing,
             velocities = VelocityFields(architecture, grid),
              pressures = PressureFields(architecture, grid),
          diffusivities = TurbulentDiffusivities(architecture, grid, tracernames(tracers), closure),
            timestepper = :AdamsBashforth,
         poisson_solver = PoissonSolver(architecture, PoissonBCs(boundary_conditions), grid)
    )

    if architecture == GPU()
        !has_cuda() && throw(ArgumentError("Cannot create a GPU model. No CUDA-enabled GPU was detected!"))
        if mod(grid.Nx, 16) != 0 || mod(grid.Ny, 16) != 0
            throw(ArgumentError("For GPU models, Nx and Ny must be multiples of 16."))
        end
    end

    timestepper = TimeStepper(timestepper, float_type, architecture, grid, tracernames(tracers))

    tracers = TracerFields(architecture, grid, tracers)
    validate_buoyancy(buoyancy, tracernames(tracers))

    # Regularize forcing, boundary conditions, and closure for given tracer fields
    forcing = ModelForcing(tracernames(tracers), forcing)
    boundary_conditions = ModelBoundaryConditions(tracernames(tracers), boundary_conditions)
    closure = with_tracers(tracernames(tracers), closure)

    return Model(architecture, grid, clock, buoyancy, coriolis, surface_waves, velocities, tracers,
                 pressures, forcing, closure, boundary_conditions, timestepper,
                 poisson_solver, diffusivities, output_writers, diagnostics, parameters)
end

"""Show the innards of a `Model` in the REPL."""
Base.show(io::IO, model::Model) =
    print(io, "Oceananigans.Model on a ", typeof(model.architecture), " architecture ",
                                          "(time = ",  prettytime(model.clock.time),
                                          ", iteration = ", model.clock.iteration, ") \n",
              "├── grid: ", short_show(model.grid), '\n',
              "├── tracers: ", tracernames(model.tracers), '\n',
              "├── closure: ", typeof(model.closure), '\n',
              "├── buoyancy: ", typeof(model.buoyancy), '\n',
              "├── coriolis: ", typeof(model.coriolis), '\n',
              "├── output writers: ", ordered_dict_show(model.output_writers, "│"), '\n',
              "└── diagnostics: ", ordered_dict_show(model.diagnostics, " "))

"""
    ChannelModel(; kwargs...)

Construct a `Model` with walls in the y-direction. This is done by imposing
`FreeSlip` boundary conditions in the y-direction instead of `Periodic`.

kwargs are passed to the regular `Model` constructor.
"""
ChannelModel(; boundary_conditions=ChannelSolutionBCs(), kwargs...) =
    Model(; boundary_conditions=boundary_conditions, kwargs...)

"""
    NonDimensionalModel(; N, L, Re, Pr=0.7, Ro=Inf, float_type=Float64, kwargs...)

Construct a "Non-dimensional" `Model` with resolution `N`, domain extent `L`,
precision `float_type`, and the four non-dimensional numbers:

    * `Re = U λ / ν` (Reynolds number)
    * `Pr = U λ / κ` (Prandtl number)
    * `Ro = U / f λ` (Rossby number)

for characteristic velocity scale `U`, length-scale `λ`, viscosity `ν`,
tracer diffusivity `κ`, and Coriolis parameter `f`. Buoyancy is scaled
with `λ U²`, so that the Richardson number is `Ri=B`, where `B` is a
non-dimensional buoyancy scale set by the user via initial conditions or
forcing.

Note that `N`, `L`, and `Re` are required.

Additional `kwargs` are passed to the regular `Model` constructor.
"""
function NonDimensionalModel(; grid, float_type=Float64, Re, Pr=0.7, Ro=Inf,
    buoyancy = BuoyancyTracer(),
    coriolis = FPlane(float_type, f=1/Ro),
     closure = ConstantIsotropicDiffusivity(float_type, ν=1/Re, κ=1/(Pr*Re)),
    kwargs...)

    return Model(; float_type=float_type, grid=grid, closure=closure,
                   coriolis=coriolis, tracers=(:b,), buoyancy=buoyancy, kwargs...)
end

#####
##### Utils
#####

float_type(m::AbstractModel) = eltype(model.grid)

"""
    VelocityFields(arch, grid)

Return a NamedTuple with fields `u`, `v`, `w` initialized on
the architecture `arch` and `grid`.
"""
function VelocityFields(arch, grid)
    u = FaceFieldX(arch, grid)
    v = FaceFieldY(arch, grid)
    w = FaceFieldZ(arch, grid)
    return (u=u, v=v, w=w)
end

"""
    TracerFields(arch, grid)

Return a NamedTuple with tracer fields initialized
as `CellField`s on the architecture `arch` and `grid`.
"""
function TracerFields(arch, grid, tracernames)
    tracerfields = Tuple(CellField(arch, grid) for c in tracernames)
    return NamedTuple{tracernames}(tracerfields)
end

TracerFields(arch, grid, ::Union{Tuple{}, Nothing}) = NamedTuple{()}(())
TracerFields(arch, grid, tracer::Symbol) = TracerFields(arch, grid, tuple(tracer))
TracerFields(arch, grid, tracers::NamedTuple) = tracers

tracernames(::Nothing) = ()
tracernames(name::Symbol) = tuple(name)
tracernames(names::NTuple{N, Symbol}) where N = :u ∈ names ? names[4:end] : names
tracernames(::NamedTuple{names}) where names = tracernames(names)

"""
    PressureFields(arch, grid)

Return a NamedTuple with pressure fields `pHY′` and `pNHS`
initialized as `CellField`s on the architecture `arch` and `grid`.
"""
function PressureFields(arch, grid)
    pHY′ = CellField(arch, grid)
    pNHS = CellField(arch, grid)
    return (pHY′=pHY′, pNHS=pNHS)
end

"""
    Tendencies(arch, grid, tracernames)

Return a NamedTuple with tendencies for all solution fields
(velocity fields and tracer fields), initialized on
the architecture `arch` and `grid`.
"""
function Tendencies(arch, grid, tracernames)

    velocities = (u = FaceFieldX(arch, grid),
                  v = FaceFieldY(arch, grid),
                  w = FaceFieldZ(arch, grid))

    tracers = TracerFields(arch, grid, tracernames)

    return merge(velocities, tracers)
end


text:
using Documenter, PastaQ

DocMeta.setdocmeta!(PastaQ, :DocTestSetup, :(using PastaQ); recursive=true)

sitename = "PastaQ.jl"

settings = Dict(
  :modules => [PastaQ],
  :pages => [
    "Introduction" => "index.md",
    "Documentation" => [
      "Circuit Simulator" => "Circuits.md",
      "Generation of measurement data" => "GetSamples.md",
      "Random states" => "RandomStates.md",
      "Quantum Tomography" => "QuantumTomography.md",
    ],
  ],
  :format => Documenter.HTML(; assets=["assets/favicon.ico"], prettyurls=false),
  :doctest => true,
  :checkdocs => :none,
)


text:
export NodalChannel,
NodalData,
info_dump,
read_nodal,
resample!,
resample


text:
function write(document::Markdown, block::Blockquote)
    text = block.text
    # TODO put a > before all \n
    write(document.iostream, ">$text")
end

function write(document::Markdown, author::Author)
    write(document.iostream, "")
end

function write(document::Markdown, header::Header)
   write(document.iostream, "#"^header.level*header.text*"\n\n")
end

function write(document::Markdown, table::Table)
    # only works on string data arrays so far, so data need to be converted
    header = ""
    underlines = ""
    dash = "-"
    rows = fill("", table.nrows)
    # reshape = shape the data to fit the nrows and ncolumns - TODO: check for validity
    data = reshape(table.data, table.nrows, table.ncolumns)
    for i = 1:length(table.header)
        # first get the maximum size for each column
        hlength = length(table.header[i])
        maxlength = hlength
        for j = 1:table.nrows
            datalength = length(string(data[j,i]))
            if datalength > maxlength
                maxlength = datalength
            end
        end
        # then write out the rows with spaces extending the cell to the maximum column size
        for j = 1:table.nrows
            cdata = string(data[j,i])
            nspaces = maxlength - length(cdata) + 1
            rows[j] = rows[j]*cdata*" "^nspaces
        end

        nspaces = maxlength - hlength + 1
        header = header*table.header[i]*" "^nspaces
        underlines = underlines*dash^maxlength*" "
    end

    write(document.iostream,header*"\n")
    write(document.iostream,underlines*"\n")
    for r in rows
        write(document.iostream,r*"\n")
    end
    write(document.iostream,underlines*"\n")
    caption = table.caption
    if caption != ""
        write(document.iostream,"Table: $caption\n")
    end
    write(document.iostream,"\n\n")
end

function write(document::Markdown, link::Link)
    write(document.iostream,"["*link.text*"]("*link.url*")")
end

function write(document::Markdown, figure::Figure)
    # to cover situations where the figure file has an absolute path
    if document.figurefolder != ""
        url = joinpath(document.figurefolder,figure.url)
    else
        url = figure.url
    end
    write(document.iostream,"!["*figure.caption*"]("*url*")\n\n")
end

function write(document::Markdown, code::Code)
    lang = code.language
    head = ""
    if code.identifier != ""
        id = code.identifier
        head = head*"#$id "
    end
    if lang != ""
        head = head*".$lang "
    end
    if code.numberLines
        startFrom = code.startFrom
        head = head*".numberLines startFrom=$startFrom"
    end

    if head != ""
        head = "{"*head*"}"
    end

    write(document.iostream, "~~~~~~~~"*head*"\n")
    write(document.iostream, code.code*"\n")
    write(document.iostream, "~~~~~~~~\n")
end

function write(document::Markdown, paragraph::Paragraph)
    write(document.iostream, paragraph.text*"\n")
end

function write(document::Markdown, text::String)
    write(document.iostream, text*"\n")
end


text:
"""
    get_pages(
        target::String;
        concatenate::Bool=false
    ) -> Vector{Int32}

Extract the pages from a passed PDF or all PDFs found recursively
in the passed directory.

# Arguments
- `target::String`: $(TARGET_PDF_ARGUMENT)

# Keywords
- `concatenate::Bool=false`: $(CONCATENATE_KEYWORD)

# Returns
- `Vector{Int32}`: the pages

# Throws
- Exceptions from: [`get_highlights_comments_pages`](@ref)

# Example
```jldoctest; output = false
using PDFHighlights

path_to_pdf_dir = joinpath(pathof(PDFHighlights) |> dirname |> dirname, "test", "pdf")
path_to_pdf = joinpath(path_to_pdf_dir, "TestPDF.pdf")

get_pages(path_to_pdf_dir) == get_pages(path_to_pdf) == Int32[1, 2, 3, 4, 5, 6, 7, 8, 9]

# output

true
```
"""
function get_pages(target::String; concatenate::Bool=false)::Vector{Int32}
    _, _, pages = get_highlights_comments_pages(target; concatenate)
    return pages
end


text:
module Phonon

using Comonicon: @cast
using Express.PhononWorkflow.Recipes: buildworkflow
using QuantumESPRESSOExpress
using SimpleWorkflows: run!
# using ..ExpressCommands: @load_plugin

"""
Run a `config` file, better with absolute path.

# Arguments
- `config`: the file to be run. Available extensions are `.json`, `.yaml`, `.yml` or `.toml`.
"""
@cast function run(config)
    # @load_plugin
    workflow = buildworkflow(config)
    run!(workflow)
    return workflow
end

end


text:
@testset "PlanarGraph.jl" begin
    rz = planar_rz()
    @test incident(rz, 1, :face) == [2, 4, 30]
    @test ϕ(rz, incident(rz, 1, :ϕ)[1]) == 1
end


text:
# startup file for simple ODE test problem:
#   q = t^2 + x^2

import PDESolver.solvePDE
#=
"""
  This function invokes the solver for the simpleODE equation, using the
  specified input file

  Inputs:
    input_file: a string containing the path to an input file, or just a file
                name.  If it is just a file name, it is taken to be in the
                users pwd.

  Outputs:
    mesh: an AbstractMesh object
    sbp: the SBP operator used in solving the equation
    eqn: the AbstractSolutionData object used to solve the equation
    opts: the options dictonary

"""
function run_simpleode(input_file::AbstractString)

  mesh, sbp, eqn, opts, pmesh = createObjects(input_file)
  solve_simpleODE(mesh, sbp, eqn, opts, pmesh)

  return mesh, sbp, eqn, opts
end
=#

"""
  This function creates and initializes the mesh, sbp, eqn, and opts objects

  **Inputs**

   * opts: options dictionary

  **Outputs**

   * mesh: an AbstractMesh.  The concrete type is determined by the options
          dictionary
   * sbp: an AbstractSBP.  The concrete type is determined by the options
         dictionary
   * eqn: an EulerData object
   * opts: the options dictionary
   * pmesh: mesh used for preconditioning, can be same object as mesh
"""
function createObjects(opts::Dict)

  dim = opts["dimensions"]
  dofpernode = 1

  sbp, mesh, pmesh, Tsol, Tres, Tmsh, mesh_time = createMeshAndOperator(opts, dofpernode)

  # Create a simpleODE equation object
  Tdim = dim
  eqn = SimpleODEData_{Tsol, Tres, Tdim, Tmsh}(mesh, sbp, opts)
  eqn.params.time.t_meshinit = mesh_time

  init(mesh, sbp, eqn, opts)

  return mesh, sbp, eqn, opts, pmesh
end

"""
  Given fully initialized mesh, sbp, eqn, opts, this function solves
  the simpleODE equation.  The 4 object should be obtained from
  createObjects().

  Specifically, it applies an initial condition and invokes a nonlinear
  solver according to the options dictionary.

  **Inputs**

   * mesh: an AbstractMesh
   * sbp: an AbstractSBP
   * eqn: an AdvectionData
   * opts: the options dictionary.  This must be the options dictionary returned
          by createObjects().  Changing values in the options dictionary after
          calling createObjects() results in undefined behavior.
   * pmesh: mesh used for preconditioning, can be same object as mesh.
           default value of mesh

"""
function solvePDE(mesh::AbstractMesh, sbp::AbstractSBP,
                  eqn::SimpleODEData, opts::Dict, pmesh::AbstractMesh=mesh)

  # timestepping parameters
  delta_t = opts["delta_t"]
  t_max = opts["t_max"]
  # flag determines whether to calculate u, dR/du, or dR/dx (1, 2, or 3)
  flag = opts["run_type"]

  myrank = mesh.myrank

  fill!(eqn.res, 0.0)
  fill!(eqn.res_vec, 0.0)

  q_vec = eqn.q_vec
  # Apply IC's
  println("\nEvaluating initial condition")
  ICfunc_name = opts["IC_name"]
  ICfunc = ICDict[ICfunc_name]
  println("ICfunc = ", ICfunc)
  ICfunc(mesh, sbp, eqn, opts, q_vec)
  println("finished initializing q")

  saveSolutionToMesh(mesh, real(eqn.q_vec))
  writeVisFiles(mesh, "solution_ic")

  writeSolutionFiles(mesh, sbp, eqn, opts, "IC")

  # because simpleODE is different in how it calls rk4, it does not use
  # call_nlsolver
#  call_nlsolver(mesh, sbp, eqn, opts, pmesh)

  if opts["solve"]

    # handle flags

    if flag == 1        # RK4

      # NOTE: needs to have a custom pde_pre_func & pde_post_func suited to ODEs, named ode_pre_func & ode_post_func

      rk4(evalResidual, opts["delta_t"], t_max, eqn.q_vec, eqn.res_vec, ode_pre_func, ode_post_func,
          (mesh, sbp, eqn), opts;
          majorIterationCallback=eqn.majorIterationCallback, res_tol=opts["res_abstol"], real_time=opts["real_time"])

    elseif flag == 20   # Crank-Nicolson
      @time crank_nicolson(evalResidual, opts["delta_t"], t_max, mesh, sbp, eqn,
                           opts, opts["res_abstol"], opts["real_time"])

    else

      throw(ErrorException("only RK4 and CN are implemented now; user selected different solve flag"))
      return nothing

    end

    array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
    evalResidual(mesh, sbp, eqn, opts, eqn.params.t)

    eqn.res_vec[:] = 0.0
    array3DTo1D(mesh, sbp, eqn, opts, eqn.res, eqn.res_vec)

    saveSolutionToMesh(mesh, real(eqn.q_vec))
    writeVisFiles(mesh, "solution_done")

    postproc(mesh, sbp, eqn, opts)
    # TODO: comparison with exact solution

  end   # end of if opts["solve"]

  MPI.Barrier(mesh.comm)

  return mesh, sbp, eqn, opts
end  # end function

"""
  This function does post processing, if requested by the input options.
  Typical post processing includes calculation of errors, norms of important
  quantities, writing of files. etc.

  Inputs:
    mesh
    sbp
    eqn
    opts
"""
function postproc(mesh, sbp, eqn, opts)
  if opts["do_postproc"] && opts["solve"]
    exfname = opts["exact_soln_func"]
    if haskey(ICDict, exfname)
      exfunc = ICDict[exfname]
      q_exact = zeros(Tsol, mesh.numDof)
      exfunc(mesh, sbp, eqn, opts, q_exact)

    end
  end     # end of if opts["do_postproc"]

  return nothing
end


text:
PAM_250 = Dict()
PAM_250["A"] = Dict("A" => 2, "R" => -2, "N" => 0, "D" => 0, "C" => -2, "Q" => 0, "E" => 0, "G" => 1,
                    "H" => -1, "I" => -1, "L" => -2, "K" => -1, "M" => -1, "F" => -3, "P" => 1, "S" => 1,
                    "T" => 1, "W" => -6, "Y" => -3, "V" => 0, "B" => 0, "J" => -1, "Z" => 0, "X" => -1)
PAM_250["R"] = Dict("A" => -2, "R" => 6, "N" => 0, "D" => -1, "C" => -4, "Q" => 1, "E" => -1, "G" => -3,
                    "H" => 2, "I" => -2, "L" => -3, "K" => 3, "M" => 0, "F" => -4, "P" => 0, "S" => 0,
                    "T" => -1, "W" => 2, "Y" => -4, "V" => -2, "B" => -1, "J" => -3, "Z" => 0, "X" => -1)
PAM_250["N"] = Dict("A" => 0, "R" => 0, "N" => 2, "D" => 2, "C" => -4, "Q" => 1, "E" => 1, "G" => 0,
                    "H" => 2, "I" => -2, "L" => -3, "K" => 1, "M" => -2, "F" => -3, "P" => 0, "S" => 1,
                    "T" => 0, "W" => -4, "Y" => -2, "V" => -2, "B" => 2, "J" => -3, "Z" => 1, "X" => -1)
PAM_250["D"] = Dict("A" => 0, "R" => -1, "N" => 2, "D" => 4, "C" => -5, "Q" => 2, "E" => 3, "G" => 1,
                    "H" => 1, "I" => -2, "L" => -4, "K" => 0, "M" => -3, "F" => -6, "P" => -1, "S" => 0,
                    "T" => 0, "W" => -7, "Y" => -4, "V" => -2, "B" => 3, "J" => -3, "Z" => 3, "X" => -1)
PAM_250["C"] = Dict("A" => -2, "R" => -4, "N" => -4, "D" => -5, "C" => 12, "Q" => -5, "E" => -5, "G" => -3,
                    "H" => -3, "I" => -2, "L" => -6, "K" => -5, "M" => -5, "F" => -4, "P" => -3, "S" => 0,
                    "T" => -2, "W" => -8, "Y" => 0, "V" => -2, "B" => -4, "J" => -5, "Z" => -5, "X" => -1)
PAM_250["Q"] = Dict("A" => 0, "R" => 1, "N" => 1, "D" => 2, "C" => -5, "Q" => 4, "E" => 2, "G" => -1,
                    "H" => 3, "I" => -2, "L" => -2, "K" => 1, "M" => -1, "F" => -5, "P" => 0, "S" => -1,
                    "T" => -1, "W" => -5, "Y" => -4, "V" => -2, "B" => 1, "J" => -2, "Z" => 3, "X" => -1)
PAM_250["E"] = Dict("A" => 0, "R" => -1, "N" => 1, "D" => 3, "C" => -5, "Q" => 2, "E" => 4, "G" => 0,
                    "H" => 1, "I" => -2, "L" => -3, "K" => 0, "M" => -2, "F" => -5, "P" => -1, "S" => 0,
                    "T" => 0, "W" => -7, "Y" => -4, "V" => -2, "B" => 3, "J" => -3, "Z" => 3, "X" => -1)
PAM_250["G"] = Dict("A" => 1, "R" => -3, "N" => 0, "D" => 1, "C" => -3, "Q" => -1, "E" => 0, "G" => 5,
                    "H" => -2, "I" => -3, "L" => -4, "K" => -2, "M" => -3, "F" => -5, "P" => 0, "S" => 1,
                    "T" => 0, "W" => -7, "Y" => -5, "V" => -1, "B" => 0, "J" => -4, "Z" => 0, "X" => -1)
PAM_250["H"] = Dict("A" => -1, "R" => 2, "N" => 2, "D" => 1, "C" => -3, "Q" => 3, "E" => 1, "G" => -2,
                    "H" => 6, "I" => -2, "L" => -2, "K" => 0, "M" => -2, "F" => -2, "P" => 0, "S" => -1,
                    "T" => -1, "W" => -3, "Y" => 0, "V" => -2, "B" => 1, "J" => -2, "Z" => 2, "X" => -1)
PAM_250["I"] = Dict("A" => -1, "R" => -2, "N" => -2, "D" => -2, "C" => -2, "Q" => -2, "E" => -2, "G" => -3,
                    "H" => -2, "I" => 5, "L" => 2, "K" => -2, "M" => 2, "F" => 1, "P" => -2, "S" => -1,
                    "T" => 0, "W" => -5, "Y" => -1, "V" => 4, "B" => -2, "J" => 3, "Z" => -2, "X" => -1)
PAM_250["L"] = Dict("A" => -2, "R" => -3, "N" => -3, "D" => -4, "C" => -6, "Q" => -2, "E" => -3, "G" => -4,
                    "H" => -2, "I" => 2, "L" => 6, "K" => -3, "M" => 4, "F" => 2, "P" => -3, "S" => -3,
                    "T" => -2, "W" => -2, "Y" => -1, "V" => 2, "B" => -3, "J" => 5, "Z" => -3, "X" => -1)
PAM_250["K"] = Dict("A" => -1, "R" => 3, "N" => 1, "D" => 0, "C" => -5, "Q" => 1, "E" => 0, "G" => -2,
                    "H" => 0, "I" => -2, "L" => -3, "K" => 5, "M" => 0, "F" => -5, "P" => -1, "S" => 0,
                    "T" => 0, "W" => -3, "Y" => -4, "V" => -2, "B" => 1, "J" => -3, "Z" => 0, "X" => -1)
PAM_250["M"] = Dict("A" => -1, "R" => 0, "N" => -2, "D" => -3, "C" => -5, "Q" => -1, "E" => -2, "G" => -3,
                    "H" => -2, "I" => 2, "L" => 4, "K" => 0, "M" => 6, "F" => 0, "P" => -2, "S" => -2,
                    "T" => -1, "W" => -4, "Y" => -2, "V" => 2, "B" => -2, "J" => 3, "Z" => -2, "X" => -1)
PAM_250["F"] = Dict("A" => -3, "R" => -4, "N" => -3, "D" => -6, "C" => -4, "Q" => -5, "E" => -5, "G" => -5,
                    "H" => -2, "I" => 1, "L" => 2, "K" => -5, "M" => 0, "F" => 9, "P" => -5, "S" => -3,
                    "T" => -3, "W" => 0, "Y" => 7, "V" => -1, "B" => -4, "J" => 2, "Z" => -5, "X" => -1)
PAM_250["P"] = Dict("A" => 1, "R" => 0, "N" => 0, "D" => -1, "C" => -3, "Q" => 0, "E" => -1, "G" => 0,
                    "H" => 0, "I" => -2, "L" => -3, "K" => -1, "M" => -2, "F" => -5, "P" => 6, "S" => 1,
                    "T" => 0, "W" => -6, "Y" => -5, "V" => -1, "B" => -1, "J" => -2, "Z" => 0, "X" => -1)
PAM_250["S"] = Dict("A" => 1, "R" => 0, "N" => 1, "D" => 0, "C" => 0, "Q" => -1, "E" => 0, "G" => 1,
                    "H" => -1, "I" => -1, "L" => -3, "K" => 0, "M" => -2, "F" => -3, "P" => 1, "S" => 2,
                    "T" => 1, "W" => -2, "Y" => -3, "V" => -1, "B" => 0, "J" => -2, "Z" => 0, "X" => -1)
PAM_250["T"] = Dict("A" => 1, "R" => -1, "N" => 0, "D" => 0, "C" => -2, "Q" => -1, "E" => 0, "G" => 0,
                    "H" => -1, "I" => 0, "L" => -2, "K" => 0, "M" => -1, "F" => -3, "P" => 0, "S" => 1,
                    "T" => 3, "W" => -5, "Y" => -3, "V" => 0, "B" => 0, "J" => -1, "Z" => -1, "X" => -1)
PAM_250["W"] = Dict("A" => -6, "R" => 2, "N" => -4, "D" => -7, "C" => -8, "Q" => -5, "E" => -7, "G" => -7,
                    "H" => -3, "I" => -5, "L" => -2, "K" => -3, "M" => -4, "F" => 0, "P" => -6, "S" => -2,
                    "T" => -5, "W" => 17, "Y" => 0, "V" => -6, "B" => -5, "J" => -3, "Z" => -6, "X" => -1)
PAM_250["Y"] = Dict("A" => -3, "R" => -4, "N" => -2, "D" => -4, "C" => 0, "Q" => -4, "E" => -4, "G" => -5,
                    "H" => 0, "I" => -1, "L" => -1, "K" => -4, "M" => -2, "F" => 7, "P" => -5, "S" => -3,
                    "T" => -3, "W" => 0, "Y" => 10, "V" => -2, "B" => -3, "J" => -1, "Z" => -4, "X" => -1)
PAM_250["V"] = Dict("A" => 0, "R" => -2, "N" => -2, "D" => -2, "C" => -2, "Q" => -2, "E" => -2, "G" => -1,
                    "H" => -2, "I" => 4, "L" => 2, "K" => -2, "M" => 2, "F" => -1, "P" => -1, "S" => -1,
                    "T" => 0, "W" => -6, "Y" => -2, "V" => 4, "B" => -2, "J" => 2, "Z" => -2, "X" => -1)
PAM_250["B"] = Dict("A" => 0, "R" => -1, "N" => 2, "D" => 3, "C" => -4, "Q" => 1, "E" => 3, "G" => 0,
                    "H" => 1, "I" => -2, "L" => -3, "K" => 1, "M" => -2, "F" => -4, "P" => -1, "S" => 0,
                    "T" => 0, "W" => -5, "Y" => -3, "V" => -2, "B" => 3, "J" => -3, "Z" => 2, "X" => -1)
PAM_250["J"] = Dict("A" => -1, "R" => -3, "N" => -3, "D" => -3, "C" => -5, "Q" => -2, "E" => -3, "G" => -4,
                    "H" => -2, "I" => 3, "L" => 5, "K" => -3, "M" => 3, "F" => 2, "P" => -2, "S" => -2,
                    "T" => -1, "W" => -3, "Y" => -1, "V" => 2, "B" => -3, "J" => 5, "Z" => -2, "X" => -1)
PAM_250["Z"] = Dict("A" => 0, "R" => 0, "N" => 1, "D" => 3, "C" => -5, "Q" => 3, "E" => 3, "G" => 0,
                    "H" => 2, "I" => -2, "L" => -3, "K" => 0, "M" => -2, "F" => -5, "P" => 0, "S" => 0,
                    "T" => -1, "W" => -6, "Y" => -4, "V" => -2, "B" => 2, "J" => -2, "Z" => 3, "X" => -1)
PAM_250["X"] = Dict("A" => -1, "R" => -1, "N" => -1, "D" => -1, "C" => -1, "Q" => -1, "E" => -1, "G" => -1,
                    "H" => -1, "I" => -1, "L" => -1, "K" => -1, "M" => -1, "F" => -1, "P" => -1, "S" => -1,
                    "T" => -1, "W" => -1, "Y" => -1, "V" => -1, "B" => -1, "J" => -1, "Z" => -1, "X" => -1)

# score same as fitness but without building the string
function score_gap_list_pam250(input_sequence_set, gap_list, sequence_length, score_list, gap_ext_penalty)
    sequence_length = Int(sequence_length)
    score = 0
    num_items = length(input_sequence_set)
    # start offsets for accessing with respect to current index,
    # ex: [AAA] w/ gap [2,3] = [A__AA] - a 5 item string
    # but instead of building string, we want to look at A, then gap of 2, 3, then AAA
    # so index will be 1-5, when needing second A, which should be index 4, offsets will be
    # +2 gap_list, -2 sequence_list, so index of 4 will give 2 into sequence_list and show
    # that we've looked past the end of the gap list (because we already had 2 gaps)
    gap_offset_list = fill(0, num_items)

    prev_gaps = fill(false, num_items)

    gap_list_begin = fill(0, num_items)
    gap_list_end = fill(sequence_length + 1, num_items)

    #println("starting")
    for i = 1:num_items
        for j = 1:length(gap_list[i])
            if gap_list[i][j] == gap_list_begin[i] + 1
                gap_list_begin[i] = gap_list[i][j]
                #println("here")
            else
                break
            end
        end
    end

    for i = 1:num_items
        for j = length(gap_list[i]):-1:1
            if gap_list[i][j] == gap_list_end[i] - 1
                gap_list_end[i] = gap_list[i][j]
                #println("here2")
            else
                break
            end
        end
    end

    # make the last item an overflow position
    for i = 1:num_items
        push!(gap_list[i], -1)
    end

    # save variable allocation later
    score_buffer = 0
    item_one = nothing
    item_two = nothing

    # go through sequence length, to index by column, all input sequences + gaps should = sequence length
    for i = 1:sequence_length
        score_buffer = 0

        # go through each input sequence, which is a single line of characters
        # last one will just pass through, either not doing a score or it was a gap
        # so that gap count and offset still needs to increase
        for j = 1:num_items
            # get either a gap or the character
            # if it is a gap, this layer is all gap gap_penalty
            # so save a lot of calculation by mult and continue
            if i == gap_list[j][1 + gap_offset_list[j]]
                gap_offset_list[j] += 1
                item_one = '_'
            else
                item_one = input_sequence_set[j][i - gap_offset_list[j]]
            end

            # go through each input seqence from the next item up to the last
            # final pattern of access = [1, 2, 3] => 12, 13, 23
            for k = j + 1:num_items
                # get either a blank space or the character
                if i == gap_list[k][1 + gap_offset_list[k]]
                    item_two = '_'
                else
                    item_two = input_sequence_set[k][i - gap_offset_list[k]]
                end

                if item_one == '_' || item_two == '_'

                    if item_one == '_' && (i <= gap_list_begin[j] || i >= gap_list_end[j])
                        #=
                        if i <= gap_list_begin[j]
                            score_buffer += min(abs(score_list[3]), i)
                        end

                        if i >= gap_list_end[j]
                            score_buffer += min(abs(score_list[3]), sequence_length - i)
                        end
                        =#

                        #score_buffer += 1
                        continue
                    end

                    if item_two == '_' && (i <= gap_list_begin[k] || i >= gap_list_end[k])
#=
                        if i <= gap_list_begin[k]
                            score_buffer += min(abs(score_list[3]), i)
                        end

                        if i >= gap_list_end[k]
                            score_buffer += min(abs(score_list[3]), sequence_length - i)
                        end
=#
                        #score_buffer += 1
                        continue
                    end

                    if item_one != item_two
                        ratio = length(gap_list[j])
                        count = 1
                        if item_one == '_' && i - 1 in gap_list[j]
                            for item = i - 2:-1;1
                                if item in gap_list[j]
                                    count += 1
                                end
                            end

                            score_buffer += (1 - (count / ratio)) * gap_ext_penalty
                            #score_buffer += gap_ext_penalty
                        elseif item_two == '_' && i - 1 in gap_list[k]
                            for item = i - 2:-1;1
                                if item in gap_list[k]
                                    count += 1
                                end
                            end

                            score_buffer += (1 - (count / ratio)) * gap_ext_penalty
                            #score_buffer += gap_ext_penalty
                        else
                            score_buffer += score_list[3]
                        end

                        #score_buffer += score_list[3]
                    end

                    continue
                end

                score_buffer += PAM_250[string(item_one)][string(item_two)]
            end
        end

        # dont add to score if all items ended up being a gap
        score += score_buffer
    end

    # remove the added overflow position
    for i = 1:num_items
        pop!(gap_list[i])
    end

    return score
end

function score_sequence_pam250(input_sequence_set, sequence_length, score_list, gap_ext_penalty)
    score = 0
    score_buffer = 0
    item_one = nothing
    item_two = nothing
    num_items = length(input_sequence_set)

    try
        # go through sequence length, to index by column, all input sequences + gaps should = sequence length
        for i = 1:sequence_length
            score_buffer = 0

            for j = 1:num_items - 1
                for k = j + 1:num_items
                    item_one = input_sequence_set[j][i]
                    item_two = input_sequence_set[k][i]

                    # get either a blank space or the character
                    if item_one == '_' || item_two == '_'
                        if item_one != item_two

                            if i != 1 && item_one == '_'
                                if input_sequence_set[j][i - 1] != '_'
                                    score_buffer += score_list[3]
                                else
                                    score_buffer += gap_ext_penalty
                                end
                            elseif i != 1 && item_two == '_'
                                if input_sequence_set[k][i - 1] != '_'
                                    score_buffer += score_list[3]
                                else
                                    score_buffer += gap_ext_penalty
                                end
                            else
                                score_buffer += score_list[3]
                            end

                            #score_buffer += score_list[3]
                        end
                    else
                        score_buffer += PAM_250[string(item_one)][string(item_two)]
                    end
                end
            end

            score += score_buffer
        end
    catch e
        println("ERROR")
        println(item_one)
        println(item_two)
    end

    return score
end

export score_gap_list_pam250
export score_sequence_pam250


text:
"""
Variational inference for Bayesian NMF
"""
module NMF
using Distributions
using StatsFuns.logsumexp, SpecialFunctions.digamma

export NMFModel
export sample_data, VI

####################
## Types
struct NMFModel
    a_t::Array{Float64, 2} # D x K
    b_t::Array{Float64, 2} # D x L
    a_v::Float64 # 1 dim
    b_v::Float64 # 1 dim
end

function sqsum(mat::Array, idx)
    return squeeze(sum(mat, idx), idx)
end

####################
## functions
function init(X::Array{Int64, 2}, model::NMFModel)
    D, N = size(X)
    K = size(model.a_t, 2)
    S = zeros(D, K, N)
    A_t = rand(D, K)
    B_t = rand(D, K)
    A_v = rand(K, N)
    B_v = rand(K, N)
    for d in 1 : D
        for k in 1 : K
            for n in 1 : N
                S[d,k,n] = X[d,n] * A_t[d,k] * B_t[d,k] * A_v[k,n] * B_v[k,n]
            end
        end
    end
    return S, A_t, B_t, A_v, B_v
end

function update_S(X::Array{Int64, 2}, A_t::Array{Float64, 2}, B_t::Array{Float64, 2}, A_v::Array{Float64, 2}, B_v::Array{Float64, 2})
    D, K = size(A_t)
    N = size(A_v, 2)
    S = zeros(D, K, N)
    for d in 1 : D
        for n in 1 : N
            # K dim
            ln_P = (digamma.(A_t[d,:]) + log.(B_t[d,:])
                    + digamma.(A_v[:,n]) + log.(B_v[:,n])
                    )
            ln_P = ln_P - logsumexp(ln_P)
            S[d,:,n] = X[d,n] * exp.(ln_P)
        end
    end
    return S
end

function update_T(S::Array{Float64, 3}, A_v::Array{Float64, 2}, B_v::Array{Float64, 2}, model::NMFModel)
    D, K, N = size(S)
    a_t = model.a_t # DxK
    b_t = model.b_t # DxK
    A_t = a_t + sqsum(S, 3)
    B_t = (a_t ./ b_t + repmat(sqsum(A_v.*B_v, 2)', D, 1)).^(-1)
    return A_t, B_t
end

function update_V(S::Array{Float64, 3}, A_t::Array{Float64, 2}, B_t::Array{Float64, 2}, model::NMFModel)
    a_v = model.a_v
    b_v = model.b_v
    D, K, N = size(S)
    A_v = a_v + sqsum(S, 1)
    B_v = (a_v / b_v + repmat(sqsum(A_t.*B_t, 1), 1, N)).^(-1)
    return A_v, B_v
end

"""
Sample data given hyperparameters.
"""
function sample_data(N::Int, model::NMFModel)
    # TODO; check b or 1/b ?
    D, K = size(model.a_t)

    T = zeros(D, K)
    for d in 1 : D
        for k in 1 : K
            T[d,k] = rand(Gamma(model.a_t[d,k], 1.0/model.b_t[d,k])) # TODO: check
        end
    end

    V = reshape(rand(Gamma(model.a_v, 1.0/model.b_v), K*N) , K, N) # TODO: check

    S = zeros(D, K, N)
    for d in 1 : D
        for k in 1 : K
            for n in 1 : N
                S[d,k,n] = T[d,k] * V[k,n]
            end
        end
    end
    #X = sqsum(S, 2) + 0.0 # zero noise
    X = sqsum(S, 2)
    return X, T, S, V
end

function update_model(A_t::Array{Float64, 2}, B_t::Array{Float64, 2}, model::NMFModel)
    return NMFModel(A_t, B_t, model.a_v, model.b_v)
end

"""
Compute variational posterior distributions.
"""
function VI(X::Array{Int64, 2}, model::NMFModel, max_iter::Int)
    K = size(model.a_t, 2)
    D, N = size(X)
    S, A_t, B_t, A_v, B_v = init(X, model)
    for iter in 1 : max_iter
        # latent
        S = update_S(X, A_t, B_t, A_v, B_v)
        A_v, B_v = update_V(S, A_t, B_t, model)

        # param
        A_t, B_t = update_T(S, A_v, B_v, model)
    end

    return update_model(A_t, B_t, model), S, A_t.*B_t, A_v.*B_v
end

end


text:
module Day3

export parse_input
export build_used_patches
export part_1
export part_2

struct Rectangle
    id::Int
    left::Int
    top::Int
    width::Int
    height::Int
end

function Rectangle(m::RegexMatch)
	id, left, top, width, height = map(x -> parse(Int, x), m.captures)
	return Rectangle(id, left, top, width, height)
end

function parse_input(input::AbstractString)::Array{Rectangle}
	re = r"#(\d+) @ (\d+),(\d+): (\d+)x(\d+)"
	return map(l -> Rectangle(match(re, l)), split(strip(input), "\n"))
end

function build_used_patches(rectangles::Array{Rectangle})

	patches_used = Dict{Tuple{Int, Int}, Int}()
	for rectangle in rectangles
		for w in 0:rectangle.width - 1
			for h in 0:rectangle.height - 1
				x = rectangle.left + w
				y = rectangle.top + h
				if !in((x, y), keys(patches_used))
					patches_used[(x, y)] = 0
				end
				patches_used[(x, y)] += 1
			end
		end
	end
	return patches_used
end

function part_1(patches_used::Dict{Tuple{Int, Int}, Int})
	return count(x -> x > 1, values(patches_used))
end

function part_2(rectangles::Array{Rectangle}, patches_used::Dict{Tuple{Int, Int}, Int})

	for rectangle in rectangles
		alone = true
		for w in 0:rectangle.width - 1
			if !alone
				break
			end
			for h in 0:rectangle.height - 1
				x = rectangle.left + w
				y = rectangle.top + h
				if patches_used[(x, y)] > 1
					alone = false
					break
				end
			end
		end
		if alone
			return rectangle.id
		end
	end
end

end


text:
function mk_original_test(x, α = 0.05)
    n = length(x)
    s = mk_score(x)
    var_s = variance_s(x,length(x))
    τ = s/(.5*n*(n-1))
    z = z_score(s, var_s)

    p, h, trend = p_value(z, α)

    return (trend=trend,  h=h, p=p, z=z, τ=τ, s=s, var_s=var_s)
end


text:
using DiffEqEnvironments
using ReinforcementLearningBase
using Test

# Test uniform_sampler on multivariate ICs
s0_lb = -[1,2,3]
s0_ub = [4,5,6]
sampler = UniformSampler(s0_lb, s0_ub)

@test sampler isa DiffEqEnvironments.ICSampler
for _ ∈ 1:10
    @test sampler() in MultiContinuousSpace(s0_lb, s0_ub)
end

# Test uniform_sampler on univariate ICs
s0_lb = -1
s0_ub = 2
sampler = UniformSampler(s0_lb, s0_ub)

@test sampler isa DiffEqEnvironments.ICSampler
for _ ∈ 1:10
    @test sampler()[] in ContinuousSpace(s0_lb, s0_ub)
end

# Sample constant ICs by using lb = ub (used in default constructor)
s0_lb = [1,2,3]
s0_ub = s0_lb
sampler = UniformSampler(s0_lb, s0_ub)

@test sampler() == s0_ub

# Test whether ArgumentError triggers
s0_lb = [1,2]
s0_ub = [1,2,3]
@test_throws ArgumentError UniformSampler(s0_lb, s0_ub)