text:
export CuArray, CuVector, CuMatrix, CuVecOrMat, cu, is_unified

## array storage

# array storage is shared by arrays that refer to the same data, while keeping track of
# the number of outstanding references

struct ArrayStorage{B}
  buffer::B

  ctx::CuContext

  # the refcount also encodes the state of the array:
  # < 0: unmanaged
  # = 0: freed
  # > 0: referenced
  refcount::Threads.Atomic{Int}
end

ArrayStorage(buf::B, ctx, state::Int) where {B} =
  ArrayStorage{B}(buf, ctx, Threads.Atomic{Int}(state))

## array type

mutable struct CuArray{T,N,B} <: AbstractGPUArray{T,N}
  storage::Union{Nothing,ArrayStorage{B}}

  maxsize::Int  # maximum data size; excluding any selector bytes
  offset::Int

  dims::Dims{N}

  function CuArray{T,N,B}(::UndefInitializer, dims::Dims{N}) where {T,N,B}
    Base.allocatedinline(T) || error("CuArray only supports element types that are stored inline")
    maxsize = prod(dims) * sizeof(T)
    bufsize = if Base.isbitsunion(T)
      # type tag array past the data
      maxsize + prod(dims)
    else
      maxsize
    end
    buf = alloc(B, bufsize)
    storage = ArrayStorage(buf, context(), 1)
    obj = new{T,N,B}(storage, maxsize, 0, dims)
    finalizer(unsafe_finalize!, obj)
  end

  function CuArray{T,N}(storage::ArrayStorage{B}, dims::Dims{N};
                        maxsize::Int=prod(dims) * sizeof(T), offset::Int=0) where {T,N,B}
    Base.allocatedinline(T) || error("CuArray only supports element types that are stored inline")
    return new{T,N,B}(storage, maxsize, offset, dims,)
  end
end

"""
    CUDA.unsafe_free!(a::CuArray, [stream::CuStream])

Release the memory of an array for reuse by future allocations. This function is
automatically called by the finalizer when an array goes out of scope, but can be called
earlier to reduce pressure on the memory allocator.

By default, the operation is performed on the task-local stream. During task or process
finalization however, that stream may be destroyed already, so be sure to specify a safe
stream (i.e. `CuDefaultStream()`, which will ensure the operation will block on other
streams) when calling this function from a finalizer. For simplicity, the `unsafe_finalize!`
function does exactly that.
"""
function unsafe_free!(xs::CuArray, stream::CuStream=stream())
  # this call should only have an effect once, because both the user and the GC can call it
  if xs.storage === nothing
    return
  elseif xs.storage.refcount[] < 0
    throw(ArgumentError("Cannot free an unmanaged buffer."))
  end

  refcount = Threads.atomic_add!(xs.storage.refcount, -1)
  if refcount == 1
    @context! skip_destroyed=true xs.storage.ctx begin
      free(xs.storage.buffer; stream)
    end
  end

  # this array object is now dead, so replace its storage by a dummy one
  xs.storage = nothing

  return
end

function unsafe_finalize!(xs::CuArray)
  # during task or process finalization, the local stream might be destroyed already, so
  # use the default stream. additionally, since we don't use per-thread APIs, this default
  # stream follows legacy semantics and will synchronize all other streams. this protects
  # against freeing resources that are still in use.
  #
  # TODO: although this is still an asynchronous operation, even when using the default
  # stream, it synchronizes "too much". we could do better, e.g., by keeping track of all
  # streams involved, or by refcounting uses and decrementing that refcount after the
  # operation using `cuLaunchHostFunc`. See CUDA.jl#778 and CUDA.jl#780 for details.
  unsafe_free!(xs, CuDefaultStream())
end

## alias detection

Base.dataids(A::CuArray) = (UInt(pointer(A.storage.buffer)),)

Base.unaliascopy(A::CuArray) = copy(A)

function Base.mightalias(A::CuArray, B::CuArray)
  rA = pointer(A):pointer(A)+sizeof(A)
  rB = pointer(B):pointer(B)+sizeof(B)
  return first(rA) <= first(rB) < last(rA) || first(rB) <= first(rA) < last(rB)
end

## convenience constructors

CuVector{T} = CuArray{T,1}
CuMatrix{T} = CuArray{T,2}
CuVecOrMat{T} = Union{CuVector{T},CuMatrix{T}}

# default to non-unified memory
CuArray{T,N}(::UndefInitializer, dims::Dims{N}) where {T,N} =
  CuArray{T,N,Mem.DeviceBuffer}(undef, dims)
is_unified(a::CuArray) = isa(a.storage.buffer, Mem.UnifiedBuffer)

# type and dimensionality specified, accepting dims as series of Ints
CuArray{T,N,B}(::UndefInitializer, dims::Integer...) where {T,N,B} =
  CuArray{T,N,B}(undef, dims)
CuArray{T,N}(::UndefInitializer, dims::Integer...) where {T,N} =
  CuArray{T,N}(undef, dims)

# type but not dimensionality specified
CuArray{T}(::UndefInitializer, dims::Dims{N}) where {T,N} =
  CuArray{T,N}(undef, dims)
CuArray{T}(::UndefInitializer, dims::Integer...) where {T} =
  CuArray{T}(undef, convert(Tuple{Vararg{Int}}, dims))

# empty vector constructor
CuArray{T,1,B}() where {T,B} = CuArray{T,1,B}(undef, 0)
CuArray{T,1}() where {T} = CuArray{T,1}(undef, 0)

# do-block constructors
for (ctor, tvars) in (:CuArray => (),
                      :(CuArray{T}) => (:T,),
                      :(CuArray{T,N}) => (:T, :N),
                      :(CuArray{T,N,B}) => (:T, :N, :B))
  @eval begin
    function $ctor(f::Function, args...) where {$(tvars...)}
      xs = $ctor(args...)
      try
        f(xs)
      finally
        unsafe_free!(xs)
      end
    end
  end
end

Base.similar(a::CuArray{T,N,B}) where {T,N,B} =
  CuArray{T,N,B}(undef, size(a))
Base.similar(a::CuArray{T,<:Any,B}, dims::Base.Dims{N}) where {T,N,B} =
  CuArray{T,N,B}(undef, dims)
Base.similar(a::CuArray{<:Any,<:Any,B}, ::Type{T}, dims::Base.Dims{N}) where {T,N,B} =
  CuArray{T,N,B}(undef, dims)

function Base.copy(a::CuArray{T,N}) where {T,N}
  b = similar(a)
  @inbounds copyto!(b, a)
end

"""
  unsafe_wrap(::CuArray, ptr::CuPtr{T}, dims; own=false, ctx=context())

Wrap a `CuArray` object around the data at the address given by `ptr`. The pointer
element type `T` determines the array element type. `dims` is either an integer (for a 1d
array) or a tuple of the array dimensions. `own` optionally specified whether Julia should
take ownership of the memory, calling `cudaFree` when the array is no longer referenced. The
`ctx` argument determines the CUDA context where the data is allocated in.
"""
function Base.unsafe_wrap(::Union{Type{CuArray},Type{CuArray{T}},Type{CuArray{T,N}}},
                          ptr::CuPtr{T}, dims::NTuple{N,Int};
                          own::Bool=false, ctx::CuContext=context()) where {T,N}
  Base.isbitstype(T) || error("Can only unsafe_wrap a pointer to a bits type")
  sz = prod(dims) * sizeof(T)

  # identify the buffer
  buf = try
    typ = memory_type(ptr)
    if is_managed(ptr)
      Mem.UnifiedBuffer(ptr, sz)
    elseif typ == CU_MEMORYTYPE_DEVICE
      # TODO: can we identify whether this pointer was allocated asynchronously?
      Mem.DeviceBuffer(ptr, sz, false)
    elseif typ == CU_MEMORYTYPE_HOST
      error("Cannot unsafe_wrap a host pointer with a CuArray")
    else
      error("Unknown memory type; please file an issue.")
    end
  catch err
      error("Could not identify the buffer type; are you passing a valid CUDA pointer to unsafe_wrap?")
  end

  storage = ArrayStorage(buf, ctx, own ? 1 : -1)
  CuArray{T, length(dims)}(storage, dims)
end

function Base.unsafe_wrap(Atype::Union{Type{CuArray},Type{CuArray{T}},Type{CuArray{T,1}}},
                          p::CuPtr{T}, dim::Integer;
                          own::Bool=false, ctx::CuContext=context()) where {T}
  unsafe_wrap(Atype, p, (dim,); own, ctx)
end

Base.unsafe_wrap(T::Type{<:CuArray}, ::Ptr, dims::NTuple{N,Int}; kwargs...) where {N} =
  throw(ArgumentError("cannot wrap a CPU pointer with a $T"))

## array interface

Base.elsize(::Type{<:CuArray{T}}) where {T} = sizeof(T)

Base.size(x::CuArray) = x.dims
Base.sizeof(x::CuArray) = Base.elsize(x) * length(x)

## derived types

export DenseCuArray, DenseCuVector, DenseCuMatrix, DenseCuVecOrMat,
       StridedCuArray, StridedCuVector, StridedCuMatrix, StridedCuVecOrMat,
       AnyCuArray, AnyCuVector, AnyCuMatrix, AnyCuVecOrMat

# dense arrays: stored contiguously in memory
#
# all common dense wrappers are currently represented as CuArray objects.
# this simplifies common use cases, and greatly improves load time.
# CUDA.jl 2.0 experimented with using ReshapedArray/ReinterpretArray/SubArray,
# but that proved much too costly. TODO: revisit when we have better Base support.
DenseCuArray{T,N} = CuArray{T,N}
DenseCuVector{T} = DenseCuArray{T,1}
DenseCuMatrix{T} = DenseCuArray{T,2}
DenseCuVecOrMat{T} = Union{DenseCuVector{T}, DenseCuMatrix{T}}

# strided arrays
StridedSubCuArray{T,N,I<:Tuple{Vararg{Union{Base.RangeIndex, Base.ReshapedUnitRange,
                                            Base.AbstractCartesianIndex}}}} =
  SubArray{T,N,<:CuArray,I}
StridedCuArray{T,N} = Union{CuArray{T,N}, StridedSubCuArray{T,N}}
StridedCuVector{T} = StridedCuArray{T,1}
StridedCuMatrix{T} = StridedCuArray{T,2}
StridedCuVecOrMat{T} = Union{StridedCuVector{T}, StridedCuMatrix{T}}

Base.pointer(x::StridedCuArray{T}) where {T} = Base.unsafe_convert(CuPtr{T}, x)
@inline function Base.pointer(x::StridedCuArray{T}, i::Integer) where T
    Base.unsafe_convert(CuPtr{T}, x) + Base._memory_offset(x, i)
end

# anything that's (secretly) backed by a CuArray
AnyCuArray{T,N} = Union{CuArray{T,N}, WrappedArray{T,N,CuArray,CuArray{T,N}}}
AnyCuVector{T} = AnyCuArray{T,1}
AnyCuMatrix{T} = AnyCuArray{T,2}
AnyCuVecOrMat{T} = Union{AnyCuVector{T}, AnyCuMatrix{T}}

## interop with other arrays

@inline function CuArray{T,N,B}(xs::AbstractArray{<:Any,N}) where {T,N,B}
  A = CuArray{T,N,B}(undef, size(xs))
  copyto!(A, convert(Array{T}, xs))
  return A
end

@inline CuArray{T,N}(xs::AbstractArray{<:Any,N}) where {T,N} =
  CuArray{T,N,Mem.Device}(xs)

@inline CuArray{T,N}(xs::CuArray{<:Any,N,B}) where {T,N,B} =
  CuArray{T,N,B}(xs)

# underspecified constructors
CuArray{T}(xs::AbstractArray{S,N}) where {T,N,S} = CuArray{T,N}(xs)
(::Type{CuArray{T,N} where T})(x::AbstractArray{S,N}) where {S,N} = CuArray{S,N}(x)
CuArray(A::AbstractArray{T,N}) where {T,N} = CuArray{T,N}(A)

# idempotency
CuArray{T,N,B}(xs::CuArray{T,N,B}) where {T,N,B} = xs
CuArray{T,N}(xs::CuArray{T,N,B}) where {T,N,B} = xs

## conversions

Base.convert(::Type{T}, x::T) where T <: CuArray = x

## interop with C libraries

Base.unsafe_convert(::Type{Ptr{T}}, x::CuArray{T}) where {T} =
  throw(ArgumentError("cannot take the CPU address of a $(typeof(x))"))
Base.unsafe_convert(::Type{CuPtr{T}}, x::CuArray{T}) where {T} =
  convert(CuPtr{T}, pointer(x.storage.buffer)) + x.offset

## interop with device arrays

function Base.unsafe_convert(::Type{CuDeviceArray{T,N,AS.Global}}, a::DenseCuArray{T,N}) where {T,N}
  CuDeviceArray{T,N,AS.Global}(size(a), reinterpret(LLVMPtr{T,AS.Global}, pointer(a)),
                               a.maxsize - a.offset)
end

## interop with CPU arrays

typetagdata(a::Array, i=1) = ccall(:jl_array_typetagdata, Ptr{UInt8}, (Any,), a) + i - 1
typetagdata(a::CuArray, i=1) =
  convert(CuPtr{UInt8}, pointer(a.storage.buffer) + a.maxsize) + a.offset÷Base.elsize(a) + i - 1

# We don't convert isbits types in `adapt`, since they are already
# considered GPU-compatible.

Adapt.adapt_storage(::Type{CuArray}, xs::AT) where {AT<:AbstractArray} =
  isbitstype(AT) ? xs : convert(CuArray, xs)

# if an element type is specified, convert to it
Adapt.adapt_storage(::Type{<:CuArray{T}}, xs::AT) where {T, AT<:AbstractArray} =
  isbitstype(AT) ? xs : convert(CuArray{T}, xs)

function Base.copyto!(dest::DenseCuArray{T}, doffs::Integer, src::Array{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::DenseCuArray{T}, src::Array{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::Array{T}, doffs::Integer, src::DenseCuArray{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::Array{T}, src::DenseCuArray{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.copyto!(dest::DenseCuArray{T}, doffs::Integer, src::DenseCuArray{T}, soffs::Integer,
                      n::Integer) where T
  n==0 && return dest
  @boundscheck checkbounds(dest, doffs)
  @boundscheck checkbounds(dest, doffs+n-1)
  @boundscheck checkbounds(src, soffs)
  @boundscheck checkbounds(src, soffs+n-1)
  unsafe_copyto!(dest, doffs, src, soffs, n)
  return dest
end

Base.copyto!(dest::DenseCuArray{T}, src::DenseCuArray{T}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

function Base.unsafe_copyto!(dest::DenseCuArray{T}, doffs, src::Array{T}, soffs, n) where T
  if !is_pinned(pointer(src))
    # operations on unpinned memory cannot be executed asynchronously, and synchronize
    # without yielding back to the Julia scheduler. prevent that by eagerly synchronizing.
    synchronize()
  end
  GC.@preserve src dest begin
    unsafe_copyto!(pointer(dest, doffs), pointer(src, soffs), n; async=true)
    if Base.isbitsunion(T)
      unsafe_copyto!(typetagdata(dest, doffs), typetagdata(src, soffs), n; async=true)
    end
  end
  return dest
end

function Base.unsafe_copyto!(dest::Array{T}, doffs, src::DenseCuArray{T}, soffs, n) where T
  if !is_pinned(pointer(dest))
    # operations on unpinned memory cannot be executed asynchronously, and synchronize
    # without yielding back to the Julia scheduler. prevent that by eagerly synchronizing.
    synchronize()
  end
  GC.@preserve src dest begin
    unsafe_copyto!(pointer(dest, doffs), pointer(src, soffs), n; async=true)
    if Base.isbitsunion(T)
      unsafe_copyto!(typetagdata(dest, doffs), typetagdata(src, soffs), n; async=true)
    end
  end
  synchronize() # users expect values to be available after this call
  return dest
end

function Base.unsafe_copyto!(dest::DenseCuArray{T}, doffs, src::DenseCuArray{T}, soffs, n) where T
  GC.@preserve src dest begin
    unsafe_copyto!(pointer(dest, doffs), pointer(src, soffs), n; async=true)
    if Base.isbitsunion(T)
      unsafe_copyto!(typetagdata(dest, doffs), typetagdata(src, soffs), n; async=true)
    end
  end
  return dest
end

# XXX: defining deepcopy_internal, as per the deepcopy documentation, results in a ton
#      of invalidations, so we redefine deepcopy itself (see JuliaGPU/CUDA.jl#632)
function Base.deepcopy(x::CuArray)
  dict = IdDict()
  haskey(dict, x) && return dict[x]::typeof(x)
  return dict[x] = copy(x)
end

## opinionated gpu array adaptor

# eagerly converts Float64 to Float32, for performance reasons

struct CuArrayAdaptor{B} end

Adapt.adapt_storage(::CuArrayAdaptor{B}, xs::AbstractArray{T,N}) where {T,N,B} =
  isbits(xs) ? xs : CuArray{T,N,B}(xs)

Adapt.adapt_storage(::CuArrayAdaptor{B}, xs::AbstractArray{T,N}) where {T<:AbstractFloat,N,B} =
  isbits(xs) ? xs : CuArray{Float32,N,B}(xs)

Adapt.adapt_storage(::CuArrayAdaptor{B}, xs::AbstractArray{T,N}) where {T<:Complex{<:AbstractFloat},N,B} =
  isbits(xs) ? xs : CuArray{ComplexF32,N,B}(xs)

# not for Float16
Adapt.adapt_storage(::CuArrayAdaptor{B}, xs::AbstractArray{T,N}) where {T<:Union{Float16,BFloat16},N,B} =
  isbits(xs) ? xs : CuArray{T,N,B}(xs)

@inline cu(xs; unified::Bool=false) = adapt(CuArrayAdaptor{unified ? Mem.UnifiedBuffer : Mem.DeviceBuffer}(), xs)
Base.getindex(::typeof(cu), xs...) = CuArray([xs...])

## utilities

zeros(T::Type, dims...) = fill!(CuArray{T}(undef, dims...), 0)
ones(T::Type, dims...) = fill!(CuArray{T}(undef, dims...), 1)
zeros(dims...) = zeros(Float32, dims...)
ones(dims...) = ones(Float32, dims...)
fill(v, dims...) = fill!(CuArray{typeof(v)}(undef, dims...), v)
fill(v, dims::Dims) = fill!(CuArray{typeof(v)}(undef, dims...), v)

# optimized implementation of `fill!` for types that are directly supported by memset
memsettype(T::Type) = T
memsettype(T::Type{<:Signed}) = unsigned(T)
memsettype(T::Type{<:AbstractFloat}) = Base.uinttype(T)
const MemsetCompatTypes = Union{UInt8, Int8,
                                UInt16, Int16, Float16,
                                UInt32, Int32, Float32}
function Base.fill!(A::DenseCuArray{T}, x) where T <: MemsetCompatTypes
  U = memsettype(T)
  y = reinterpret(U, convert(T, x))
  Mem.set!(convert(CuPtr{U}, pointer(A)), y, length(A))
  A
end

## views

# optimize view to return a CuArray when contiguous

struct Contiguous end
struct NonContiguous end

# NOTE: this covers more cases than the I<:... in Base.FastContiguousSubArray
CuIndexStyle() = Contiguous()
CuIndexStyle(I...) = NonContiguous()
CuIndexStyle(::Union{Base.ScalarIndex, CartesianIndex}...) = Contiguous()
CuIndexStyle(i1::Colon, ::Union{Base.ScalarIndex, CartesianIndex}...) = Contiguous()
CuIndexStyle(i1::AbstractUnitRange, ::Union{Base.ScalarIndex, CartesianIndex}...) = Contiguous()
CuIndexStyle(i1::Colon, I...) = CuIndexStyle(I...)

cuviewlength() = ()
@inline cuviewlength(::Real, I...) = cuviewlength(I...) # skip scalar

if VERSION >= v"1.8.0-DEV.120"
@inline cuviewlength(i1::AbstractUnitRange, I...) = (Base.length(i1), cuviewlength(I...)...)
@inline cuviewlength(i1::AbstractUnitRange, ::Base.ScalarIndex...) = (Base.length(i1),)
else
@inline cuviewlength(i1::AbstractUnitRange, I...) = (Base.unsafe_length(i1), cuviewlength(I...)...)
@inline cuviewlength(i1::AbstractUnitRange, ::Base.ScalarIndex...) = (Base.unsafe_length(i1),)
end

# we don't really want an array, so don't call `adapt(Array, ...)`,
# but just want CuArray indices to get downloaded back to the CPU.
# this makes sure we preserve array-like containers, like Base.Slice.
struct BackToCPU end
Adapt.adapt_storage(::BackToCPU, xs::CuArray) = convert(Array, xs)

@inline function Base.view(A::CuArray, I::Vararg{Any,N}) where {N}
    J = to_indices(A, I)
    @boundscheck begin
        # Base's boundscheck accesses the indices, so make sure they reside on the CPU.
        # this is expensive, but it's a bounds check after all.
        J_cpu = map(j->adapt(BackToCPU(), j), J)
        checkbounds(A, J_cpu...)
    end
    J_gpu = map(j->adapt(CuArray, j), J)
    unsafe_view(A, J_gpu, CuIndexStyle(I...))
end

@inline function unsafe_view(A, I, ::Contiguous)
    unsafe_contiguous_view(Base._maybe_reshape_parent(A, Base.index_ndims(I...)), I, cuviewlength(I...))
end
@inline function unsafe_contiguous_view(a::CuArray{T}, I::NTuple{N,Base.ViewIndex}, dims::NTuple{M,Integer}) where {T,N,M}
    offset = Base.compute_offset1(a, 1, I) * sizeof(T)

    refcount = a.storage.refcount[]
    @assert refcount != 0
    if refcount > 0
      Threads.atomic_add!(a.storage.refcount, 1)
    end

    b = CuArray{T,M}(a.storage, dims; a.maxsize, offset=a.offset+offset)
    if refcount > 0
        finalizer(unsafe_finalize!, b)
    end
    return b
end

@inline function unsafe_view(A, I, ::NonContiguous)
    Base.unsafe_view(Base._maybe_reshape_parent(A, Base.index_ndims(I...)), I...)
end

# pointer conversions
## contiguous
function Base.unsafe_convert(::Type{CuPtr{T}}, V::SubArray{T,N,P,<:Tuple{Vararg{Base.RangeIndex}}}) where {T,N,P}
    return Base.unsafe_convert(CuPtr{T}, parent(V)) +
           Base._memory_offset(V.parent, map(first, V.indices)...)
end
## reshaped
function Base.unsafe_convert(::Type{CuPtr{T}}, V::SubArray{T,N,P,<:Tuple{Vararg{Union{Base.RangeIndex,Base.ReshapedUnitRange}}}}) where {T,N,P}
   return Base.unsafe_convert(CuPtr{T}, parent(V)) +
          (Base.first_index(V)-1)*sizeof(T)
end

## PermutedDimsArray

Base.unsafe_convert(::Type{CuPtr{T}}, A::PermutedDimsArray) where {T} =
    Base.unsafe_convert(CuPtr{T}, parent(A))

## reshape

# optimize reshape to return a CuArray

function Base.reshape(a::CuArray{T,M}, dims::NTuple{N,Int}) where {T,N,M}
  if prod(dims) != length(a)
      throw(DimensionMismatch("new dimensions $(dims) must be consistent with array size $(size(a))"))
  end

  if N == M && dims == size(a)
      return a
  end

  refcount = a.storage.refcount[]
  @assert refcount != 0
  if refcount > 0
    Threads.atomic_add!(a.storage.refcount, 1)
  end

  b = CuArray{T,N}(a.storage, dims; a.maxsize, a.offset)
  if refcount > 0
      finalizer(unsafe_finalize!, b)
  end
  return b
end

## reinterpret

# optimize reshape to return a CuArray

struct _CuReinterpretBitsTypeError{T,A} <: Exception end
function Base.showerror(io::IO, ::_CuReinterpretBitsTypeError{T, <:AbstractArray{S}}) where {T, S}
  print(io, "cannot reinterpret an `$(S)` array to `$(T)`, because not all types are bitstypes")
end

struct _CuReinterpretZeroDimError{T,A} <: Exception end
function Base.showerror(io::IO, ::_CuReinterpretZeroDimError{T, <:AbstractArray{S,N}}) where {T, S, N}
  print(io, "cannot reinterpret a zero-dimensional `$(S)` array to `$(T)` which is of a different size")
end

struct _CuReinterpretDivisibilityError{T,A} <: Exception
  dim::Int
end
function Base.showerror(io::IO, err::_CuReinterpretDivisibilityError{T, <:AbstractArray{S,N}}) where {T, S, N}
  dim = err.dim
  print(io, """
      cannot reinterpret an `$(S)` array to `$(T)` whose first dimension has size `$(dim)`.
      The resulting array would have non-integral first dimension.
      """)
end

struct _CuReinterpretFirstIndexError{T,A,Ax1} <: Exception
  ax1::Ax1
end
function Base.showerror(io::IO, err::_CuReinterpretFirstIndexError{T, <:AbstractArray{S,N}}) where {T, S, N}
  ax1 = err.ax1
  print(io, "cannot reinterpret a `$(S)` array to `$(T)` when the first axis is $ax1. Try reshaping first.")
end

function _reinterpret_exception(::Type{T}, a::AbstractArray{S,N}) where {T,S,N}
  if !isbitstype(T) || !isbitstype(S)
    return _CuReinterpretBitsTypeError{T,typeof(a)}()
  end
  if N == 0 && sizeof(T) != sizeof(S)
    return _CuReinterpretZeroDimError{T,typeof(a)}()
  end
  if N != 0 && sizeof(S) != sizeof(T)
      ax1 = axes(a)[1]
      dim = length(ax1)
      if Base.rem(dim*sizeof(S),sizeof(T)) != 0
        return _CuReinterpretDivisibilityError{T,typeof(a)}(dim)
      end
      if first(ax1) != 1
        return _CuReinterpretFirstIndexError{T,typeof(a),typeof(ax1)}(ax1)
      end
  end
  return nothing
end

function Base.reinterpret(::Type{T}, a::CuArray{S,N}) where {T,S,N}
  err = _reinterpret_exception(T, a)
  err === nothing || throw(err)

  if sizeof(T) == sizeof(S) # for N == 0
    osize = size(a)
  else
    isize = size(a)
    size1 = div(isize[1]*sizeof(S), sizeof(T))
    osize = tuple(size1, Base.tail(isize)...)
  end

  refcount = a.storage.refcount[]
  @assert refcount != 0
  if refcount > 0
    Threads.atomic_add!(a.storage.refcount, 1)
  end

  b = CuArray{T,N}(a.storage, osize; a.maxsize, a.offset)
  if refcount > 0
      finalizer(unsafe_finalize!, b)
  end
  return b
end

## resizing

"""
  resize!(a::CuVector, n::Int)

Resize `a` to contain `n` elements. If `n` is smaller than the current collection length,
the first `n` elements will be retained. If `n` is larger, the new elements are not
guaranteed to be initialized.

Note that this operation is only supported on managed buffers, i.e., not on arrays that are
created by `unsafe_wrap` with `own=false`.
"""
function Base.resize!(A::CuVector{T}, n::Int) where T
  # TODO: add additional space to allow for quicker resizing
  maxsize = n * sizeof(T)
  bufsize = if Base.isbitsunion(T)
    # type tag array past the data
    maxsize + n
  else
    maxsize
  end

  new_storage = @context! A.storage.ctx begin
    buf = alloc(typeof(A.storage.buffer), bufsize)
    ptr = convert(CuPtr{T}, buf)
    m = Base.min(length(A), n)
    unsafe_copyto!(ptr, pointer(A), m)
    ArrayStorage(buf, A.storage.ctx, 1)
  end

  unsafe_free!(A)
  A.storage = new_storage
  A.dims = (n,)
  A.maxsize = maxsize
  A.offset = 0

  A
end


text:
using JuMP, Gurobi, OneClassSampling, SVDD, OneClassActiveLearning
include(joinpath(@__DIR__, "config_syn.jl"))

experiment_name="synthetic"

solver = with_optimizer(Gurobi.Optimizer; OutputFlag=0, Threads=1)
quality_metrics = Dict(:mcc => matthews_corr, :kappa => cohens_kappa, :f1 => f1_score)

threshold_strategies = [
    GroundTruthThresholdStrategy(),
]

model_init_strategy = Dict{Symbol, Any}(
                            :svdd => Dict{Symbol, Any}(
                                :MMC_C1 => SVDD.SimpleCombinedStrategy(SVDD.ModifiedMeanCriterion(), SVDD.FixedCStrategy(1)),
                            ),
                            :sample_qe => Dict{Symbol, Any}(
                                :threshold_strategies => threshold_strategies
                        ))

KDE_GAMMA = :mod_mean_crit

competitors = [
    BPS(nothing, 0.005),
    DAEDS(30, 0.1, 0.3),
    DBSRSVDD(7, 0.3),
    HSR(20, 0.01),
    IESRSVDD(KDE_GAMMA, 0.5),
    FBPE(360),
    KFNCBD(100, 0.2),
    NDPSR(20, 7)
]

density_sampler = [
    RAPID(OutlierPercentageThresholdStrategy(0.0), KDE_GAMMA),
]

sampling_strategies = Vector{Dict{Symbol, Any}}()

for c in vcat(competitors, density_sampler)
    push!(sampling_strategies, Dict{Symbol, Any}(
        :method => PreFilteringWrapper(c, KDE_GAMMA, GroundTruthThresholdStrategy()),
        :seed => 0
    ))
end

ext: jl

max_stars_repo_name: englhardt/ocs-evaluation

max_issues_repo_name: englhardt/ocs-evaluation

max_forks_repo_name: englhardt/ocs-evaluation


text:
using Flux3D, BenchmarkTools, Printf

function setup_benchmark_record(names)
    benchmarks = Dict{String,Vector{Float64}}()
    for name in names
        benchmarks[name] = []
    end
    return benchmarks
end

function generate_point_cloud(npoints::Int)
    points = ones(3, npoints)
    points = cumsum(points, dims = 2)
    return PointCloud(points / npoints)
end

function generate_trimesh(npoints::Int)
    v = ones(3, npoints)
    v = cumsum(v, dims = 2) / npoints
    f = reshape(collect(1:npoints*3), 3, npoints)
    return TriMesh([v], [f])
end

function run_benchmarks!(benchmarks, x, npoints, generate_func, benchmark_func, device)
    for (transform, name) in x
        transform = transform |> device
        pc = generate_func(npoints) |> device
        trial = @benchmark $benchmark_func($transform, $pc)
        time = minimum(trial.times) * 1.0e-6
        println("$name: $time ms")
        push!(benchmarks[name], time)
    end
end

function realign_point_cloud(npoints)
    pc = generate_point_cloud(npoints)
    rot = RotatePointCloud(-ROT_MATRIX)
    return rot(pc)
end

function realign_trimesh(_npoints)
    pc = generate_trimesh(_npoints)
    rot = RotateTriMesh(-ROT_MATRIX)
    return rot(pc)
end

ROT_MATRIX = [
    1.0 2.0 3.0
    0.2 0.5 0.9
    3.0 2.0 1.0
]

npoint_arr = 2 .^ [12, 14, 16, 18, 20]

names = [
    "ScalePointCloud",
    "RotatePointCloud",
    "ReAlignPointCloud",
    "NormalizePointCloud",
    "Chain",
]
cpu_bm_pcloud = setup_benchmark_record(names)
gpu_bm_pcloud = setup_benchmark_record(names)

names = ["ScaleTriMesh", "RotateTriMesh", "ReAlignTriMesh", "NormalizeTriMesh", "Chain"]
cpu_bm_trimesh = setup_benchmark_record(names)
gpu_bm_trimesh = setup_benchmark_record(names)

println("DEVICE: CPU")
for _npoints in npoint_arr
    pcloud_arr = [
        (ScalePointCloud(0.5; inplace = false), "ScalePointCloud"),
        (RotatePointCloud(ROT_MATRIX; inplace = false), "RotatePointCloud"),
        (
            ReAlignPointCloud(realign_point_cloud(_npoints); inplace = false),
            "ReAlignPointCloud",
        ),
        (NormalizePointCloud(inplace = false), "NormalizePointCloud"),
        (
            Chain(
                ScalePointCloud(0.5; inplace = false),
                RotatePointCloud(ROT_MATRIX; inplace = false),
                ReAlignPointCloud(realign_point_cloud(_npoints); inplace = false),
                NormalizePointCloud(),
            ),
            "Chain",
        ),
    ]

    trimesh_arr = [
        (ScaleTriMesh(0.5; inplace = true), "ScaleTriMesh"),
        (RotateTriMesh(ROT_MATRIX; inplace = true), "RotateTriMesh"),
        (ReAlignTriMesh(realign_trimesh(_npoints); inplace = true), "ReAlignTriMesh"),
        (NormalizeTriMesh(inplace = true), "NormalizeTriMesh"),
        (
            Chain(
                ScaleTriMesh(0.5; inplace = true),
                RotateTriMesh(ROT_MATRIX; inplace = true),
                ReAlignTriMesh(realign_trimesh(_npoints); inplace = true),
                NormalizeTriMesh(inplace = true),
            ),
            "Chain",
        ),
    ]
    println("Running benchmarks for npoints = $_npoints")
    run_benchmarks!(
        cpu_bm_pcloud,
        pcloud_arr,
        _npoints,
        generate_point_cloud,
        (op, pc) -> op(pc),
        cpu,
    )
    run_benchmarks!(
        cpu_bm_trimesh,
        trimesh_arr,
        _npoints,
        generate_trimesh,
        (op, pc) -> op(pc),
        cpu,
    )
    println()
end

using CUDA
if has_cuda()
    println("CUDA is on. Running GPU Benchmarks")
    println("DEVICE: GPU")
    for _npoints in npoint_arr
        pcloud_arr = [
            (ScalePointCloud(0.5; inplace = false), "ScalePointCloud"),
            (RotatePointCloud(ROT_MATRIX; inplace = false), "RotatePointCloud"),
            (
                ReAlignPointCloud(realign_point_cloud(_npoints); inplace = false),
                "ReAlignPointCloud",
            ),
            (NormalizePointCloud(inplace = false), "NormalizePointCloud"),
            (
                Chain(
                    ScalePointCloud(0.5; inplace = false),
                    RotatePointCloud(ROT_MATRIX; inplace = false),
                    ReAlignPointCloud(realign_point_cloud(_npoints); inplace = false),
                    NormalizePointCloud(),
                ),
                "Chain",
            ),
        ]

        trimesh_arr = [
            (ScaleTriMesh(0.5; inplace = true), "ScaleTriMesh"),
            (RotateTriMesh(ROT_MATRIX; inplace = true), "RotateTriMesh"),
            (ReAlignTriMesh(realign_trimesh(_npoints); inplace = true), "ReAlignTriMesh"),
            (NormalizeTriMesh(inplace = true), "NormalizeTriMesh"),
            (
                Chain(
                    ScaleTriMesh(0.5; inplace = true),
                    RotateTriMesh(ROT_MATRIX; inplace = true),
                    ReAlignTriMesh(realign_trimesh(_npoints); inplace = true),
                    NormalizeTriMesh(inplace = true),
                ),
                "Chain",
            ),
        ]
        println("Running benchmarks for npoints = $_npoints")
        run_benchmarks!(
            gpu_bm_pcloud,
            pcloud_arr,
            _npoints,
            generate_point_cloud,
            (op, pc) -> (CUDA.@sync op(pc)),
            gpu,
        )
        run_benchmarks!(
            gpu_bm_trimesh,
            trimesh_arr,
            _npoints,
            generate_trimesh,
            (op, pc) -> (CUDA.@sync op(pc)),
            gpu,
        )
        println()
    end
end

function save_bm(fname, rep, cpu_benchmarks, gpu_benchmarks)
    open(fname, "w") do io
        device = "cpu"
        for (key, values) in cpu_benchmarks
            for (p, v) in zip(npoint_arr, values)
                Printf.@printf(io, "%s %s %s %d %f ms\n", rep, device, key, p, v)
            end
        end

        device = "gpu"
        for (key, values) in gpu_benchmarks
            for (p, v) in zip(npoint_arr, values)
                Printf.@printf(io, "%s %s %s %d %f ms\n", rep, device, key, p, v)
            end
        end
    end
end

fname = joinpath(@__DIR__, "bm_flux3d.txt")
save_bm(fname, "PointCloud", cpu_bm_pcloud, gpu_bm_pcloud)
save_bm(fname, "TriMesh", cpu_bm_trimesh, gpu_bm_trimesh)
@info "Benchmarks have been saved at $fname"


text:
function get_nass(args...; format="json")
    key = ENV["USDA_QUICK_SURVEY_KEY"]

    header = string(usda_url, "/api/api_GET/?key=", key, "&format=$format")

    query = ""
    for i in args
        arg = string("&", i)
        query *= arg
    end

    request("GET", string(header, query))
end


text:
@testset "894.all-possible-full-binary-trees.jl" begin
    @test all_possible_fbt(7) == [TreeNode{Int64}(0, TreeNode{Int64}(0, nothing, nothing),
                           TreeNode{Int64}(0, TreeNode{Int64}(0, nothing, nothing),
                                           TreeNode{Int64}(0,
                                                           TreeNode{Int64}(0, nothing,
                                                                           nothing),
                                                           TreeNode{Int64}(0, nothing,
                                                                           nothing)))),
           TreeNode{Int64}(0, TreeNode{Int64}(0, nothing, nothing),
                           TreeNode{Int64}(0,
                                           TreeNode{Int64}(0,
                                                           TreeNode{Int64}(0, nothing,
                                                                           nothing),
                                                           TreeNode{Int64}(0, nothing,
                                                                           nothing)),
                                           TreeNode{Int64}(0, nothing, nothing))),
           TreeNode{Int64}(0,
                           TreeNode{Int64}(0, TreeNode{Int64}(0, nothing, nothing),
                                           TreeNode{Int64}(0, nothing, nothing)),
                           TreeNode{Int64}(0, TreeNode{Int64}(0, nothing, nothing),
                                           TreeNode{Int64}(0, nothing, nothing))),
           TreeNode{Int64}(0,
                           TreeNode{Int64}(0, TreeNode{Int64}(0, nothing, nothing),
                                           TreeNode{Int64}(0,
                                                           TreeNode{Int64}(0, nothing,
                                                                           nothing),
                                                           TreeNode{Int64}(0, nothing,
                                                                           nothing))),
                           TreeNode{Int64}(0, nothing, nothing)),
           TreeNode{Int64}(0,
                           TreeNode{Int64}(0,
                                           TreeNode{Int64}(0,
                                                           TreeNode{Int64}(0, nothing,
                                                                           nothing),
                                                           TreeNode{Int64}(0, nothing,
                                                                           nothing)),
                                           TreeNode{Int64}(0, nothing, nothing)),
                           TreeNode{Int64}(0, nothing, nothing))]
    @test all_possible_fbt(3) == [TreeNode{Int64}(0, TreeNode{Int64}(0, nothing, nothing),
                           TreeNode{Int64}(0, nothing, nothing))]
    @test all_possible_fbt(2) == TreeNode{Int}[]
end


text:
# rotate quaternion
quaternion_rotate(q1::Quaternion,q2::Quaternion) = q2 * q1 / q2

# rotate vector
vector_rotate(v::AbstractVector,q::Quaternion) = Vmat(quaternion_rotate(Quaternion(v), q))
# vector_rotate(v::StaticVector,q::Quaternion) = q*v
∂vector_rotate∂p(p::AbstractVector, q::Quaternion) = VRᵀmat(q) * LVᵀmat(q)
∂vector_rotate∂q(p::AbstractVector, q::Quaternion) = VLmat(q) * Lmat(Quaternion(p)) * Tmat() + VRᵀmat(q) * Rmat(Quaternion(p))

# rotation matrix
rotation_matrix(q::Quaternion) = VRᵀmat(q) * LVᵀmat(q)
∂rotation_matrix∂q(q::Quaternion, p::AbstractVector) = ∂VRᵀmat∂q(LVᵀmat(q) * p) + VRᵀmat(q) * ∂LVᵀmat∂q(p)
∂rotation_matrix_inv∂q(q::Quaternion, p::AbstractVector) = ∂rotation_matrix∂q(inv(q), p) * Tmat()


text:
using AbstractMCMC
using AdvancedMH
using Clustering
using Distributions
using DistributionsAD
using FiniteDifferences
using ForwardDiff
using MCMCChains
using Memoization
using NamedArrays
using Optim
using PDMats
using ReverseDiff
using SpecialFunctions
using StatsBase
using StatsFuns
using Tracker
using Turing
using Turing.Inference
using Turing.RandomMeasures
using Zygote

using LinearAlgebra
using Pkg
using Random
using Test

using AdvancedPS: ResampleWithESSThreshold, resample_systematic, resample_multinomial
using AdvancedVI: TruncatedADAGrad, DecayedADAGrad, apply!
using Distributions: Binomial, logpdf
using DynamicPPL: getval, getlogp
using ForwardDiff: Dual
using MCMCChains: Chains
using StatsFuns: binomlogpdf, logistic, logsumexp
using Turing: BinomialLogit, ForwardDiffAD, Sampler, SampleFromPrior, NUTS, TrackerAD,
                Variational, ZygoteAD, getspace, gradient_logp
using Turing.Core: TuringDenseMvNormal, TuringDiagMvNormal
using Turing.Variational: TruncatedADAGrad, DecayedADAGrad, AdvancedVI

setprogress!(false)

include("test_utils/AllUtils.jl")

@testset "Turing" begin
    @testset "core" begin
        include("core/ad.jl")
    end

    Turing.setrdcache(false)
    for adbackend in (:forwarddiff, :tracker, :reversediff)
        Turing.setadbackend(adbackend)
        @testset "inference: $adbackend" begin
            @testset "samplers" begin
                include("inference/gibbs.jl")
                include("inference/gibbs_conditional.jl")
                include("inference/hmc.jl")
                include("inference/is.jl")
                include("inference/mh.jl")
                include("inference/ess.jl")
                include("inference/emcee.jl")
                include("inference/AdvancedSMC.jl")
                include("inference/Inference.jl")
                include("contrib/inference/dynamichmc.jl")
                include("contrib/inference/sghmc.jl")
            end
        end

        @testset "variational algorithms : $adbackend" begin
            include("variational/advi.jl")
        end

        @testset "modes" begin
            include("modes/ModeEstimation.jl")
        end
    end
    @testset "variational optimisers" begin
        include("variational/optimisers.jl")
    end

    Turing.setadbackend(:forwarddiff)
    @testset "stdlib" begin
        include("stdlib/distributions.jl")
        include("stdlib/RandomMeasures.jl")
    end

    @testset "utilities" begin
        # include("utilities/stan-interface.jl")
        include("inference/utilities.jl")
    end
end


text:
# These functions are taken from algorithms.jl in
# https://github.com/JuliaImages/Images.jl/. We don't want to depend on the
# entire Images package just for these functions hence I have just copied and
# pasted them here for now. However, as we are refactoring the Images.jl package
# we might want to move these (and some other functions) into a separate
# package. Once that is completed we can import that package and remove this
# file.

"""
`m = minfinite(A)` calculates the minimum value in `A`, ignoring any values that are not finite (Inf or NaN).
"""
function minfinite(A::AbstractArray{T}) where T
    ret = sentinel_min(T)
    for a in A
        ret = minfinite_scalar(a, ret)
    end
    ret
end
function minfinite(f, A::AbstractArray)
    ret = sentinel_min(typeof(f(first(A))))
    for a in A
        ret = minfinite_scalar(f(a), ret)
    end
    ret
end

"""
`m = maxfinite(A)` calculates the maximum value in `A`, ignoring any values that are not finite (Inf or NaN).
"""
function maxfinite(A::AbstractArray{T}) where T
    ret = sentinel_max(T)
    for a in A
        ret = maxfinite_scalar(a, ret)
    end
    ret
end
function maxfinite(f, A::AbstractArray)
    ret = sentinel_max(typeof(f(first(A))))
    for a in A
        ret = maxfinite_scalar(f(a), ret)
    end
    ret
end

minfinite_scalar(a::T, b::T) where {T} = isfinite(a) ? (b < a ? b : a) : b
maxfinite_scalar(a::T, b::T) where {T} = isfinite(a) ? (b > a ? b : a) : b
minfinite_scalar(a::T, b::T) where {T<:Union{Integer,FixedPoint}} = b < a ? b : a
maxfinite_scalar(a::T, b::T) where {T<:Union{Integer,FixedPoint}} = b > a ? b : a
minfinite_scalar(a, b) = minfinite_scalar(promote(a, b)...)
maxfinite_scalar(a, b) = maxfinite_scalar(promote(a, b)...)

function minfinite_scalar(c1::C, c2::C) where C<:AbstractRGB
    C(minfinite_scalar(c1.r, c2.r),
      minfinite_scalar(c1.g, c2.g),
      minfinite_scalar(c1.b, c2.b))
end
function maxfinite_scalar(c1::C, c2::C) where C<:AbstractRGB
    C(maxfinite_scalar(c1.r, c2.r),
      maxfinite_scalar(c1.g, c2.g),
      maxfinite_scalar(c1.b, c2.b))
end

sentinel_min(::Type{T}) where {T<:Union{Integer,FixedPoint}} = typemax(T)
sentinel_max(::Type{T}) where {T<:Union{Integer,FixedPoint}} = typemin(T)
sentinel_min(::Type{T}) where {T<:AbstractFloat} = convert(T, NaN)
sentinel_max(::Type{T}) where {T<:AbstractFloat} = convert(T, NaN)
sentinel_min(::Type{C}) where {C<:AbstractRGB} = _sentinel_min(C, eltype(C))
_sentinel_min(::Type{C},::Type{T}) where {C<:AbstractRGB,T} = (s = sentinel_min(T); C(s,s,s))
sentinel_max(::Type{C}) where {C<:AbstractRGB} = _sentinel_max(C, eltype(C))
_sentinel_max(::Type{C},::Type{T}) where {C<:AbstractRGB,T} = (s = sentinel_max(T); C(s,s,s))
sentinel_min(::Type{C}) where {C<:AbstractGray} = _sentinel_min(C, eltype(C))
_sentinel_min(::Type{C},::Type{T}) where {C<:AbstractGray,T} = C(sentinel_min(T))
sentinel_max(::Type{C}) where {C<:AbstractGray} = _sentinel_max(C, eltype(C))
_sentinel_max(::Type{C},::Type{T}) where {C<:AbstractGray,T} = C(sentinel_max(T))


text:
include("src/base.jl")
include("src/envs/toytext_mdp.jl")
include("src/policies/value_table_policy.jl")
include("src/solvers/sarsa_policy_solver.jl")

gym_env = ToyTextMDP(gym.make("Taxi-v2"))
solver = SarsaPolicySolver(5000; max_frame_iterations = 100, epsilon = 0.25, epsilon_discount = 0.999, learning_rate = 0.15, discount = 0.99)
policy = solve(solver, gym_env, verbose = true)

# println(policy.value_map)

max_stars_repo_name: dmitrijsc/practical-rl

max_issues_repo_name: dmitrijsc/practical-rl

max_forks_repo_name: dmitrijsc/practical-rl


text:
using Camoforge
using Documenter

makedocs(;
    modules=[Camoforge],
    authors="bradlindblad <bradley.lindblad@gmail.com> and contributors",
    repo="https://github.com/bradlindblad/Camoforge.jl/blob/{commit}{path}#L{line}",
    sitename="Camoforge.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://bradlindblad.github.io/Camoforge.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/bradlindblad/Camoforge.jl",
    devbranch="main"
)


text:
using PyCall

@pyimport serial

s = serial.Serial("/dev/ttyACM1", 250000)
s[:flushOutput]()
s[:flushInput]()
s[:setDTR](0)
sleep(1)
s[:setDTR](1)
sleep(2)
s[:write]("G91\n")
s[:write]("G1 X100 F1000\n")


text:
export Token

struct Token
  id::Integer
  name::String
end

function is_sent_start(token::Token)
  token.id == 1
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
@json struct Commit
    id::String
    short_id::String
    title::String
    author_name::String
    author_email::String
    committer_name::String
    committer_email::String
    created_at::DateTime
    message::String
    committed_date::DateTime
    authored_date::DateTime
    parent_ids::Vector{String}
    last_pipeline::Pipeline
    status::String
end

endpoint(::GitLabAPI, ::typeof(get_commit), owner::AStr, repo::AStr, ref::AStr) =
    Endpoint(:GET, "/projects/$(encode(owner, repo))/repository/commits/$ref")
endpoint(::GitLabAPI, ::typeof(get_commit), project::Integer, ref::AStr) =
    Endpoint(:GET, "/projects/$project/repository/commits/$ref")
into(::GitLabAPI, ::typeof(get_commit)) = Commit


text:
# packing and unpacking of Graph related types

using KernelDensityEstimate
using IncrementalInference

using Base.Test

fg = emptyFactorGraph()

N=100

doors = zeros(1,4)
doors[1,1:4] = [-100.0;0.0;100.0;300.0]
pd = kde!(doors,[3.0])
pd = resample(pd,N);
bws = getBW(pd)[:,1]
bws2 = Array{Float64}(length(bws),1)
bws2[:,1] = bws[:]
doors2 = getPoints(pd);
v1 = addNode!(fg,:x1, ContinuousScalar,N=N)
f1  = addFactor!(fg,[v1],Obsv2( doors2, bws2, [1.0]))

tem = 2.0*randn(1,N)+getVal(v1)+50.0
v2 = addNode!(fg,:x2, ContinuousScalar, N=N)
f2 = addFactor!(fg, [:x1, :x2], Odo(50.0*ones(1,1),2.0*ones(1,1),[1.0]))

@testset "Testing conversion to packed function node data structure and back" begin
  topack = getData(f1)
  dd = convert(PackedFunctionNodeData{PackedObsv2},topack)
  upd = convert(FunctionNodeData{CommonConvWrapper{Obsv2}}, dd)

  @test compare(topack, upd)

  topack = getData(f2) #fg.f[4].attributes["data"]
  dd =  convert(IncrementalInference.PackedFunctionNodeData{PackedOdo},topack)
  upd = convert(IncrementalInference.FunctionNodeData{CommonConvWrapper{Odo}}, dd)

  @test compare(topack, upd)

  # packedv4data = FNDencode(IncrementalInference.PackedFunctionNodeData{PackedOdo}, getData(f2))
  # upv4data = FNDdecode(IncrementalInference.FunctionNodeData{CommonConvWrapper{Odo}}, packedv4data)
  # @test compare(getData(f2), upv4data)
end

@testset "Testing conversion to packed variable node data structure and back" begin
  dat = getData(fg,:x1) # 1
  pd = convert(IncrementalInference.PackedVariableNodeData, dat)
  unpckd = convert(IncrementalInference.VariableNodeData, pd)

  @test IncrementalInference.compare(getData(fg, :x1), unpckd)
  @test getData(fg,:x1) == unpckd

  # packedv1data = VNDencoder(IncrementalInference.PackedVariableNodeData, getData(v1))
  # upv1data = VNDdecoder(IncrementalInference.VariableNodeData, packedv1data)
  # @test compare(getData(v1), upv1data)
end

#


text:
function ergodic_distribution(model, dr, exo_grid:: UnstructuredGrid, endo_grid:: CartesianGrid, dp)
    N_m = n_nodes(exo_grid)
    N_s = n_nodes(endo_grid)
    N = N_m*N_s
    Π = zeros(N_m, N_s, N_m, endo_grid.n...)
    s = nodes(endo_grid)
    p = SVector(model.calibration[:parameters]...)
    a = SVector(endo_grid.min...)
    b = SVector(endo_grid.max...)
    for i_m in 1:n_nodes(exo_grid)
        x = dr(i_m, s)
        m = node(exo_grid, i_m)
        for i_M in 1:n_inodes(dp, i_m)
            M = inode(Point, dp, i_m, i_M)
            w = iweight(dp, i_m, i_M)
            S = transition(model, m, s, x, M, p)
            S = [(S[n]-a)./(b-a) for n=1:length(S)]
            trembling_hand!(view(Π,i_m, :, i_M, :), S, w)
        end
    end
    Π0 = (reshape(Π,N,N))' - Matrix(I, N, N)
    Π0[end,:] .= 1
    B = zeros(N)
    B[end] = 1.0
    μ = Π0\B

    μ = reshape(μ, N_m, N_s)
    return Π, μ
end

function trembling_hand!(A, x, w)
    N,n0 = size(A)
    δ0 = 1.0./(n0-1.0)
    for n in 1:N
        x0 = x[n][1]
        x0 = min.(max.(x0, 0.0),1.0)
        q0 = div.(x0, δ0)
        q0 = max.(0, q0)
        q0 = min.(q0, n0-2)
        λ0 = (x0./δ0-q0) # ∈[0,1[ by construction
        q0_ = round.(Int,q0) + 1
        A[n, q0_]   += (1-λ0)*w
        A[n, q0_+1] += λ0*w
    end
end


text:
using Test

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

for testscen in 1:2
    valdir, scenario, use_permafrost, use_seaice = get_scenario(testscen)
    println(scenario)

    m = climatemodel(scenario, use_permafrost, use_seaice)

    rt_g = m[:ClimateTemperature, :rt_g_globaltemperature]
    rt_g_compare = readpagedata(m, "test/validationdata/$valdir/rt_g_globaltemperature.csv")

    @test rt_g ≈ rt_g_compare rtol = 1e-3
end


text:
module SSSS

using IJulia, Pkg

export notebooks

# function tutorial()
#     cmd = IJulia.find_jupyter_subcommand("notebook")
#     push!(cmd.exec, joinpath(@__DIR__, "..", "notebooks", "tutorial.ipynb"))
#     return IJulia.launch(cmd, joinpath(@__DIR__, "..", "notebooks"), false)
# end

function notebooks()
    return IJulia.notebook(dir=joinpath(@__DIR__, ".."))
end

# REQUIRE = [
#     "GR",
#     "PyCall",
#     "IJulia",
#     "Revise",
#     "Plots",
#     "Latexify",
#     "FFTW",
#     PackageSpec(name="Flux", rev="master"),
#     "BitBasis",
#     "KrylovKit",
#     PackageSpec(url="https://github.com/QuantumBFS/QuAlgorithmZoo.jl.git", rev="master"),
#     PackageSpec(name="IRTools", rev="master"),
#     PackageSpec(name="NNlib", rev="master"),
#     PackageSpec(name="Zygote", rev="master"),
#     PackageSpec(name="Yao", rev="master"),
#     PackageSpec(name="YaoBlocks", rev="master"),
#     PackageSpec(name="YaoArrayRegister", rev="master"),
# ]

# function __init__()
#     for each in REQUIRE
#         if each in keys(Pkg.installed())
#             continue
#         else
#             Pkg.add(each)
#         end
#     end
# end

end # module


text:
# cuStateVec types

## custatevec compute type

function Base.convert(::Type{custatevecComputeType_t}, T::DataType)
    if T == Float16
        return CUSTATEVEC_COMPUTE_16F
    elseif T == Float32
        return CUSTATEVEC_COMPUTE_32F
    elseif T == Float64
        return CUSTATEVEC_COMPUTE_64F
    elseif T == UInt8
        return CUSTATEVEC_COMPUTE_8U
    elseif T == Int8
        return CUSTATEVEC_COMPUTE_8I
    elseif T == UInt32
        return CUSTATEVEC_COMPUTE_32U
    elseif T == Int32
        return CUSTATEVEC_COMPUTE_32I
    else
        throw(ArgumentError("CUSTATEVEC type equivalent for compute type $T does not exist!"))
    end
end

function Base.convert(::Type{Type}, T::custatevecComputeType_t)
    if T == CUSTATEVEC_COMPUTE_16F
        return Float16
    elseif T == CUSTATEVEC_COMPUTE_32F
        return Float32
    elseif T == CUSTATEVEC_COMPUTE_64F
        return Float64
    elseif T == CUSTATEVEC_COMPUTE_8U
        return UInt8
    elseif T == CUSTATEVEC_COMPUTE_32U
        return UInt32
    elseif T == CUSTATEVEC_COMPUTE_8I
        return Int8
    elseif T == CUSTATEVEC_COMPUTE_32I
        return Int32
    else
        throw(ArgumentError("Julia type equivalent for compute type $T does not exist!"))
    end
end

function compute_type(sv_type::DataType, mat_type::DataType)
    if sv_type == ComplexF64 && mat_type == ComplexF64
        return Float64
    elseif sv_type == ComplexF32 && mat_type == ComplexF64
        return Float32
    elseif sv_type == ComplexF32 && mat_type == ComplexF32
        return Float32
    end
end

abstract type Pauli end
struct PauliX <: Pauli end
struct PauliY <: Pauli end
struct PauliZ <: Pauli end
struct PauliI <: Pauli end
CuStateVecPauli(pauli::PauliX) = CUSTATEVEC_PAULI_X
CuStateVecPauli(pauli::PauliY) = CUSTATEVEC_PAULI_Y
CuStateVecPauli(pauli::PauliZ) = CUSTATEVEC_PAULI_Z
CuStateVecPauli(pauli::PauliI) = CUSTATEVEC_PAULI_I

mutable struct CuStateVec{T}
    data::CuVector{T}
    nbits::UInt32
end
function CuStateVec(T, n_qubits::Int)
    data = CUDA.zeros(T, 2^n_qubits)
    # in most cases, taking the hit here for setting one element
    # is cheaper than building the entire thing on the CPU and
    # copying it over
    CUDA.@allowscalar data[1] = one(T)
    CuStateVec{T}(data, n_qubits)
end

Base.eltype(sv::CuStateVec{T}) where T = T

mutable struct CuStateVecSampler
    handle::custatevecSamplerDescriptor_t
    ws_size::Csize_t
    function CuStateVecSampler(sv::CuStateVec, shot_count::UInt32)
        desc_ref   = Ref{custatevecSamplerDescriptor_t}()
        extra_size = Ref{Csize_t}()
        custatevecSampler_create(handle(), pointer(sv.data), eltype(sv), sv.nbits, desc_ref, shot_count, extra_size)
        obj = new(desc_ref[], extra_size[])
        #finalizer(custatevecSampler_destroy, obj) # apparently there is no such method?
        obj
    end
end

Base.unsafe_convert(::Type{custatevecSamplerDescriptor_t}, desc::CuStateVecSampler) = desc.handle


text:
# wrappers of low-level functionality

function cufftGetProperty(property::CUDAapi.libraryPropertyType)
  value_ref = Ref{Cint}()
  cufftGetProperty(property, value_ref)
  value_ref[]
end

version() = VersionNumber(cufftGetProperty(CUDAapi.MAJOR_VERSION),
                          cufftGetProperty(CUDAapi.MINOR_VERSION),
                          cufftGetProperty(CUDAapi.PATCH_LEVEL))


text:
##############################################################################
##
## VcovData stores data you need to compute errors
##
##############################################################################

struct VcovData{T, N}
    regressors::Matrix{Float64}       # X
    crossmatrix::T                    # X'X in the simplest case. Can be Matrix but preferably Factorization
    residuals::Array{Float64, N}      # vector or matrix of residuals (matrix in the case of IV, residuals of Xendo on (Z, Xexo))
    dof_residual::Int
end

VcovDataVector{T} = VcovData{T, 1}
VcovDataMatrix{T} = VcovData{T, 2}

##############################################################################
##
## Any Method to COmpute Standard Errors must Define two Types
##
##############################################################################

abstract type AbstractVcovFormula end
# this type is used for syntax
# the argument vcov = v(x) calls the function VcovFormula(::Type{Val{:v}}, x) which must result into a type that inherits from AbstractVcovFormula
# Moreover, this type must define
# (i) a function allvars that gives all variables needed to compute standard errors and a Vcovmethod
allvars(x::AbstractVcovFormula) = Symbol[]
# (ii) a VcovMethod that transform a dataframe and a object <: ABstractVcovFormula into a AbstractVcovMethod

# this type must implement
# vcov!(::VcovSimpleMethod, x::VcovData)
# shat!(::VcovSimpleMethod, x::VcovData)
abstract type AbstractVcovMethod end
sandwich(H, S::Matrix{Float64}) = H \ S * inv(H)
df_FStat(::AbstractVcovMethod, x::VcovData, hasintercept::Bool) = x.dof_residual - hasintercept


text:
#Write a program that asks the user for their name and greets them with their name.

println("Please Enter your name:")
name = readline()

function greeter(x)
    println("Hello $x, it's nice to meet you.")
end

greeter(name)

max_stars_repo_name: atott/Simple-Programming-Problems-Julia

max_issues_repo_name: atott/Simple-Programming-Problems-Julia

max_forks_repo_name: atott/Simple-Programming-Problems-Julia


text:
# setting working directory
cd("/Volumes/SSD Hans/Github/MacroFall2020/AdvMacro")

#Defining vars
a = 1/3
b = 0.6
k_star = (a*b)^(1/(1-a))
r = 1/b

# Initializing array for loop
m, n = 100,2
K = fill(0.0, (m, n))

#defining policy fn
g = x -> a*b*(x^a)

# loop for capital
k = 0.8*k_star
for i in 1:m
    K[i, 1] = k
    K[i, 2] = g(K[i,1])
    global k = K[i, 2]
end

# defining fn's for the rest of the variables
w = x -> (1-a)x^a
y = x -> x^a
c = (p,q) -> y(p)-q
r = x -> a*(x^(a-1))
int = r.(K[1:100,2])
wage = w.(K[1:100,1])
output = y.(K[1:100,1])
cons = c.(K[1:100,1],K[1:100,2])

# importing packages for plotting
import Pkg; Pkg.add("Plots")
using Plots
gr()

# plotting capital
plot(1:10,fill(k_star,10), label="capital SS")
plot!(1:10,K[1:10,1], label="capital path", linestyle = :dot, legend = :bottomright)
xlabel!("t")
title!("Capital")
savefig("graphs/capital")

# plotting output

plot(1:10,fill(y(k_star),10), label="output SS")
plot!(1:10,output[1:10], label="output path", linestyle = :dot, legend = :bottomright)
xlabel!("t")
title!("Output")
savefig("graphs/output")

# plotting consumption

plot(1:10,fill(c(k_star,k_star),10), label="consumption SS")
plot!(1:10,cons[1:10], label="consumption path", linestyle = :dot, legend = :bottomright)
xlabel!("t")
title!("Consumption")
savefig("graphs/consumption")

# plotting wage

plot(1:10,fill(w(k_star),10), label="wage SS")
plot!(1:10,wage[1:10], label="wage path", linestyle = :dot, legend = :bottomright)
xlabel!("t")
title!("Wage")
savefig("graphs/wage")

# plotting capital rental rate
plot(1:10,fill(r(k_star),10), label="capital rental rate SS")
plot!(1:10,int[1:10], label="rate path", linestyle = :dot, legend = :topright)
xlabel!("t")
title!("Capital rental rate")
savefig("graphs/capitalr")

# Increase in productivity permanently by 5%

k_starp = (a*b*1.5)^(1/(1-a))
#defining policy fn with productivity increases
g1 = x -> a*b*1.05*(x^a)

# Initializing array for loop
K1 = fill(0.0, (m, n))

# loop for capital
k1 = k_starp
for i in 1:m
    K1[i, 1] = k1
    K1[i, 2] = g(K[i,1])
    global k1 = K[i, 2]
end

# defining fn's for the rest of the variables
w1 = x -> (1-a)*(x^a) #does not change
y1 = x -> 1.5*(x^a)
c1 = (p,q) -> y(p)-q #does not change
r1 = x -> a*(1.5*x^(a-1))
int_p = r1.(K1[1:100,2])
wage_p = w1.(K1[1:100,1])
output_p = y1.(K1[1:100,1])
cons_p = c1.(K1[1:100,1],K1[1:100,2])

# plotting capital
plot(1:10,fill(k_star,10), label="capital SS")
plot!(1:10,K1[1:10,1], label="capital path", linestyle = :dot, legend = :bottomright)
xlabel!("t")
title!("Capital (productivity increase)")
savefig("graphs/capitalp")

# plotting output

plot(1:10,fill(y1(k_star),10), label="output SS")
plot!(1:10,output_p[1:10], label="output path", linestyle = :dot, legend = :bottomright)
xlabel!("t")
title!("Output (productivity increase)")
savefig("graphs/outputp")

# plotting consumption

plot(1:10,fill(c1(k_star,k_star),10), label="consumption SS")
plot!(1:10,cons_p[1:10], label="consumption path", linestyle = :dot, legend = :bottomright)
xlabel!("t")
title!("Consumption (productivity increase)")
savefig("graphs/consumptionp")

# plotting wage

plot(1:10,fill(w1(k_star),10), label="wage SS")
plot!(1:10,wage_p[1:10], label="wage path", linestyle = :dot, legend = :bottomright)
xlabel!("t")
title!("Wage (productivity increase)")
savefig("graphs/wagep")

# plotting capital rental rate
plot(1:10,fill(r1(k_star),10), label="capital rental rate SS")
plot!(1:10,int_p[1:10], label="rate path", linestyle = :dot, legend = :topright)
xlabel!("t")
title!("Capital rental rate (productivity increase)")
savefig("graphs/capitalrp")

max_stars_repo_name: py-r-hans/AdvMacro

max_issues_repo_name: py-r-hans/AdvMacro

max_forks_repo_name: py-r-hans/AdvMacro


text:
"""
Data type defining an airfoil

* `polar` should be an array `[a, b, c]` such that `CD = CL ^ 2 * a + CL * b + c`
"""
struct Airfoil{Fg <: Real}
	CL0::Fg
	CLα::Fg
	Cm0::Fg
	Cmα::Fg
	polar::Vector{Fg}
end
Airfoil(CL0::Fg, args...) where Fg = Airfoil{Fg}(CL0, args...)

Base.:+(afl1::Airfoil, afl2::Airfoil) = Airfoil(
    (afl1.CL0 + afl2.CL0) / 2,
    (afl1.CLα + afl2.CLα) / 2,
    (afl1.Cm0 + afl2.Cm0) / 2,
    (afl1.Cmα + afl2.Cmα) / 2,
    (afl1.polar .+ afl2.polar) ./ 2
)
Base.:*(η::Number, afl::Airfoil) = Airfoil(
    η * afl.CL0,
    η * afl.CLα,
    η * afl.Cm0,
    η * afl.Cmα,
    η .* afl.polar
)

"""
Default, inviscid flat plate airfoil
"""
flat_plate = Airfoil{Float64}(0.0, 2 * π, 0.0, 0.0, [0.0, 0.0, 0.0])

max_stars_repo_name: Equipe-AeroDesign-ITA/PropBiology

max_issues_repo_name: Equipe-AeroDesign-ITA/PropBiology

max_forks_repo_name: Equipe-AeroDesign-ITA/PropBiology


text:
#!/usr/bin/env julia

using Luxor

using Test

using Random
Random.seed!(42)

function bezierpathtest(fname)
    Random.seed!(3)
    currentwidth = 1200
    currentheight = 1200
    Drawing(currentwidth, currentheight, fname)
    origin()
    background(setgray(0.8)...)
    tiles = Tiler(currentwidth, currentheight, 4, 4, margin=20)
    for (pos, n) in tiles
        @layer begin
            translate(pos)
            pts = randompointarray(Point(-tiles.tilewidth/2, -tiles.tilewidth/2), Point(tiles.tilewidth/2, tiles.tilewidth/2), rand(3:5))
            setopacity(0.7)
            sethue("black")
            prettypoly(pts, :stroke, close=true)
            randomhue()
            drawbezierpath(makebezierpath(pts, smoothing=rand(0:0.1:1)), :fill)
        end
    end
    @test finish() == true
end

function pathtobezierpath(fname)
    Random.seed!(3)
    currentwidth = 800
    currentheight = 800
    Drawing(currentwidth, currentheight, fname)
    origin()
    background("white")
    st = "julia"
    setopacity(0.7)
    thefontsize = 200
    sethue("red")
    fontsize(thefontsize)
    move(Point(-350, 0))
    textpath(st)
    nbps = pathtobezierpaths()
    for nbp in nbps
        setline(.25)
        sethue("grey50")
        drawbezierpath(nbp, :stroke)
        # exaggerate position of Bezier control points
        for p in nbp
            sethue("magenta")
            exthandle1 = between(p[2], p[1], -3)
            exthandle2 = between(p[3], p[4], -3)
            circle(exthandle1, 1, :fill)
            circle(exthandle2, 1, :fill)
            line(p[1], exthandle1, :stroke)
            line(exthandle2, p[4], :stroke)
            if p[1] != p[4]
                sethue("black")
                circle(p[1], 1, :fill)
                circle(p[4], 1, :fill)
            end
        end
    end
    @test finish() == true
end

fname = "bezierpath.png"
bezierpathtest(fname)

println("...bezier test: output in $(fname)")

fname = "pathtobezierpath.svg"
pathtobezierpath(fname)

println("...finished bezier tests: output in $(fname)")


text:
displays = split.(readlines("day8/displays.txt"), " | ")
displays = [map(split, display) for display in displays]

is_easy_digit(pattern) = length(pattern) in [2, 3, 4, 7]

p1 = sum(sum(is_easy_digit.(outputs)) for (_, outputs) in displays)

@show p1

charsort(string) = join(sort(collect(string)))

contains_pattern(p1, p2) = all(p in p1 for p in p2)

function find_superset(patterns, pattern)
    findfirst(p -> contains_pattern(p, pattern), patterns)
end

function find_subset(patterns, pattern)
    findfirst(p -> contains_pattern(pattern, p), patterns)
end

function decode_display(display)
    signals, outputs = display
    signals = [charsort(signal) for signal in signals]
    signals = sort(signals, by = length)

    pattern = Dict([
        (1, signals[1]),  # only signal of length 2
        (7, signals[2]),  # only signal of length 3
        (4, signals[3]),  # only signal of length 4
        (8, signals[10])  # only signal of length 7
    ])

    signals_5 = signals[4:6]  # three signals (2, 3, 5) of length 5
    signals_6 = signals[7:9]  # three signals (0, 6, 9) of length 6

    pattern[3] = popat!(signals_5, find_superset(signals_5, pattern[7]))
    pattern[9] = popat!(signals_6, find_superset(signals_6, pattern[4]))
    pattern[0] = popat!(signals_6, find_superset(signals_6, pattern[1]))
    pattern[5] = popat!(signals_5, find_subset(signals_5, pattern[9]))
    pattern[2] = pop!(signals_5)
    pattern[6] = pop!(signals_6)

    digit = Dict((v, k) for (k, v) in pattern)
    output_digits = [digit[charsort(output)] for output in outputs]
    parse(Int, join(output_digits))
end

p2 = sum(decode_display.(displays))

@show p2


text:
using Subsequences
using Base.Test

a = "abc"
b = "abc"
seq, a_range, b_range = longest_common_subsequence(a, b)
@test seq == "abc"
@test a_range == b_range == 1:3

a = "-abc-"
b = "_abc_"
seq, a_range, b_range = longest_common_subsequence(a, b)
@test seq == "abc"
@test a_range == b_range == 2:4

a = "abc-abc----"
b = "___ab_c__a_bc"
seq, a_range, b_range = longest_common_subsequence(a, b)
@test seq == "abcabc"
@test a_range == 1:7
@test b_range == 4:13

# Test for type-stability with identical element types
a = [:v, :x, :x, :y, :x, :z, :y, :z]
b = [:x, :z, :x, :y]
@inferred longest_common_subsequence(a, b)
seq, a_range, b_range = longest_common_subsequence(a, b)
@test seq == [:x, :x, :y]
@test typeof(seq) == Vector{Symbol}
@test a_range == 2:4
@test b_range == 1:4

# Test type-stability for different element types
a = [1.0, 2.0, 3.0, 4.0]
b = [1, 3]
@inferred longest_common_subsequence(a, b)
seq, a_range, b_range = longest_common_subsequence(a, b)
@test typeof(seq) == Vector{promote_type(Float64, Int)}
@test seq == [1.0, 3.0]

# -------

a = "abc"
b = "abc"
seq, a_range, b_range = longest_contiguous_subsequence(a, b)
@test a_range == b_range == 1:3

a = "-abc-"
b = "_abc_"
seq, a_range, b_range = longest_contiguous_subsequence(a, b)
@test seq == "abc"
@test a_range == b_range == 2:4

a = "abc-abc----"
b = "___ab_c__abc"
seq, a_range, b_range = longest_contiguous_subsequence(a, b)
@test seq == "abc"
@test a_range == 1:3
@test b_range == 10:12


text:
module GridPotentialFlow

using Reexport
using StructArrays
using DocStringExtensions
using RecipesBase

@reexport using CartesianGrids
@reexport using RigidBodyTools

include("suctionparameter.jl")
include("bodies.jl")

include("solver/righthandside.jl")
include("solver/solution.jl")
include("solver/systems.jl")

include("vortex.jl")
include("vortexmodel.jl")

end


text:
# Hola Mundo Julia

#=
My first code in Julia
=#
println("Hello World")

# Variable

s = ("Hello World")

println(s)

# Data types

anInteger = 763543254
aFloat    = .2345
aBoolean  = true
aString   = "A string"
aChar     = '+'

println(typeof(anInteger))
println(typeof(aFloat))
println(typeof(aBoolean))
println(typeof(aString))
println(typeof(aChar))

# Casting

i2c = Char(120)
f2i = UInt8(trunc(3.14))
s2f = parse(Float64 , "1")
s2i = parse(Int8    , "1")

println("i2c: ", i2c)
println("f21: ", f2i)
println("s2f: ", s2f)
println("s2i: ", s2i)

# String

t = "This is a string"

println(s)
println("Size : ", length(s)) # Size of string
println("Index# : ", t[1] )   # The 1st index is at 1
println("Last : ", t[end] )   # Get last character
println("Range : ", t[1:4])   # Get a range

println("Index of x : ", findfirst(isequal('i'), t) ) # Find index
println("Presence check : ", occursin("is", t) ) # Find a substring
println("And : ", t * " ,cool" ) # And
println("Concatenation : ", string(t, " ,cool") ) # Concatenation


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using UnitCommitment, Test, LinearAlgebra

@testset "Screening" begin
    @testset "Violation filter" begin
        instance = UnitCommitment.read_benchmark("test/case14")
        filter = ViolationFilter(max_per_line=1, max_total=2)

        offer(filter, Violation(time=1,
                                monitored_line=instance.lines[1],
                                outage_line=nothing,
                                amount=100.))

        offer(filter, Violation(time=1,
                                monitored_line=instance.lines[1],
                                outage_line=instance.lines[1],
                                amount=300.))

        offer(filter, Violation(time=1,
                                monitored_line=instance.lines[1],
                                outage_line=instance.lines[5],
                                amount=500.))

        offer(filter, Violation(time=1,
                                monitored_line=instance.lines[1],
                                outage_line=instance.lines[4],
                                amount=400.))

        offer(filter, Violation(time=1,
                                monitored_line=instance.lines[2],
                                outage_line=instance.lines[1],
                                amount=200.))

        offer(filter, Violation(time=1,
                                monitored_line=instance.lines[2],
                                outage_line=instance.lines[8],
                                amount=100.))

        actual = query(filter)
        expected = [Violation(time=1,
                              monitored_line=instance.lines[2],
                              outage_line=instance.lines[1],
                              amount=200.),
                    Violation(time=1,
                              monitored_line=instance.lines[1],
                              outage_line=instance.lines[5],
                              amount=500.)]
        @test actual == expected
    end

    @testset "find_violations" begin
        instance = UnitCommitment.read_benchmark("test/case14")
        for line in instance.lines, t in 1:instance.time
            line.normal_flow_limit[t] = 1.0
            line.emergency_flow_limit[t] = 1.0
        end
        isf = UnitCommitment.injection_shift_factors(lines=instance.lines,
                                                     buses=instance.buses)
        lodf = UnitCommitment.line_outage_factors(lines=instance.lines,
                                                  buses=instance.buses,
                                                  isf=isf)
        inj = [1000.0 for b in 1:13, t in 1:instance.time]
        overflow = [0.0 for l in instance.lines, t in 1:instance.time]
        violations = UnitCommitment.find_violations(instance=instance,
                                                    net_injections=inj,
                                                    overflow=overflow,
                                                    isf=isf,
                                                    lodf=lodf)

        @test length(violations) == 20
    end
end


text:
using PkgTemplates

DirPath=joinpath(ENV["HOME"],".julia/dev")
PackageName="githubactions4julia.jl"
user="wakakusa"
cd(DirPath)

t = Template(;
                  :julia=>v"1",
                  :user=>user,
                  :authors=>[user],
                  :dir=>DirPath,

                  plugins=[
                      License(; name="MIT"),
                      Git(; manifest=false, ssh=true),
                      GitHubActions(; x86=true),
                      Documenter{GitHubActions}(),
                  ],
              )

t(PackageName)

using DocumenterTools
using githubactions4julia

repo=string("git@github.com:",user,"/",PackageName,".git")
# 鍵が２つ出力されるので見落とさないこと
DocumenterTools.genkeys(;user=user, repo=repo)

## ドキュメントローカル作成
cd $DirPath/$PackageName
julia --project=docs -e 'using Pkg;Pkg.develop(PackageSpec(path=pwd()));Pkg.instantiate();include("docs/make.jl");'


text:
# ---
# title: 1423. Maximum Points You Can Obtain from Cards
# id: problem1423
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array, Dynamic Programming, Sliding Window
# link: <https://leetcode.com/problems/maximum-points-you-can-obtain-from-cards/description/>
# hidden: true
# ---
#
# There are several cards  **arranged in a row** , and each card has an
# associated number of points The points are given in the integer array
# `cardPoints`.
#
# In one step, you can take one card from the beginning or from the end of the
# row. You have to take exactly `k` cards.
#
# Your score is the sum of the points of the cards you have taken.
#
# Given the integer array `cardPoints` and the integer `k`, return the _maximum
# score_ you can obtain.
#
#
#
# **Example 1:**
#
#
#
#     Input: cardPoints = [1,2,3,4,5,6,1], k = 3
#     Output: 12
#     Explanation: After the first step, your score will always be 1. However, choosing the rightmost card first will maximize your total score. The optimal strategy is to take the three cards on the right, giving a final score of 1 + 6 + 5 = 12.
#
#
# **Example 2:**
#
#
#
#     Input: cardPoints = [2,2,2], k = 2
#     Output: 4
#     Explanation: Regardless of which two cards you take, your score will always be 4.
#
#
# **Example 3:**
#
#
#
#     Input: cardPoints = [9,7,7,9,7,7,9], k = 7
#     Output: 55
#     Explanation: You have to take all the cards. Your score is the sum of points of all cards.
#
#
# **Example 4:**
#
#
#
#     Input: cardPoints = [1,1000,1], k = 1
#     Output: 1
#     Explanation: You cannot take the card in the middle. Your best score is 1.
#
#
# **Example 5:**
#
#
#
#     Input: cardPoints = [1,79,80,1,1,1,200,1], k = 3
#     Output: 202
#
#
#
#
# **Constraints:**
#
#   * `1 <= cardPoints.length <= 10^5`
#   * `1 <= cardPoints[i] <= 10^4`
#   * `1 <= k <= cardPoints.length`
#
#
## @lc code=start
using LeetCode

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


text:
# by lazarusA # HIDE
# using GLMakie # HIDE
using GLMakie
GLMakie.activate!() # HIDE
let
    x = -2:0.005:2
    y = -2:0.005:2
    f(z) = (z^2 + 1)/(z^2 - 1)
    fvals = [f(u+1im*v) for u in x, v in y]
    fvalues = abs.(fvals)
    fargs = angle.(fvals)
    indxCut = fvalues .> 3
    fvalues[indxCut] .= 3.01
    fig, ax, pltobj = surface(x, y, fvalues, color = fargs,
        colormap = :diverging_rainbow_bgymr_45_85_c67_n256,
        colorrange = (-π,π), backlight = 1f0, highclip = :black)
    fig
    save(joinpath(@__DIR__, "output", "surfaceComplexF.png"), fig, px_per_unit = 2.0) # HIDE
end
using Pkg # HIDE
Pkg.status(["GLMakie"]) # HIDE

max_stars_repo_name: MisterBiggs/BeautifulMakie

max_issues_repo_name: MisterBiggs/BeautifulMakie

max_forks_repo_name: MisterBiggs/BeautifulMakie


text:
using Revise

if  pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using Minvar, Random, Revise

#Revise.revise(Minvar)

const WORKERS=2
const CLEAN_WORKERS=false

#set up for parallel operations
initializepar(WORKERS, reset=CLEAN_WORKERS)

Random.seed!(11)
priortype= :uninformative
nburn = 10_000
prefix = priortype == :informative ? "inf_$nburn" : "uninf_$nburn"

#testportfolio(iter=20)
#testcov(iter=2_000)
#testtrisectedportfolio(iter=200, increment=199, displayiter = 0)
#print(mathematica2juliaEq())
#testω(iter=1000)
#testposteriors()
@time E = estimateminvariance(nburn=nburn, nsamples=100_000, parallel=(nburn≠0), recordweights=false,
  saveestimate=true, priortype=priortype, prefix=prefix)
@time makeconvergencegraphs(E, prefix=prefix)

if CLEAN_WORKERS
  cleanup()
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using MPI
using Trilinos
using Test

if isinteractive()
  using Plots
end

const USE_LOCAL_INDICES = !isempty(ARGS) && ARGS[1] == "local"

@static if USE_LOCAL_INDICES
  const IdxT = Int32
else
  const IdxT = Int
end

function replace_values(A, gid, indices, values)
  @static if USE_LOCAL_INDICES
    return Tpetra.replaceLocalValues(A, gid, indices, values)
  else
    return Tpetra.replaceGlobalValues(A, gid, indices, values)
  end
end

function global_element(rowmap,i)
  @static if USE_LOCAL_INDICES
    return Int(i)
  else
    return Tpetra.getGlobalElement(rowmap,i)
  end
end

function local_element(rowmap,gid)
  @static if USE_LOCAL_INDICES
    return gid
  else
    return Tpetra.getLocalElement(rowmap,gid)
  end
end

function is_node_global_element(rowmap,gid)
  @static if USE_LOCAL_INDICES
    return true
  else
    return Tpetra.isNodeGlobalElement(rowmap,gid)
  end
end

struct CartesianGrid
  nx::Int # Number of points in the X direction
  ny::Int # Number of points in the Y direction
  h::Float64 # Spacing
end

nnodes(g::CartesianGrid) = g.nx*g.ny

function xyindices(g::CartesianGrid, i)
  ix = i % g.nx
  iy = (i-ix) ÷ g.nx
  return (ix,iy)
end

function coordinates(g::CartesianGrid, i)
  (ix,iy) = xyindices(g,i)
  x_center = g.h*(g.nx-1)/2
  y_center = g.h*(g.ny-1)/2
  return (ix*g.h - x_center, iy*g.h - y_center)
end

"""
Matrix global indices for the 2D laplace stencil (0-based!)
"""
function laplace2d_indices!(inds_array,i,g::CartesianGrid)
  (ix,iy) = xyindices(g,i)

  n_inds = 1
  inds_array[n_inds] = i
  if iy != g.ny-1
      n_inds += 1
      inds_array[n_inds] = i+g.nx
  end
  if ix != 0
      n_inds += 1
      inds_array[n_inds] = i-1
  end
  if iy != 0
      n_inds += 1
      inds_array[n_inds] = i-g.nx
  end
  if ix != g.nx-1
      n_inds += 1
      inds_array[n_inds] = i+1
  end
  return n_inds
end

function set_source_term!(b, g::CartesianGrid)
  rowmap = Tpetra.getMap(b)
  b_view = Tpetra.device_view(b)
  n_my_elms = Tpetra.getNodeNumElements(rowmap)
  @assert n_my_elms == length(LinearIndices(b_view))
  for i in LinearIndices(b_view)
    gid = global_element(rowmap,i)
    (x,y) = coordinates(g,gid)
    b_view[i] = 2*g.h^2*((1-x^2)+(1-y^2))
  end
end

function graph_laplace2d!(crsgraph, g::CartesianGrid)
  rowmap = Tpetra.getRowMap(crsgraph)
  n_my_elms = Tpetra.getNodeNumElements(rowmap)

  # storage for the per-row indices
  row_indices = [0,0,0,0,0]

  for i in 0:n_my_elms-1
    global_row = global_element(rowmap,i)
    row_n_elems = laplace2d_indices!(row_indices, global_row, g)
    Tpetra.insertGlobalIndices(crsgraph, global_row, Teuchos.ArrayView(row_indices,row_n_elems))
  end
end

function fill_laplace2d!(A, g::CartesianGrid)
  rowmap = Tpetra.getRowMap(A)
  n_my_elms = Tpetra.getNodeNumElements(rowmap)

  # storage for the per-row values
  row_indices = IdxT[0,0,0,0,0]
  row_values = [4.0,-1.0,-1.0,-1.0,-1.0]

  for i in 0:n_my_elms-1
    global_row = global_element(rowmap,i)
    row_n_elems = laplace2d_indices!(row_indices, global_row, g)
    row_values[1] = 4.0 - (5-row_n_elems)
    replace_values(A, global_row, Teuchos.ArrayView(row_indices,row_n_elems), Teuchos.ArrayView(row_values,row_n_elems))
  end
end

# Ensure a dirichlet condition (1 at the diagonal) for the boundaries
function set_dirichlet!(A, b, g::CartesianGrid)
  rowmap = Tpetra.getRowMap(A)
  row_indices = IdxT[0,0,0,0,0]
  row_values = [1.0,0.0,0.0,0.0,0.0]

  b_view = Tpetra.device_view(b)

  # Collect the boundary nodes
  boundary_gids = zeros(Int,2*(g.nx+g.ny))
  gid_idx = 1
  # left and right
  for iy in 0:(g.ny-1)
    boundary_gids[gid_idx] = iy*g.nx
    boundary_gids[gid_idx+1] = boundary_gids[gid_idx] + g.nx-1
    gid_idx += 2
  end

  # top and bottom
  for ix in 0:(g.nx-1)
    boundary_gids[gid_idx] = ix
    boundary_gids[gid_idx+1] = ix + (g.ny-1)*g.nx
    gid_idx += 2
  end

  # apply
  for gid in boundary_gids
    (x,y) = coordinates(g,gid)
    if is_node_global_element(rowmap, gid)
      row_n_elems = laplace2d_indices!(row_indices, gid, g)
      replace_values(A, gid, Teuchos.ArrayView(row_indices,row_n_elems), Teuchos.ArrayView(row_values,row_n_elems))
      b_view[local_element(rowmap, gid)] = (1-x^2)*(1-y^2)
    end
  end
end

function check_solution(sol, g::CartesianGrid)
  solmap = Tpetra.getMap(sol)
  solview = Tpetra.device_view(sol)
  result = 0
  for i in LinearIndices(solview)
    gid = global_element(solmap, i)
    (x,y) = coordinates(g,gid)
    result +=  abs(solview[i] - (1-x^2)*(1-y^2)) > 1e-10
  end
  return result == 0
end

"""
Assemble the 2D Poisson problem on a structured grid
"""
function laplace2d(comm, g::CartesianGrid)
  # Construct map
  rowmap = Tpetra.Map(nnodes(g), 0, comm)

  # Construct graph (i.e. the sparsity pattern)
  matrix_graph = Tpetra.CrsGraph(rowmap, 0)
  println("-----------------------------------")
  println("Graph construction time:")
  @time graph_laplace2d!(matrix_graph, g)
  Tpetra.fillComplete(matrix_graph)

  # Construct the vectors
  b = Tpetra.Vector(Tpetra.getRangeMap(matrix_graph))
  println("Source term time:")
  @time set_source_term!(b,g)

  # Matrix construction
  A = Tpetra.CrsMatrix(matrix_graph)
  Tpetra.resumeFill(A)
  @test Tpetra.isLocallyIndexed(A)
  println("Matrix fill time:")
  @time fill_laplace2d!(A,g)
  println("Dirichlet time:")
  @time set_dirichlet!(A,b,g)
  Tpetra.fillComplete(A)
  @test Tpetra.isLocallyIndexed(A)

  # This prints the matrix if uncommented
  # Tpetra.describe(A, Teuchos.VERB_EXTREME)

  params = Trilinos.default_parameters()
  maxiter = 1000
  solver_params = params["Linear Solver Types"]["Belos"]["Solver Types"]["BLOCK GMRES"]
  solver_params["Convergence Tolerance"] = 1e-12
  solver_params["Verbosity"] = Belos.StatusTestDetails + Belos.FinalSummary + Belos.TimingDetails
  solver_params["Num Blocks"] = Int32(maxiter)
  solver_params["Maximum Iterations"] = Int32(maxiter)

  params["Preconditioner Type"] = "MueLu"

  solver = TpetraSolver(A, params)

  return (solver,b)
end

function solve_laplace2d(comm)

  nx = 1001
  h = 2/(nx-1)
  grid = CartesianGrid(nx,nx,h)

  (lows,b) = laplace2d(comm, grid)
  (lows,b) = laplace2d(comm, grid)
  (lows,b) = laplace2d(comm, grid)

  # Compute the solution
  sol = lows \ b

  @test check_solution(sol, grid)
  println("-----------------------------------")
  println("test time:")
  @time check_solution(sol, grid)
  @time check_solution(sol, grid)
  @time check_solution(sol, grid)

  x = linspace(-1.0,1.0,grid.nx)
  y = linspace(-1.0,1.0,grid.ny)*grid.h*(grid.ny-1)/2
  return (x,y,Tpetra.device_view(sol),grid)
end

# MPI setup
if !MPI.Initialized()
  MPI.Init()
  MPI.finalize_atexit()
end
comm = Teuchos.MpiComm(MPI.CComm(MPI.COMM_WORLD))

my_rank = Teuchos.getRank(comm)
if my_rank != 0
  redirect_stdout(open("/dev/null", "w"))
  redirect_stderr(open("/dev/null", "w"))
end

nb_procs = Teuchos.getSize(comm)
if USE_LOCAL_INDICES && nb_procs > 1
  error("Local indices may only be used on single core execution")
end

(x,y,f) = solve_laplace2d(comm)

if isinteractive()
  plot(x, y, f, aspect_ratio=1, seriestype=:heatmap)
end


text:
#Stationary iterative methods
#Templates, section 2.2
export jacobi, jacobi!, gauss_seidel, gauss_seidel!, sor, sor!, ssor, ssor!

jacobi(A::AbstractMatrix, b;
       tol=size(A,2)^3*eps(typeof(real(b[1]))), maxiter=size(A,2)^2) =
    jacobi!(zerox(A, b), A, b; tol=tol, maxiter=maxiter)

function jacobi!(x, A::AbstractMatrix, b;
        tol=size(A,2)^3*eps(typeof(real(b[1]))), maxiter=size(A,2)^2)
    history = ConvergenceHistory()
    history[:tol] = tol
    reserve!(history,:resnorm, maxiter)
    jacobi_method!(history, x, A, b; tol=tol, maxiter=maxiter)
    x, history
end

function jacobi_method!(log::ConvergenceHistory, x, A::AbstractMatrix, b;
        tol=size(A,2)^3*eps(typeof(real(b[1]))), maxiter=size(A,2)^2)
    iter=0
    n = size(A,2)
    xold = copy(x)
    z = zero(Amultype(A, x))
    tol = tol * norm(b)
	for iter=1:maxiter
        nextiter!(log,mvps=1)
		for i=1:n
			xi = z
			for j=[1:i-1;i+1:n]
				xi += A[i,j]*xold[j]
			end
			A[i,i]==0 && throw(SingularError())
			x[i]=(b[i]-xi)/A[i,i]
		end
		#check convergence
        resnorm = norm(A*x-b)
        push!(log,:resnorm,resnorm)
		resnorm < tol && (setconv(log, resnorm>=0); break)
		copy!(xold, x)
	end
    shrink!(log)
	x
end

gauss_seidel(A::AbstractMatrix, b;
        tol=size(A,2)^3*eps(typeof(real(b[1]))), maxiter=size(A,2)^2) =
    gauss_seidel!(zerox(A, b), A, b; tol=tol, maxiter=maxiter)

function gauss_seidel!(x, A::AbstractMatrix, b;
        tol=size(A,2)^3*eps(typeof(real(b[1]))), maxiter=size(A,2)^2)
    history = ConvergenceHistory()
    history[:tol] = tol
    reserve!(history,:resnorm, maxiter)
    gauss_seidel_method!(history, x, A, b; tol=tol, maxiter=maxiter)
    x, history
end

function gauss_seidel_method!(log::ConvergenceHistory, x, A::AbstractMatrix, b;
        tol=size(A,2)^3*eps(typeof(real(b[1]))), maxiter=size(A,2)^2)
    iter=0
    n = size(A,2)
    xold = copy(x)
    z = zero(Amultype(A, x))
    tol = tol * norm(b)
	for iter=1:maxiter
        nextiter!(log,mvps=1)
		for i=1:n
			σ=z
			for j=1:i-1
				σ+=A[i,j]*x[j]
			end
			for j=i+1:n
				σ+=A[i,j]*xold[j]
			end
			A[i,i]==0 && throw(SingularError())
			x[i]=(b[i]-σ)/A[i,i]
		end
		#check convergence
		resnorm = norm(A*x-b)
        push!(log,:resnorm,resnorm)
        resnorm < tol && (setconv(log, resnorm>=0); break)
		copy!(xold, x)
	end
    shrink!(log)
	x
end

#Successive overrelaxation
sor(A::AbstractMatrix, b, ω::Real;
    tol=size(A,2)^3*eps(typeof(real(b[1]))), maxiter=size(A,2)^2) =
    sor!(zerox(A, b), A, b, ω; tol=tol, maxiter=maxiter)

function sor!(x, A::AbstractMatrix, b, ω::Real;
        tol=size(A,2)^3*eps(typeof(real(b[1]))), maxiter=size(A,2)^2)
    history = ConvergenceHistory()
    history[:tol] = tol
    reserve!(history,:resnorm, maxiter)
    sor_method!(history, x, A, b, ω; tol=tol, maxiter=maxiter)
    x, history
end

function sor_method!(log::ConvergenceHistory, x, A::AbstractMatrix, b, ω::Real;
        tol=size(A,2)^3*eps(typeof(real(b[1]))), maxiter=size(A,2)^2)
	0 < ω < 2 || warn("ω = $ω lies outside the range 0<ω<2 which is required for convergence")

    iter=0
	n = size(A,2)
    xold = copy(x)
    z = zero(Amultype(A, x))
    tol = tol * norm(b)
	for iter=1:maxiter
        nextiter!(log,mvps=1)
		for i=1:n
			σ=z
			for j=1:i-1
				σ+=A[i,j]*x[j]
			end
			for j=i+1:n
				σ+=A[i,j]*xold[j]
			end
			A[i,i]==0 && throw(SingularError())
			σ=(b[i]-σ)/A[i,i]
			x[i]=xold[i]+ω*(σ-xold[i])
		end
		#check convergence
		resnorm = norm(A*x-b)
        push!(log,:resnorm,resnorm)
        resnorm < tol && (setconv(log, resnorm>=0); break)
		copy!(xold, x)
	end
    shrink!(log)
	x
end

#Symmetric successive overrelaxation
#A must be symmetric
ssor(A::AbstractMatrix, b, ω::Real;
     tol=size(A,2)^3*eps(typeof(real(b[1]))), maxiter=size(A,2)) =
    ssor!(zerox(A, b), A, b, ω; tol=tol, maxiter=maxiter)

function ssor!(x, A::AbstractMatrix, b, ω::Real;
        tol=size(A,2)^3*eps(typeof(real(b[1]))), maxiter=size(A,2))
    history = ConvergenceHistory()
    history[:tol] = tol
    reserve!(history,:resnorm, maxiter)
    ssor_method!(history, x, A, b, ω; tol=tol, maxiter=maxiter)
    x, history
end

function ssor_method!(log::ConvergenceHistory, x, A::AbstractMatrix, b, ω::Real;
        tol=size(A,2)^3*eps(typeof(real(b[1]))), maxiter=size(A,2))
	0 < ω < 2 || warn("ω = $ω lies outside the range 0<ω<2 which is required for convergence")

    iter=0
	n = size(A,2)
    xold = copy(x)
    z = zero(Amultype(A, x))
    tol = tol * norm(b)
	for iter=1:maxiter
        nextiter!(log,mvps=1)
		for i=1:n #Do a SOR sweep
			σ=z
			for j=1:i-1
				σ+=A[i,j]*x[j]
			end
			for j=i+1:n
				σ+=A[i,j]*xold[j]
			end
			A[i,i]==0 && throw(SingularError())
			σ=(b[i]-σ)/A[i,i]
			x[i]=xold[i]+ω*(σ-xold[i])
		end
		copy!(xold, x)
		for i=n:-1:1 #Do a backward SOR sweep
			σ=z
			for j=1:i-1
				σ+=A[i,j]*xold[j]
			end
			for j=i+1:n
				σ+=A[i,j]*x[j]
			end
			A[i,i]==0 && throw(SingularError())
			σ=(b[i]-σ)/A[i,i] #This line is missing in the Templates reference
			x[i]=xold[i]+ω*(σ-xold[i])
		end
		#check convergence
		resnorm = norm(A*x-b)
        push!(log,:resnorm,resnorm)
        resnorm < tol && (setconv(log, resnorm>=0); break)
		copy!(xold, x)
	end
    shrink!(log)
	x
end


text:
mutable struct Position
    x::Int
    y::Int
end

struct Size
    width::Int
    height::Int
end

# ----
@enum Weapon Laser Missile

# A thing should have a position and size
abstract type Thing end
position(t::Thing) = t.position
size(t::Thing) = t.size

struct Spaceship <: Thing
    position::Position
    size::Size
    weapon::Weapon
end

Base.show(io::IO, s::Spaceship) =
    print(io, "Spaceship (", s.position.x, ",", s.position.y, ") " ,
        s.size.width, "x", s.size.height,
        "/", s.weapon)

struct Asteroid <: Thing
    position::Position
    size::Size
end

Base.show(io::IO, s::Asteroid) =
    print(io, "Asteroid (", s.position.x, ",", s.position.y,
        ") ", s.size.width, "x", s.size.height)

# Test
s1 = Spaceship(Position(0,0), Size(30,5), Missile)
s2 = Spaceship(Position(10,0), Size(30,5), Laser)
a1 = Asteroid(Position(20,0), Size(20,20))
a2 = Asteroid(Position(0,20), Size(20,20))

# ----------------------------------------------------------------------
# parametric methods

# explode an array of objects
function explode(things::AbstractVector{Any})
    for t in things
        println("Exploding ", t)
    end
end

# explode an array of objects (parametric version)
function explode(things::AbstractVector{T}) where {T}
    for t in things
        println("Exploding ", t)
    end
end

explode([a1, a2])
#= REPL
julia> explode([a1, a2])
Exploding Asteroid (20,0) 20x20
Exploding Asteroid (0,20) 20x20
=#

explode([:building, :hill])
#= REPL
julia> explode([:building, :hill])
Exploding building
Exploding hill

julia> Vector{Asteroid} <: AbstractVector{Asteroid}
true
=#

# Same function with a more narrow type
function explode(things::AbstractVector{T}) where {T <: Thing}
    for t in things
        println("Exploding thing => ", t)
    end
end

explode([a1, a2])
#= REPL
julia> explode([a1, a2])
Exploding thing => Asteroid (20,0) 20x20
Exploding thing => Asteroid (0,20) 20x20
=#

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

# specifying abstract/concrete types in method signature
function tow(A::Spaceship, B::Thing)
    "tow 1"
end

methods(tow)
#=
julia> methods(tow)
# 1 method for generic function "tow":
[1] tow(A::Spaceship, B::Thing) in Main at REPL[17]:8
=#

# equivalent of parametric type
function tow(A::Spaceship, B::T) where {T <: Thing}
    "tow 2"
end

methods(tow)
#= REPL
julia> methods(tow)
# 1 method for generic function "tow":
[1] tow(A::Spaceship, B::T) where T<:Thing in Main at REPL[20]:2
=#

# --- ensure type consistency

function group_anything(A::Thing, B::Thing)
    println("Grouped ", A, " and ", B)
end

group_anything(s1, s2)
group_anything(a1, a2)
group_anything(s1, a1)
group_anything(a1, s1)

#= REPL
julia> group_anything(s1, s2)
Grouped Spaceship((0,0), 30 x 5, Missile::Weapon = 1) and Spaceship((10,0), 30 x 5, Laser::Weapon = 0)

julia> group_anything(a1, a2)
Grouped Asteroid((20,0), 20 x 20) and Asteroid((0,20), 20 x 20)

julia> group_anything(s1, a1)
Grouped Spaceship((0,0), 30 x 5, Missile::Weapon = 1) and Asteroid((20,0), 20 x 20)

julia> group_anything(a1, s1)
Grouped Asteroid((20,0), 20 x 20) and Spaceship((0,0), 30 x 5, Missile::Weapon = 1)
=#

function group_same_things(A::T, B::T) where {T <: Thing}
    println("Grouped ", A, " and ", B)
end

group_same_things(s1, s2)
group_same_things(a1, a2)
group_same_things(s1, a1)
group_same_things(a1, s1)

#=
julia> group_same_things(s1, s2)
Grouped Spaceship((0,0), 30 x 5, Missile::Weapon = 1) and Spaceship((10,0), 30 x 5, Laser::Weapon = 0)

julia> group_same_things(a1, a2)
Grouped Asteroid((20,0), 20 x 20) and Asteroid((0,20), 20 x 20)

julia> group_same_things(s1, a1)
ERROR: MethodError: no method matching group_same_things(::Spaceship, ::Asteroid)
Closest candidates are:
  group_same_things(::T<:Thing, ::T<:Thing) where T<:Thing at REPL[88]:2

julia> group_same_things(a1, s1)
ERROR: MethodError: no method matching group_same_things(::Asteroid, ::Spaceship)
Closest candidates are:
  group_same_things(::T<:Thing, ::T<:Thing) where T<:Thing at REPL[88]:2
=#

# ----------------------------------------------------------------------
# Extracting type info

eltype(things::AbstractVector{T}) where {T <: Thing} = T

eltype([s1, s2])
eltype([a1, a2])
eltype([s1, s2, a1, a2])

typeof([s1, s2])
typeof([a1, a2])
typeof([s1, s2, a1, a2])

# julia> eltype([s1, s2])
# Spaceship

# julia> eltype([a1, a2])
# Asteroid

# julia> eltype([s1, s2, a1, a2])
# Thing

# julia> typeof([s1, s2])
# Array{Spaceship,1}

# julia> typeof([a1, a2])
# Array{Asteroid,1}

# julia> typeof([s1, s2, a1, a2])
# Array{Thing,1}


text:
using AbstractDAQ
using Test

@testset "AbstractDAQ.jl" begin
    # Write your tests here.
end


text:
"""
    PubMedMiner

Utilities to mine results from a PubMed/Medline search
Authors: Isabel Restrepo, Mary McGrath
BCBI - Brown University
"""
module PubMedMiner

using MySQL
using DataFrames

export get_semantic_occurrences_df,
    mesh_stats
include("Stats.jl")
include("MySQL_DB.jl")

# include("API_DF.jl")
# include("Plotting.jl")

end #PubMedMiner


text:
#=
To do :
- verify support -Inf Inf for truncated variable bounds (marche pas avec uniform)
=#

using Distributions, Test
using FastElitismGaussianEDA
using FastElitismGaussianEDA: AbstractVariable, rand, fit_model, update_variable!

@testset "Variable constructors and type consistency" begin
    v = ContinuousVariable(1.0, 2.0)
    @test isa(rand(v), Float64)
    v = ContinuousVariable(1f0, 2f0)
    @test isa(rand(v), Float32)
    v = ContinuousVariable(1f0, 2.0)
    @test isa(rand(v), Float64)
    v = ContinuousVariable(1, 2)
    @test isa(rand(v), Float64)
    v = ContinuousVariable(Float16(1), Float16(2))
    @test isa(rand(v), Float16)

    v = ContinuousBoundedVariable(1.0, 2.0)
    @test isa(rand(v), Float64)
    v = ContinuousBoundedVariable(1f0, 2f0)
    @test isa(rand(v), Float32)
    v = ContinuousBoundedVariable(1f0, 2.0)
    @test isa(rand(v), Float64)
    v = ContinuousBoundedVariable(1, 2)
    @test isa(rand(v), Float64)
    v = ContinuousBoundedVariable(Float16(1), Float16(2))
    @test isa(rand(v), Float16)

    rng = 1.0:0.01:2.0
    v = DiscreteVariable(rng)
    @test isa(rand(v), eltype(rng))
    rng = 1f0:1f-1:2f0
    v = DiscreteVariable(rng)
    @test isa(rand(v), eltype(rng))
    rng = 1:1:5
    v = DiscreteVariable(rng)
    @test isa(rand(v), eltype(rng))

    rng = 1.0:0.01:2.0
    v = DiscreteBoundedVariable(rng)
    @test isa(rand(v), eltype(rng))
    rng = 1f0:1f-1:2f0
    v = DiscreteBoundedVariable(rng)
    @test isa(rand(v), eltype(rng))
    rng = 1:1:5
    v = DiscreteBoundedVariable(rng)
    @test isa(rand(v), eltype(rng))
end

@testset "VariableSet" begin
    v1 = ContinuousVariable(-10, 10)
    v2 = ContinuousBoundedVariable(-10, 10)
    v3 = DiscreteVariable(-10:0.1:10)
    v4 = DiscreteBoundedVariable(-10:0.1:10)

    s = VariableSet(v1, v2, v3, v4)
    @test isa(s.vars, AbstractArray)
    @test first(s) === v1
    @test last(s) === v4
    @test s[2] === v2
    tmp = DiscreteBoundedVariable(-10:0.1:10)
    s[3] = tmp
    @test s[3] === tmp
    @test length(values.(s)) == 4

    s = VariableSet(:x => ContinuousVariable(-10, 10),
                    :y => DiscreteVariable(-10:0.1:10),
                    :z => DiscreteBoundedVariable(-10:0.1:10))
    @test isa(s[1], ContinuousVariable)
    @test isa(s[2], DiscreteVariable)
    @test isa(s[3], DiscreteBoundedVariable)
    @test isa(s[:x], ContinuousVariable)
    @test isa(s[:y], DiscreteVariable)
    @test isa(s[:z], DiscreteBoundedVariable)
    tmp = DiscreteBoundedVariable(-10:0.1:10)
    s[:z] = tmp
    @test s[3] === tmp
    for v ∈ s
        @test isa(v, AbstractVariable)
    end
    push!(s, ContinuousVariable(-10, 10))
    @test length(s) == 4
    push!(s, ContinuousVariable(-10, 10), ContinuousBoundedVariable(-10, 10))
    @test length(s) == 6
end

@testset "Variable sampling" begin
    v = ContinuousVariable(1.0, 2.0)
    @test 1.0 <= rand(v) <= 2.0

    v = ContinuousBoundedVariable(1.0, 2.0)
    @test 1.0 <= rand(v) <= 2.0

    v = DiscreteVariable(1.0:0.01:2.0)
    @test 1.0 <= rand(v) <= 2.0

    v = DiscreteBoundedVariable(1.0:0.01:2.0)
    @test 1.0 <= rand(v) <= 2.0
end

@testset "Population constructor and initialisation" begin
    f(x) = -x^2

    v = ContinuousVariable(-10, 10)
    p = Population(f, v, ngen = 100, psze = 30, nslc = 10, nelt = 3)
    @test length(p.vars) == 1

    g(x, y) = -x^2 - y^2
    v1 = ContinuousVariable(-10, 10)
    v2 = ContinuousBoundedVariable(-10, 10)
    p = Population(g, (v1, v2), ngen = 100, psze = 30, nslc = 10, nelt = 3)
    @test length(p.vars) == 2

    @test_throws ArgumentError Population(f, (v1, v2); ngen = 100, psze = 30, nslc = 100, nelt = 3)

    @test_throws ArgumentError Population(f, (v1, v2); ngen = 100, psze = 30, nslc = 10, nelt = 30)

    p = Population(f, (v, v1, v2), ngen = 100, psze = 30, nslc = 10, nelt = 3)
    @test length(p.vars) == 3
end

@testset "fit_model test" begin
    v = ContinuousVariable(-1f0, 1f0)
    v.vpop = randn(10)
    m = fit_model(v, 1:10)
    @test isa(m, Normal)
    @test isa(rand(v), Float32)

    v = ContinuousBoundedVariable(-1f0, 1f0)
    v.vpop = randn(10)
    m = fit_model(v, 1:10)
    @test isa(m, Truncated)
    @test isa(rand(v), Float32)
end

@testset "update_model test" begin
    v = ContinuousVariable(-1f0, 1f0)
    v.vpop = randn(300)
    update_variable!(v, 5, 1:100, 1:20)
    @test length(v.vpop) == 20
end

@testset "Optimisation" begin
    f(x, y) = -x^2 - y^2

    phst = []
    cb1(p) = push!(phst, values(p))
    bhst = []
    cb2(p) = push!(bhst, optimum(p))
    v1 = ContinuousVariable(-10, 10)
    v2 = ContinuousBoundedVariable(-10, 10)
    p = Population(f, (v1, v2), ngen = 100, psze = 50, nslc = 25, nelt = 5)
    optimise!(p, cbbs = cb1, cbas = cb2)
    (x̂, ŷ), ff = optimum(p)
    @test isapprox(x̂, zero(x̂), atol = eps(typeof(x̂)))
    @test isapprox(ŷ, zero(ŷ), atol = eps(typeof(ŷ)))
    @test isapprox(ff, zero(ff), atol = eps(typeof(ff)))
    @test length(phst) == 100
    @test length(bhst) == 100
    @test neval(p) == 99 * 45 + 50

    v1 = DiscreteVariable(-10:0.1:10)
    v2 = DiscreteBoundedVariable(-10:0.1:10)
    p = Population(f, (v1, v2), ngen = 100, psze = 50, nslc = 25, nelt = 5)
    optimise!(p)
    (x̂, ŷ), ff = optimum(p)
    @test isapprox(x̂, zero(x̂), atol = eps(typeof(x̂)))
    @test isapprox(ŷ, zero(ŷ), atol = eps(typeof(ŷ)))
    @test isapprox(ff, zero(ff), atol = eps(typeof(ff)))

    g(x, y) = x^2 + y^2

    v1 = ContinuousVariable(-10, 10)
    v2 = ContinuousBoundedVariable(-10, 10)
    p = Population(g, (v1, v2), ngen = 100, psze = 50, nslc = 25, nelt = 5, fscl = -1)
    optimise!(p)
    (x̂, ŷ), ff = optimum(p)
    @test isapprox(x̂, zero(x̂), atol = eps(typeof(x̂)))
    @test isapprox(ŷ, zero(ŷ), atol = eps(typeof(ŷ)))
    @test isapprox(ff, zero(ff), atol = eps(typeof(ff)))

    v1 = DiscreteVariable(-10:0.1:10)
    v2 = DiscreteBoundedVariable(-10:0.1:10)
    p = Population(g, (v1, v2), ngen = 100, psze = 50, nslc = 25, nelt = 5, fscl = -1)
    optimise!(p, mt = true)
    (x̂, ŷ), ff = optimum(p)
    @test isapprox(x̂, zero(x̂), atol = eps(typeof(x̂)))
    @test isapprox(ŷ, zero(ŷ), atol = eps(typeof(ŷ)))
    @test isapprox(ff, zero(ff), atol = eps(typeof(ff)))

    f1(x, y) = 0.5 - (sin(sqrt(x^2 + y^2))^2 - 0.5) / (1 + 0.001 * (x^2 + y^2))^2
    v = VariableSet(:x => ContinuousBoundedVariable(-5.12, 5.12),
                    :y => ContinuousBoundedVariable(-5.12, 5.12))
    p = Population(f1, v, ngen = 10, psze = 50, nslc = 25, nelt = 5)
    optimise!(p)
    (x̂, ŷ), ff = optimum(p)
    @test_throws ArgumentError optimise!(p)
    s = state(p)
    @test s[:fitness][:mean] == mean(p.fpop)
    @test s[:fitness][:std] == std(p.fpop)
    @test s[:x][:mean] == mean(v[:x].modl)
    @test s[:x][:std] == std(v[:x].modl)
    @test s[:y][:mean] == mean(v[:y].modl)
    @test s[:y][:std] == std(v[:y].modl)
end

#rosenbrock(x, y) = (1 - x)^2 + 100 * (y - x^2)^2

#=
A tester :
x Constructeur Population
x build_model
x update_model!
x rand(o::Population)
-

pyplot()
rosenbrock(x, y) = (1 - x)^2 + 100 * (y - x^2)^2

x = -1.5:0.01:2
y = -0.5:0.01:3

z = Surface(rosenbrock, x, y)
surface(x, y, z, color = :jet1, camera = (-30, 30))

p1 = contour(x, y, rosenbrock, fill = true)

X = repeat(reshape(x, 1, :), length(y), 1)
Y = repeat(y, 1, length(x))
Z = map(rosenbrock, X, Y)
heatmap(x, y, Z, color=:jet1)

function rosenbrock(x::Vector)
  return (1.0 - x[1])^2 + 100.0 * (x[2] - x[1]^2)^2
end

default(size=(600,600), fc=:heat)
x, y = -1.5:0.1:1.5, -1.5:0.1:1.5
z = Surface((x,y)->rosenbrock([x,y]), x, y)
surface(x,y,z, linealpha = 0.3)

=#


text:
###############################################################################
#
# Layer of soil and air to run ODE
#
###############################################################################
"""
    mutable struct VerticalLayers{FT}

Struct that store trace gas information along vertical layers.

# Fields
$(DocStringExtensions.FIELDS)
"""
Base.@kwdef mutable struct VerticalLayers{FT}
    "Number of layers"
    n_layer::Int = 20
    "Thickness per layer"
    d_layer::Array{FT,1} = FT[1 for i in 1:n_layer]
    "Vertical distances among layers (0.5 for boundaries)"
    x_layer::Array{FT,1} = FT[0.5; [1 for i in 1:n_layer-2]; 0.5]

    # Temperature
    "Temperature at each layer"
    T::Array{FT,1} = ones(FT,n_layer) .* WaterPhysics.K_25(FT)

    # CO₂ partial pressure
    "Diffusion coefficient array"
    d_CO₂::Array{FT,1} = diffusive_coefficient.(T,
                                                [TraceGasCO₂()],
                                                [TraceGasAir()])
    "Vertical CO₂ partial pressure"
    p_CO₂::Array{FT,1} = ones(FT,n_layer) .* 41
    "derivative operator for CO₂ partial pressure"
    Δ_CO₂::DerivativeOperator = CenteredDifference(2, 2, FT(1), n_layer)
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land