text:
using AMRVW
using Test

const  A = AMRVW

T = Float64
S = Complex{T}

## Test roots function
p4 = Float64[24, -50, 35, -10, 1] #a0 a1 a2 a4 a4=1
p5 = Float64[-120.0, 274.0, -225.0, 85.0, -15.0, 1.0] #a0,a1,a2,a4,a4,a5=1
p6 = Float64[720, -1764, 1624, -735, 175, -21, 1]
pc = Complex{Float64}[40.0 + 10.0im, -76.0 + 60.0im, 10.0 - 95.0im, 25.0 + 40.0im, -10.0 - 5.0im, 1.0 + 0.0im]
pc_rts = S[im, 1+im,2+im, 3.0+im, 4+im]

@testset "simple cases" begin

    ## constants
    T = Float64; S = Complex{T}
    ps = [1.0]
    rts = T[]
    @test A.roots(ps) == rts
    @test A.roots(vcat(zeros(2), ps)) == vcat(rts, zeros(T,2))
    @test A.roots(vcat(ps, zeros(2))) == rts
    @test A.roots(vcat(zeros(2), ps, zeros(2))) == vcat(rts, zeros(T,2))

    ps = [1.0, 1.0]
    rts = [-1.0]
    @test A.roots(ps) == rts
    @test A.roots(vcat(zeros(2), ps)) == vcat(rts, zeros(T,2))
    @test A.roots(vcat(ps, zeros(2))) == rts
    @test A.roots(vcat(zeros(2), ps, zeros(2))) == vcat(rts, zeros(T,2))

    ps = [1.0, 2, 1]
    rts = [-1.0, -1.0]
    @test A.roots(ps) == rts
    @test A.roots(vcat(zeros(2), ps)) == vcat(rts, zeros(T,2))
    @test A.roots(vcat(ps, zeros(2))) == rts
    @test A.roots(vcat(zeros(2), ps, zeros(2))) == vcat(rts, zeros(T,2))

    ###

    ps = S[1.0]
    rts = S[]
    @test A.roots(ps) == rts
    @test A.roots(vcat(zeros(2), ps)) == vcat(rts, zeros(S,2))
    @test A.roots(vcat(ps, zeros(2))) == rts
    @test A.roots(vcat(zeros(2), ps, zeros(2))) == vcat(rts, zeros(S,2))

    ps = S[1.0, 1.0]
    rts = S[-1.0]
    @test A.roots(ps) == rts
    @test A.roots(vcat(zeros(2), ps)) == vcat(rts, zeros(S,2))
    @test A.roots(vcat(ps, zeros(2))) == rts
    @test A.roots(vcat(zeros(2), ps, zeros(2))) == vcat(rts, zeros(S,2))

    ps = S[1.0, 2, 1]
    rts = S[-1.0, -1.0]
    @test A.roots(ps) == rts
    @test A.roots(vcat(zeros(2), ps)) == vcat(rts, zeros(S,2))
    @test A.roots(vcat(ps, zeros(2))) == rts
    @test A.roots(vcat(zeros(2), ps, zeros(2))) == vcat(rts, zeros(S,2))

    ## Pencil
    ps = S[1.0, 1]
    rts = S[-1.0]

    vs, ws = A.basic_pencil(ps)
    @test A.roots(vs, ws) == rts
    vs, ws = A.basic_pencil(vcat(zeros(2), ps))
    @test A.roots(vs, ws) == vcat(rts, zeros(S,2))
    vs, ws = A.basic_pencil(vcat(ps, zeros(2)))
    @test A.roots(vs, ws) == rts
    vs, ws = A.basic_pencil(vcat(zeros(2), ps, zeros(2)))
    @test A.roots(vs, ws) == vcat(rts, zeros(T,2))

    ps = S[1.0, 2, 1]
    rts = S[-1.0, -1.0]

    vs, ws = A.basic_pencil(ps)
    @test A.roots(vs, ws) == rts
    vs, ws = A.basic_pencil(vcat(zeros(2), ps))
    @test A.roots(vs, ws) == vcat(rts, zeros(S,2))
    vs, ws = A.basic_pencil(vcat(ps, zeros(2)))
    @test A.roots(vs, ws) == rts
    vs, ws = A.basic_pencil(vcat(zeros(2), ps, zeros(2)))
    @test A.roots(vs, ws) == vcat(rts, zeros(T,2))

end

@testset "RDS" begin

    @test all(A.roots(p4) .≈ 1.0:4)
    @test all(A.roots(p5) .≈ 1.0:5)
    @test all(A.roots(p6) .≈ 1.0:6)

    rs = rand(T, 500)
    rts = A.roots(rs)
    @test all(!iszero(rts))

end

@testset "CSS" begin

    T = Float64
    @test all(isapprox.(A.roots(pc), pc_rts, atol=sqrt(eps(T))))

    rs = rand(S, 500)
    rts = A.roots(rs)
    @test all(!iszero(rts))

end

@testset "RDS, QZ" begin

    vs, ws = A.basic_pencil(p4)
    @test all(A.roots(vs, ws) .≈ 1.0:4)

    vs, ws = A.basic_pencil(p5)
    @test all(A.roots(vs, ws) .≈ 1.0:5)

    vs, ws = A.basic_pencil(p6)
    @test all(A.roots(vs, ws) .≈ 1.0:6)

    rs = rand(T, 500)
    vs, ws = A.basic_pencil(rs)
    rts = A.roots(vs, ws)
    @test all(!iszero(rts))

end

@testset "CSS, QZ" begin

    vs, ws = A.basic_pencil(pc)
    @test all(round.(A.roots(vs, ws), digits=12) .≈  pc_rts)

    rs = rand(S, 500)
    vs, ws = A.basic_pencil(rs)
    rts = A.roots(vs, ws)
    @test all(!iszero(rts))

    function pencil_split(ps, n)
        vs = zeros(eltype(ps), length(ps)-1)
        ws = zeros(eltype(ps), length(ps)-1)

        vs[1:n] = ps[1:n]
        ws[n:end] = ps[n+1:end]

        vs, ws
    end

    ps = [24.0, -50.0, 35.0, -10.0, 1.0]
    for i in 1:4
        resid = eigvals(A.amrvw(pencil_split(ps, i)...)) - [1,2,3,4]
        @test norm(resid) <= sqrt(eps())
    end

end

@testset "qr_factorization" begin

    # unitary=false
    n = 15
    M = triu(rand(T, n, n), -1)
    e1 = eigvals(M)
    e2 = eigvals(A.qr_factorization(M, unitary=false))
    @test norm(e1 - e2) <= sqrt(eps())

    n = 15
    M = triu(rand(S, n, n), -1)
    e1 = eigvals(M)
    e2 = eigvals(A.qr_factorization(M, unitary=false))
    @test norm(e1 - e2) <= sqrt(eps())

    # unitary = true
    Qs = A.random_rotator.(T, 1:n-1)
    D = diagm(0 => rand([-1.0, 1.0], n))
    M = Qs *  D
     e1 = eigvals(M)
    e2 = eigvals(A.qr_factorization(M, unitary=false))
    @test norm(e1 - e2) <= sqrt(eps())

    Qs = A.random_rotator.(S, 1:n-1)
    D = diagm(0 => [complex(sincos(t)...) for t in rand(n)])
    M = Qs *  D
    e1 = eigvals(M)
    e2 = eigvals(A.qr_factorization(M, unitary=false))
    @test norm(e1 - e2) <= sqrt(eps())

end


text:
import .API: CustomRuleType
import LLVM: refcheck

LLVM.@checked struct TypeAnalysis
    ref::API.EnzymeTypeAnalysisRef
end
Base.unsafe_convert(::Type{API.EnzymeTypeAnalysisRef}, ta::TypeAnalysis) = ta.ref
LLVM.dispose(ta::TypeAnalysis) = API.FreeTypeAnalysis(ta)

function TypeAnalysis(logic, typerules::Dict{String, CustomRuleType}=Dict{String,CustomRuleType}())
    rulenames = String[]
    rules = CustomRuleType[]
    for (rulename, rule) in typerules
        push!(rulenames, rulename)
        push!(rules, rule)
    end
    ref = API.CreateTypeAnalysis(logic, rulenames, rules)
    TypeAnalysis(ref)
end

# typedef bool (*CustomRuleType)(int /*direction*/, CTypeTree * /*return*/,
#                                CTypeTree * /*args*/, size_t /*numArgs*/,
#                                LLVMValueRef)=T

ext: jl


text:
function parse_instructions(instructions::Vector{AbstractBFInstruction})
    expressions = Vector{Expr}()
    loops = Vector{Expr}()
    for (index, instruction) in enumerate(instructions)
        if instruction isa LoopEnd
            isempty(loops) && error("Loop close at $index does not have matching loop opening")
            expression = pop!(loops)
        else
            expression = parse(instruction)
        end
        if instruction isa LoopStart
            push!(loops, expression)
        elseif isempty(loops)
            push!(expressions, expression)
        else
            push!(last(last(loops).args).args, expression)
        end
    end
    return expressions
end

Base.parse(::Increment) = :(memory[index] += one(UInt8))
Base.parse(::Decrement) = :(memory[index] -= one(UInt8))
Base.parse(::RightShift) = :(index += 1)
Base.parse(::LeftShift) = :(index -= 1)
Base.parse(::Write) = :(print(output, Char(memory[index])))
Base.parse(::Read) = :(memory[index] = eof(characters) ? memory[index] : first(read(characters, 1)))
Base.parse(::LoopStart) = :(while !iszero(memory[index]) end)


text:
"""
# Module CpuHints

The module `CpuHints` enables to emit low-level CPU instructions that aim at
optimising or manipulating how data caches are operated and where instruction
re-ordering related to memory access is undesirable.

Primary targeted is to improve benchmarking, but potentially also for
low-level optimisations of tight memory or bandwidth bound algorithms.
"""
module CpuHints

export lfence, sfence, mfence, prefetch, prefetcht0, prefetcht1, prefetcht2,
       prefetchnta, prefetchw, clflush, clflushopt, clwb, reorder_barrier,
       elimination_barrier

using Base:  llvmcall, @_inline_meta
using CpuId: cpufeature, CLWB, CLFSH, CLFLUSH, PREFETCHW, SSE, SSE2

"""
    lfence()

Load memory fence: The LFENCE instruction prevents the CPU from re-ordering
speculative loads from memory to pass the fence.  That is, loads will not be
performed prior to the fence if the load instruction is located after the
fence, and vice versa.

The LFENCE instruction requires SSE2 extensions.
"""
function lfence end

__lfence() = (@_inline_meta; llvmcall(
      raw"""
        tail call void asm sideeffect "lfence", "~{memory},~{dirflag},~{fpsr},~{flags}"()
        ret void
      """
    , Void, Tuple{}
   ))

"""
    sfence()

Store memory fence: The SFENCE instruction prevents the CPU from re-ordering
deferred stores to memory to pass the fence.  That is, stores will not be
performed after the fence if the store instruction is located before the
fence, and vice versa.

The SFENCE instruction requires SSE extensions.
"""
function sfence end

__sfence() = (@_inline_meta; llvmcall(
      raw"""
        tail call void asm sideeffect "sfence", "~{memory},~{dirflag},~{fpsr},~{flags}"()
        ret void
      """
    , Void, Tuple{}
   ))

"""
    mfence()

Load and store memory fence: The MFENCE instruction prevents the CPU from
re-ordering deferred stores and speculative loads to and from the memory to
pass the fence.  That is, stores will not be performed after the fence if the
store instruction is located before the fence, and loads will not be performed
prior to the fence if the load instruction is located after the fence, and vice
versa.  The MFENCE instruction is a such the combination of an LFENCE and
SFENCE.

The MFENCE instruction requires SSE2 extensions.
"""
function mfence end

__mfence() = (@_inline_meta; llvmcall(
      raw"""
        tail call void asm sideeffect "mfence", "~{memory},~{dirflag},~{fpsr},~{flags}"()
        ret void
      """
    , Void, Tuple{}
   ))

"""
    prefetch(ptr::Ptr)

Load the memory block of size of one cache line that associated with the given
address into a higher cache level.  The target cache level is CPU-wise
implementation dependent.
"""
prefetch(ptr::Ptr) = (@_inline_meta; llvmcall(
      raw"""
        call void asm sideeffect "prefetch $0", "*m,~{dirflag},~{fpsr},~{flags}"(i8* nonnull %0)
        ret void
      """
    , Void, Tuple{Ptr{Void}}
    , Ptr{Void}(ptr)
   ))

"""
    prefetcht0(ptr::Ptr)

Load the memory block of size of one cache line that is associated with the
given address into the lowest cache hierarchy (L1 cache).  The CPU is free to
ignore this request.
"""
prefetcht0(ptr::Ptr) = (@_inline_meta; llvmcall(
      raw"""
        call void asm sideeffect "prefetcht0 $0", "*m,~{dirflag},~{fpsr},~{flags}"(i8* nonnull %0)
        ret void
      """
    , Void, Tuple{Ptr{Void}}
    , Ptr{Void}(ptr)
   ))

"""
    prefetcht1(ptr::Ptr)

Load the memory block of size of one cache line that is associated with the
given address into the second lowest cache hierarchy (L2 cache).  The CPU is free to
ignore this request.
"""
prefetcht1(ptr::Ptr) = (@_inline_meta; llvmcall(
      raw"""
        call void asm sideeffect "prefetcht1 $0", "*m,~{dirflag},~{fpsr},~{flags}"(i8* nonnull %0)
        ret void
      """
    , Void, Tuple{Ptr{Void}}
    , Ptr{Void}(ptr)
   ))

"""
    prefetcht2(ptr::Ptr)

Load the memory block of size of one cache line that is associated with the
given address into the third lowest cache hierarchy (L3 cache).  The CPU is free to
ignore this request.
"""
prefetcht2(ptr::Ptr) = (@_inline_meta; llvmcall(
      raw"""
        call void asm sideeffect "prefetcht2 $0", "*m,~{dirflag},~{fpsr},~{flags}"(i8* nonnull %0)
        ret void
      """
    , Void, Tuple{Ptr{Void}}
    , Ptr{Void}(ptr)
   ))

"""
    prefetchnta(ptr::Ptr)

Prepare to load the memory block of size of one cache line that is associated
with the given address into the CPU, bypassing the cache hierarchy
(non-temporal load). The CPU is free to ignore this request.
"""
prefetchnta(ptr::Ptr) = (@_inline_meta; llvmcall(
      raw"""
        call void asm sideeffect "prefetchnta $0", "*m,~{dirflag},~{fpsr},~{flags}"(i8* nonnull %0)
        ret void
      """
    , Void, Tuple{Ptr{Void}}
    , Ptr{Void}(ptr)
   ))

"""
    prefetchw(ptr::Ptr)

Load the memory block of size of one cache line that is associated
with the given address into the CPU, in anticipation of a subsequent write to
same address.  This also invalidates the cache line for other cores the cache
line might be shared with.  The CPU is free to ignore this request.

Requires a CPU having the cpuid feature flag PREFETCHW set.
"""
function prefetchw end

__prefetchw(ptr::Ptr) = (@_inline_meta; llvmcall(
      raw"""
        call void asm sideeffect "prefetchw $0", "*m,~{dirflag},~{fpsr},~{flags}"(i8* nonnull %0)
        ret void
      """
    , Void, Tuple{Ptr{Void}}
    , Ptr{Void}(ptr)
   ))

"""
    clflush(ptr::Ptr)

Write and invalidate a cache line from all levels of the CPU's cache
hierarchy, if the associated cache line is present.

Requires a CPU having the cpuid feature flag CLFSH set, leaf 0x01, EDX[19].
"""
function clflush end

__clflush(ptr::Ptr) = (@_inline_meta; llvmcall(
      raw"""
        call void asm sideeffect "clflush $0", "*m,~{dirflag},~{fpsr},~{flags}"(i8* nonnull %0)
        ret void
      """
    , Void, Tuple{Ptr{Void}}
    , Ptr{Void}(ptr)
   ))

"""
    clflushopt(ptr::Ptr)

Write and invalidate a cache line from all levels of the CPU's cache
hierarchy, if the associated cache line is present, with optimised memory
throughput.

Requires a CPU having the cpuid feature flag CLFLUSH set, leaf 0x07, EBX[23].
"""
function clflushopt end

__clflushopt(ptr::Ptr) = (@_inline_meta; llvmcall(
      raw"""
        call void asm sideeffect "clflush $0", "*m,~{dirflag},~{fpsr},~{flags}"(i8* nonnull %0)
        ret void
      """
    , Void, Tuple{Ptr{Void}}
    , Ptr{Void}(ptr)
   ))

"""
    clwb(ptr)

Write a cache line back to main memory, if modified, but may retain the data
in the cache hierarchy for future use.

Requires a CPU having the cpuid feature flag CLWB set, leaf 0x07, EBX[24].
"""
function clwb end

__clwb(ptr::Ptr) = (@_inline_meta; llvmcall(
      raw"""
        call void asm sideeffect "clwb $0", "*m,~{dirflag},~{fpsr},~{flags}"(i8* nonnull %0)
        ret void
      """
    , Void, Tuple{Ptr{Void}}
    , Ptr{Void}(ptr)
   ))

#
#  Compiler-related functions
#

"""
    reorder_barrier()

Prevent LLVM from reordering any instructions from before to after or vice
versa relative to this barrier.  This is attained by emitting a fake
instruction that pretends to manipulate all of the memory, and LLVM couldn't
possibly infer what has happened, and, hence, must keep everything in
sequence.
"""
reorder_barrier() = (@_inline_meta; llvmcall(
      raw"""
        tail call void asm sideeffect "", "~{memory},~{dirflag},~{fpsr},~{flags}"()
        ret void
      """
    , Void, Tuple{}
   ))

"""
    elimination_barrier(ptr::Ptr)

Prevent LLVM from eliminating function calls related to data associated
with memory address `ptr` by emitting a fake instruction that pretends to
manipulate that piece of memory.  This is, so to say, a 'touch ptr' function.
"""
elimination_barrier(ptr::Ptr) = (@_inline_meta; llvmcall(
      raw"""
        tail call void asm sideeffect "", "imr,~{memory},~{dirflag},~{fpsr},~{flags}"(i8* %0)
        ret void
      """
    , Void, Tuple{Ptr{Void}}
    , Ptr{Void}(ptr)
   ))

"""
    __noop()
    __noop(ptr)

A do-nothing function for unsupported CPU instructions.  Chosen at module
initialisation time to prevent undesired process termination signals from the
operating system due to illegal instruction errors.
"""
function __noop end

__noop() = info("This instruction is not supported by this CPU.")
__noop(ptr::Ptr) = info("This instruction is not supported by this CPU.")

#
# Module initialisation:
# Enable/Disable all functions for which we have no CPU support
#

function __init__()
    eval( :(sfence()      = cpufeature(:SSE)       ? __sfence()      : __noop() ) )
    eval( :(lfence()      = cpufeature(:SSE2)      ? __lfence()      : __noop() ) )
    eval( :(mfence()      = cpufeature(:SSE2)      ? __mfence()      : __noop() ) )
    eval( :(prefetchw(p)  = cpufeature(:PREFETCHW) ? __prefetchw(p)  : __noop(p) ) )
    eval( :(clflush(p)    = cpufeature(:CLFSH)     ? __clflush(p)    : __noop(p) ) )
    eval( :(clflushopt(p) = cpufeature(:CLFLUSH)   ? __clflushopt(p) : __noop(p) ) )
    eval( :(clwb(p)       = cpufeature(:CLWB)      ? __clwb(p)       : __noop(p) ) )
end

end # module CpuHints


text:
to_alongs(::Val{L}, alongs::NTuple{L,TypedBool}) where {L} = alongs

const TypedBoolLike = Union{Colon, typeof(*), True, False}
to_alongs(::Val{L}, alongs::NTuple{L,TypedBoolLike}) where {L} = _to_alongs(alongs)
@inline function _to_alongs(alongs::TupleN{TypedBoolLike})
    (_to_along(first(alongs)), _to_alongs(tail(alongs))...)
end
@inline _to_alongs(alongs::Tuple{}) = ()
@inline _to_along(::Colon) = True()
@inline _to_along(::typeof(*)) = False()
@inline _to_along(b::TypedBool) = b

@pure function to_alongs(::Val{L}, dims::TupleN{Integer}) where {L}
    ntuple(dim -> (@_inline_meta; dim in dims ? True() : False()), Val(L))
end

to_alongs(::Val{L}, alongs::Tuple{}) where {L} = ntuple(_ -> False(), Val(L))

to_alongs(::Val{L}, alongs::Tuple) where {L} = throw_invalid_alongs(L, alongs)

function to_alongs(::Val{L}, ::Val{M}) where {M,L}
    (ntuple(_ -> True(), Val(M))..., ntuple(_ -> False(), Val(L-M))...)
end

@inline to_alongs(::Val{L}, alongs...) where {L} = to_alongs(Val(L), alongs)

@noinline throw_invalid_alongs(l::Integer, alongs) = argerror("Expected $l alongs. Got: $alongs.")
@noinline function throw_invalid_alongs(m::Integer, n::Integer, alongs)
    argerror("Expected $(m+n) alongs with $m sliced dimensions. Got: $alongs.")
end

# returns true iff any number of True's followed by any number of False's
iscontiguous(alongs::Tuple{}) = true
iscontiguous(alongs::Tuple{True}) = true
iscontiguous(alongs::Tuple{False}) = true
iscontiguous(alongs::Tuple{True, Vararg{False}}) = true
iscontiguous(alongs::Tuple{True, Vararg{True}}) = true
iscontiguous(alongs::Tuple{False, Vararg{False}}) = true
iscontiguous(alongs::Tuple{False, Vararg{TypedBool}}) = false
@inline iscontiguous(alongs::Tuple{True, Vararg{TypedBool}}) = iscontiguous(tail(alongs))

@pure sliceaxes(A::AbstractArray) = sliceaxes(axes(A))
@pure sliceaxes(axes::Tuple) = tuple_map(Base.Slice, axes)

####
#### Misc
####

@noinline argerror(msg::AbstractString) = throw(ArgumentError(msg))
@noinline dimerror(msg::AbstractString) = throw(DimensionMismatch(msg))

# modified from Base (https://github.com/JuliaLang/julia/blob/04cf2d5f26b7b1ce58ece7c076ff97561db01722/base/indices.jl#L189)
@inline function setindex_shape_check(szdest, szsrc)
    szdest === szsrc || throw_setindex_mismatch(szdest, szsrc)
end

@noinline function throw_setindex_mismatch(szdest, szsrc)
    if length(szdest) == 1
        dimerror("tried to assign $(szsrc[1]) elements to $(szdest[1]) destinations")
    else
        dimerror("tried to assign $(Base.dims2string(szsrc)) array to $(Base.dims2string(szdest)) destination")
    end
end

@inline function viewtype(A::AA, I::II) where {AA<:AbstractArray, II<:Tuple}
    T = Core.Compiler.return_type(view, Tuple{AA,II.parameters...})
    return isconcretetype(T) ? T : _viewtype(A, I)
end

@inline viewtype(A::AbstractArray, I...) = viewtype(A, I)

@noinline function _viewtype(A::AA, I::II) where {AA<:AbstractArray, II<:Tuple}
    try
        return @inbounds typeof(view(A, I...))
    catch e
        msg = """
        Unable to infer the return type of $(_call2str(view, AA, II.parameters...))
        and typeof($(_call2str(view, typeof(A), I...))) threw the below error.
        Try passing in valid indices or implement:
            $(_call2str(viewtype, typeof(A), typeof(I)))

        """
        printstyled(stderr, msg, color = :light_red)
        rethrow(e)
    end
end

function _call2str(f, args...)
    argstrings = map(a -> a isa Type ? "::$a" : string(a), args)
    return length(args) > 0 ? "$f($(join(argstrings, ", ")))" : string(f)
end


text:
module Kernels

using Test
using Gridap
using Gridap.Helpers
using Base.Cartesian: @nloops, @nexprs, @nref
using TensorValues

export NumberKernel
export ArrayKernel
export NumberKernelFromFunction
export ArrayKernelFromBroadcastedFunction
export compute_type
export compute_value
export compute_value!
export compute_size
export compute_ndim
export test_number_kernel
export test_array_kernel

# Interfaces

abstract type NumberKernel end

function compute_type(::NumberKernel,::Vararg{<:Type})::Type{<:NumberLike}
  @abstractmethod
end

function compute_value(::NumberKernel,::Vararg)::NumberLike
  @abstractmethod
end

abstract type ArrayKernel end

function compute_type(::ArrayKernel,::Vararg{<:Type})::Type{<:NumberLike}
  @abstractmethod
end

function compute_ndim(::ArrayKernel,::Vararg{Int})::Int
  @abstractmethod
end

function compute_size(::ArrayKernel,::Vararg{<:NTuple})::NTuple
  @abstractmethod
end

function compute_value!(::AbstractArray,::ArrayKernel,::Vararg)
  @abstractmethod
end

function compute_value(k::ArrayKernel,i::Vararg)
  s = [ _size_for_broadcast(ii) for ii in i ]
  T = _compute_T(k,i)
  N = _compute_N(k,i)
  si = compute_size(k,s...)
  r = Array{T,N}(undef,si)
  compute_value!(r,k,i...)
  r
end

function _compute_T(k,v)
  t = _compute_eltype(v...)
  T = compute_type(k,t...)
  @assert T <: NumberLike
  T
end

function _compute_N(k,v)
  d = _compute_ndims(v...)
  compute_ndim(k,d...)
end

_nd(v::CellNumber) = 0

_nd(v::CellArray{T,N}) where {T,N} = N

_nd(v::CellMap{S,M,T,N}) where {S,M,T,N} = N

_nd(v::Type{<:Map{S,M,T,N}}) where {S,M,T,N} = N

_nd(v::Type{<:AbstractArray{T,N}}) where {T,N} = N

_nd(v::Map{S,M,T,N}) where {S,M,T,N} = N

_nd(v::AbstractArray{T,N}) where {T,N} = N

_nd(v::T) where T<:NumberLike = 0

_nd(v::Type{T}) where T<:NumberLike = 0

# TODO use a generated function here
_compute_ndims(v...) = @notimplemented
_compute_ndims(v1) = (_nd(v1),)
_compute_ndims(v1,v2) = (_nd(v1),_nd(v2))
_compute_ndims(v1,v2,v3) = (_nd(v1),_nd(v2),_nd(v3))
_compute_ndims(v1,v2,v3,v4) = (_nd(v1),_nd(v2),_nd(v3),_nd(v4))
_compute_ndims(v1,v2,v3,v4,v5) = (_nd(v1),_nd(v2),_nd(v3),_nd(v4),_nd(v5))
_compute_ndims(v1,v2,v3,v4,v5,v6) = (_nd(v1),_nd(v2),_nd(v3),_nd(v4),_nd(v5),_nd(v6))

_eltype(v::Map{S,M,T}) where {S,M,T} = T

_eltype(v::AbstractArray{T}) where T = T

_eltype(v::T) where T<:NumberLike = T

_eltype(v::CellNumber{T}) where T = T

_eltype(v::CellArray{T}) where T = T

_eltype(v::CellMap{S,M,T}) where {S,M,T} = T

_eltype(v::Type{<:Map{S,M,T}}) where {S,M,T} = T

_eltype(v::Type{<:AbstractArray{T}}) where T = T

_eltype(v::Type{T}) where T<:NumberLike = T

const _et = _eltype

# TODO use a generated function here
_compute_eltype(v...) = @notimplemented
_compute_eltype(v1) = (_et(v1),)
_compute_eltype(v1,v2) = (_et(v1),_et(v2))
_compute_eltype(v1,v2,v3) = (_et(v1),_et(v2),_et(v3))
_compute_eltype(v1,v2,v3,v4) = (_et(v1),_et(v2),_et(v3),_et(v4))
_compute_eltype(v1,v2,v3,v4,v5) = (_et(v1),_et(v2),_et(v3),_et(v4),_et(v5))
_compute_eltype(v1,v2,v3,v4,v5,v6) = (_et(v1),_et(v2),_et(v3),_et(v4),_et(v5),_et(v6))

# Testers

function test_number_kernel(k::NumberKernel,o::T,i::Vararg) where T
  t = [ typeof(ii) for ii in i ]
  S = compute_type(k,t...)
  @test S == T
  r = compute_value(k,i...)
  @test r == o
end

function test_array_kernel(
  k::ArrayKernel,o::AbstractArray{T,N},i::Vararg) where {T,N}
  t = [ eltype(ii) for ii in i ]
  S = compute_type(k,t...)
  @test S == T
  s = [ _size_for_broadcast(ii) for ii in i ]
  si = compute_size(k,s...)
  @test size(o) == si
  d = [length(si) for si in s]
  n = compute_ndim(k,d...)
  @test n == length(si)
  r = Array{T,N}(undef,si)
  compute_value!(r,k,i...)
  @test r == o
  r = compute_value(k,i...)
  @test r == o
end

@inline _size_for_broadcast(a) = size(a)

@inline _size_for_broadcast(a::NumberLike) = ()

# Implementations

struct NumberKernelFromFunction{F<:Function} <: NumberKernel
  fun::F
end

function compute_type(
  self::NumberKernelFromFunction,types::Vararg{<:Type})
  T = Base._return_type(self.fun,types)
  @assert T <: NumberLike
  T
end

function compute_value(self::NumberKernelFromFunction,args::Vararg)
  self.fun(args...)
end

struct ArrayKernelFromBroadcastedFunction{F<:Function} <: ArrayKernel
  fun::F
end

function compute_type(
  self::ArrayKernelFromBroadcastedFunction,etypes::Vararg{<:Type})
  Base._return_type(self.fun,etypes)
end

function compute_size(::ArrayKernelFromBroadcastedFunction,s::Vararg{<:NTuple})
  Base.Broadcast.broadcast_shape(s...)
end

function compute_value!(
  v::AbstractArray, k::ArrayKernelFromBroadcastedFunction, a::Vararg)
  broadcast!(k.fun,v,a...)
end

function compute_ndim(::ArrayKernelFromBroadcastedFunction,d::Vararg{Int})
  n = 0
  for di in d
    n = max(n,di)
  end
  n
end

struct CellSumKernel{D} <: ArrayKernel end

compute_type(k::CellSumKernel,T::Type) = T

compute_ndim(k::CellSumKernel,nd::Int) = nd-1

@generated function compute_size(
  k::CellSumKernel{D},asize::NTuple{N,Int}) where {D,N}

  @assert N > 0
  @assert D <= N
  str = join([ "asize[$i]," for i in 1:N if i !=D ])
  Meta.parse("($str)")
end

@generated function compute_value!(
  B::AbstractArray,
  ::CellSumKernel{D},
  A::AbstractArray{T,N}) where {T,N,D}

  @assert N > 0
  @assert D <= N
  quote
    @nloops $(N-1) b B begin
      (@nref $(N-1) B b) = zero(T)
    end
    @nloops $N a A begin
      @nexprs $(N-1) j->(b_j = a_{ j < $D ? j : j+1  } )
      (@nref $(N-1) B b) += @nref $N A a
    end
  end
end

struct LinCombKernel <: ArrayKernel end

function compute_type(::LinCombKernel,::Type{T},::Type{S}) where {T,S}
  Base._return_type(outer,Tuple{T,S})
end

function compute_ndim(::LinCombKernel,n::Int,m::Int)
  @assert n == 2
  @assert m == 1
  1
end

function compute_size(::LinCombKernel,sa::NTuple{2,Int},sb::NTuple{1,Int})
  ndofs, npoints = sa
  ndofsb, = sb
  @assert ndofsb == ndofs
  (npoints,)
end

function compute_value!(
  v::AbstractArray{T},::LinCombKernel,a::AbstractArray,b::AbstractArray) where T
  ndofs, npoints = size(a)
  for i in eachindex(v)
    @inbounds v[i] = zero(T)
  end
  for j in 1:npoints
    for i in 1:ndofs
      @inbounds v[j] += outer(a[i,j],b[i])
    end
  end
end

struct VarinnerKernel <: ArrayKernel end

function compute_type(::VarinnerKernel,::Type{T},::Type{S}) where {T,S}
  Base._return_type(inner,Tuple{T,S})
end

function compute_ndim(::VarinnerKernel,n::Int,m::Int)
  (n-1) + (m-1) + 1
end

function compute_size(::VarinnerKernel,sa::NTuple{2,Int},sb::NTuple{1,Int})
  ndofs, npoints = sa
  npointsb, = sb
  @assert npoints == npointsb
  sa
end

function compute_size(::VarinnerKernel,sa::NTuple{2,Int},sb::NTuple{2,Int})
  ndofs, npoints = sa
  ndofsb, npointsb = sb
  @assert npoints == npointsb
  (ndofs,ndofsb,npoints)
end

function compute_value!(
  v::AbstractArray,::VarinnerKernel,a::AbstractArray{T,2},b::AbstractArray{S,1}) where {T,S}
  ndofs, npoints = size(a)
  for j in 1:npoints
    for i in 1:ndofs
      @inbounds v[i,j] = inner(a[i,j],b[j])
    end
  end
end

function compute_value!(
  v::AbstractArray,::VarinnerKernel,a::AbstractArray{T,2},b::AbstractArray{S,2}) where {T,S}
  ndofsa, npoints = size(a)
  ndofsb, _ = size(b)
  for k in 1:npoints
    for j in 1:ndofsb
      for i in 1:ndofsa
        @inbounds v[i,j,k] = inner(a[i,k],b[j,k])
      end
    end
  end
end

struct PhysGradKernel <: ArrayKernel end

function compute_type(::PhysGradKernel,::Type{T},::Type{S}) where {T,S}
  Base._return_type(*,Tuple{T,S})
end

function compute_ndim(::PhysGradKernel,n::Int,m::Int)
  @assert n == 1
  @assert m == 2
  2
end

function compute_size(::PhysGradKernel,sa::NTuple{1,Int},sb::NTuple{2,Int})
  npointsa, = sa
  ndofs, npoints = sb
  @assert npoints == npointsa
  sb
end

function compute_value!(
  v::AbstractArray,::PhysGradKernel,a::AbstractArray,b::AbstractArray)
  ndofs, npoints = size(b)
  for j in 1:npoints
    @inbounds invJ = inv(a[j])
    for i in 1:ndofs
      @inbounds v[i,j] = invJ*b[i,j]
    end
  end
end

struct IntegrateNumberKernel{T} <: NumberKernel end

function compute_type(
  ::IntegrateNumberKernel{T},
  ::Type{<:AbstractVector{F}},
  ::Type{<:AbstractVector{J}},
  ::Type{<:AbstractVector{W}}) where {T,F,J,W}
  S = _compute_type_ink(F,J,W)
  @assert S == T
  T
end

function compute_value(k::IntegrateNumberKernel{T},f,j,w) where T
  _compute_value_ink(f,j,w,T)
end

@inline function _compute_type_ink(::Type{F},::Type{J},::Type{W}) where {F,J,W}
  Base._return_type(*,(F,J,W))
end

@inline function _compute_type_ink(
  ::AbstractVector{F},
  ::AbstractVector{J},
  ::AbstractVector{W}) where {F,J,W}
  _compute_type_ink(F,J,W)
end

function _compute_value_ink(f,j,w,::Type{T}) where T
  r = zero(T)
  for i in 1:length(f)
    r += f[i]*j[i]*w[i]
  end
  r
end

struct IntegrateArrayKernel <: ArrayKernel end

function compute_type(
  ::IntegrateArrayKernel, ::Type{F}, ::Type{J}, ::Type{W}) where {F,J,W}
  _compute_type_ink(F,J,W)
end

function compute_ndim(::IntegrateArrayKernel,nf::Int,nj::Int,nw::Int)
  @assert nf > 1
  @assert nj == 1
  @assert nw == 1
  nf-1
end

@generated function compute_size(
  ::IntegrateArrayKernel,sf::NTuple{N,Int},sj::NTuple,sw::NTuple) where N
  M = N-1
  str = join([ "sf[$i]," for i in 1:M ])
  Meta.parse("($str)")
end

@generated function compute_value!(
  B::AbstractArray{T,N},::IntegrateArrayKernel,F,J,W) where {T,N}
  quote
    z = zero(T)
    for i in eachindex(B)
      B[i] = z
    end
    @nloops $(N+1) f F begin
      @nexprs $N k->(b_k = f_{ k } )
      @nexprs 1 k->(w_k = f_{$(N+1)} )
      (@nref $N B b) += (@nref $(N+1) F f) * ( @nref 1 J w ) * ( @nref 1 W w )
    end
  end
end

IntegrateKernel(::Val{N},::Type{T}) where {N,T} = IntegrateArrayKernel()

IntegrateKernel(::Val{1},::Type{T}) where T = IntegrateNumberKernel{T}()

IntegrateKernel(i::Int,T) = IntegrateKernel(Val(i),T)

# CellNewAxisKernel

struct CellNewAxisKernel{D} <: ArrayKernel end

compute_type(k::CellNewAxisKernel,T::Type) = T

compute_ndim(k::CellNewAxisKernel,nd::Int) = nd+1

function compute_size(k::CellNewAxisKernel{D},asize::NTuple{N,Int}) where {D,N}
  newaxis_size(Val(D),asize)
end

function compute_value!(
  v::AbstractArray,
  ::CellNewAxisKernel{D},
  a::AbstractArray{T,N}) where {T,N,D}
  newaxis_kernel!(Val(D),a,v)
end

@generated function newaxis_size(::Val{A},asize::NTuple{M,Int}) where {A,M}
  @assert A <= M+1
  str = ["asize[$i]," for i in 1:M]
  insert!(str,A,"1,")
  Meta.parse("($(join(str)))")
end

@generated function newaxis_kernel!(::Val{D}, A::AbstractArray{T,M}, B::AbstractArray{T,N}) where {D,T,M,N}
  @assert N == M + 1
  @assert D <= N
  quote
    @nloops $M a A begin
      @nexprs $N j->(b_j = j == $D ? 1 : a_{ j < $D ? j : j-1 } )
      (@nref $N B b) = @nref $M A a
    end
  end
end

end # module


text:
abstract Graph

type StaticGraph <: Graph
	fun::Function
	points::Dict{Integer, Vector2f}
	accuracy::Integer
	color::Color
	thickness::Real
end

type AnimatedGraph <: Graph
	fun::Function
	points::Dict{Integer, Vector2f}
	accuracy::Integer
	xval::Real
	startx::Real
	speed::Real
	clock::Clock
	rtime::Real
	color::Color
	thickness::Real
end

function StaticGraph(fun::Function; thickness=2, color=SFML.red)
	points = Dict{Integer, Vector2f}()
	StaticGraph(fun, points, 1, color, thickness)
end

function AnimatedGraph(fun::Function; thickness = 2, color = SFML.red, speed = 2, startx = 0)
	points = Dict{Integer, Vector2f}()
	clock = Clock()
	AnimatedGraph(fun, points, 1, startx, startx, speed, clock, 0, color, thickness)
end

function advance(graph::AnimatedGraph, ppu::Real)
	if graph.xval == graph.startx
		graph.startx *= ppu
		graph.xval *= ppu
		restart(graph.clock)
		advancex(graph)
	end
	elapsedtime = as_seconds(get_elapsed_time(graph.clock)) + graph.rtime
	interval = 1/graph.speed * graph.accuracy / ppu
	while elapsedtime >= interval
		try
			restart(graph.clock)
			result = graph.fun(graph.xval/ppu)
			if result == nothing
				advancex(graph)
				elapsedtime -= interval
				continue
			end
			pos = Vector2f(0, 0)
			if typeof(result) <: Real
				pos = Vector2f(graph.xval, ppu*result)
			else
				pos = Vector2f(ppu*result[1], ppu*result[2])
			end
			graph.points[graph.xval] = pos
		catch exception
			println(exception)
		end
		advancex(graph)
		elapsedtime -= interval
	end
	graph.rtime = elapsedtime
end

function advancex(graph::AnimatedGraph)
	graph.xval += graph.accuracy
end

function add_point(graph::Graph, index::Integer, ppu::Real)
	try
		pos = Vector2f(0, 0)
		result = graph.fun(index/ppu)
		if typeof(result) <: Real
			pos = Vector2f(index, ppu*result)
		else
			pos = Vector2f(ppu*result[1], ppu*result[2])
		end
		graph.points[index] = pos
	catch exception
	end
end

function add_point(graph::Graph, index::Integer, pos::Vector2f)
	graph.points[index] = pos
end

function draw(window::RenderWindow, graph::Graph)
	last_pos = 0
	points = graph.points
	sortedkeys = sort(collect(keys(graph.points)))
	for key in sortedkeys
		if last_pos != 0
			l = Line(points[key], last_pos, graph.thickness)
			set_fillcolor(l, graph.color)
			SFML.draw(window, l)
		end
		last_pos = points[key]
	end
end

function remove_points_outside(graph::Graph, size::Vector2f)
	for keyval in graph.points
		if !(size.x < keyval[1] < size.y)
			delete!(graph.points, keyval[1])
		end
	end
end

export Graph, StaticGraph, AnimatedGraph


text:
### Default APIs
## Environments
function State(env::AbstractEnv)
    error("Define the structure of state for your environment")
end

function Params(env::AbstractEnv)
    error("Define the structure of parameters for your environment")
end

function Dynamics!(env::AbstractEnv)
    error("Define in-place dynamics")
end

function Dynamics(env::AbstractEnv)
    @warn "It is recommended to use `Dynamics!`; see https://github.com/JinraeKim/FlightSims.jl/issues/16"
    error("Define out-of-place dynamics")
end

function Command(controller::AbstractController)
    error("Define command rule for your controller")
end


text:
module UndirectedGraphFunctions
    include("./UndirectedGraphGates.jl")
    using Reexport
    @reexport using .UndirectedGraphGates
    export create_buckets, contract_graph, gate_parser

    mutable struct StateTensor
        item::Array{ComplexF64, 1}
        target::Int
        worldline1::Int
    end

    struct Tensor
        indices::Array{Int, 1}
        item::Array{ComplexF64}
    end

    function create_buckets(_n_qregs::Int, gates::Array{UG_Gate,1}, input_state::Array{Int,1}, output_state::Array{Int,1})
        worldlines = ones(Int64, _n_qregs)
        elements = []
        for i in 1:length(gates)
            obj, worldlines = get_tensor(gates[i], worldlines)
            push!(elements, obj)
        end

        num_variables = sum(worldlines)
        buckets = []
        for i in 1:sum(num_variables)
            push!(buckets, [])
        end

        for element in elements
            if typeof(element) == OneQubitGateTensor
                element.worldline1 += sum(worldlines[1:(element.target-1)])
                element.worldline2 += sum(worldlines[1:(element.target-1)])
            elseif typeof(element) == TwoQubitGateTensor
                element.worldline1 += sum(worldlines[1:(element.control-1)])
                element.worldline2 += sum(worldlines[1:(element.target-1)])
            end
        end

        for qreg in 1:_n_qregs
            input_val = convert(Array{ComplexF64, 1}, [1 - input_state[qreg], input_state[qreg]])
            _input_state = StateTensor(input_val, qreg, 1 + sum(worldlines[1:(qreg-1)]))
            insert!(elements, qreg, _input_state)
            output_val = convert(Array{ComplexF64, 1}, [1 - output_state[qreg], output_state[qreg]])
            _output_state = StateTensor(output_val, qreg, sum(worldlines[1:qreg]))
            push!(elements, _output_state)
        end

        for element in elements
            tensor = make_tensor(element)
            first_index = element.worldline1
            push!(buckets[first_index], tensor)
        end
        return buckets
    end

    function get_tensor(gate::T, worldlines::Array{Int, 1}) where T <: UG_OneQubitGate
        mat = gate_mat(gate)
        worldline = worldlines[gate._target]
        obj = OneQubitGateTensor(mat, gate._target, worldline, worldline + 1)
        worldlines[gate._target] = worldline + 1
        return obj, worldlines
    end

    function get_tensor(gate::T, worldlines::Array{Int, 1}) where T <: UG_TwoQubitGate
        mat = gate_mat(gate)
        target1 = min(gate._control, gate._target)
        target2 = max(gate._control, gate._target)
        worldline1 = worldlines[target1]
        worldline2 = worldlines[target2]
        obj = TwoQubitGateTensor(mat, target1, target2, worldline1, worldline2)
        return obj, worldlines
    end

    function make_tensor(arr::T) where T <: UndirectedGraphTensor
        index1 = arr.worldline1
        index2 = arr.worldline2
        A = Tensor([index1, index2], arr.item)
        return A
    end

    function make_tensor(arr::StateTensor)
        index = arr.worldline1
        A = Tensor([index], arr.item)
        return A
    end

    function process_bucket(bucket)
        res = bucket[1]
        l = length(bucket)
        for tensor in bucket[2:l]
            res = broadcast_product(res, tensor)
        end
        return reduce_sum(res)
    end

    # element-wise product with broadcast, in the case size of each dim = 2
    function broadcast_product(A::Tensor, B::Tensor)
        indA = A.indices
        indB = B.indices
        indOut = sort(collect(union(Set(indA), Set(indB))))
        _indA = [Int(indOut[i] in Set(indA)) for i in 1:length(indOut)]
        _indB = [Int(indOut[i] in Set(indB)) for i in 1:length(indOut)]
        reshapeA = Tuple(_indA .+ 1)
        reshapeB = Tuple(_indB .+ 1)
        expandA = Tuple((_indA .- 1) .* -1 .+ 1)
        expandB = Tuple((_indB .- 1) .* -1 .+ 1)
        _A = repeat(reshape(A.item, reshapeA), outer=expandA)
        _B = repeat(reshape(B.item, reshapeB), outer=expandB)
        return Tensor(indOut, _A .* _B)
    end

    function reduce_sum(A::Tensor)
        arr = A.item
        l = length(A.indices)
        _arr = sum(A.item, dims=1)
        if length(size(_arr)) > 1
            _arr = dropdims(_arr; dims=1)
        end
        return Tensor(A.indices[2:l], _arr)
    end

    function contract_graph(buckets)
        result = 0
        for bucket in buckets
            if length(bucket) > 0
                tensor = process_bucket(bucket)
                if max(size(tensor.item)...) > 1
                    first_index = tensor.indices[1]
                    to_bucket = buckets[first_index]
                    push!(to_bucket, tensor)
                else
                    if result != 0
                        result = result .* tensor.item
                    else
                        result = tensor.item
                    end
                end
            end
        end
        return result[1]
    end

    function gate_mat(gate::X)
        return [0 1; 1 0]
    end

    function gate_mat(gate::Y)
        return [0 im; -im 0]
    end

    function gate_mat(gate::Z)
        return [1 0; 0 -1]
    end

    function gate_mat(gate::H)
        return 1 / sqrt(2) * [1 1; 1 -1]
    end

    function gate_mat(gate::S)
        return [1 0; 0 _exp(im*pi/2)]
    end

    function gate_mat(gate::T)
        return [1 0; 0 _exp(im*pi/4)]
    end

    function gate_mat(gate::CZ)
        return [1 1; 1 -1]
    end

    function gate_mat(gate::RX)
        theta = gate._theta
        thetapi = theta / pi
        return [cospi(thetapi/2) -im*sinpi(thetapi/2); -im*sinpi(thetapi/2) cospi(thetapi/2)]
    end

    function gate_mat(gate::RY)
        theta = gate._theta
        thetapi = theta / pi
        return [cospi(thetapi/2) sinpi(thetapi/2); -sinpi(thetapi/2) cospi(thetapi/2)]
    end

    function gate_mat(gate::RZ)
        theta = gate._theta
        return [_exp(-im*theta/2) 0; 0 _exp(im*theta/2)]
    end

    function gate_parser(gate::CX)
        parsed_gates = []
        push!(parsed_gates, H(gate._target))
        push!(parsed_gates, CZ(gate._control, gate._target))
        push!(parsed_gates, H(gate._target))
        return parsed_gates
    end

    function gate_parser(gate::T) where T <: UG_Gate
        return [gate]
    end

    function _exp(theta::ComplexF64)
        thetapi = theta / pi / im
        if thetapi == conj(thetapi)
            thetapi = convert(Float64, thetapi)
        else
            error("Invalid gate parameters.")
        end
        res = cospi(thetapi) + im * sinpi(thetapi)
        return res
    end

end

max_stars_repo_name: ryuNagai/blueqat-julia

max_issues_repo_name: ryuNagai/blueqat-julia

max_forks_repo_name: ryuNagai/blueqat-julia


text:
print_with_color(:magenta, "
 _____         _
|   __|___ ___|_|___
|  |  | -_|   | | -_|
|_____|___|_|_|_|___|

")


text:
using LinearAlgebra
using Test
using RecipesBase
using DynamicPolynomials
using SetProg, SetProg.Sets

RecipesBase.is_key_supported(k::Symbol) = true # Plots normally defines this
function recipe(set, npoints)
    result = RecipesBase.apply_recipe(Dict{Symbol, Any}(:npoints => npoints), set)
    return result[1].args
end
function recipe_test(set, exp_x, exp_y, npoints=4)
    x, y = recipe(set, npoints)
    @test x ≈ exp_x
    @test y ≈ exp_y
end

@testset "Recipe" begin
    Q = [1.0 0.0; 0.0 1.0]
    @testset "Ellipsoid" begin
        @testset "Circle" begin
            for circle in [Sets.HyperSphere(2),
                           Sets.Ellipsoid(Symmetric(Q))]
                recipe_test(circle,
                            [1.0, 0.0, -1.0, 0.0], [0.0, 1.0, 0.0, -1.0])
                recipe_test(Sets.polar(circle),
                            [1.0, 0.0, -1.0, 0.0], [0.0, 1.0, 0.0, -1.0])
            end
        end
        @testset "Shifted Circle" begin
            for circle in [Sets.HyperSphere(2),
                           Sets.Ellipsoid(Symmetric(Q))]
                shifted = Sets.Translation(circle, [1.0, 2.0])
                recipe_test(shifted, [2.0, 1.0, 0.0, 1.0], [2.0, 3.0, 2.0, 1.0])
            end
        end
        @testset "Scaled circle" begin
            scaled_circle = Sets.Ellipsoid(Symmetric(2Q))
            recipe_test(scaled_circle,
                        [1/√2, 0.0, -1/√2, 0.0], [0.0, 1/√2, 0.0, -1/√2])
            recipe_test(Sets.polar(scaled_circle),
                        [√2, 0.0, -√2, 0.0], [0.0, √2, 0.0, -√2])
        end
    end
    @testset "Polynomial" begin
        @polyvar z x y
        @testset "Circle" begin
            p = SetProg.GramMatrix{Float64}((i, j) -> convert(Float64, i == j),
                                               monovec([x, y]))
            circle = Sets.PolySet(2, p)
            recipe_test(circle,
                        [1.0, 0.0, -1.0, 0.0], [0.0, 1.0, 0.0, -1.0])
        end
        @testset "Convex Circle" begin
            p = SetProg.GramMatrix{Float64}((i, j) -> convert(Float64, i == j),
                                               monovec([x, y]))
            circle = Sets.ConvexPolySet(2, p, nothing)
            recipe_test(circle,
                        [1.0, 0.0, -1.0, 0.0], [0.0, 1.0, 0.0, -1.0])
            recipe_test(Sets.polar(circle),
                        [-1.0, -1.0, 1.0,  1.0, -1.0],
                        [-1.0,  1.0, 1.0, -1.0, -1.0])
        end
        @testset "Scaled circle" begin
            p = SetProg.GramMatrix{Float64}((i, j) -> 2convert(Float64, i == j),
                                               monovec([x, y]))
            circle = Sets.ConvexPolySet(2, p, nothing)
            recipe_test(circle,
                        [1/√2, 0.0, -1/√2, 0.0], [0.0, 1/√2, 0.0, -1/√2])
            recipe_test(Sets.polar(circle),
                        [-√2, -√2, √2,  √2, -√2],
                        [-√2,  √2, √2, -√2, -√2])
        end
        @testset "Non-homogeneous Circle" begin
            @testset "Basic" begin
                q = SetProg.GramMatrix(Float64[0 0 0
                                               0 1 0
                                               0 0 1], monovec([z, x, y]))
                shifted_circle = SetProg.perspective_dual_polyset(2, q, SetProg.InteriorPoint(zeros(2)), z, [x, y])
                recipe_test(shifted_circle,
                            [-1.0, -1.0, 1.0,  1.0, -1.0],
                            [-1.0,  1.0, 1.0, -1.0, -1.0])
            end
            @testset "Scaled" begin
                q = SetProg.GramMatrix(Float64[0 0 0
                                               0 2 0
                                               0 0 2], monovec([z, x, y]))
                shifted_circle = SetProg.perspective_dual_polyset(2, q, SetProg.InteriorPoint(zeros(2)), z, [x, y])
                recipe_test(shifted_circle,
							[-√2, -√2, √2,  √2, -√2],
							[-√2,  √2, √2, -√2, -√2])
            end
            @testset "z-Scaled" begin
                # z: -1/2 + 1 = 1/2
                q = SetProg.GramMatrix([1/2 0 0
                                        0   1 0
                                        0   0 1], monovec([z, x, y]))
                shifted_circle = SetProg.perspective_dual_polyset(2, q, SetProg.InteriorPoint(zeros(2)), z, [x, y])
				recipe_test(shifted_circle,
							[-√2, -√2, √2,  √2, -√2],
							[-√2,  √2, √2, -√2, -√2])
            end
        end
    end
    @testset "Piecewise" begin
        polytope = polyhedron(HalfSpace([1, 0], 1) ∩ HalfSpace([0, 1], 1) ∩ HalfSpace([-1, 1], 1) ∩ HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1), lib)
        Q_1 = [1.0 0.0
               0.0 0.0]
        Q_2 = [0.0 0.0
               0.0 1.0]
        Q_3 = [1.0 0.0
               0.0 1.0]
        Q_4 = [1.0 1.0
               1.0 1.0]
        Q_5 = [ 1.0 -0.5
               -0.5  1.0]
        set = Sets.Piecewise(Sets.Ellipsoid.(Symmetric.([Q_1, Q_2, Q_3, Q_4, Q_5])), polytope)
        recipe_test(
            set,
            [1.0, 0.0, -1.0, 0.0],
            [0.0, 1.0, 0.0, -1.0]
        )
        recipe_test(
            set,
            [1.0, 1.0, 0.0, -1/√2, -1.0, -0.5, -0.0,  1/√3],
            [0.0, 1.0, 1.0,  1/√2,  0.0, -0.5, -1.0, -1/√3],
            8
        )
        α = 0.4142135623730951
        β = 0.732050807568877
        recipe_test(
            Sets.polar(set),
            [-1.0, -1.0, -α, 0.0, 1.0, 1.0, β, -1.0],
            [-1.0, α, 1.0, 1.0, 0.0, -β, -1.0, -1.0],
            8)
    end
end


text:
function quibble(arr::Array)
    if isempty(arr) rst = "" else rst = "$(arr[end])" end
    if length(arr) > 1 rst = join(arr[1:end-1], ", ") * " and " * rst end
    return "{" * rst * "}"
end

@show quibble([])
@show quibble(["ABC"])
@show quibble(["ABC", "DEF"])
@show quibble(["ABC", "DEF", "G", "H"])

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module PyMNE

#####
##### Dependencies
#####

using PyCall

#####
##### Exports
#####

export get_info

######
###### Actual functionality
######

const mne = PyNULL()

# TODO: examine how to do wrappers for subsubmodules. for now it's not a huge deal
#       because the wrapper just gets put in the top-level package namespace and so
#       is still accessible

include("wrappers.jl")

function __init__()
    # all of this is __init__() so that it plays nice with precompilation
    # see https://github.com/JuliaPy/PyCall.jl/#using-pycall-from-julia-modules
    copy!(mne, pyimport("mne"))
    # don't eval into the module while precompiling; this breaks precompilation
    # of downstream modules (see #4)
    if ccall(:jl_generating_output, Cint, ()) == 0
        # delegate everything else to mne
        for pn in propertynames(mne)
            isdefined(@__MODULE__, pn) && continue
            prop = getproperty(mne, pn)
            @eval $pn = $prop
        end
    end
    return nothing
end

"""
    install_sklearn(ver="")

Install scikit-learn using the specified version.

The default version is the latest stable version.
"""
function install_sklearn(version="latest"; verbose=false)
    verbose && @info "Installing scikit-learn"
    pip = pyimport("pip")
    flags = split(get(ENV, "PIPFLAGS", ""))
    packages = ["scikit-learn" * (version == "latest" ? "" : "==$version")]
    if verbose
        @info "Package requirements:" packages
        @info "Flags for pip install:" flags
        @info "scikit-learn version:" version
    end
    pip.main(["install"; flags; packages])
    return nothing
end

end # module


text:
module Pidfile

export mkpidlock

using Base:
    IOError, UV_EEXIST, UV_ESRCH,
    Process

using Base.Libc: rand

using Base.Filesystem:
    File, open, JL_O_CREAT, JL_O_RDWR, JL_O_RDONLY, JL_O_EXCL,
    rename, samefile, path_separator

using ..FileWatching: watch_file
using Base.Sys: iswindows

"""
    mkpidlock([f::Function], at::String, [pid::Cint, proc::Process]; kwopts...)

Create a pidfile lock for the path "at" for the current process
or the process identified by pid or proc. Can take a function to execute once locked,
for usage in `do` blocks, after which the lock will be automatically closed. If the lock fails
and `wait` is false, then an error is thrown.

The lock will be released by either `close`, a `finalizer`, or shortly after `proc` exits.
Make sure the return value is live through the end of the critical section of
your program, so the `finalizer` does not reclaim it early.

Optional keyword arguments:
 - `mode`: file access mode (modified by the process umask). Defaults to world-readable.
 - `poll_interval`: Specify the maximum time to between attempts (if `watch_file` doesn't work)
 - `stale_age`: Delete an existing pidfile (ignoring the lock) if its mtime is older than this.
     The file won't be deleted until 25x longer than this if the pid in the file appears that it may be valid.
     By default this is disabled (`stale_age` = 0), but a typical recommended value would be about 3-5x an
     estimated normal completion time.
 - `refresh`: Keeps a lock from becoming stale by updating the mtime every interval of time that passes.
     By default, this is set to `stale_age/2`, which is the recommended value.
 - `wait`: If true, block until we get the lock, if false, raise error if lock fails.
"""
function mkpidlock end

# mutable only because we want to add a finalizer
mutable struct LockMonitor
    const path::String
    const fd::File
    const update::Union{Nothing,Timer}

    global function mkpidlock(at::String, pid::Cint; stale_age::Real=0, refresh::Real=stale_age/2, kwopts...)
        local lock
        atdir, atname = splitdir(at)
        isempty(atdir) && (atdir = pwd())
        at = realpath(atdir) * path_separator * atname
        fd = open_exclusive(at; stale_age=stale_age, kwopts...)
        update = nothing
        try
            write_pidfile(fd, pid)
            if refresh > 0
                # N.b.: to ensure our finalizer works we are careful to capture
                # `fd` here instead of `lock`.
                update = Timer(t -> isopen(t) && touch(fd), refresh; interval=refresh)
            end
            lock = new(at, fd, update)
            finalizer(close, lock)
        catch ex
            tryrmopenfile(at)
            close(fd)
            rethrow(ex)
        end
        return lock
    end
end

mkpidlock(at::String; kwopts...) = mkpidlock(at, getpid(); kwopts...)
mkpidlock(f::Function, at::String; kwopts...) = mkpidlock(f, at, getpid(); kwopts...)

function mkpidlock(f::Function, at::String, pid::Cint; kwopts...)
    lock = mkpidlock(at, pid; kwopts...)
    try
        return f()
    finally
        close(lock)
    end
end

function mkpidlock(at::String, proc::Process; kwopts...)
    lock = mkpidlock(at, getpid(proc); kwopts...)
    closer = @async begin
        wait(proc)
        close(lock)
    end
    isdefined(Base, :errormonitor) && Base.errormonitor(closer)
    return lock
end

"""
    Base.touch(::Pidfile.LockMonitor)

Update the `mtime` on the lock, to indicate it is still fresh.

See also the `refresh` keyword in the [`mkpidlock`](@ref) constructor.
"""
Base.touch(lock::LockMonitor) = (touch(lock.fd); lock)

"""
    write_pidfile(io, pid)

Write our pidfile format to an open IO descriptor.
"""
function write_pidfile(io::IO, pid::Cint)
    print(io, "$pid $(gethostname())")
end

"""
    parse_pidfile(file::Union{IO, String}) => (pid, hostname, age)

Attempt to parse our pidfile format,
replaced an element with (0, "", 0.0), respectively, for any read that failed.
"""
function parse_pidfile(io::IO)
    fields = split(read(io, String), ' ', limit = 2)
    pid = tryparse(Cuint, fields[1])
    pid === nothing && (pid = Cuint(0))
    hostname = (length(fields) == 2) ? fields[2] : ""
    when = mtime(io)
    age = time() - when
    return (pid, hostname, age)
end

function parse_pidfile(path::String)
    try
        existing = open(path, JL_O_RDONLY)
        try
            return parse_pidfile(existing)
        finally
            close(existing)
        end
    catch ex
        isa(ex, EOFError) || isa(ex, IOError) || rethrow(ex)
        return (Cuint(0), "", 0.0)
    end
end

"""
    isvalidpid(hostname::String, pid::Cuint) :: Bool

Attempt to conservatively estimate whether pid is a valid process id.
"""
function isvalidpid(hostname::AbstractString, pid::Cuint)
    # can't inspect remote hosts
    (hostname == "" || hostname == gethostname()) || return true
    # pid < 0 is never valid (must be a parser error or different OS),
    # and would have a completely different meaning when passed to kill
    !iswindows() && pid > typemax(Cint) && return false
    # (similarly for pid 0)
    pid == 0 && return false
    # see if the process id exists by querying kill without sending a signal
    # and checking if it returned ESRCH (no such process)
    return ccall(:uv_kill, Cint, (Cuint, Cint), pid, 0) != UV_ESRCH
end

"""
    stale_pidfile(path::String, stale_age::Real) :: Bool

Helper function for `open_exclusive` for deciding if a pidfile is stale.
"""
function stale_pidfile(path::String, stale_age::Real)
    pid, hostname, age = parse_pidfile(path)
    age < -stale_age && @warn "filesystem time skew detected" path=path
    if age > stale_age
        if (age > stale_age * 25) || !isvalidpid(hostname, pid)
            return true
        end
    end
    return false
end

"""
    tryopen_exclusive(path::String, mode::Integer = 0o444) :: Union{Void, File}

Try to create a new file for read-write advisory-exclusive access,
return nothing if it already exists.
"""
function tryopen_exclusive(path::String, mode::Integer = 0o444)
    try
        return open(path, JL_O_RDWR | JL_O_CREAT | JL_O_EXCL, mode)
    catch ex
        (isa(ex, IOError) && ex.code == UV_EEXIST) || rethrow(ex)
    end
    return nothing
end

"""
    open_exclusive(path::String; mode, poll_interval, stale_age) :: File

Create a new a file for read-write advisory-exclusive access.
If `wait` is `false` then error out if the lock files exist
otherwise block until we get the lock.

For a description of the keyword arguments, see [`mkpidlock`](@ref).
"""
function open_exclusive(path::String;
                        mode::Integer = 0o444 #= read-only =#,
                        poll_interval::Real = 10 #= seconds =#,
                        wait::Bool = true #= return on failure if false =#,
                        stale_age::Real = 0 #= disabled =#)
    # fast-path: just try to open it
    file = tryopen_exclusive(path, mode)
    file === nothing || return file
    if !wait
        if file === nothing && stale_age > 0
            if stale_age > 0 && stale_pidfile(path, stale_age)
                @warn "attempting to remove probably stale pidfile" path=path
                tryrmopenfile(path)
            end
            file = tryopen_exclusive(path, mode)
        end
        if file === nothing
            error("Failed to get pidfile lock for $(repr(path)).")
        else
            return file
        end
    end
    # fall-back: wait for the lock

    while true
        # start the file-watcher prior to checking for the pidfile existence
        t = @async try
            watch_file(path, poll_interval)
        catch ex
            isa(ex, IOError) || rethrow(ex)
            sleep(poll_interval) # if the watch failed, convert to just doing a sleep
        end
        # now try again to create it
        file = tryopen_exclusive(path, mode)
        file === nothing || return file
        Base.wait(t) # sleep for a bit before trying again
        if stale_age > 0 && stale_pidfile(path, stale_age)
            # if the file seems stale, try to remove it before attempting again
            # set stale_age to zero so we won't attempt again, even if the attempt fails
            stale_age -= stale_age
            @warn "attempting to remove probably stale pidfile" path=path
            tryrmopenfile(path)
        end
    end
end

function _rand_filename(len::Int=4) # modified from Base.Libc
    slug = Base.StringVector(len)
    chars = b"0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ"
    for i = 1:len
        slug[i] = chars[(Libc.rand() % length(chars)) + 1]
    end
    return String(slug)
end

function tryrmopenfile(path::String)
    # Deleting open file on Windows is a bit hard
    # if we want to reuse the name immediately after:
    # we need to first rename it, then delete it.
    if Sys.iswindows()
        try
            local rmpath
            rmdir, rmname = splitdir(path)
            while true
                rmpath = string(rmdir, isempty(rmdir) ? "" : path_separator,
                    "\$", _rand_filename(), rmname, ".deleted")
                ispath(rmpath) || break
            end
            rename(path, rmpath)
            path = rmpath
        catch ex
            isa(ex, IOError) || rethrow(ex)
        end
    end
    return try
        rm(path)
        true
    catch ex
        isa(ex, IOError) || rethrow(ex)
        ex
    end
end

"""
    close(lock::LockMonitor)

Release a pidfile lock.
"""
function Base.close(lock::LockMonitor)
    update = lock.update
    update === nothing || close(update)
    isopen(lock.fd) || return false
    removed = false
    path = lock.path
    pathstat = try
            # Windows sometimes likes to return EACCES here,
            # if the path is in the process of being deleted
            stat(path)
        catch ex
            ex isa IOError || rethrow()
            removed = ex
            nothing
        end
    if pathstat !== nothing && samefile(stat(lock.fd), pathstat)
        # try not to delete someone else's lock
        removed = tryrmopenfile(path)
    end
    close(lock.fd)
    havelock = removed === true
    havelock || @warn "failed to remove pidfile on close" path=path removed=removed
    return havelock
end

end # module

max_stars_repo_name: SmiteWindows/julia

max_issues_repo_name: SmiteWindows/julia

max_forks_repo_name: SmiteWindows/julia


text:
abstract type los_direction end

struct los_radial <: los_direction end
struct los_earth <: los_direction end

Base.broadcastable(l::los_direction) = Ref(l)

#######################################################################################################

# covariant helicity components of eₓ
function line_of_sight_covariant((θ,ϕ)::Tuple{Real,Real},
	los::los_earth = los_earth(), eulerangles::Nothing = nothing)

	m1 = +1/√2 * (cos(θ)cos(ϕ) + im*sin(ϕ))
	z  = sin(θ)cos(ϕ)
	p1 = -1/√2 * (cos(θ)cos(ϕ) - im*sin(ϕ))
	OffsetVector(SVector{3}(m1, z, p1), -1:1)
end

function rotatedpt((θ,ϕ), eulerangles::NTuple{3,Real})
	α, β, γ = eulerangles
	invrot = inv(RotZYZ(α, β, γ)) # active equivalent
	v = invrot * radiusvector(θ, ϕ)
	p = Point2D(v)
	p.θ, p.ϕ
end

function line_of_sight_covariant((θ,ϕ)::NTuple{2,Real},
	los::los_earth, eulerangles::NTuple{3,Real})

	θ′,ϕ′ = rotatedpt((θ,ϕ), eulerangles)

	p1 = -1/√2 * (cos(θ′)cos(ϕ′) - im*sin(ϕ′))
	m1 = +1/√2 * (cos(θ′)cos(ϕ′) + im*sin(ϕ′))
	z  = sin(θ′)cos(ϕ′)
	OffsetVector(SVector{3}(m1, z, p1),-1:1)
end

function line_of_sight_covariant(n::SphericalPoint,
	los::los_earth = los_earth(), eulerangles = nothing)

	line_of_sight_covariant((n.θ,n.ϕ), los, eulerangles)
end

function line_of_sight_covariant(n, ::los_radial, args...)
	OffsetArray(SVector{3,Float64}(0,1,0), -1:1)
end

function polcoords(v)
    # This is not unique at the poles, as (θ=0, ϕ) map to the north pole for all ϕ
	# We choose to retutn ϕ = 0 at the poles
	vunit = normalize(v)
	if vunit[3] == 1
		# north pole
		return SVector{2}(promote(oftype(float(vunit[3]), 0), 0))
	elseif vunit[3] == -1
		# south pole
		return SVector{2}(promote(oftype(float(vunit[3]), pi), 0))
	end
	θ = acos(vunit[3])
	ϕ = mod2pi(atan(vunit[2], vunit[1]))
	SVector{2}(promote(θ, ϕ))
end
rotn1n1′(θ1, ϕ1, θ1′, ϕ1′) = RotZYZ(ϕ1′, θ1′ - θ1, -ϕ1)
isapproxdefault(x, y) = isapprox(x, y, atol = 1e-14, rtol = 1e-8)
function rotn1′n2n2′(r1′, r2′′, r2′)
	r1′x, r1′y, r1′z = r1′
	a = r2′′ - r1′
	b = r2′ - r1′
	# Avoid errors due to machine precision
	if (isapproxdefault(cross(r1′, a), zero(SVector{3})) || isapproxdefault(cross(r1′, b), zero(SVector{3})))
		return RotZYZ(one(SMatrix{3,3,eltype(r1′),9}))
	end
	sinω = b ⋅ cross(r1′, a) / norm(cross(r1′, a))^2
	cosω = b ⋅ (a - (r1′⋅a)r1′) / norm(a - (r1′⋅a)r1′)^2
	r1′C = SMatrix{3,3}(0, r1′z, -r1′y, -r1′z, 0, r1′x, r1′y, -r1′x, 0)
	r1′r1′ = r1′ * r1′'
	RotZYZ(cosω*I + sinω * r1′C + (1 - cosω)*r1′r1′)
end

function _rotation_points(::Cartesian, n1, n2, n1′, n2′)
	(θ1, ϕ1), (θ1′, ϕ1′) = map(Point2D, (n1, n1′))
	r1, r2, r1′, r2′ = map(radiusvector∘Point2D	, (n1, n2, n1′, n2′))
	if !isapproxdefault(r1 ⋅ r2, r1′ ⋅ r2′)
		throw(ArgumentError("points are not related by a rotation"))
	end
	R11′ = rotn1n1′(θ1, ϕ1, θ1′, ϕ1′)
	r2′′ = R11′ * r2
	R2′′2′ = rotn1′n2n2′(r1′, r2′′, r2′)
	R = RotZYZ(R2′′2′ * R11′)
	R * r1 ≈ r1′ || error("rotation does not map r1 to r1′")
	R * r2 ≈ r2′ || error("rotation does not map r2 to r2′")
	return R
end

function rotation_points_12(Basis, n1, n2, n1′, n2′)
	R = _rotation_points(Cartesian(), Point2D(n1), Point2D(n2), Point2D(n1′), Point2D(n2′))
	Un1, Un1′, Un2, Un2′ = map((n1, n1′, n2, n2′)) do n
		VectorSphericalHarmonics.basisconversionmatrix(Cartesian(), Basis, Point2D(n)...)
	end
	M′1 = Un1′ * R * Un1'
	M′2 = Un2′ * R * Un2'
	return R, M′1, M′2
end
function _rotation_points(Basis, n1, n2, n1′, n2′)
	_, M′1, M′2 = rotation_points_12(Basis, n1, n2, n1′, n2′)
	LinearAlgebra.kron(M′1, M′2)
end
rotation_points(Basis, n1, n2, n1′, n2′) = _rotation_points(Basis, n1, n2, n1′, n2′)
rotation_points(Basis::HelicityCovariant, n1, n2, n1′, n2′) = Diagonal(_rotation_points(Basis, n1, n2, n1′, n2′))


text:
exec '/Applications/Julia-0.6.app/Contents/Resources/julia/bin/julia'
alias julia=/Applications/Julia-0.6.app/Contents/Resources/julia/bin/julia

Pkg.update()
Pkg.add("Watershed")
Pkg.checkout("Watershed")

### test
using Watershed
aff = rand(Float32, 124,124,12,3);
# watershed(aff)
#println("watershed ...")
#@time watershed(aff)
# @profile watershed(aff)
# Profile.print()
low = 0.1
high = 0.8
thresholds = []
dust_size = 1
watershed(aff; is_threshold_relative=true)

### test assets
chmod +x ~/.julia/v0.6/Watershed/test/test_assets.jl
julia ~/.julia/v0.6/Watershed/test/test_assets.jl

Pkg.add("EMIRT")
Pkg.add("HDF5")

using Watershed
using EMIRT
using HDF5
using Base.Test

aff = h5read(joinpath("/Users/michielk/.julia/v0.6/Watershed", "assets/piriform.aff.h5"), "main")

#seg = atomicseg(aff)
seg, rg = watershed(aff; is_threshold_relative = true)
#h5write(joinpath(Pkg.dir(), "Watershed/assets/piriform.seg.h5", "seg", seg)
#h5write("seg.h5", "seg", seg)
# compare with segmentation

seg0 = readseg(joinpath(("/Users/michielk/.julia/v0.6/Watershed", "assets/piriform.seg.h5"))

err = segerror(seg0, seg)
@show err
@test_approx_eq err[:re]  0

### test MKdata

using Watershed
using EMIRT
using HDF5
using Base.Test

"""
# in conda root env
import os
import h5py
import numpy as np
from wmem import utils
datadir = "/Users/michielk/oxdata/P01/EM/Myrf_01/SET-B/B-NT-S10-2f_ROI_00/blocks_0500";
pred_file = "B-NT-S10-2f_ROI_00_00480-01020_00480-01020_00000-00184_probs_vol01+vol06";
h5path_in = os.path.join(datadir, pred_file + '.h5', 'data')
h5file_in, ds_in, elsize, axlab = utils.h5_load(h5path_in)
grad = np.array(np.absolute(np.gradient(ds_in[:], 1)))
h5file_in.close()
h5path_out = os.path.join(datadir, pred_file + '_absgrad.h5', 'main')
h5file_out, ds_out = utils.h5_write(None, grad.shape, grad.dtype,
                                    h5path_out,
                                    element_size_um=np.insert(elsize, 0, 1),
                                    axislabels=np.insert(axlab, 0, 'c'))
ds_out[:] = grad
h5file_out.close()
"""

datadir = "/Users/michielk/oxdata/P01/EM/Myrf_01/SET-B/B-NT-S10-2f_ROI_00/blocks_0500";
pred_file = "B-NT-S10-2f_ROI_00_00480-01020_00480-01020_00000-00184_probs_vol01+vol06";
aff = h5read(joinpath(datadir, pred_file * "_absgrad.h5"), "main");

low = 0.1;
high = 0.8;
thresholds = Tuple{Int64,Float32}[(8000, 0.0)];
dust_size = 600;
seg, rg = watershed(aff; low = low, high = high, thresholds = thresholds, dust_size = dust_size, is_threshold_relative = true);
h5write(joinpath(datadir, pred_file * "_absgrad_seg.h5"), "seg", seg);

"""
datadir="/Users/michielk/oxdata/P01/EM/Myrf_01/SET-B/B-NT-S10-2f_ROI_00"
datastem="B-NT-S10-2f_ROI_00_00480-01020_00480-01020_00000-00184"
ipf=; ids=data; ods=$ids; args='-i zyx -o xyz';
ipf=_probs_vol01+vol06; ids=data; ods=$ids; args='-i zyx -o xyz';
ipf=_probs_vol01+vol06_absgrad; ids=main; ods=${ids}_vol00; args='-i czyx -o xyzc -D 0 1 1 0 0 1 0 0 1 0 0 1';
pf=_probs_vol01+vol06_absgrad_seg; ids=seg; ods=$ids; args='-i zyx -o xyz -e 0.007 0.007 0.1';
python $scriptdir/wmem/stack2stack.py \
    $datadir/blocks_0500/$datastem$ipf.h5/$ids \
    $datadir/blocks_0500/$datastem${ipf}_$ods.nii.gz $args
"""

datadir = "/Users/michielk/oxdata/P01/EM/Myrf_01/SET-B/B-NT-S10-2f_ROI_00/zws";
aff = h5read(joinpath(datadir, "B-NT-S10-2f_ROI_00ds7_probs_main_vol00_grad.h5"), "main");

# watershed(aff::Array{Float32,4}; low, high, thresholds, dust_size, is_threshold_relative)
# thresholds=[(merge_size, merge_threshold)]

seg = atomicseg(aff);
h5write(joinpath(datadir, "B-NT-S10-2f_ROI_00ds7_probs_main_vol00_grad.h5"), "atomseg", seg);
seg, rg = watershed(aff; is_threshold_relative = true);
# low 0.1, high 0.8, thresholds Tuple{Int64,Float32}[(800, 0.2)]
# INFO: absolute watershed threshold: low: -0.100817, high: 0.019489, thresholds: Tuple{Int64,Float32}[(800, -0.020919)], dust: 600
# steepestascent...
# divideplateaus...
# findbasins!
# found: 21239023 components
# regiongraph...
# Region graph size: 98651500
# use region graph construction code > julia 0.4
# mergeregions...
# Done with merging
# Done with remapping, total: 5931 regions
# Done with updating the region graph, size: 41675
h5write(joinpath(datadir, "B-NT-S10-2f_ROI_00ds7_probs_main_vol00_grad.h5"), "seg", seg);

seg, rg = watershed(aff; low = 0.0, high = 1.0, is_threshold_relative = true);
# INFO: use percentage threshold: low 0.0, high 1.0, thresholds Tuple{Int64,Float32}[(800, -0.020919)]
# INFO: absolute watershed threshold: low: -0.49285802, high: 0.489795, thresholds: Tuple{Int64,Float32}[(800, -0.492858)], dust: 600
# steepestascent...
# divideplateaus...
# findbasins!
# found: 41816208 components
# regiongraph...
# Region graph size: 218970938
# use region graph construction code > julia 0.4
# mergeregions...
# Done with merging
# Done with remapping, total: 27812 regions
# Done with updating the region graph, size: 307919
h5write(joinpath(datadir, "B-NT-S10-2f_ROI_00ds7_probs_main_vol00_grad.h5"), "segall", seg);

seg, rg = watershed(aff; low = 0.0, high = 1.0, dust_size = 10, is_threshold_relative = true);
# INFO: use percentage threshold: low 0.0, high 1.0, thresholds Tuple{Int64,Float32}[(800, -0.492858)]
# INFO: absolute watershed threshold: low: -0.49285802, high: 0.489795, thresholds: Tuple{Int64,Float32}[(800, -0.492858)], dust: 10
# steepestascent...
# divideplateaus...
# findbasins!
# found: 41816208 components
# regiongraph...
# Region graph size: 218970938
# use region graph construction code > julia 0.4
# mergeregions...
# Done with merging
# Done with remapping, total: 27812 regions
# Done with updating the region graph, size: 307919
h5write(joinpath(datadir, "B-NT-S10-2f_ROI_00ds7_probs_main_vol00_grad.h5"), "segdust", seg);

seg, rg = watershed(aff; low = 0.0, high = 1.0, thresholds = Tuple{Int64,Float32}[(800, 0.2)], dust_size = 10, is_threshold_relative = true);
# INFO: use percentage threshold: low 0.0, high 1.0, thresholds Tuple{Int64,Float32}[(800, 0.2)]
# INFO: absolute watershed threshold: low: -0.49285802, high: 0.489795, thresholds: Tuple{Int64,Float32}[(800, -0.020919)], dust: 10
# steepestascent...
# divideplateaus...
# findbasins!
# found: 41816208 components
# regiongraph...
# Region graph size: 218970938
# use region graph construction code > julia 0.4
# mergeregions...
# Done with merging
# Done with remapping, total: 27812 regions
# Done with updating the region graph, size: 307919
h5write(joinpath(datadir, "B-NT-S10-2f_ROI_00ds7_probs_main_vol00_grad.h5"), "segmthr", seg);

low = 0.0;
high = 1.0;
thresholds = Tuple{Int64,Float32}[(100, 0.2)];
dust_size = 10;
seg, rg = watershed(aff; low = low, high = high, thresholds = thresholds, dust_size = dust_size, is_threshold_relative = true);
# INFO: use percentage threshold: low 0.0, high 1.0, thresholds Tuple{Int64,Float32}[(100, 0.2)]
# INFO: absolute watershed threshold: low: -0.49285802, high: 0.489795, thresholds: Tuple{Int64,Float32}[(100, -0.020919)], dust: 10
# steepestascent...
# divideplateaus...
# findbasins!
# found: 41816208 components
# regiongraph...
# Region graph size: 218970938
# use region graph construction code > julia 0.4
# mergeregions...
# Done with merging
# Done with remapping, total: 389363 regions
# Done with updating the region graph, size: 4798320
h5write(joinpath(datadir, "B-NT-S10-2f_ROI_00ds7_probs_main_vol00_grad.h5"), "segmsize", seg);

# new_rg = mergeregions!(seg, rg, counts, thresholds, dust_size)

"""
pf=_probs_main_vol00; ids=main; ods=$ids; args='-i czyx -o xyzc';
pf=_probs_main_vol00_grad; ids=main; ods=${ids}_vol00; args='-i czyx -o xyzc -D 0 1 1 0 0 1 0 0 1 0 0 1';
pf=_probs_main_vol00_grad; ids=atomseg; ods=$ids; args='-i zyx -o xyz -e 0.049 0.049 0.1';
pf=_probs_main_vol00_grad; ids=seg; ods=$ids; args='-i zyx -o xyz -e 0.049 0.049 0.1';
pf=_probs_main_vol00_grad; ids=segall; ods=$ids; args='-i zyx -o xyz -e 0.049 0.049 0.1';
pf=_probs_main_vol00_grad; ids=segdust; ods=$ids; args='-i zyx -o xyz -e 0.049 0.049 0.1';
pf=_probs_main_vol00_grad; ids=segmthr; ods=$ids; args='-i zyx -o xyz -e 0.049 0.049 0.1';
pf=_probs_main_vol00_grad; ids=segmsize; ods=$ids; args='-i zyx -o xyz -e 0.049 0.049 0.1';
python $scriptdir/wmem/stack2stack.py \
    $datadir/zws/$dataset_ds$pf.h5/$ids \
    $datadir/zws/$dataset_ds${pf}_$ods.nii.gz $args
"""

max_stars_repo_name: michielkleinnijenhuis/EM

max_issues_repo_name: michielkleinnijenhuis/EM

max_forks_repo_name: michielkleinnijenhuis/EM


text:
using Test

using Random
Random.seed!(42)

ENV["JULIA_DEBUG"] = Main

function run_all_tests()

    @testset "Basics" begin
        include("basic-test.jl")
    end

    @testset "animation (Unix only)" begin
        if ! Sys.iswindows()
            include("animation-test.jl")
        end
        include("easingfunctions.jl")
    end

    @testset "polygons" begin
        include("boxtests.jl")
        include("offset-poly-tests.jl")
        include("path-to-poly.jl")
        include("point-arithmetic.jl")
        include("point-inside-polygon.jl")
        include("point-intersection.jl")
        include("polycross-test.jl")
        include("polyfit-test.jl")
        include("polysample.jl")
        include("polygon-centroid-sort-test.jl")
        include("polygon-test.jl")
        include("polysplit.jl")
        include("polysplit-2.jl")
        include("polysmooth-tests.jl")
        include("pretty-poly-test.jl")
        include("simplify-polygons.jl")
        include("spirals.jl")
        include("star-test-1.jl")
        include("star-test.jl")
        include("test-holes.jl")
        include("morepolytests.jl")
        include("polytriangulate.jl")
        include("pointcircletangent-test.jl")
    end

    @testset "text" begin
        include("text-alignment.jl")
        include("text-rotations.jl")
        include("pro-text-test.jl")
        include("text-curve-centered.jl")
        include("text-path-clipping.jl")
        include("text-wrapping.jl")
        include("text-boxes.jl")
        include("label-tests.jl")
    end

    @testset "curves" begin
        include("arc-twopoints.jl")
        include("circletests.jl")
        include("circleintersectiontests.jl")
        include("circle-as-path.jl")
        include("randomsinecurves.jl")
        include("sector-test.jl")
        include("sector-rounded.jl")
        include("ellipse-test.jl")
        include("hypertrochoid-test.jl")
        include("epitrochoid-test.jl")
        include("pie-test.jl")
        include("bezierpath.jl")
        include("bezierpathtopoly.jl")
        include("bezierstroke-test.jl")
        include("circletangent-test.jl")
        include("arc-sagitta-test.jl")
        include("circlecircletangent-test.jl")
        include("pointcircletangent-test.jl")
    end

    @testset "color" begin
        include("color-blend-test.jl")
        include("blendmodes.jl")
        include("palette_test.jl")
        include("meshtest.jl")
    end

    @testset "images" begin
        include("images-with-alpha.jl")
        include("test-image.jl")
#        include("tiling-images.jl")
    end

    @testset "imagematrix" begin
        include("imagematrix.jl")
    end

    @testset "arrows and dimensions" begin
        include("arrow-arc-test.jl")
        include("arrow-bezier-test.jl")
        include("arrow-line-test.jl")
        include("dimensioning.jl")
    end

    @testset "getpath" begin
        include("get_path.jl")
        include("get_path_flat.jl")
    end

    @testset "julia logos" begin
        include("heart-julia.jl")
        include("julia-logo-draw-eps.jl")
        include("julia-logo-draw.jl")
        include("julia-logo-testing.jl")
    end

    @testset "rulers" begin
        include("rulers-test.jl")
    end

    @testset "boundingbox" begin
        include("boundingboxtests.jl")
        include("boundingboxpointtests.jl")
    end

    @testset "clipping" begin
        include("clipping-test.jl")
    end
    @testset "tilerstables" begin
        include("pagetiler-test.jl")
        include("table-tests.jl")
    end
    @testset "matrix" begin
        include("matrix-tests.jl")
    end
    @testset "sierpinski" begin
        include("sierpinski-svg.jl")
        include("sierpinski.jl")
    end
    @testset "turtle" begin
        include("turtle-test.jl")
    end
    @testset "intersection" begin
        include("line-intersection-options.jl")
        include("intersection-line-circle.jl")
        include("line-intersections.jl")
        include("polygonintersect-test.jl")
    end
    @testset "misc" begin
        include("luxor-test1.jl")
        include("rules.jl")
        include("barstest.jl")
        include("unit-conversions.jl")
        include("hex-grid-test.jl")
        include("cropmarkstest.jl")
        include("boxmaptest.jl")
        include("noise-test.jl")
        include("dashtests.jl")
        include("triangles.jl")
    end
end

if get(ENV, "LUXOR_KEEP_TEST_RESULTS", false) == "true"
        # they changed mktempdir in v1.3
        if VERSION <= v"1.2"
            cd(mktempdir())
        else
            cd(mktempdir(cleanup=false))
        end
        @info("...Keeping the results in: $(pwd())")
        run_all_tests()
        @info("Test images were saved in: $(pwd())")
else
    mktempdir() do tmpdir
        cd(tmpdir) do
            @info("running tests in: $(pwd())")
            @info("but not keeping the results")
            @info("because you didn't do: ENV[\"LUXOR_KEEP_TEST_RESULTS\"] = \"true\"")
            run_all_tests()
            @info("Test images weren't saved. To see the test images, next time do this before running:")
            @info(" ENV[\"LUXOR_KEEP_TEST_RESULTS\"] = \"true\"")
        end
    end
end


text:
""" List appender: holds events, messages & data in memory

    This appender is primarily used for testing. Use in a real environment
    is discouraged as the it could eventually grow out of memory.
"""
type List <: Appender
    name::AbstractString
    layout::LAYOUT
    filter::FILTER
    state::LifeCycle.State

    events::Vector{Event}
    messages::Vector{Message}
    data::Vector{UInt8}

    raw::Bool
    newLine::Bool

    List(name::AbstractString) = new(name, LAYOUT(), FILTER(), LifeCycle.INITIALIZED, Event[], Message[], UInt8[], false, false)
    function List(name::AbstractString, lyt::LAYOUT, flt::FILTER; raw=false, newline=false)
        apndr = new(name, lyt, flt, LifeCycle.INITIALIZED, Event[], Message[], UInt8[], raw, newline)
        if !isnull(apndr.layout)
            hdr = header(layout)
            if length(hdr) > 0
                write(apndr, hdr)
            end
        end
        apndr
    end
end
function List(config::Dict{AbstractString,Any})
    nm  = get(config, "name", "List")
    lyt = get(config, "layout", nothing)
    flt = get(config, "filter", nothing)
    raw = get(config, "raw", false)
    nl  = get(config, "newline", false)
    List(nm, LAYOUT(lyt), FILTER(flt), raw=raw, newline=nl)
end
List(;kwargs...) = map(e->(string(e[1]),e[2]), kwargs) |> Dict{AbstractString,Any} |> List

show(io::IO, apnd::List) = print(io, "List["* (isnull(apnd.layout) ? "evnts=$(length(apnd.events))" : "msgs=$(length(apnd.messages))")*", state=$(apnd.state))")
name(apnd::List) = isempty(apnd.name) ? string(typeof(apnd)) : apnd.name

function append!(apnd::List, evnt::Event)
    if isnull(apnd.layout)
        push!(apnd.events, evnt)
    else
        write(serialize(apnd.layout, evnt))
    end
end

function write(apnd::List, data::Vector{UInt8})
    if apnd.raw
        push!(apnd.data, data)
    else
        msg = bytestring(data)
        if apnd.newline
            for part in split(msg,['\n','\r'],keep=false)
                push!(apnd.messages, part)
            end
        else
            push!(apnd.messages, msg)
        end
    end
end

function empty!(apnd::List)
    empty!(apnd.events)
    empty!(apnd.messages)
    empty!(apnd.data)
end


text:
using Phylo
using Phylo.API
using Printf

function show(io::IO, object::Tuple{<: AbstractTree, <: AbstractNode},
              n::String = "")
    node = "node"
    if !isempty(n)
        node *= " $n"
    end
    if !hasinbound(object[1], object[2])
        if outdegree(object[1], object[2]) > 0
            blank = repeat(" ", length(" [root $node]") + (isempty(n) ? 0 : 1))
            for (i, bn) in zip(Base.OneTo(_outdegree(object[1], object[2])),
                               getoutbounds(object[1], object[2]))
                b = typeof(bn) <: Number ? "$bn" : "\"$bn\""
                if outdegree(object[1], object[2]) == 1
                    print(io, "[root $node]-->[branch $b]")
                elseif get(io, :compact, false)
                    if i == 1
                        print(io, "[root $node]-->[branches $b")
                    elseif i < outdegree(object[1], object[2])
                        print(io, ", $b")
                    else
                        print(io, " and $b]")
                    end
                else # multiline view
                    if i == 1
                        println(io, "[root $node]-->[branch $b]")
                    elseif i < outdegree(object[1], object[2])
                        println(io, "$blank-->[branch $b]")
                    else
                        print(io, "$blank-->[branch $b]")
                    end
                end
            end
        else
            print(io, "[unattached $node]")
        end
    else # hasinbound
        inb = typeof(getinbound(object[1], object[2])) <: Number ?
            "$(_getinbound(object[1], object[2]))" :
            "\"$(_getinbound(object[1], object[2]))\""
        if outdegree(object[1], object[2]) == 0
            print(io, "[branch $inb]-->[leaf $node]")
        elseif hasinbound(object[1], object[2])
            blank = repeat(" ",
                           length(" [branch $inb]-->[internal $node]") +
                           (isempty(n) ? 0 : 1))
            for (i, bn) in zip(Base.OneTo(_outdegree(object[1], object[2])),
                               getoutbounds(object[1], object[2]))
                b = typeof(bn) <: Number ? "$bn" : "\"$bn\""
                if outdegree(object[1], object[2]) == 1
                    print(io, "[branch $inb]-->[internal $node]-->[branch $b]")
                elseif get(io, :compact, false)
                    if i == 1
                        print(io, "[branch $inb]-->[internal $node]-->" *
                              "[branches $b")
                    elseif i < outdegree(object[1], object[2])
                        print(io, ", $b")
                    else
                        print(io, " and $b]")
                    end
                else # multiline view
                    if i == 1
                        println(io, "[branch $inb]-->[internal $node]-->" *
                                "[branch $b]")
                    elseif i < outdegree(object[1], object[2])
                        println(io, "$blank-->[branch $b]")
                    else
                        print(io, "$blank-->[branch $b]")
                    end
                end
            end
        end
    end
end

function show(io::IO, p::Pair{NT, N}) where {N <: AbstractNode, NT}
    n = NT <: Number ? "$(p[1])" : "\"$(p[1])\""
    show(io, p[2], "$n")
end

function show(io::IO, object::Tuple{<: AbstractTree{OneTree, RT, NL},
                                    <: AbstractBranch{RT, NL}}) where {RT, NL}
    source = NL <: Number ? "$(_src(object[1], object[2]))" :
        "\"$(_src(object[1], object[2]))\""
    destination = NL <: Number ? "$(_dst(object[1], object[2]))" :
        "\"$(_dst(object[1], object[2]))\""
    println(io, "[node $source]-->[$(_getlength(object[1], object[2])) " *
        "length branch]--> [node $destination]")
end

function show(io::IO, p::Pair{BT, B}) where {BT, B <: AbstractBranch}
    NT = typeof(_src(p[2]))
    source = NT <: Number ? "$(_src(p[2]))" : "\"$(_src(p[2]))\""
    destination = NT <: Number ? "$(_dst(p[2]))" : "\"$(_dst(p[2]))\""
    branch = BT <: Number ? "$(p[1])" : "\"$(p[1])\""
    println(io, "[node $source]-->[$(_getlength(p[2])) length branch $branch]-->[node $destination]")
end

function showsimple(io::IO, object::TreeSet)
    @printf(io, "TreeSet with %d trees, each with %d tips.\n",
            ntrees(object), nleaves(object))
    tn = collect(gettreenames(object))
    if length(tn) < 10
        println(io, "Tree names are " * join(tn, ", ", " and "))
    else
        println(io, "Tree names are " * join(tn[1:5], ", ") *
                " ... $(length(tn) - 6) missing ... $(tn[end])")
    end
end

function show(io::IO, object::TreeSet)
    showsimple(io, object)
    tn =  sort(collect(_gettreenames(object)))
    index = 0
    for name in tn
        index += 1
        if index ≤ 5 || index == length(tn)
            @printf(io, "\n%s: ", name)
            showsimple(io, object[name])
        elseif index == 6
            println("\n[$(length(tn)-6) trees omitted]\n")
        end
    end
end

function showsimple(io::IO, object::TREE) where {TREE <: AbstractTree}
    println(io, "$TREE with $(_nleaves(object)) tips, " *
            "$(_nnodes(object)) nodes and " *
            "$(_nbranches(object)) branches.")
    ln = getleafnames(object)
    if length(ln) < 10
        println(io, "Leaf names are " * join(ln, ", ", " and "))
    else
        println(io, "Leaf names are " * join(ln[1:5], ", ") *
                ", ... [$(length(ln) - 6) omitted] ... and $(ln[end])")
    end
end

function show(io::IO, object::AbstractTree)
    showsimple(io, object)
    if !get(io, :compact, true)
        println(io, "Nodes:")
        println(io, [(object, node) for node in _getnodes(object)])
        println(io, "Branches:")
        println(io, [(object, branch) for branch in _getbranches(object)])
        if _nodedatatype(TREE) !== Nothing
            println(io, "Node records:")
            println(io, Dict(name => _getnodedata(object, name)
                             for name in _getnodenames(object)))
        end
    end
end


text:
#Bitterness is measured in  International Bitterness Units
type Bitter{T<:Real}
    IBU :: T
end

isbitter(A) = false
isbitter(A::Bitter) = A.IBU > 0


text:
const ARGS_BLOCKS = Dict{Symbol,BlockTemplate}(e.opening => e for e in [
   BlockTemplate(:STRING,  :SINGLE_QUOTE, :SINGLE_QUOTE ),
   BlockTemplate(:STRING,  :TRIPLE_QUOTE, :TRIPLE_QUOTE ),
   ])


text:
module MatDeforNeohookeanADModule

using FinEtools.FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import FinEtoolsDeforLinear.DeforModelRedModule: AbstractDeforModelRed, DeforModelRed3D, DeforModelRed2DStrain, DeforModelRed2DStress, DeforModelRed2DAxisymm, DeforModelRed1D, nstressstrain, nthermstrain
import FinEtoolsDeforLinear.MatDeforModule: AbstractMatDefor, stressvtot!, stressttov!, strainttov!, strainvdet, strainvtr
import ..MatDeforNonlinearModule: AbstractMatDeforNonlinear, totlag2currsymm!
using LinearAlgebra: Transpose, Diagonal, mul!
At_mul_B!(C, A, B) = mul!(C, Transpose(A), B)
A_mul_B!(C, A, B) = mul!(C, A, B)
using LinearAlgebra: eigen, eigvals, norm, cholesky, cross, dot, log, diagm, det
using ForwardDiff: gradient, hessian

"""
	MatDeforNeohookeanAD{MR<:AbstractDeforModelRed, MTAN<:Function, MUPD<:Function} <: AbstractMatDeforNonlinear

Type for triaxial neohookean hyperelastic material.

"""
struct  MatDeforNeohookeanAD{MR<:AbstractDeforModelRed, MTAN<:Function, MUPD<:Function} <: AbstractMatDeforNonlinear
	mr::Type{MR} # model reduction type
	mass_density::FFlt # mass density
	E::FFlt # Young's modulus
	nu::FFlt # Poisson ratio
	_lambda::FFlt
	_mu::FFlt
	tangentmoduli!::MTAN # Function to return the tangent moduli matrix
	update!::MUPD # Function to update the material state
end

"""
    MatDeforNeohookeanAD(mr::Type{MR}, E::FFlt, nu::FFlt) where {MR}

Create neohookean isotropic elastic material.

The mass density is the default value.
"""
function MatDeforNeohookeanAD(mr::Type{MR}, E::FFlt, nu::FFlt) where {MR}
	mass_density = 1.0
	return MatDeforNeohookeanAD(mr, mass_density, E, nu)
end

const D3D = DeforModelRed3D

# Note that the strain energy is defined in terms of the right Cauchy-Green tensor.
function strainenergy(Cv, lambda, mu)
	trC = strainvtr(D3D, Cv)
	J = sqrt(strainvdet(D3D, Cv))
	lJ = log(J)
	return mu/2*(trC-3) - mu*lJ + lambda/2*(lJ)^2;
end

################################################################################
# 3-D solid model
################################################################################

"""
    MatDeforNeohookeanAD(mr::Type{DeforModelRed3D}, mass_density::FFlt, E::FFlt, nu::FFlt)

Create triaxial neohookean hyperelastic material.
"""
function MatDeforNeohookeanAD(mr::Type{DeforModelRed3D}, mass_density::FFlt, E::FFlt, nu::FFlt)
	function tangentmoduli3d!(self::MatDeforNeohookeanAD, D::FFltMat, statev::FFltVec, Fn1::FFltMat, Fn::FFltMat, tn::FFlt, dtn::FFlt, loc::FFltMat, label::FInt)
		C = Fn1'*Fn1;
		Cv = fill(0.0, 6)
		strainttov!(mr, Cv, C)
		Dtotal = 4 .* hessian(Cv -> strainenergy(Cv, self._lambda, self._mu), Cv);
		return totlag2currsymm!(D, Dtotal, Fn1)
	end
	function update3d!(self::MatDeforNeohookeanAD, statev::FFltVec, cauchy::FFltVec, output::FFltVec, Fn1::FFltMat, Fn::FFltMat, tn::FFlt, dtn::FFlt, loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)
		@assert length(cauchy) == nstressstrain(self.mr)
		C = Fn1'*Fn1;
		Cv = fill(0.0, 6)
		strainttov!(mr, Cv, C)
		Sv = 2 * gradient(Cv -> strainenergy(Cv, self._lambda, self._mu), Cv);
		S = fill(0.0, 3, 3)
		stressvtot!(mr, S, Sv)
		cauchyt = Fn1*(S/det(Fn1))*Fn1'; # Cauchy stress
		stressttov!(mr, cauchy, cauchyt)
		if quantity == :nothing
			#Nothing to be copied to the output array
		elseif quantity == :cauchy || quantity == :Cauchy
			(length(output) >= 6) || (output = zeros(6)) # make sure we can store it
			copyto!(output, cauchy);
		elseif quantity == :pressure || quantity == :Pressure
			output[1]  =  -sum(cauchy[1:3])/3.
		elseif quantity == :princCauchy || quantity == :princcauchy
			ep = eigen(sigma);
			(length(output) >= 3) || (output = zeros(3)) # make sure we can store it
			copyto!(output,  sort(ep.values, rev=true));
		elseif quantity==:vonMises || quantity==:vonmises || quantity==:von_mises || quantity==:vm
			s1=cauchy[1]; s2=cauchy[2]; s3=cauchy[3];
			s4=cauchy[4]; s5=cauchy[5]; s6=cauchy[6];
			(length(output) >= 1) || (output = zeros(1)) # make sure we can store it
			output[1] = sqrt(1.0/2*((s1-s2)^2+(s1-s3)^2+(s2-s3)^2+6*(s4^2+s5^2+s6^2)))
		end
		return output
	end
	_lambda = E * nu / (1 + nu) / (1 - 2*(nu));
	_mu = E / 2. / (1+nu);
	return MatDeforNeohookeanAD(mr, mass_density, E, nu, _lambda, _mu,
		tangentmoduli3d!, update3d!)
end

################################################################################
# 2-D plane strain model
################################################################################

"""
    MatDeforNeohookeanAD(mr::Type{DeforModelRed2DStrain}, mass_density::FFlt, E::FFlt, nu::FFlt)

Create triaxial neohookean hyperelastic material.
"""
function MatDeforNeohookeanAD(mr::Type{DeforModelRed2DStrain}, mass_density::FFlt, E::FFlt, nu::FFlt)
	function tangentmoduli2dstrain!(self::MatDeforNeohookeanAD, D::FFltMat, statev::FFltVec, Fn1::FFltMat, Fn::FFltMat, tn::FFlt, dtn::FFlt, loc::FFltMat, label::FInt)
		Fn13d = fill(0.0, 3, 3)
		Fn13d[1:2, 1:2] = Fn1
		Fn13d[3, 3] = 1.0
		C = Fn13d'*Fn13d;
		Cv = fill(0.0, 6)
		strainttov!(D3D, Cv, C)
		Dtotal = 4 .* hessian(Cv -> strainenergy(Cv, self._lambda, self._mu), Cv);
		Dcurrent = fill(0.0, size(Dtotal))
		totlag2currsymm!(Dcurrent, Dtotal, Fn13d)
		fill!(D, 0.0)
		D[1:2, 1:2] = Dcurrent[1:2, 1:2]
		D[3, 3] = Dcurrent[4, 4]
		return D
	end
	function update2dstrain!(self::MatDeforNeohookeanAD, statev::FFltVec, cauchy::FFltVec, output::FFltVec, Fn1::FFltMat, Fn::FFltMat, tn::FFlt, dtn::FFlt, loc::FFltMat=zeros(3,1), label::FInt=0, quantity=:nothing)
		@assert length(cauchy) == nstressstrain(self.mr)
		Fn13d = fill(0.0, 3, 3)
		Fn13d[1:2, 1:2] = Fn1
		Fn13d[3, 3] = 1.0
		C = Fn13d'*Fn13d;
		Cv = fill(0.0, 6)
		strainttov!(D3D, Cv, C)
		Sv = 2 * gradient(Cv -> strainenergy(Cv, self._lambda, self._mu), Cv);
		S = fill(0.0, 3, 3)
		stressvtot!(D3D, S, Sv)
		cauchyt = Fn13d*(S/det(Fn1))*Fn13d'; # Cauchy stress, 3D
		stressttov!(mr, cauchy, cauchyt) # 2D stress vector
		if quantity == :nothing
			#Nothing to be copied to the output array
		elseif quantity == :cauchy || quantity == :Cauchy
			(length(output) >= 3) || (output = zeros(3)) # make sure we can store it
			copyto!(output, cauchy);
		elseif quantity == :pressure || quantity == :Pressure
			output[1]  =  -sum(tr(cauchyt))/3.
		elseif quantity == :princCauchy || quantity == :princcauchy
			ep = eigen(cauchyt);
			(length(output) >= 3) || (output = zeros(3)) # make sure we can store it
			copyto!(output,  sort(ep.values, rev=true));
		elseif quantity==:vonMises || quantity==:vonmises || quantity==:von_mises || quantity==:vm
			s1=cauchyt[1, 1]; s2=cauchyt[2, 2]; s3=cauchyt[3, 3];
			s4=cauchyt[1, 2]; s5=cauchyt[1, 3]; s6=cauchyt[2, 3];
			(length(output) >= 1) || (output = zeros(1)) # make sure we can store it
			output[1] = sqrt(1.0/2*((s1-s2)^2+(s1-s3)^2+(s2-s3)^2+6*(s4^2+s5^2+s6^2)))
		end
		return output
	end
	_lambda = E * nu / (1 + nu) / (1 - 2*(nu));
	_mu = E / 2. / (1+nu);
	return MatDeforNeohookeanAD(mr, mass_density, E, nu, _lambda, _mu,
		tangentmoduli2dstrain!, update2dstrain!)
end

end


text:
#==============================================================================#
# AWSException.jl
#
# Copyright Sam O'Connor 2014 - All rights reserved
#==============================================================================#

import JSON: json

export AWSException

abstract AWSException <: Exception

function Base.show(io::IO,e::AWSException)
    println(io, string(e.code,
                       e.message == "" ? "" : (" -- " * e.message), "\n",
                       e.cause))
end

function AWSException(e::HTTPException)

    code = string(http_status(e))
    message = "AWSException"
    info = Dict()

    # Extract API error code from Lambda-style JSON error message...
    if ismatch(r"json$", content_type(e))
        info = JSON.parse(http_message(e))
        message = get(info, "message", message)
    end

    # Extract API error code from JSON error message...
    if content_type(e) == "application/x-amz-json-1.0"
        json = JSON.parse(http_message(e))
        if haskey(json, "__type")
            code = split(json["__type"], "#")[2]
        end
    end

    # Extract API error code from XML error message...
    if (content_type(e) in ["", "application/xml", "text/xml"]
    &&  length(http_message(e)) > 0)
        info = parse_xml(http_message(e))
    end

    info = get(info, "Errors", info)
    info = get(info, "Error", info)
    code = get(info, "Code", code)
    message = get(info, "Message", message)

    # Create specialised exception object based on "code"...
    etype = Symbol(code)
    @repeat 2 try
        e = eval(:($etype($code, $message, $info, $e)))
    catch x
        @retry if isa(x, UndefVarError)
            eval(:(type $etype <: AWSException code; message; info; cause end))
        end
    end
end

AWSException(e) = e

#==============================================================================#
# End of file.
#==============================================================================#