text:
module Ast

import ..IntPred, ..Attribute

using ..OrderedDict

abstract LLVMAstNode

Base.(:(==))(node1::LLVMAstNode, node2::LLVMAstNode) = begin
    typeof(node1) === typeof(node2) || return false
    for field in names(node1)
        getfield(node1, field) == getfield(node2, field) || return false
    end
    return true
end

abstract LLVMLinkage <: LLVMAstNode

immutable Linkage{T} <: LLVMLinkage end

#TODO: get rid of this
Base.convert(::Type{LLVMLinkage}, enum::Integer) =
    enum == 0 ? Linkage{:External}() : error("unknown")

# -----------------------------------------------------------------------------
# Ast Names
# -----------------------------------------------------------------------------
abstract LLVMName <: LLVMAstNode

immutable Name <: LLVMName
    val::String
end

immutable UnName <: LLVMName
    val::Int
end

# -----------------------------------------------------------------------------
# Parameter Attributes
# -----------------------------------------------------------------------------
abstract ParamAttr <: LLVMAstNode

for typ in [:SignExt, :InReg, :SRet, :NoAlias, :ByVal, :NoCapture, :Nest]
    @eval begin
        immutable $typ <: ParamAttr
        end
    end
end

const param_attr_map = (Int => ParamAttr)[
    Attribute.s_ext => SignExt(),
    Attribute.in_reg => InReg(),
    Attribute.no_alias => NoAlias(),
    Attribute.by_val => ByVal(),
    Attribute.nest => Nest(),
    Attribute.no_capture => NoCapture(),
]

# -----------------------------------------------------------------------------
# Function Attributes
# -----------------------------------------------------------------------------
abstract FuncAttr <: LLVMAstNode

immutable Alignment <: FuncAttr
    val::Int
end

immutable StackAlignment <: FuncAttr
    val::Int
end

for typ in [:NoReturn, :NoUnwind, :ReadNone, :ReadOnly, :NoInline,
            :AlwaysInline, :OptimizeForSize, :StackProtect, :StackProtectReq,
            :NoRedZone, :NoImplicitFloat, :Naked, :InlineHint, :ReturnsTwice,
            :UWTable, :NonLazyBind]
@eval begin
        immutable $typ <: FuncAttr
        end
    end
end

const func_attr_map = (Int => FuncAttr)[
    Attribute.no_return => NoReturn(),
    Attribute.no_unwind => NoUnwind(),
    Attribute.read_none => ReadNone(),
    Attribute.read_only => ReadOnly(),
    Attribute.no_inline => NoInline(),
    Attribute.always_inline => AlwaysInline(),
    Attribute.optimize_for_size => OptimizeForSize(),
    Attribute.stack_protect => StackProtect(),
    Attribute.stack_protect_req => StackProtectReq(),
    Attribute.no_red_zone => NoRedZone(),
    Attribute.no_implicit_float => NoImplicitFloat(),
    Attribute.naked => Naked(),
    Attribute.inline_hint => InlineHint(),
    Attribute.returns_twice => ReturnsTwice(),
    Attribute.uw_table => UWTable()
]

# -----------------------------------------------------------------------------
# Calling Conventions
# -----------------------------------------------------------------------------
abstract CallingConvention <: LLVMAstNode

immutable CConvention <: CallingConvention
end

immutable FastConvention <: CallingConvention
end

immutable ColdConvention <: CallingConvention
end

immutable NumConvention <: CallingConvention
    val::Int
end

# -----------------------------------------------------------------------------
# Pointer Address Space
# -----------------------------------------------------------------------------
immutable AddrSpace <: LLVMAstNode
    val::Int
end

# -----------------------------------------------------------------------------
# Floating Point Format
# -----------------------------------------------------------------------------
abstract FloatFormat <: LLVMAstNode

immutable IEEE <: FloatFormat
end

immutable DoubleExtended <: FloatFormat
end

immutable PairOfFloats <: FloatFormat
end

# -----------------------------------------------------------------------------
# Memory Ordering
# http://llvm.org/docs/LangRef.html#atomic-memory-ordering-constraints
# http://llvm.org/docs/Atomics.html
# -----------------------------------------------------------------------------
abstract MemoryOrdering <: LLVMAstNode

for typ in [:Unordered, :Monotonic, :Acquire, :Release,
            :AcquireRelease, :SequentiallyConsistent]
    @eval begin
        immutable $typ <: MemoryOrdering
        end
    end
end

# -----------------------------------------------------------------------------
# Atomicity
# http://llvm.org/docs/Atomics.html
# http://llvm.org/docs/LangRef.html#singlethread
# -----------------------------------------------------------------------------
immutable Atomicity <: LLVMAstNode
    cross_thread::Bool
    memory_ordering::MemoryOrdering
end

# -----------------------------------------------------------------------------
# Atomic Ops
# http://llvm.org/docs/LangRef.html#atomicrmw-instruction
# -----------------------------------------------------------------------------
abstract AtomicRMWOperation  <: LLVMAstNode

for op in [:Xchg, :Add, :Sub, :And, :Nand, :Or, :Xor, :Max, :Min, :UMax, :Umin]
    typ = symbol(string("Atomic", op))
    @eval begin
        immutable $typ <: AtomicRMWOperation
        end
    end
end

# -----------------------------------------------------------------------------
# Comparison Predicates
# -----------------------------------------------------------------------------
abstract CmpPredicate <: LLVMAstNode

for typ in [:UNE,:SGT,:ORD,:UGT,:SLE,:OGT,:OLT,:OEQ,:ULT,:SGE,
            :ONE,:True,:OGE,:UEQ,:UGE,:NE,:ULE,:UNO,:SLT,:False,:OLE,:EQ]
    @eval begin
        immutable $typ <: CmpPredicate
        end
    end
end

# http://llvm.org/docs/LangRef.html#icmp-instruction
typealias IntCmpPredicate Union(EQ, NE, UGT, UGE, ULT, ULE, SGT, SGE, SLT, SLE)

const int_cmp_pred_map = (Int=>IntCmpPredicate)[
    IntPred.eq  => EQ(),
    IntPred.ne  => NE(),
    IntPred.ugt => UGT(),
    IntPred.uge => UGE(),
    IntPred.ult => ULT(),
    IntPred.ule => ULE(),
    IntPred.sgt => SGT(),
    IntPred.sge => SGE(),
    IntPred.slt => SLT(),
    IntPred.sle => SLE(),
]

# http://llvm.org/docs/LangRef.html#fcmp-instruction
typealias FloatCmpPredicate Union(False, OEQ, OGT, OGE, OLT, OLE, ONE, ORD,
                                  UNO, UEQ, UGT, UGE, ULT, ULE, UNE, True)

# -----------------------------------------------------------------------------
# LLVM Types
# http://llvm.org/docs/LangRef.html#type-system
# -----------------------------------------------------------------------------
abstract LLVMType <: LLVMAstNode

typealias LLVMFloat Union(Float16, Float32, Float64)

# http://llvm.org/docs/LangRef.html#void-type
immutable VoidType <: LLVMType
end

# http://llvm.org/docs/LangRef.html#integer-type
immutable IntType <: LLVMType
    nbits::Int
end

Base.convert{T<:Integer}(::Type{LLVMType}, ::Type{T}) =
    IntType(sizeof(T) * 8)

Base.convert(::Type{LLVMType}, ::Type{BigInt}) =
    error("cannot convert BigInt to LLVMType")

# http://llvm.org/docs/LangRef.html#pointer-type
immutable PtrType <: LLVMType
    typ::LLVMType
    addrspace::AddrSpace
end
PtrType(typ) = PtrType(typ, Ast.AddrSpace(0))
Base.convert{T}(::Type{LLVMType}, ::Type{Ptr{T}}) = PtrType(T, AddrSpace(0))

# http://llvm.org/docs/LangRef.html#floating-point-types
immutable FloatType <: LLVMType
    nbits::Int
    fmt::FloatFormat
end

Base.convert(::Type{LLVMType}, ::Type{Float16}) = FloatType(16, IEEE())
Base.convert(::Type{LLVMType}, ::Type{Float32}) = FloatType(32, IEEE())
Base.convert(::Type{LLVMType}, ::Type{Float64}) = FloatType(64, IEEE())

# http://llvm.org/docs/LangRef.html#function-type
immutable FuncType <: LLVMType
    rettyp::LLVMType
    argtyps::Vector{LLVMType}
    vaargs::Bool
end

# http://llvm.org/docs/LangRef.html#vector-type
immutable VecType <: LLVMType
    typ::LLVMType
    len::Int
end

# http://llvm.org/docs/LangRef.html#structure-type
immutable StructType <: LLVMType
    packed::Bool
    typs::Vector{LLVMType}
end

# http://llvm.org/docs/LangRef.html#array-type
immutable ArrayType <: LLVMType
    typ::LLVMType
    len::Int
end

# http://llvm.org/docs/LangRef.html#opaque-structure-types
immutable NamedTypeRef <: LLVMType
    val::LLVMName
end

# http://llvm.org/docs/LangRef.html#opaque-structure-types
abstract MetadataType <: LLVMType

# -----------------------------------------------------------------------------
# Inline Assembly
# http://llvm.org/docs/LangRef.html#inline-assembler-expressions
# -----------------------------------------------------------------------------
abstract Dialect <: LLVMAstNode

# dialect of assembly used in an inline assembly string
immutable ATTDialect <: Dialect end
immutable IntelDialect <: Dialect end

# Used with 'Call' instruction
type InlineAssembly <: LLVMAstNode
    typ::LLVMType
    assembly::String
    constraints::String
    has_side_effects::Bool
    alignstack::Bool
    dialect::Dialect
end

# -----------------------------------------------------------------------------
# Operand
# -----------------------------------------------------------------------------
# an 'Operand' is roughly an argument to a 'LLVM.Ast.Instruction'
abstract Operand  <: LLVMAstNode
abstract Constant <: LLVMAstNode

# 'MetadataNodeId' is a unique number for identifying a metadata node
immutable MetadataNodeID
    val::Int
end

# http://llvm.org/docs/LangRef.html#metadata
abstract LLVMMetadata <: LLVMAstNode

immutable MetadataNode <: LLVMMetadata
    val::Vector{Union(Nothing, Operand)}
end

immutable MetadataNodeRef <: LLVMMetadata
    val::MetadataNodeID
end

# %foo
type LocalRef <: Operand
    typ::LLVMType
    name::LLVMName
end

type ConstOperand <: Operand
    val::Constant
end

type MetadataStrOperand <: Operand
    val::String
end

type MetadataNodeOperand <: Operand
    val::MetadataNode
end

# the 'LLVM.Ast.Call' instruction is special.  As the callee can be inline asm
typealias CallableOperand Union(Operand, InlineAssembly)

# -----------------------------------------------------------------------------
# Constant
# http://llvm.org/docs/LangRef.html#constants
# http://llvm.org/docs/LangRef.html#constant-expressions
# -----------------------------------------------------------------------------
typealias MaybeConstant Union(Nothing, Constant)

immutable ConstInt <: Constant
    nbits::Int
    val::Integer
end

Base.(:(==))(c1::ConstInt, c2::ConstInt) =
    c1.nbits == c2.nbits && ((v1, v2) = promote(c1.val, c2.val); v1 == v2)

Base.convert{T<:Integer}(::Type{Constant}, v::T) = ConstInt(sizeof(T)*8, v)
Base.convert{T<:Integer}(::Type{MaybeConstant}, v::T) = convert(Constant, v)
Base.convert(::Type{MaybeConstant}, n::Nothing) = nothing

immutable ConstFloat <: Constant
    val::LLVMFloat
end

Base.convert{T<:LLVMFloat}(::Type{Constant}, v::T) = ConstFloat(v)

immutable ConstNull <: Constant
    typ::LLVMType
end

immutable ConstStruct <: Constant
    name::Union(Nothing, LLVMName)
    packed::Bool
    vals::Vector{Constant}
end

immutable ConstArray <: Constant
    typ::LLVMType
    vals::Vector{Constant}
end

immutable ConstVector{T<:Constant} <: Constant
    vals::Vector{T}
end

immutable ConstUndef <: Constant
    typ::LLVMType
end

immutable ConstBlockAddress <: Constant
    func::Name
    block::Name
end

immutable ConstGlobalRef <: Constant
    typ::LLVMType
    name::LLVMName
end

type ConstGetElementPtr <: Constant
    inbounds::Bool
    addr::Constant
    idxs::Vector{Constant}
end

for op in [:Add, :Sub, :Mul, :Shl]
    typ = symbol(string("Const", op))
    @eval begin
        type $typ <: Constant
            nsw::Bool
            nuw::Bool
            op1::Constant
            op2::Constant
        end
    end
end

for op in [:FAdd, :FSub, :FMul, :FDiv, :URem, :SRem, :FRem, :And, :Or, :Xor]
    typ = symbol(string("Const", op))
    @eval begin
        type $typ <: Constant
            op1::Constant
            op2::Constant
        end
    end
end

for op in [:UDiv, :SDiv, :LShr, :AShr]
    typ = symbol(string("Const", op))
    @eval begin
        type $typ <: Constant
            exact::Bool
            op1::Constant
            op2::Constant
        end
    end
end

for op in [:Trunc, :ZExt, :Sext, :FPToUI, :FPToSI, :UIToFP, :SIToFP,
           :FPTrunc, :FPExt, :PtrToInt, :IntToPtr, :BitCast]
    typ = symbol(string("Const", op))
    @eval begin
        type $typ <: Constant
            op::Constant
            typ::LLVMType
        end
    end
end

type ConstICmp <: Constant
    pred::IntCmpPredicate
    op1::Constant
    op2::Constant
end

type ConstFCmp <: Constant
    pred::FloatCmpPredicate
    op1::Constant
    op2::Constant
end

type ConstSelect <: Constant
    cond::Constant
    tval::Operand
    fval::Operand
end

type ConstExtractElement <: Constant
    vec::Constant
    idx::Constant
end

type ConstInsertElement <: Constant
    vec::Constant
    elem::Constant
    idx::Constant
end

type ConstShuffleVector <: Constant
    op1::Constant
    op2::Constant
    mask::Constant
end

type ConstExtractValue <: Constant
    aggregate::Constant
    idxs::Vector{Int}
end

type ConstInsertValue <: Constant
    aggregate::Constant
    elem::Operand
    idxs::Vector{Int}
end

# -----------------------------------------------------------------------------
# Terminal Instructions
# http://llvm.org/docs/LangRef.html#terminators
# -----------------------------------------------------------------------------
abstract Terminator <: LLVMAstNode

# http://llvm.org/docs/LangRef.html#metadata-nodes-and-metadata-strings
typealias InstructionMetadata Vector{(String, MetadataNode)}

type Ret <: Terminator
    op::Union(Nothing, Operand)
    metadata::InstructionMetadata
end
Ret(op) = Ret(op, InstructionMetadata[])

type CondBr <: Terminator
    cond::Operand
    tdest::Name
    fdest::Name
    metadata::InstructionMetadata
end

type Br <: Terminator
    dest::Name
    metadata::InstructionMetadata
end

type Switch <: Terminator
    op::Operand
    defaultdest::Name
    dests::Vector{(Constant, Name)}
    metadata::InstructionMetadata
end

type IndirectBr <: Terminator
    op::Operand
    dests::Vector{Name}
    metadata::InstructionMetadata
end

type Invoke <: Terminator
    convention::CallingConvention
    retattrs::Vector{ParamAttr}
    func::CallableOperand
    args::Vector{(Operand, Vector{ParamAttr})}
    funcattrs::Vector{FuncAttr}
    retdest::Name
    exdest::Name
    metadata::InstructionMetadata
end

type Resume <: Terminator
    op::Operand
    metadata::InstructionMetadata
end

type Unreachable <: Terminator
    metadata::InstructionMetadata
end

# -----------------------------------------------------------------------------
# Non-Terminal Instructions
# http://llvm.org/docs/LangRef.html#binaryops
# http://llvm.org/docs/LangRef.html#bitwiseops
# http://llvm.org/docs/LangRef.html#memoryops
# http://llvm.org/docs/LangRef.html#otherops
# -----------------------------------------------------------------------------
abstract Instruction <: LLVMAstNode
abstract FastMathFlags <: LLVMAstNode

abstract LandingPadClause <: LLVMAstNode

immutable Catch <: LandingPadClause
    val::Constant
end

immutable Filter <: LandingPadClause
    val::Constant
end

type Add <: Instruction
    nsw::Bool
    nuw::Bool
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type FAdd <: Instruction
    fmflags::FastMathFlags
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type Sub <: Instruction
    nsw::Bool
    nuw::Bool
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type FSub <: Instruction
    fmflags::FastMathFlags
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type Mul <: Instruction
    nsw::Bool
    nuw::Bool
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type FMul <: Instruction
    fmflags::FastMathFlags
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type UDiv <: Instruction
    exact::Bool
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type SDiv <: Instruction
    exact::Bool
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type FDiv <: Instruction
    fmflags::FastMathFlags
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type URem <: Instruction
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type SRem <: Instruction
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type FRem <: Instruction
    fmflags::FastMathFlags
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type Sh1 <: Instruction
    nsw::Bool
    nuw::Bool
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type LShr <: Instruction
    exact::Bool
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type AShr <: Instruction
    exact::Bool
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type And <: Instruction
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type Or <: Instruction
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type XOr <: Instruction
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type Alloca <: Instruction
    typ::LLVMType
    len::Union(Nothing, Operand)
    align::Int
    metadata::InstructionMetadata
end

type Load <: Instruction
    volatile::Bool
    addr::Operand
    val::Operand
    atomicity::Union(Nothing, Atomicity)
    align::Int
    metadata::InstructionMetadata
end

type Store <: Instruction
    volatile::Bool
    addr::Operand
    val::Operand
    atomicity::Union(Nothing, Atomicity)
    align::Int
    metadata::InstructionMetadata
end

type GetElemPtr <: Instruction
    inbounds::Bool
    addr::Operand
    idxs::Vector{Operand}
    metadata::InstructionMetadata
end

type Fence <: Instruction
    atomicity::Atomicity
    metadata::InstructionMetadata
end

type CmpXchg <: Instruction
    volatile::Bool
    addr::Operand
    expect::Operand
    replace::Operand
    atomicity::Atomicity
    metadata::InstructionMetadata
end

type AtomicRMW <: Instruction
    volatile::Bool
    op::AtomicRMWOperation
    addr::Operand
    val::Operand
    atomicity::Atomicity
    metadata::InstructionMetadata
end

for typ in [:Trunc, :ZExt, :SExt, :FPToUI, :FPToSI, :UIToFP, :SIToFP,
            :FPTrunc, :FPExt, :PtrToInt, :IntToPtr, :BitCast, :AddrSpaceCast]
    @eval begin
        type $typ <: Instruction
            op::Operand
            typ::LLVMType
            metadata::InstructionMetadata
        end
    end
end

type ICmp <: Instruction
    pred::IntCmpPredicate
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type FCmp <: Instruction
    pred::FloatCmpPredicate
    op1::Operand
    op2::Operand
    metadata::InstructionMetadata
end

type Phi <: Instruction
    typ::LLVMType
    incoming::Vector{(Operand, Name)}
    metadata::InstructionMetadata
end

typealias CallArgs Vector{(Operand, Vector{ParamAttr})}

type Call <: Instruction
    istailcall::Bool
    callcov::CallingConvention
    retattrs::Vector{ParamAttr}
    func::CallableOperand
    args::CallArgs
    funcattrs::Vector{FuncAttr}
    metadata::InstructionMetadata
end

type Select <: Instruction
    cond::Operand
    tval::Operand
    fval::Operand
    metadata::InstructionMetadata
end

type VAArg <: Instruction
    args::Operand
    typ::LLVMType
    metadata::InstructionMetadata
end

type ExtractElement <: Instruction
    vec::Operand
    idx::Operand
    metadata::InstructionMetadata
end

type InsertElement <: Instruction
    vec::Operand
    elem::Operand
    idx::Operand
    metadata::Operand
end

type ShuffleVector <: Instruction
    op1::Operand
    op2::Operand
    mask::Constant
    metadata::InstructionMetadata
end

type ExtractValue <: Instruction
    aggregate::Operand
    idxs::Vector{Int}
    metadata::InstructionMetadata
end

type InsertValue <: Instruction
    aggregate::Operand
    elem::Operand
    idxs::Vector{Int}
    metadata::InstructionMetadata
end

type LandingPad <: Instruction
    typ::LLVMType
    func::Operand
    cleanup::Bool
    clauses::Vector{LandingPadClause}
    metadata::InstructionMetadata
end

# -----------------------------------------------------------------------------
# Visibility
# http://llvm.org/docs/LangRef.html#visibility
# -----------------------------------------------------------------------------
abstract Visibility <: LLVMAstNode

immutable DefaultVisibility <: Visibility
end

immutable HiddenVisibility <: Visibility
end

immutable ProtectedVisibility <: Visibility
end

# -----------------------------------------------------------------------------
# Globals
# http://llvm.org/doxygen/classllvm_1_1GlobalValue.html
# -----------------------------------------------------------------------------
abstract LLVMGlobal <: LLVMAstNode

# Parameters for 'LLVM.Ast.Function'
type Param <: LLVMAstNode
    typ::LLVMType
    name::LLVMName
    attrs::Vector{ParamAttr}
end
Param(typ, name) = Param(typ, name, ParamAttr[])

# LLVM code in a function is a sequence of 'LLVM.Ast.BasicBlocks' with a
# label, some instructions, and terminator
# http://llvm.org/doxygen/classllvm_1_1BasicBlock.html>
type BasicBlock <: LLVMAstNode
    name::LLVMName
    insts::Vector{Instruction}
    term::Terminator
end
BasicBlock(name::LLVMName, term::Terminator) = BasicBlock(name, Instruction[], term)

# http://llvm.org/docs/LangRef.html#global-variables
type GlobalVar <: LLVMGlobal
    name::LLVMName
    linkage::LLVMLinkage
    visibility::Visibility
    threadlocal::Bool
    addrspace::AddrSpace
    unamedaddr::Bool
    isconst::Bool
    typ::LLVMType
    init::Union(Nothing, Constant)
    section::Union(Nothing, String)
    alignment::Int
end

#TODO: get rid of this
Base.convert(::Type{Visibility}, enum::Integer) =
    enum == 0 ? DefaultVisibility() : error("unhandled visibility")

GlobalVar(;name=error("global variable name not defined"),
           linkage=Linkage{:External}(),
           visibility=DefaultVisibility(),
           threadlocal=false,
           addrspace=AddrSpace(0),
           has_unnamed_addr=false,
           isconst=false,
           typ=error("global variable typ not defined"),
           init=nothing,
           section=nothing,
           alignment=0) = begin
    GlobalVar(name,
              linkage,
              visibility,
              threadlocal,
              addrspace,
              has_unnamed_addr,
              isconst,
              typ,
              init,
              section,
              alignment)
end

# http://llvm.org/docs/LangRef.html#aliases
type GlobalAlias <: LLVMGlobal
    name::LLVMName
    linkage::LLVMLinkage
    visibility::Visibility
    typ::LLVMType
    aliasee::Constant
end

# http://llvm.org/docs/LangRef.html#functions
type Func <: LLVMGlobal
    linkage::LLVMLinkage
    visibility::Visibility
    callingcov::CallingConvention
    retattrs::Vector{ParamAttr}
    rettyp::LLVMType
    name::LLVMName
    params::Vector{Param}
    vaargs::Bool
    attrs::Vector{FuncAttr}
    section::Union(Nothing, String)
    alignment::Int
    gcname::Union(Nothing, String)
    blocks::Vector{BasicBlock}
end

Func(;linkage=Linkage{:External}(),
      visibility=DefaultVisibility(),
      callingcov=CConvention(),
      retattrs=ParamAttr[],
      rettyp=error("return type not defined"),
      name=error("function name not defined"),
      params=Param[],
      vaargs=false,
      attrs=FuncAttr[],
      section=nothing,
      alignment=0,
      gcname=nothing,
      blocks=BasicBlock[]) = Func(linkage,
                                  visibility,
                                  callingcov,
                                  retattrs,
                                  rettyp,
                                  name,
                                  params,
                                  vaargs,
                                  attrs,
                                  section,
                                  alignment,
                                  gcname,
                                  blocks)

# -----------------------------------------------------------------------------
# DataLayout
# http://llvm.org/docs/LangRef.html#data-layout>
# -----------------------------------------------------------------------------
abstract Endianness <: LLVMAstNode

immutable LittleEndian <: Endianness
end

immutable BigEndian <: Endianness
end

# describes how a given type must be aligned
immutable AlignmentInfo <: LLVMAstNode
    abi::Int
    perferred::Union(Nothing, Int)
end

# typeo of type for which 'AlignmentInfo' may be specified
abstract AlignType <: LLVMAstNode

immutable IntegerAlign <: AlignType end
immutable VectorAlign <: AlignType end
immutable FloatAlign <: AlignType end
immutable AggregateAlign <: AlignType end

# style of name mangling
abstract Mangling <: LLVMAstNode

immutable ELFMangling <: Mangling end
immutable MIPSMangling <: Mangling end
immutable MachOMangling <: Mangling end
immutable WindowsCOFFMangling <: Mangling end

typealias PointerLayoutMap Associative{AddrSpace, (Int, AlignmentInfo)}
typealias PointerLayoutDict OrderedDict{AddrSpace, (Int, AlignmentInfo)}

typealias TypeLayoutMap Associative{(AlignType, Int), AlignmentInfo}
typealias TypeLayoutDict OrderedDict{(AlignType, Int), AlignmentInfo}

# description of various data layout properties which may be used during optimization
type DataLayout <: LLVMAstNode
    endianness::Union(Nothing, Endianness)
    mangling::Union(Nothing, Mangling)
    stackalignment::Union(Nothing, Int)
    pointerlayouts::PointerLayoutMap
    typelayouts::TypeLayoutMap
    nativesizes::Union(Nothing, Set{Int})
end

DataLayout(;endianness=nothing,
            mangling=nothing,
            stackalignment=nothing,
            pointerlayouts=PointerLayoutDict(),
            typelayouts=TypeLayoutDict(),
            nativesizes=nothing) = DataLayout(endianness,
                                              mangling,
                                              stackalignment,
                                              pointerlayouts,
                                              typelayouts,
                                              nativesizes)

# -----------------------------------------------------------------------------
# Module
# -----------------------------------------------------------------------------
# Anything that can be at the top level of a 'Module'
abstract Definition <: LLVMAstNode

type GlobalDef <: Definition
    val::LLVMGlobal
end

#TODO: revist this
type TypeDef <: Definition
    name::LLVMName
    typs#::Vector{LLVMType}
end

type MetadataNodeDef <: Definition
    id::MetadataNodeID
    ops::Vector{Operand}
end

type NamedMetadataDef <: Definition
    name::String
    ids::Vector{MetadataNodeID}
end

type ModInlineAsm <: Definition
    source::String
end

type Mod <: LLVMAstNode
    name::String
    layout::Union(Nothing, DataLayout)
    target::Union(Nothing, String)
    defs::Vector{Definition}
end

end


text:
@generated function tri2int(x::NTuple{N, Int}) where N
    ex = :(x[1])
    for i = 2:N
        base = 3^(i-1)
        ex = :($ex + x[$i] * $base)
    end
    return ex
end

tri2int(x::Int...) = tri2int(x)

# 1 + 2 * 3 + 3 * 9
@show tri2int(1, 2, 3)

ext: jl


text:
if haskey(ENV, "CI") && ENV["CI"] == "true"
    ENV["PYTHON"] = ""
    using Pkg
    Pkg.build("PyCall")
end

using Test
using CUDA
using PyPlot
using Oceananigans

using ConvergenceTests
using ConvergenceTests.OneDimensionalCosineAdvectionDiffusion: run_test
using ConvergenceTests.OneDimensionalUtils: plot_solutions!, plot_error_convergence!, unpack_errors

""" Run advection-diffusion test for all Nx in resolutions. """
function run_convergence_test(κ, U, resolutions, arch)

    # Determine safe time-step
             Lx = 2π
      stop_time = 0.01
              h = Lx / maximum(resolutions)
    proposal_Δt = 1e-3 * min(h / U, h^2 / κ)

    # Adjust time-step
    stop_iteration = round(Int, stop_time / proposal_Δt)
                Δt = stop_time / stop_iteration

    # Run the tests
    results = [run_test(architecture=arch, Nx=Nx, Δt=Δt, stop_iteration=stop_iteration, U=U, κ=κ)
               for Nx in resolutions]

    return results
end

#####
##### Run test
#####

arch = CUDA.has_cuda() ? GPU() : CPU()

Nx = 2 .^ (3:7) # N = 8 through N = 256
diffusion_results = run_convergence_test(1e-1, 0, Nx, arch)
advection_results = run_convergence_test(1e-6, 3, Nx, arch)
advection_diffusion_results = run_convergence_test(1e-2, 1, Nx, arch)

#####
##### Plot solution and error profile
#####

all_results = (diffusion_results, advection_results, advection_diffusion_results)
names = ("diffusion only", "advection only", "advection-diffusion")
linestyles = ("-", "--", ":")
specialcolors = ("xkcd:black", "xkcd:indigo", "xkcd:wine red")

# Solutions
close("all")

fig, axs = subplots(nrows=2, figsize=(12, 6), sharex=true)

legends = plot_solutions!(axs, all_results, names, linestyles, specialcolors)

filename = "cosine_advection_diffusion_solutions_$(typeof(arch)).png"
filepath = joinpath(@__DIR__, "figs", filename)
mkpath(dirname(filepath))
savefig(filepath, dpi=480, bbox_extra_artists=legends, bbox_inches="tight")

# Error profile
fig, axs = subplots()

legend = plot_error_convergence!(axs, Nx, all_results, names)

filename = "cosine_advection_diffusion_error_convergence_$(typeof(arch)).png"
filepath = joinpath(@__DIR__, "figs", filename)
mkpath(dirname(filepath))
savefig(filepath, dpi=480, bbox_extra_artists=(legend,), bbox_inches="tight")

# Test rate of convergence
for (results, name) in zip(all_results, names)
    name = "1D cosine " * name
    @info "Testing rate of convergence for $name..."

    cx_L₁, cy_L₁, cz_L₁,
           uy_L₁, uz_L₁,
    vx_L₁,        vz_L₁,
    wx_L₁, wy_L₁,
    cx_L∞, cy_L∞, cz_L∞,
           uy_L∞, uz_L∞,
    vx_L∞,        vz_L∞,
    wx_L∞, wy_L∞         = unpack_errors(results)

    test_rate_of_convergence(uy_L₁, Nx, expected=-2.0, atol=0.01, name=name*" ux_L₁")
    test_rate_of_convergence(uz_L₁, Nx, expected=-2.0, atol=0.01, name=name*" uz_L₁")
    test_rate_of_convergence(vx_L₁, Nx, expected=-2.0, atol=0.01, name=name*" vx_L₁")
    test_rate_of_convergence(vz_L₁, Nx, expected=-2.0, atol=0.01, name=name*" vz_L₁")
    test_rate_of_convergence(wx_L₁, Nx, expected=-2.0, atol=0.01, name=name*" wx_L₁")
    test_rate_of_convergence(wy_L₁, Nx, expected=-2.0, atol=0.01, name=name*" wy_L₁")
    test_rate_of_convergence(cx_L₁, Nx, expected=-2.0, atol=0.01, name=name*" cx_L₁")
    test_rate_of_convergence(cy_L₁, Nx, expected=-2.0, atol=0.01, name=name*" cy_L₁")
    test_rate_of_convergence(cz_L₁, Nx, expected=-2.0, atol=0.01, name=name*" cz_L₁")

    test_rate_of_convergence(uy_L∞, Nx, expected=-2.0, atol=0.01, name=name*" ux_L∞")
    test_rate_of_convergence(uz_L∞, Nx, expected=-2.0, atol=0.01, name=name*" uz_L∞")
    test_rate_of_convergence(vx_L∞, Nx, expected=-2.0, atol=0.01, name=name*" vx_L∞")
    test_rate_of_convergence(vz_L∞, Nx, expected=-2.0, atol=0.01, name=name*" vz_L∞")
    test_rate_of_convergence(wx_L∞, Nx, expected=-2.0, atol=0.01, name=name*" wx_L∞")
    test_rate_of_convergence(wy_L∞, Nx, expected=-2.0, atol=0.01, name=name*" wy_L∞")
    test_rate_of_convergence(cx_L∞, Nx, expected=-2.0, atol=0.01, name=name*" cx_L∞")
    test_rate_of_convergence(cy_L∞, Nx, expected=-2.0, atol=0.01, name=name*" cy_L∞")
    test_rate_of_convergence(cz_L∞, Nx, expected=-2.0, atol=0.01, name=name*" cz_L∞")

    @test uy_L₁ ≈ uz_L₁ ≈ vx_L₁ ≈ vz_L₁ ≈ wx_L₁ ≈ wy_L₁ ≈ cx_L₁ ≈ cy_L₁ ≈ cz_L₁
    @test uy_L∞ ≈ uz_L∞ ≈ vx_L∞ ≈ vz_L∞ ≈ wx_L∞ ≈ wy_L∞ ≈ cx_L∞ ≈ cy_L∞ ≈ cz_L∞
end


text:
using Pkg
using RetroCap
using Test

@testset "compress_versions()" begin
    # Test exact version matching
    vs = [v"1.1.0", v"1.1.1", v"1.1.2"]
    @test RetroCap.Compress.compress_versions(vs, [vs[2]]) == Pkg.Types.VersionSpec("1.1.1")

    # Test holes
    vs = [v"1.1.0", v"1.1.1", v"1.1.4"]
    @test RetroCap.Compress.compress_versions(vs, [vs[2]]) == Pkg.Types.VersionSpec("1.1.1")

    # Test patch variation with length(subset) > 1
    vs = [v"1.1.0", v"1.1.1", v"1.1.2", v"1.1.3", v"1.2.0"]
    @test RetroCap.Compress.compress_versions(vs, [vs[2], vs[3]]) == Pkg.Types.VersionSpec("1.1.1-1.1.2")

    # Test minor variation
    vs = [v"1.1.0", v"1.1.1", v"1.2.0"]
    @test RetroCap.Compress.compress_versions(vs, [vs[2]]) == Pkg.Types.VersionSpec("1.1.1-1.1")

    # Test major variation
    vs = [v"1.1.0", v"1.1.1", v"1.2.0", v"2.0.0"]
    @test RetroCap.Compress.compress_versions(vs, [vs[2], vs[3]]) == Pkg.Types.VersionSpec("1.1.1-1")

    # Test build numbers and prerelease values are ignored
    vs = [v"1.1.0-alpha", v"1.1.0+0", v"1.1.0+1"]
    @test RetroCap.Compress.compress_versions(vs, [vs[2]]) == Pkg.Types.VersionSpec("1")
end

@testset "_my_string()" begin
    @test RetroCap.Compress._my_string("foo") == "foo"
    @test RetroCap.Compress._my_string(v"1.2.3") == "1.2.3"
    @test RetroCap.Compress._my_string(VersionNumber(1, 2, 3)) == "1.2.3"
end

@testset "_make_keys_strings" begin
    a = Dict{Any, Any}()
    a[v"1.2.3"] = "foo"
    @test eltype([collect(keys(a))...]) <: VersionNumber
    b = RetroCap.Compress._make_keys_strings(a)
    @test eltype([collect(keys(b))...]) <: AbstractString
end

@testset "_save_uncompressed" begin
    RetroCap.with_temp_dir() do tmp_dir
        cd(tmp_dir)
        run(`git clone https://github.com/JuliaRegistries/General.git`)
        cd(joinpath(tmp_dir, "General", "P", "PredictMD"))
        compat = RetroCap.Compress.load("Compat.toml")
        rm("Compat.toml"; force = true, recursive = true)
        RetroCap.Compress._save_uncompressed("Compat.toml", compat)
        deps = RetroCap.Compress.load("Deps.toml")
        rm("Deps.toml"; force = true, recursive = true)
        RetroCap.Compress._save_uncompressed("Deps.toml", deps)
        compat = RetroCap.Compress.load("Compat.toml")
        rm("Compat.toml"; force = true, recursive = true)
        RetroCap.Compress.save("Compat.toml", compat)
        deps = RetroCap.Compress.load("Deps.toml")
        rm("Deps.toml"; force = true, recursive = true)
        RetroCap.Compress.save("Deps.toml", deps)
    end
end


text:
# Copyright (c) 2022 Code Komali
#
# This software is released under the MIT License.
# https://opensource.org/licenses/MIT

using Test, Printf, Memoize

# A global variable just for simple benchmarking
called = 0

# Using the Memoize package
@memoize Dict function knapsack(values, weights, MaxWieght)
    global called += 1
    if MaxWieght < 0
        return -Inf
    elseif isempty(values)
        return 0
    else
        max(values[end] + knapsack(values[1:end-1], weights[1:end-1], MaxWieght - weights[end]),
            knapsack(values[1:end-1], weights[1:end-1], MaxWieght))
    end
end

# Non memoized bruteforce approach
function knapsack_bruteforce(values, weights, MaxWieght)
    global called += 1
    if MaxWieght < 0
        return -Inf
    elseif isempty(values)
        return 0
    else
        max(values[end] + knapsack_bruteforce(values[1:end-1], weights[1:end-1], MaxWieght - weights[end]),
            knapsack_bruteforce(values[1:end-1], weights[1:end-1], MaxWieght))
    end
end

function compute_knapsack(values, weights, MaxWieght; memoized = true)
    empty!(memoize_cache(knapsack))
    global called = 0
    result = memoized ? knapsack(values, weights, MaxWieght) : knapsack_bruteforce(values, weights, MaxWieght)
    #uncomment the below code for call count
    #@printf "\n Function called: %d \n" called
    return result
end

@testset "0-1 knapsack problem tests" begin
    @test compute_knapsack([20, 5, 10, 40, 15, 25], [1, 2, 3, 8, 7, 4], 10, memoized = false) == 60
    @test compute_knapsack([60, 100, 120], [10, 20, 30], 50, memoized = false) == 220
    @test compute_knapsack([20, 5, 10, 40, 15, 25], [1, 2, 3, 8, 7, 4], 10) == 60
    @test compute_knapsack([60, 100, 120], [10, 20, 30], 50) == 220
    @test compute_knapsack([92, 57, 49, 68, 60, 43, 67, 84, 87, 72],
        [23, 31, 29, 44, 53, 38, 63, 85, 89, 82], 165) == 309
end

@time compute_knapsack([92, 57, 49, 68, 60, 43, 67, 84, 87, 72],
    [23, 31, 29, 44, 53, 38, 63, 85, 89, 82], 165)
@time compute_knapsack([92, 57, 49, 68, 60, 43, 67, 84, 87, 72],
    [23, 31, 29, 44, 53, 38, 63, 85, 89, 82], 165, memoized = false)

#Note:
# The memoized version is being called less number of times as expected.
# But, it is taking more time to compute. Need to investigate more.

max_stars_repo_name: codekomali/JuliaPractice

max_issues_repo_name: codekomali/JuliaPractice

max_forks_repo_name: codekomali/JuliaPractice


text:
using Test
using NQCDynamics
using NQCDynamics.InitialConditions.MetropolisHastings
using Unitful
using StatsBase
using Random
Random.seed!(10)

atoms = Atoms{Float64}([:H, :H, :C])
cell = PeriodicCell(hcat(1))
model = NQCModels.Harmonic()

Δ = Dict([(:H, 0.1), (:C, 0.1)])
monte_carlo = MetropolisHastings.MonteCarlo{Float64}(Δ, length(atoms), 100, Int[], x->true)
sim = Simulation(atoms, model; cell=cell, temperature=100u"K")

@testset "propose_move!" begin
    Rᵢ = zeros(size(sim))
    Rₚ = zero(Rᵢ)
    @test Rᵢ == Rₚ
    MetropolisHastings.propose_move!(sim, monte_carlo, Rᵢ, Rₚ, 1)
    @test Rᵢ != Rₚ
end

@testset "assess_proposal!" begin
    Rᵢ = randn(size(sim))
    Rₚ = zero(Rᵢ)
    # Accept the move
    monte_carlo.Eᵢ = 1
    output = MetropolisHastings.MonteCarloOutput(Rᵢ, sim.atoms)
    MetropolisHastings.assess_proposal!(sim, monte_carlo, Rᵢ, Rₚ, output, 1)
    @test output.R[1] == Rₚ
    # Reject the move
    monte_carlo.Eᵢ = -1
    MetropolisHastings.assess_proposal!(sim, monte_carlo, Rᵢ, Rₚ, output, 1)
    @test output.R[1] == Rᵢ
end

@testset "write_output!" begin
    Rₚ = fill(0.1, size(sim))
    output = MetropolisHastings.MonteCarloOutput(Rₚ, sim.atoms)
    MetropolisHastings.write_output!(output, Rₚ, 1.0)
    Rₚ .+= 1
    MetropolisHastings.write_output!(output, Rₚ, 1.1)
    @test all(output.R[1] .== 0.1)
    @test output.energy[1] == 1.0
    @test all(output.R[2] .== 1.1)
    @test output.energy[2] == 1.1
end

@testset "acceptance_probability" begin
    e1 = 0.0
    e2 = -0.001
    @test MetropolisHastings.acceptance_probability(e2, e1, 1.0) == 1
    e1 = 0.0
    e2 = 0.001
    @test 1 > MetropolisHastings.acceptance_probability(e2, e1, 1.0) > 0
    e1 = 0.0
    e2 = Inf
    @test MetropolisHastings.acceptance_probability(e2, e1, 1.0) == 0
end

@testset "run_monte_carlo_sampling" begin
    R0 = zeros(size(sim))
    out = MetropolisHastings.run_monte_carlo_sampling(sim, R0, Δ, 10)
    @test !(out.R[1] ≈ out.R[10])
    @test !(out.energy[1] ≈ out.energy[20])
end

@testset "run_monte_carlo_sampling" begin
    sim = RingPolymerSimulation(atoms, model, 10; cell=cell, temperature=100u"K")
    R0 = zeros(size(sim))
    out = MetropolisHastings.run_monte_carlo_sampling(sim, R0, Δ, 10)
    @test !(out.R[1] ≈ out.R[10])
    @test !(out.energy[1] ≈ out.energy[20])
end

@testset "propose_centroid_move!" begin
    monte_carlo = MetropolisHastings.PathIntegralMonteCarlo{Float64}(Δ, length(atoms), 100, [1], 1.0, 10, x->true)
    sim = RingPolymerSimulation(atoms, model, 10; cell=cell, temperature=100u"K", quantum_nuclei=[:H])
    Rᵢ = randn(3, length(sim.atoms), 10)
    Rₚ = copy(Rᵢ)
    MetropolisHastings.propose_centroid_move!(sim, monte_carlo, Rᵢ, Rₚ, 2)
    @test mean(Rₚ[:,1,:]) ≈ mean(Rᵢ[:,1,:]) # Check fixed atom does not move

    RingPolymers.transform_to_normal_modes!(sim.beads, Rₚ)
    RingPolymers.transform_to_normal_modes!(sim.beads, Rᵢ)
    @test !(Rₚ[:,:,1] ≈ Rᵢ[:,:,1]) # Test centroid has moved
    @test Rₚ[:,1:2,2:end] ≈ Rᵢ[:,1:2,2:end] # Test only centroid moves for quantum atoms
end


text:
# Recursively search for type T and perform an operation on it based on julia/base/deepcopy.jl

function deepmap(f, T, x)
    @assert !isbitstype(T) "deepmap does not support isbitstype(T)"
    deepmap_internal(f, T, Set{UInt}(), x)
end

function deepmap_internal(f, T, s::Set{UInt}, x::Union{Module,Symbol,Core.MethodInstance,Method,GlobalRef,DataType,Union,UnionAll,Task,Regex})
    if x isa T && objectid(x) ∉ s
        push!(s, objectid(x))
        f(x)
    end
end

function deepmap_internal(f, T, s::Set{UInt}, x::Union{Tuple,Core.SimpleVector})
    objectid(x) ∈ s && return
    push!(s, objectid(x))
    x isa T && f(x)
    for xi in x
        !isbits(xi) && deepmap_internal(f, T, s, xi)
    end
end

function deepmap_internal(f, T, s::Set{UInt}, x::String)
    objectid(x) ∈ s && return
    push!(s, objectid(x))
    x isa T && f(x)
end

function deepmap_internal(f, T, s::Set{UInt}, @nospecialize(x))
    X = typeof(x)::DataType
    isbitstype(X) && return
    objectid(x) ∈ s && return
    push!(s, objectid(x))
    x isa T && f(x)
    for i in 1:nfields(x)
        if isdefined(x, i)
            xi = getfield(x,i)
            !isbits(xi) && deepmap_internal(f, T, s, xi)
        end
    end
end

function deepmap_internal(f, T, s::Set{UInt}, x::Array)
    objectid(x) ∈ s && return
    push!(s, objectid(x))
    x isa T && f(x)
    isbitstype(eltype(x)) && return
    for i = 1:length(x)
        if ccall(:jl_array_isassigned, Cint, (Any, Csize_t), x, i-1) != 0
            xi = ccall(:jl_arrayref, Any, (Any, Csize_t), x, i-1)
            !isbits(xi) && deepmap_internal(f, T, s, xi)
        end
    end
end

function deepmap_internal(f, T, s::Set{UInt}, x::Union{Dict,IdDict})
    objectid(x) ∈ s && return
    push!(s, objectid(x))
    x isa T && f(x)
    isbitstype(eltype(x)) && return
    for (k, v) in x
        !isbits(k) && deepmap_internal(f, T, s, k)
        !isbits(v) && deepmap_internal(f, T, s, v)
    end
end


text:
module mmmvvoxelmm1
using FinEtools
using Compat.Test
function test()
    V = VoxelBoxVolume(Int, 5*[5,6,7], [4.0, 4.0, 5.0])

    s1 = solidsphere((1.5, 2.0, 2.0), 1.3)
    s2 = solidsphere((1.5, 1.0, 2.0), 1.3)
    fillsolid!(V, differenceop(s1, s2), 1)

    File = "mmmvvoxelmm1.vtk"
    vtkexport(File, V)
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    @test (V.data[15,13,12]==0)
    @test (V.data[15,16,15]==1)
end
end
using .mmmvvoxelmm1
mmmvvoxelmm1.test()

module mmmvvoxelmm2
using FinEtools
using Compat.Test
function test()
    V = VoxelBoxVolume(Int, 5*[5,6,7], [4.0, 4.0, 5.0])

    b1 = solidbox((0.0, 0.0, 0.0), (1.0, 4.0, 5.0))
    b2 = solidbox((0.0, 0.0, 0.0), (4.0, 1.0, 5.0))
    fillsolid!(V, unionop(b1, b2), 1)

    File = "mmmvvoxelmm2.vtk"
    vtkexport(File, V)
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    # println("$(V.data[15,13,12])")
    # println("$(V.data[1,1,12])")
    @test (V.data[15,13,12]==0)
    @test (V.data[1,1,12]==1)
end
end
using .mmmvvoxelmm2
mmmvvoxelmm2.test()

module mmmvvoxelmm3
using FinEtools
using Compat.Test
function test()
    V = VoxelBoxVolume(Int, 15*[5,6,7], [4.0, 4.0, 5.0])

    b1 = solidbox((0.0, 0.0, 0.0), (1.0, 4.0, 5.0))
    b2 = solidbox((0.0, 0.0, 0.0), (4.0, 1.0, 5.0))
    fillsolid!(V, unionop(b1, b2), 1)
    s2 = solidsphere((2.0, 2.0, 2.5), 1.0)
    fillsolid!(V, s2, 2)

    File = "mmmvvoxelmm3.vtk"
    vtkexport(File, V)
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    # println("$(V.data[15,13,12])")
    # println("$(V.data[1,1,12])")
    # @test (V.data[15,13,12]==0)
    # @test (V.data[1,1,12]==1)
end
end
using .mmmvvoxelmm3
mmmvvoxelmm3.test()

module mmmvvoxelmm4
using FinEtools
using Compat.Test
function test()
    V = VoxelBoxVolume(Int, 5*[5,6,7], [4.0, 4.0, 5.0])

    h1 = solidcylinder((2.0, 2.0, 2.5), (1.0, 0.0, 0.0), 0.5)
    fillsolid!(V, h1, 1)

    File = "mmmvvoxelmm4.vtk"
    vtkexport(File, V)
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    # println("$(V.data[15,13,12])")
    # println("$(V.data[1,15,12])")
    # println("$(V.data[13,15,20])")
    @test (V.data[15,13,12]==0)
    @test (V.data[1,15,12]==0)
    @test (V.data[13,15,20]==1)
end
end
using .mmmvvoxelmm4
mmmvvoxelmm4.test()

module mmmvvoxelmm5
using FinEtools
using Compat.Test
function test()
    V = VoxelBoxVolume(Int, 15*[5,6,7], [4.0, 4.0, 5.0])

    b1 = solidbox((0.0, 0.0, 0.0), (1.0, 4.0, 5.0))
    b2 = solidbox((0.0, 0.0, 0.0), (4.0, 1.0, 5.0))
    h1 = solidcylinder((2.0, 3.0, 2.5), (1.0, 0.0, 0.0), 0.5)
    fillsolid!(V, differenceop(unionop(b1, b2), h1), 1)

    s2 = solidsphere((2.0, 2.0, 2.5), 1.0)
    fillsolid!(V, s2, 2)

    File = "mmmvvoxelmm5.vtk"
    vtkexport(File, V)
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    # println("$(V.data[15,13,12])")
    # println("$(V.data[1,35,12])")
    # println("$(V.data[35,45,50])")
    @test (V.data[15,13,12]==1)
    @test (V.data[1,35,12]==1)
    @test (V.data[35,45,50]==2)
end
end
using .mmmvvoxelmm5
mmmvvoxelmm5.test()

module mmmvvoxelmm6
using FinEtools
using Compat.Test
function test()
    raw = zeros(UInt8, 13, 14, 15)
    raw[5:7, 2:13, 12:14] .= 1
    V = VoxelBoxVolume(raw, [4.0, 4.0, 5.0])

    File = "mmmvvoxelmm6.vtk"
    vtkexport(File, V)
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    # println("$(V.data[7,13,12])")
    # println("$(V.data[1,12,12])")
    # println("$(V.data[12,10,5])")
    @test (V.data[7,13,12]==1)
    @test (V.data[1,12,12]==0)
    @test (V.data[12,10,5]==0)
end
end
using .mmmvvoxelmm6
mmmvvoxelmm6.test()

module mmmvvoxelmm7
using FinEtools
using Compat.Test
function test()
    raw = zeros(UInt8, 13, 14, 15)
    V = VoxelBoxVolume(raw, [4.0, 4.0, 5.0])
    fillvolume!(V, 2)
    V.data[5:7, 2:13, 12:14] .= 1

    File = "mmmvvoxelmm7.vtk"
    vtkexport(File, V)
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    # println("$(V.data[7,13,12])")
    # println("$(V.data[1,12,12])")
    # println("$(V.data[12,10,5])")
    @test (V.data[7,13,12]==1)
    @test (V.data[1,12,12]==2)
    @test (V.data[12,10,5]==2)
end
end
using .mmmvvoxelmm7
mmmvvoxelmm7.test()

module mmmvvoxelmm8
using FinEtools
using Compat.Test
function test()
    raw = zeros(UInt8, 13, 14, 15)
    V = VoxelBoxVolume(raw, [4.0, 4.0, 5.0])
    V.data[5:7, 2:13, 12:14] .= 1
    Vt = trim(V, 0)
    @test size(Vt) == (3, 12, 3)
    @test size(V, 3) == 15
    Vtt = trim(Vt, 0)
    @test size(Vtt) == (3, 12, 3)

    File = "mmmvvoxelmm8.vtk"
    vtkexport(File, V)
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end
end
end
using .mmmvvoxelmm8
mmmvvoxelmm8.test()

module mmmvvoxelmm9
using FinEtools
using Compat.Test
function test()
    raw = zeros(UInt8, 13, 14, 15)
    fill!(raw, 2)
    V = VoxelBoxVolume(raw, [4.0, 4.0, 5.0])
    V.data[5:7, 2:13, 12:14] .= 1
    Vt = trim(V, 2)
    Vp = pad(Vt, (4, 6), (1, 1), (11, 1), 0)
    @test size(V) == size(Vp)
    @test size(Vt) == (3, 12, 3)
    @test size(V, 3) == 15
    Vtt = trim(Vt, 0)
    @test size(Vtt) == (3, 12, 3)

    # File = "mmmvvoxelmm9.vtk"
    # vtkexport(File, V)
    # # @async run(`"paraview.exe" $File`)
    # try rm(File) catch end

    # println("$(V.data[7,13,12])")
    # println("$(V.data[1,12,12])")
    # println("$(V.data[12,10,5])")
    # @test (V.data[7,13,12]==1)
    # @test (V.data[1,12,12]==2)
    # @test (V.data[12,10,5]==2)
end
end
using .mmmvvoxelmm9
mmmvvoxelmm9.test()

module mvoxelmm12
using FinEtools
using Compat.Test
function test()
    V = VoxelBoxVolume(Int, 5*[5,6,7], [4.0, 4.0, 5.0])

    b1 = solidbox((0.0, 0.0, 0.0), (1.0, 4.0, 5.0))
    b2 = solidbox((0.0, 0.0, 0.0), (4.0, 1.0, 5.0))
    fillsolid!(V, unionop(b1, b2), 2)
    fillsolid!(V, complementop(unionop(b1, b2)), 1)

    File = "mvoxelmm12.vtk"
    vtkexport(File, V)
    # @async run(`"paraview.exe" $File`)
    try rm(File) catch end

    # println("$(V.data[15,13,12])")
    # println("$(V.data[1,1,12])")
    @test (V.data[15,13,12]==1)
    @test (V.data[1,1,12]==2)
end
end
using .mvoxelmm12
mvoxelmm12.test()

module mvoxelmm13
using FinEtools
using Compat.Test
function test()
    V = VoxelBoxVolume(Int, 7*[5,6,7], [4.0, 4.0, 5.0])

    for k = 1:size(V, 3)
        for j = 1:size(V, 2)
            for i = 1:size(V, 1)
                V.data[i, j, k] = 35^2-(i^2+j^2+k^2)
            end
        end
    end

    threshold_value, voxel_below, voxel_above = 20, 1, 0
    V = threshold(V, threshold_value, voxel_below, voxel_above)

    # File = "mvoxelmm13.vtk"
    # vtkexport(File, V)
    # @async run(`"paraview.exe" $File`)
    # try rm(File) catch end

    # println("$(V.data[15,13,12])")
    # println("$(V.data[1,1,12])")
    # @test (V.data[15,13,12]==1)
    @test (V.data[1,1,1]==0)
    @test (V.data[end,end,end] != 0)
end
end
using .mvoxelmm13
mvoxelmm13.test()


text:
#! /usr/bin/julia

# Rosetta Code, AKS test for primes

function polycoefs(n::Int64)
    pc = typeof(n)[]
    if n < 0
        return pc
    end
    sgn = one(n)
    for k in n:-1:0
        push!(pc, sgn*binomial(n, k))
        sgn = -sgn
    end
    return pc
end

function stringpoly(n::Int64)
    if n < 0
        return ""
    end
    st = @sprintf "(x - 1)^{%d} & = & " n
    for (i, c) in enumerate(polycoefs(n))
        if i == 1
            op = ""
            ac = c
        elseif c < 0
            op = "-"
            ac = abs(c)
        else
            op = "+"
            ac = abs(c)
        end
        p = n + 1 - i
        if p == 0
            st *= @sprintf " %s %d\\\\" op ac
        elseif ac == 1
            st *= @sprintf " %s x^{%d}" op p
        else
            st *= @sprintf " %s %dx^{%d}" op ac p
        end
    end
    return st
end

function isaksprime(n::Int64)
    if n < 2
        return false
    end
    for c in polycoefs(n)[2:(end-1)]
        if c%n != 0
            return false
        end
    end
    return true
end

println("<math>")
println("\\begin{array}{lcl}")
for i in 0:10
    println(stringpoly(i))
end
println("\\end{array}")
println("</math>\n")

L = 50
print("AKS primes less than ", L, ":  ")
sep = ""
for i in 1:L
    if isaksprime(i)
        print(sep, i)
        sep = ", "
    end
end
println()

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
using Tullio
cd(joinpath(dirname(pathof(Tullio)), "..", "benchmarks", "01"))
using Pkg; pkg"activate ."

# or

using Pkg; pkg"add LoopVectorization"
using Pkg; pkg"add TensorCast NNlib ForwardDiff Zygote"

# This example is a nice usage:
# https://discourse.julialang.org/t/help-to-improve-performance-of-gradient-calculation-on-tensor-operations/37773/3

using TensorCast

function cosinesim(a, b)
    @reduce similarity[i, j, k] := sum(s) a[s, j, k] * b[s, i, k] /
        sqrt( @reduce [_, j, k] := sum(s') a[s', j, k]^2) /
        sqrt( @reduce [_, i, k] := sum(s'') b[s'', i, k]^2)
end

function mysoftmax(a)
    @cast submax[i, j, k] := a[i, j, k] - @reduce [_, j, k] := maximum(i) a[i, j, k]
    @cast r[i, j, k] := exp(submax[i, j, k]) / @reduce [_, j, k] := sum(i) exp(submax[i, j, k])
end

pairwise2(a) = mysoftmax(cosinesim(a))

N, W, R, B = 16, 64, 4, 16;
a = rand(Float32, W, R, B);
b = rand(Float32, W, N, B);

#=
First look just at the numerator. This should be easy with `NNlib.batched_mul`,
but there is [this issue](https://github.com/FluxML/Zygote.jl/issues/552)
so I'm using [my PR](https://github.com/FluxML/NNlib.jl/pull/191)...
but the end result is a gradient 4x faster.

Almost as fast as Tullio! (which is faster without multi-threading today, something is broken)
=#

using Zygote, NNlib

reducemul(a, b) = @reduce z[i, j, k] := sum(s) a[s, j, k] * b[s, i, k]
reducemul(a, b) ≈ batched_mul(batched_transpose(b), a) # true

@btime reducemul($a, $b); # 98.262 μs (52 allocations: 342.41 KiB)
@btime batched_mul(batched_transpose($b), $a); # 13.605 μs (5 allocations: 4.20 KiB) -> 6.447 μs (2 allocations: 4.14 KiB) with MKL

using Tullio, LoopVectorization, ForwardDiff

tullio_mul(a, b) = @tullio z[i, j, k] := a[s, j, k] * b[s, i, k]
tullio_mul2(a, b) = @tullio z[i, j, k] := a[s, j, k] * b[s, i, k]  grad=Dual
reducemul(a, b) ≈ tullio_mul(a, b) # true

@btime tullio_mul($a, $b); # 2.696 μs (1 allocation: 4.13 KiB) -- improved!

#=
For gradients, I've inserted `(x->x).()` into these because sum creates a `FillArray`
which some methods don't like, but any useful calculation would have something else after this.

Now not needed for @tullio!
=#

grad_a = gradient((a,b) -> sum(reducemul(a, b)), a, b)[1];
grad_a ≈ gradient((a,b) -> sum((x->x), batched_mul(PermutedDimsArray(b, (2,1,3)), a)), a, b)[1] # true
grad_a ≈ gradient((a,b) -> sum(tullio_mul(a, b)), a, b)[1] # true
grad_a ≈ gradient((a,b) -> sum(tullio_mul2(a, b)), a, b)[1] # true

@btime gradient((a,b) -> sum(reducemul(a, b)), $a, $b);       # 307.836 μs (193 allocations: 1.17 MiB)
@btime gradient((a,b) -> sum(x->x, reducemul(a, b)), $a, $b); # 378.264 μs (3316 allocations: 1.23 MiB)

@btime gradient((a,b) -> sum(x->x, batched_mul(PermutedDimsArray(b, (2,1,3)), a)), $a, $b); # 101.782 μs (3159 allocations: 151.52 KiB) -> 218.180 μs (3334 allocations: 161.70 KiB) MKL?

@btime gradient((a,b) -> sum(tullio_mul(a, b)), $a, $b);       # 21.190 μs (41 allocations: 85.89 KiB)
@btime gradient((a,b) -> sum(x->x, tullio_mul(a, b)), $a, $b); # 79.291 μs (3156 allocations: 151.19 KiB)

# version with ForwardDiff:
@btime gradient((a,b) -> sum(tullio_mul2(a, b)), $a, $b);      # 25.168 μs (37 allocations: 85.70 KiB)

#=
Tullio ought to make it easy to fuse that multiplication with the division, perhaps?
But this is slow!
=#

function cosine_nnlib(a, b)
    @reduce den1[j, k] := sum(s) a[s, j, k]^2
    @reduce den2[i, k] := sum(s) b[s, i, k]^2
    bmm = batched_mul(PermutedDimsArray(b, (2,1,3)), a)
    @cast similarity[i, j, k] := bmm[i, j, k] / sqrt(den1[j, k] * den2[i, k])
end

function cosine_fused(a, b)
    @tullio den1[j, k] := a[s, j, k]^2
    @tullio den2[i, k] := b[s, i, k]^2
    @tullio similarity[i, j, k] := a[s, j, k] * b[s, i, k] / sqrt(den1[j, k] * den2[i, k])
end

function cosine_separated(a, b)
    @tullio bmm[i, j, k] := a[s, j, k] * b[s, i, k]

    @tullio den1[j, k] := a[s, j, k]^2
    @tullio den2[i, k] := b[s, i, k]^2

    @tullio f1[j, k] := 1/sqrt(den1[j, k])
    @tullio f2[i, k] := 1/sqrt(den2[i, k])
    @tullio similarity[i, j, k] := bmm[i, j, k] * f1[j, k] * f2[i, k]
end

cosinesim(a, b) ≈ cosine_nnlib(a, b) # true
cosinesim(a, b) ≈ cosine_fused(a, b) # true
cosinesim(a, b) ≈ cosine_separated(a, b) # true

@btime cosinesim($a, $b);    # 318.341 μs (87 allocations: 425.14 KiB)
@btime cosine_nnlib($a, $b); #  34.407 μs (42 allocations: 91.20 KiB) -> 112.308 μs (179 allocations: 99.05 KiB) MKL
@btime cosine_fused($a, $b);   #  7.678 μs (3 allocations: 5.59 KiB)
@btime cosine_separated($a, $b);   #  4.672 μs (6 allocations: 11.19 KiB)

grad_a2 = gradient((a,b) -> sum(cosinesim(a, b)), a, b)[1];
grad_a2 ≈ gradient((a,b) -> sum(cosine_nnlib(a, b)), a, b)[1] # true
grad_a2 ≈ gradient((a,b) -> sum(cosine_fused(a, b)), a, b)[1] # true
grad_a2 ≈ gradient((a,b) -> sum(cosine_separated(a, b)), a, b)[1] # true

@btime gradient((a,b) -> sum(cosinesim(a, b)), $a, $b);    #  984.351 μs (1414 allocations: 3.51 MiB)
@btime gradient((a,b) -> sum(cosine_nnlib(a, b)), $a, $b); #  232.991 μs (6375 allocations: 579.78 KiB) -> 280.876 μs (3442 allocations: 520.61 KiB) MKL
@btime gradient((a,b) -> sum(cosine_fused(a, b)), $a, $b);     # 80.207 μs (129 allocations: 252.86 KiB)
@btime gradient((a,b) -> sum(cosine_separated(a, b)), $a, $b); # 80.367 μs (207 allocations: 267.08 KiB)

# cosine_fused might end up computing / sqrt N^3 times, instead of 2N^2 times,
# so it's not obviously a great idea.
# But I like LoopVectorization is smart enough to avoid that? Gradient is as fast now.

# Its gradient now looks like this:

Tullio.@printgrad x/sqrt(y*z)   x y z

#=
Aside, things about softmax:
=#

mysoftmax(a) ≈ softmax(a, dims=1)

@btime mysoftmax($a);       # 77.517 μs (37 allocations: 50.17 KiB)
@btime softmax($a, dims=1); # 45.405 μs (8 allocations: 33.03 KiB)

@btime gradient(a -> sum(mysoftmax(a)), $a);      # 640.619 μs
@btime gradient(a -> sum(softmax(a, dims=1)), $a); # 95.622 μs
@btime gradient(a -> sum(x->x, mysoftmax(a)), $a);      #
@btime gradient(a -> sum(x->x, softmax(a, dims=1)), $a); # 322.103 μs (12358 allocations: 356.17 KiB)

# add avx to versions from https://github.com/FluxML/NNlib.jl/pull/135/
function NNlib.softmax(xs::Array; dims=1)
    max_ = maximum(xs, dims=dims)
    out = @avx exp.(xs .- max_)
    @avx out .= out ./ sum!(max_, out)
end

@btime softmax($a, dims=1); # 16.084 μs -- 3x quicker

function NNlib.∇softmax(Δ::Array, xs::Array; dims=1)
    sf = softmax(xs, dims=dims)
    @avx sf .* (Δ .- sum(Δ .* sf; dims=dims))
end

@btime gradient(a -> sum(x->x, softmax(a, dims=1)), $a); # 260.441 μs (12358 allocations: 323.36 KiB)
# also 3x quicker, once you subtract the big cost of broadcasting (x->x)!

using Tracker # quicker at broadcasting, not sure it has NNlib gradient at all

@btime Tracker.gradient(a -> sum(mysoftmax(a)), $a);       # 324.806 μs
@btime Tracker.gradient(a -> sum(softmax(a, dims=1)), $a); # 252.686 μs

using ReverseDiff

@btime ReverseDiff.gradient(a -> sum(mysoftmax(a)), $a); # 6.359 ms


text:
module SeisMain
    using Printf
    using Statistics
    using Distributed
    include("ReadWrite/ReadWrite.jl")
    include("Utils/Utils.jl")
end


text:
# this is a hack to allow Mass2 to depend on unregistered Julia Packages
# it should be removed if/when there is an official way of doing this other
# than maintaining and using my own copy of METADATA
if !isdir(Pkg.dir("ReferenceMicrocalFiles"))
	Pkg.clone("https://github.com/ggggggggg/ReferenceMicrocalFiles.jl")
end
isdir(expanduser("~/.daq")) || mkdir(expanduser("~/.daq"))
isfile(expanduser("~/.daq/latest_ljh_pulse.cur")) || touch(expanduser("~/.daq/latest_ljh_pulse.cur"))


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   ModiaMath.NonlinearEquations (ModiaMath/NonlinearEquations/_module.jl)
#

"""
    module KINSOL - Solve nonlinear equation system with Sundials KINSOL

The goal is to solve the same system several times with KINSOL.
"""
module KINSOL

import Sundials
import ModiaMath
using LinearAlgebra

@noinline function old_cfunction(f, r, a)
    ccall(:jl_function_ptr, Ptr{Cvoid}, (Any, Any, Any), f, r, a)
end

mutable struct NonlinearEquationsInfo
    extraInfo                # Model-specific extra information
    name::String             # Name of equation system
    ny::Int                  # Number of equations (length of y-vector)
    getResidues!::Function   # Function of the nonlinear equation system
    y0::Vector{Float64}      # The initial y vector (to print in error messages)
    lastNorm_r::Float64
    lastrScaledNorm_r::Float64
    kin_mem::Ptr{Cvoid}      # KINSOL pointer (to access all KINgetXXX functions)

    function NonlinearEquationsInfo(name::String, ny::Int, getResidues!::Function; extraInfo=nothing)
        @assert(ny >= 0)
        new(extraInfo, name, ny, getResidues!, zeros(ny), 1.0, 1.0)
    end
end

function kinsol_f(y::Sundials.N_Vector, r::Sundials.N_Vector, eqInfo::NonlinearEquationsInfo)
    eqInfo.getResidues!(eqInfo, Sundials.asarray(y), Sundials.asarray(r))
    return Cint(0)   # indicates normal return
end

@static if VERSION < v"0.7.0-DEV.2005"
    const kinsol_fc = cfunction(kinsol_f, Cint, (Sundials.N_Vector, Sundials.N_Vector, Ref{NonlinearEquationsInfo}))
end

function kinsol_ErrHandlerFn(error_code::Cint, KINmodule::Cstring, KINfunction::Cstring,
                             message::Cstring, eqInfo::NonlinearEquationsInfo)
    if error_code > 0
        # Print information text
        println("\n\n!!! Warning from ModiaMath.NonlinearEquations.KINSOL: ", unsafe_string(KINfunction),
               "(...) returned with a [", unsafe_string(KINmodule), "] error_code  = ",
               error_code, ":\n", unsafe_string(message), "\n")
    else
        simulationModel = eqInfo.extraInfo
        tables = ModiaMath.getVariableAndResidueValues(simulationModel)

        if tables != nothing
            println("\n\nLast used values in model:\n\n",
                  tables[1], "\n\n",
                  tables[2], "\n\n")
        end

        if error_code == -11
            str1 = "\nIt might be that the Jacobian is singular (= there are redundant equations).\n"
        else
            str1 = "\n"
        end

        if typeof(simulationModel) <: ModiaMath.AbstractSimulationModel
            simState = simulationModel.simulationState
            x_names     = String[simState.getVariableName(simState.model, ModiaMath.DAE.Category_X, i)  for i = 1:simState.nx]

            error("\n\n!!! Error from ModiaMath.NonlinearEquations.KINSOL: ", unsafe_string(KINfunction),
                  "(...) returned with a [", unsafe_string(KINmodule), "] error:\n    ",
                  unsafe_string(message), "\nModiaMath info:\nlastNorm(r) = ", eqInfo.lastNorm_r,
                  ", lastNorm(rScaled*r) = ", eqInfo.lastrScaledNorm_r, ".", str1,
                  string(simState.name), ": time = " , string(simState.time) ,
                 ", stepsize of implicit Euler step = " , string(simState.hev) ,
                 ", scaleConstraintsAtEvents = " , string(simState.scaleConstraintsAtEvents) ,
                 "\nx_names     = " , string(x_names) ,
                 "\nx_start     = " , string(eqInfo.y0) ,
                 "\nx_fixed     = " , string(simState.x_fixed) ,
                 "\nx_nominal   = " , string(simState.x_nominal) ,
                 "\nx_yScale    = " , string(simState.yScale) ,
                 "\nx_rScale    = " , string(simState.rScale) ,
                 "\nx           = " , string(simState.xev) ,
                 "\nderx        = " , string(simState.derxev) ,
                 "\nresidues    = " , string(simState.residues) ,
                 "\nnx          = " , string(simState.nx) ,
                 "\nnd          = " , string(simState.nd) ,
                 "\nnc          = " , string(simState.nc) ,
                 "\nnw          = " , string(simState.nw) ,
                 "\nnz          = " , string(simState.nz) ,
                 "\ntolerance   = " , string(simState.tolerance) ,
                 "\nFTOL        = " , string(simState.FTOL))

        else
            error("\n\n!!! Error from ModiaMath.NonlinearEquations.KINSOL: ", unsafe_string(KINfunction),
                  "(...) returned with a [", unsafe_string(KINmodule), "] error:\n    ",
                  unsafe_string(message), "\nModiaMath info:\nlastNorm(r) = ", eqInfo.lastNorm_r,
                  ", lastNorm(rScaled*r) = ", eqInfo.lastrScaledNorm_r, ".", str1)
        end
    end
    return nothing
end

@static if VERSION < v"0.7.0-DEV.2005"
    const kinsol_ErrHandlerFnc = cfunction(kinsol_ErrHandlerFn, Void, (Cint, Cstring, Cstring, Cstring, Ref{NonlinearEquationsInfo}))
end

function solveNonlinearEquations!(eqInfo::NonlinearEquationsInfo, y::Vector{Float64};
                                  FTOL::Float64=eps(Float64)^(1 / 3),
                                  yScale::Vector{Float64}=ones(length(y)),
                                  rScale::Vector{Float64}=ones(length(y)))
    # Create KINSOL and info structure
    @assert(length(y) == eqInfo.ny)
    kmem = Sundials.KINCreate()
    if kmem == C_NULL
        error("ModiaMath.NonlinearEquations.KINSOL.solveNonlinearEquations!: Failed to allocate KINSOL solver object")
    end
    eqInfo.kin_mem = kmem
    eqInfo.y0 .= y

    # Run KINSOL
    try
        # Set error handler function
        @static if VERSION >= v"0.7.0-DEV.2005"
            Sundials.KINSetErrHandlerFn(kmem, old_cfunction(kinsol_ErrHandlerFn, Nothing, Tuple{Cint,Cstring,Cstring,Cstring,Ref{typeof(NonlinearEquationsInfo)}}), pointer_from_objref(eqInfo))
        else
            Sundials.KINSetErrHandlerFn(kmem, kinsol_ErrHandlerFnc, pointer_from_objref(eqInfo))
        end

        # Initialize KINSOL
        @static if VERSION >= v"0.7.0-DEV.2005"
            Sundials.KINInit(kmem, old_cfunction(kinsol_f, Cint, Tuple{Sundials.N_Vector,Sundials.N_Vector,Ref{typeof(NonlinearEquationsInfo)}}), Sundials.NVector(y))
        else
            Sundials.KINInit(kmem, kinsol_fc, y)
        end
        Sundials.KINSetUserData(kmem, eqInfo)
        Sundials.KINSetFuncNormTol(kmem, FTOL)
        Sundials.KINSetScaledStepTol(kmem, FTOL * FTOL)

        # Set maximum allowable scaled length mxnewtstep of the Newton step
        # KINSOL defines a default of mxnewtstep = 1000*norm(y.*yScale).
        # This fails if y is a zero vector. Therefore, mxnewtstep is changed to
        # take only yScale into account.
        Sundials.KINSetMaxNewtonStep(kmem, 1000.0 * norm(yScale, Inf))

        # Set linear solver
        #   nnz_jac = nnz(sim.jac)
        #   IDAKLU(mem, ny, nnz_jac)
        #   IDAKLUSetOrdering(mem, KLUorderingChoice)
        #   IDASlsSetSparseJacFn(mem);
        #else
        A = Sundials.SUNDenseMatrix(length(y), length(y))
        LS = Sundials.SUNDenseLinearSolver(Sundials.NVector(y), A)
        Sundials.KINDlsSetLinearSolver(kmem, LS, A)
        #end
        strategy = Sundials.KIN_LINESEARCH

        @static if VERSION >= v"0.7.0-DEV.2005"
            Sundials.KINSol(kmem, Sundials.NVector(y), strategy, Sundials.NVector(yScale), Sundials.NVector(rScale))
        else
            Sundials.KINSol(kmem, y, strategy, yScale, rScale)
        end

    finally
        # Free allocated memory
        Sundials.KINFree([kmem])
    end

    return nothing
end

end


text:
# Robin Hood hash table, linear probing, linear search, backshift deletion

abstract type AbstractIndex{T, V} <: AbstractVector{T} end
mutable struct Index{T, V} <: AbstractIndex{T, V}
    vals::Vector{T}
    indices::Vector{V}
    probepositions::Vector{V}

    longestprobe::V
    initiallongestprobe::V
    deletions::UInt16

    function Index{T}(nelements::Integer) where {T}
        size = ceil(nelements / 0.9) # 0.9 load factor

        mintype = minimum_unsigned_type_for_n(nelements)
        return new{T, mintype}(
            Vector{T}(undef, nelements),
            zeros(mintype, mintype(size)),
            zeros(mintype, mintype(size)),
            0,
            0,
            0,
        )
    end
end

_length(idx::Index) = length(idx.indices)

function _setindex!(idx::Index{T, V}, elem::T, position::Unsigned) where {T, V}
    idx.vals[position] = elem
    probeposition = 0x0
    k = position
    location = hash(elem) % _length(idx)
    prevlongestprobe = idx.longestprobe
    @inbounds while k != 0x0
        probeposition += 0x1
        location = location == _length(idx) ? UInt64(0x1) : location + 0x1
        recordpos = idx.probepositions[location]
        if probeposition > recordpos
            k, idx.indices[location] = idx.indices[location], k
            idx.probepositions[location] = probeposition
            idx.longestprobe = max(idx.longestprobe, probeposition)
            probeposition = recordpos
        end
    end

    # reset longestprobe after some deletions deletions, this should maintain stable performance for _getindex and _getindex_array
    if idx.longestprobe == prevlongestprobe &&
       idx.longestprobe > idx.initiallongestprobe &&
       idx.deletions >= length(idx) ÷ 2
        idx.longestprobe = maximum(idx.probepositions)
        idx.deletions = 0x0
    end
end

function _getindex(idx::Index{T}, elem::T) where {T}
    ilength = _length(idx)
    if ilength > 0x0
        location = hash(elem) % ilength + 0x1
        for probeposition in 0x1:(idx.longestprobe)
            pos = idx.indices[location]
            if pos > 0x0 && isequal(idx.vals[pos], elem)
                return location
            elseif pos == 0x0
                return 0x0
            end
            location += 0x1
            if location > _length(idx)
                location = location % ilength
            end
        end
    end
    return 0x0
end

function _getindex_array(idx::Index{T, V}, elem::T) where {T, V}
    locations = Vector{V}()
    ilength = _length(idx)
    if ilength > 0x0
        location = hash(elem) % ilength + 0x1
        for probeposition in 0x1:(idx.longestprobe)
            pos = idx.indices[location]
            if pos > 0x0 && isequal(idx.vals[pos], elem)
                push!(locations, location)
            elseif pos == 0x0
                return locations
            end
            location += 0x1
            if location > ilength
                location = location % ilength
            end
        end
    end
    return locations
end

function _getindex_byposition(idx::Index{T}, i::Integer) where {T}
    elem = idx.vals[i]
    ilength = _length(idx)
    if ilength > 0x0
        location = hash(elem) % ilength + 0x1
        for probeposition in 0x1:(idx.longestprobe)
            pos = idx.indices[location]
            if pos > 0x0 && pos == i
                return location
            elseif pos == 0x0
                break
            end
            location += 0x1
            if location > ilength
                location = location % ilength
            end
        end
    end
    error("Element not found. This should never happen.")
end

function _delete!(idx::Index, oldkeyindex::Integer)
    lastidx = findnext(x -> x <= 0x1, idx.probepositions, oldkeyindex + 0x1)
    @inbounds if !isnothing(lastidx)
        idx.indices[oldkeyindex:(lastidx - 0x2)] .=
            @view idx.indices[(oldkeyindex + 0x1):(lastidx - 0x1)]
        idx.probepositions[oldkeyindex:(lastidx - 0x2)] .=
            @view(idx.probepositions[(oldkeyindex + 0x1):(lastidx - 0x1)]) .- 0x1
        idx.indices[lastidx - 0x1] = idx.probepositions[lastidx - 0x1] = 0x0
    else
        if oldkeyindex < _length(idx)
            idx.indices[oldkeyindex:(end - 0x1)] .= @view idx.indices[(oldkeyindex + 0x1):end]
            idx.probepositions[oldkeyindex:(end - 0x1)] .=
                @view(idx.probepositions[(oldkeyindex + 0x1):end]) .- 0x1
        end
        if idx.probepositions[1] > 0x1
            idx.indices[end] = idx.indices[1]
            idx.probepositions[end] = idx.probepositions[1] - 0x1

            lastidx = findnext(x -> x <= 0x1, idx.probepositions, 2)
            if isnothing(lastidx)
                lastidx = oldkeyindex
            end
            idx.indices[0x1:(lastidx - 0x2)] .= @view idx.indices[0x2:(lastidx - 0x1)]
            idx.probepositions[0x1:(lastidx - 0x2)] .=
                @view(idx.probepositions[0x2:(lastidx - 0x1)]) .- 0x1
            idx.indices[lastidx - 0x1] = idx.probepositions[lastidx - 0x1] = 0x0
        else
            idx.indices[end] = idx.probepositions[end] = 0x0
        end
    end

    idx.deletions += 0x1
end

function Index(elems::AbstractVector{T}) where {T}
    idx = Index{T}(length(elems))
    for (i, e) in enumerate(elems)
        _setindex!(idx, e, UInt(i))
    end
    idx.initiallongestprobe = idx.longestprobe
    return idx
end

Base.getindex(idx::AbstractIndex{T}, elem::T) where {T} = getindex(idx, elem, Val(false))
Base.getindex(idx::AbstractIndex{T}, elem::T, x::Bool) where {T} = getindex(idx, elem, Val(x))
Base.getindex(idx::AbstractIndex{T}, elem::T, x::Bool, y::Bool) where {T} =
    getindex(idx, elem, Val(x), Val(y))
Base.getindex(idx::AbstractIndex{T}, elems::AbstractVector{T}) where {T} =
    getindex(idx, elems, Val(false))
Base.getindex(idx::AbstractIndex{T}, elems::AbstractVector{T}, x::Bool) where {T} =
    getindex(idx, elems, Val(x))
Base.getindex(
    idx::AbstractIndex{T, V},
    elems::AbstractVector{T},
    x::Union{Val{true}, Val{false}},
) where {T, V} = isempty(elems) ? V[] : reduce(vcat, (getindex(idx, elem, x) for elem in elems))
Base.in(elem::T, idx::AbstractIndex{T}) where {T} =
    getindex(idx, elem, Val(false), Val(false)) != 0x0

Base.getindex(idx::Index{T}, elem::T, ::Val{true}, ::Val{false}) where {T} =
    @inbounds idx.indices[_getindex_array(idx, elem)] # exceptions may be undesirable in high-performance scenarios
function Base.getindex(idx::Index{T}, elem::T, ::Val{true}) where {T}
    d = _getindex_array(idx, elem)
    if length(d) == 0
        throw(KeyError(elem))
    end
    @inbounds return idx.indices[d]
end

function Base.getindex(idx::Index{T, V}, elem::T, ::Val{false}, ::Val{false}) where {T, V}
    i = _getindex(idx, elem)
    return i == 0x0 ? V(0x0) : @inbounds idx.indices[i]
end
function Base.getindex(idx::Index{T}, elem::T, ::Val{false}) where {T}
    i = getindex(idx, elem, Val(false), Val(false))
    if i == 0x0
        throw(KeyError(elem))
    end
    return i
end

function Base.getindex(idx::Index, i::Integer)
    return idx.vals[i]
end

function Base.setindex!(idx::Index{T}, newval::T, i::Integer) where {T}
    @boundscheck checkbounds(idx.vals, i)
    oldkeyindex = _getindex_byposition(idx, i)
    validx = idx.indices[oldkeyindex]
    _delete!(idx, oldkeyindex)
    _setindex!(idx, newval, validx)
    return idx
end

function Base.setindex!(idx::Index{T}, newval::T, oldval::T) where {T}
    oldkeyindex = _getindex(idx, oldval)
    if oldkeyindex == 0x0
        throw(KeyError(oldval))
    end
    validx = idx.indices[oldkeyindex]
    _delete!(idx, oldkeyindex)
    _setindex!(idx, newval, validx)
    return idx
end

Base.length(idx::Index) = length(idx.vals)
Base.size(idx::Index) = (length(idx),)
Base.values(idx::Index) = idx.vals
Base.convert(::Type{<:Index}, x::AbstractArray) = Index(x)
Base.convert(::Type{<:Index}, x::Index) = x

struct SubIndex{T, V, I} <: AbstractIndex{T, V}
    parent::Index{T, V}
    indices::I
    revmapping::Union{Nothing, Index}
end
SubIndex(idx::Index{T, V}, indices::I) where {T, V, I} = SubIndex{T, V, I}(idx, indices, nothing)

function Base.view(idx::Index, I::Union{AbstractRange, Colon})
    @boundscheck checkbounds(idx, I)
    return SubIndex(idx, I)
end
function Base.view(idx::Index{T, V}, I::AbstractArray{<:Integer}) where {T, V}
    @boundscheck checkbounds(idx, I)
    return SubIndex(idx, I, Index(V.(I)))
end
function Base.view(idx::Index{T, V}, I::AbstractArray{Bool}) where {T, V}
    @boundscheck checkbounds(idx, I)
    return SubIndex(idx, I, Index(V.(findall(I))))
end
function Base.view(idx::Index, I::Integer)
    @boundscheck checkbounds(idx, I)
    return @inbounds view(idx, I:I)
end
function Base.view(idx::SubIndex, I)
    @boundscheck checkbounds(idx, I)
    return @inbounds view(parent(idx), Base.reindex((parentindices(idx),), (I,))[1])
end

Base.copy(si::SubIndex) = Index(si)
Base.parent(si::SubIndex) = si.parent
Base.parentindices(si::SubIndex) = si.indices
Base.length(si::SubIndex) = length(parentindices(si))
Base.length(si::SubIndex{T, V, Colon}) where {T, V} = length(parent(si))
Base.length(si::SubIndex{T, V, I}) where {T, V, I <: AbstractArray{Bool}} = length(si.revmapping)
Base.size(si::SubIndex) = (length(si),)
Base.values(si::SubIndex) = parent(si)[parentindices(si)]
function Base.getindex(si::SubIndex{T}, elem::T, ::Val{true}, ::Val{false}) where {T}
    res = getindex(parent(si), elem, Val(true), Val(false))
    i = 0
    @inbounds for r in res
        position = findfirst(isequal(r), parentindices(si))
        if !isnothing(position)
            i += 1
            res[i] = position
        end
    end
    resize!(res, i)
    return res
end
function Base.getindex(
    si::SubIndex{T, V, I},
    elem::T,
    ::Val{true},
    ::Val{false},
) where {T, V, I <: AbstractArray{Bool}}
    res = getindex(parent(si), elem, Val(true), Val(false))
    res = res[parentindices(si)[res]]
    @inbounds for (i, r) in enumerate(res)
        res[i] = si.revmapping[r, false, false]
    end

    return res
end
function Base.getindex(
    si::SubIndex{T, V, I},
    elem::T,
    ::Val{true},
    ::Val{false},
) where {T, V, I <: AbstractArray{<:Integer}}
    res = getindex(parent(si), elem, Val(true), Val(false))
    i = 0
    @inbounds for r in res
        position = si.revmapping[r, false, false]
        if position > 0x0
            i += 1
            res[i] = position
        end
    end
    return resize!(res, i)
end
function Base.getindex(
    si::SubIndex{T, V, I},
    elem::T,
    ::Val{true},
    ::Val{false},
) where {T, V, I <: AbstractRange}
    res = getindex(parent(si), elem, Val(true), Val(false))
    j = 0
    for i in 1:length(res)
        pos = findfirst(==(res[i]), parentindices(si))
        if !isnothing(pos)
            j += 1
            res[j] = pos
        end
    end
    return resize!(res, j)
end
function Base.getindex(si::SubIndex{T}, elem::T, ::Val{true}) where {T}
    res = getindex(si, elem, Val(true), Val(false))
    if length(res) == 0
        throw(KeyError(elem))
    end
    return res
end

function Base.getindex(si::SubIndex{T, V}, elem::T, ::Val{false}, ::Val{false}) where {T, V}
    res = getindex(parent(si), elem, Val(false), Val(false))
    position = findfirst(isequal(res), parentindices(si))
    return isnothing(position) ? V(0x0) : V(position)
end
function Base.getindex(
    si::SubIndex{T, V, I},
    elem::T,
    ::Val{false},
    ::Val{false},
) where {T, V, I <: AbstractArray{Bool}}
    res = getindex(parent(si, elem, Val(false), Val(false)))
    if res == 0x0 || res > 0x0 && !(@inbounds parentindices(si)[res])
        return V(0x0)
    end
    return si.revmapping[res, false, false]
end
function Base.getindex(
    si::SubIndex{T, V, I},
    elem::T,
    ::Val{false},
    ::Val{false},
) where {T, V, I <: AbstractArray{<:Integer}}
    res = getindex(parent(si), elem, Val(false), Val(false))
    return res > 0x0 ? si.revmapping[res, false, false] : res
end
function Base.getindex(
    si::SubIndex{T, V, I},
    elem::T,
    ::Val{false},
    ::Val{false},
) where {T, V, I <: AbstractRange}
    res = getindex(parent(si), elem, Val(false), Val(false))
    if res > 0x0
        i = findfirst(==(res), parentindices(si))
        return isnothing(i) ? 0x0 : i
    else
        return 0x0
    end
end
function Base.getindex(si::SubIndex{T}, elem::T, ::Val{false}) where {T}
    res = getindex(si, elem, Val(false), Val(false))
    if res == 0x0
        throw(KeyError(elem))
    end
    return res
end
function Base.getindex(si::SubIndex, i::Union{Integer, AbstractVector{<:Integer}})
    @boundscheck checkbounds(si, i)
    return @inbounds parent(si)[Base.reindex((parentindices(si),), (i,))[1]]
end
function Base.getindex(si::SubIndex{T, V, Colon}, i::Union{Integer, AbstractVector{<:Integer}}) where {T, V}
    @boundscheck checkbounds(si, i)
    return @inbounds parent(si)[i]
end

function Base.setindex!(si::SubIndex{T}, newval::T, i::Integer) where {T}
    @boundscheck checkbounds(si, i)
    @inbounds setindex!(parent(si), newval, Base.reindex((parentindices(si),), (i,))[1])
    return si
end
function Base.setindex!(si::SubIndex{T, V, Colon}, newval::T, i::Integer) where {T, V}
    @boundscheck checkbounds(si, i)
    @inbounds setindex!(parent(si), newval, i)
end
function Base.setindex!(si::SubIndex{T}, newval::T, oldval::T) where {T}
    oldidx = parent(si)[oldval, true]
    foldidx = findfirst(in(parentindices(si)), oldidx)
    if isnothing(foldidx)
        throw(KeyError(oldval))
    end
    @inbounds parent(si)[oldidx[foldidx]] = newval
    return si
end
function Base.setindex!(
    si::SubIndex{T, V, I},
    newval::T,
    oldval::T,
) where {T, V, I <: AbstractArray{Bool}}
    oldidx = parent(si)[oldval, true]
    oldidx = @inbounds oldidx[parentindices(si)[oldix]]
    if length(oldidx) == 0
        throw(KeyError(oldval))
    end
    @inbounds parent(si)[oldidx[1]] = newval
    return si
end
function Base.setindex!(
    si::SubIndex{T, V, I},
    newval::T,
    oldval::T,
) where {T, V, I <: AbstractArray{<:Integer}}
    oldidx = parent(si)[oldval, true]
    foldidx = findfirst(in(si.revmapping), oldidx)
    if isnothing(foldidx)
        throw(KeyError(oldval))
    end
    @inbounds parent(si)[oldidx[foldidx]] = newval
    return si
end


text:
module CircleFinding

using CircleFit
using DelimitedFiles
using Plots
using StatsBase
using Random
using Statistics
using Distributions
using Optim
using ForwardDiff

function read_profiles(csv_file)
    a = readdlm(csv_file,',')
    b = map(v -> typeof(v) == Int64 ? v / 1000.0 : NaN,a)[:,2:end]
    x = b[:,1:2:end]
    y = b[:,2:2:end]
    return (x,y)
end

function read_apex_profile(csv_file)
    a = readdlm(csv_file,',')
    x = (a[:,1]./1000.0)'
    y = (a[:,2]./1000.0)'
    return (x,y)
end

function read_single_profile(csv_file)
    a = readdlm(csv_file,',')
    x = a[:,2:2:end]./1000.0
    y = a[:,3:2:end]./1000.0
end

"""

"""
function circle_ransac(data, radius; threshold = 0.4, max_iters = 2000, number_points = 8)
    data_length = size(data)[2]
    best_xc = 0
    best_yc = 0
    best_rfit = 0;
    best_err = 10000000;
    for i = 1:max_iters
        ri = randperm(data_length)
        maybe_in = data[:,ri[1:3]]
        xc,yc,rfit = circfit(maybe_in)
        also_in = hcat(filter(xy -> abs(sqrt((xy[1]-xc)^2 + (xy[2]-yc)^2) - radius) < threshold, collect(eachcol(data[:,ri[4:end]])))...)

        if size(also_in)[1] == 2 && size(also_in)[2] > number_points && abs(rfit - radius) < threshold
            current_data = hcat(maybe_in,also_in)
            better_xc, better_yc, better_rfit = circfit(current_data)

            err = circle_error(current_data,better_xc,better_yc,radius)

            if err < best_err
                best_xc = better_xc
                best_yc = better_yc
                best_rfit = better_rfit
                best_err = err
            end
        end
    end
    return (best_xc, best_yc, best_rfit, best_err)
end

circfit(xy) = circfit(xy[1,:],xy[2,:])

circle_error(data, xc, yc, radius) = sqrt(mean((radius .- sqrt.(sum((data .- [xc;yc]).^2,dims=1))).^2))

function filterdata(x,y)
    data = vcat(x[:]',y[:]')
    return hcat(filter(xy -> isfinite(xy[1]) && isfinite(xy[2]),collect(eachcol(data)))...)
end

function fit_circle(x1,y1,x2,y2,radius)
    x3 = (x1+x2) ./ 2
    y3 = (y1+y2) ./ 2
    q = sqrt((x1-x2)^2 + (y1-y2)^2)

    xc = x3 - sqrt(radius^2 - (q/2)^2) * (y1 - y2)/q
    yc = y3 - sqrt(radius^2 - (q/2)^2) * (x2 - x2)/q
    return(xc,yc)
end

function joe_circle(data,radius; xbins = (-15,15,500), ybins = (-30, 30, 1000))
    bins = zeros(ybins[3],xbins[3])

    d1 = data[:,1:end-1]
    d2 = data[:,2:end]

end

function circle_hough_map(data,radius; xlims = (-15,15), ylims = (-30, 30), step_size=0.05)

    bx = range(xlims...,step=step_size)
    len_x = length(bx)
    by = range(ylims...,step=step_size)
    len_y = length(by)

    bins = zeros(len_y,len_x)

    #d = Normal(radius, step_size)
    d = SymTriangularDist(radius,step_size)
    for xy = eachcol(data)
        xstart = findindex(xy[1] - radius - step_size, xlims[1], step_size, len_x)
        xend = findindex(xy[1] + radius + step_size, xlims[1], step_size, len_x)
        ystart = findindex(xy[2] - radius - step_size, ylims[1], step_size, len_y)
        yend = findindex(xy[2] + step_size, ylims[1], step_size, len_y)
        bins[ystart:yend,xstart:xend] .+= [pdf(d, sqrt((xy[1] - x)^2 + (xy[2] - y)^2)) for y = by[ystart:yend], x = bx[xstart:xend]]
    end
    return bins
end

function circle_hough(data,radius; xlims = (-15,15), ylims = (-30, 30), step_size = 0.05)

    bx = range(xlims...,step=step_size)
    by = range(ylims...,step=step_size)

    bins = circle_hough_map(data, radius, xlims=xlims, ylims=ylims, step_size=step_size)
    max_tuple = findmax(bins)
    return (bx[max_tuple[2][2]], by[max_tuple[2][1]],max_tuple[1])
end

function circle_hybrid(data,radius)
    x,y,w = circle_hough(data,radius)

    good_data = hcat(filter(xy -> radius-0.1 < sqrt( (xy[1] - x)^2 + (xy[2] - y)^2) < radius+0.1,collect(eachcol(data)))...)'
    cir_fit_err(p, x, y) = sum((sqrt.((x .- p[1]).^2 .+ (y .- p[2]).^2) .- p[3]).^2)

    res = optimize(b -> cir_fit_err(b,good_data[:,1],good_data[:,2]),[x,y,radius], LBFGS(), autodiff=:forward)
    return vcat(Optim.minimizer(res),length(good_data),Optim.minimum(res))
end

function findindex(x, initial, step_size, count)
    index = round(Integer, (x - initial) / step_size) + 1
    return clamp(index, 1, count)
end

end

max_stars_repo_name: joeynelson/CircleFinding

max_issues_repo_name: joeynelson/CircleFinding

max_forks_repo_name: joeynelson/CircleFinding


text:
import Statistics: mean
import Base.Threads: @threads, @sync, @spawn, nthreads, threadid

export CloudInCell

struct CloudInCell end

export compute_rho

function compute_rho(p, kernel::CloudInCell, m :: TwoDGrid)

    nx, ny = m.nx, m.ny
    dx, dy = m.dx, m.dy
    rho = zeros(nx+1, ny+1)
    nbpart = size(p.array)[2]

    @inbounds for ipart = 1:nbpart
        xp = p.array[1, ipart] / dx
        yp = p.array[2, ipart] / dy

        i = trunc(Int, xp) + 1
        j = trunc(Int, yp) + 1

        dxp = xp - i + 1
        dyp = yp - j + 1

        a1 = (1 - dxp) * (1 - dyp)
        a2 = dxp * (1 - dyp)
        a3 = dxp * dyp
        a4 = (1 - dxp) * dyp

        ip1 = mod1(i + 1, nx)
        jp1 = mod1(j + 1, ny)

        w = p.array[5, ipart]

        rho[i, j] += a1 * w
        rho[ip1, j] += a2 * w
        rho[ip1, jp1] += a3 * w
        rho[i, jp1] += a4 * w

    end

    for i = 1:nx
       rho[i,1] += rho[i,ny+1]
    end
    for j = 1:ny
       rho[1,j] += rho[nx+1,j]
    end

    rho ./= ( dx * dy )

    return rho[1:nx,1:ny] .- mean(rho[1:nx,1:ny])

end

export compute_current!

function compute_current!(jx, jy, m::TwoDGrid, kernel::CloudInCell, p)

    nbpart = size(p.array)[2]
    nx, ny = m.nx, m.ny
    dx, dy = m.dx, m.dy

    fill!(jx, 0)
    fill!(jy, 0)

    scaling = 1 / ( dx * dy )

    ntid = nthreads()
    jxloc = [zero(jx) for _ in 1:ntid]
    jyloc = [zero(jy) for _ in 1:ntid]

    chunks = Iterators.partition(1:nbpart, nbpart÷ntid)

    @sync for chunk in chunks
        @spawn begin
            tid = threadid()
            @inbounds for ipart in chunk
                xp = p.array[1, ipart] / dx
                yp = p.array[2, ipart] / dy

                i = floor(Int, xp) + 1
                j = floor(Int, yp) + 1

                dxp = xp - i + 1
                dyp = yp - j + 1

                a1 = (1 - dxp) * (1 - dyp)
                a2 = dxp * (1 - dyp)
                a3 = dxp * dyp
                a4 = (1 - dxp) * dyp

                w = p.array[5, ipart]

                w1 = p.array[3, ipart] * scaling * w
                w2 = p.array[4, ipart] * scaling * w

                jxloc[tid][i, j] += a1 * w1
                jyloc[tid][i, j] += a1 * w2

                jxloc[tid][i+1, j] += a2 * w1
                jyloc[tid][i+1, j] += a2 * w2

                jxloc[tid][i+1, j+1] += a3 * w1
                jyloc[tid][i+1, j+1] += a3 * w2

                jxloc[tid][i, j+1] += a4 * w1
                jyloc[tid][i, j+1] += a4 * w2
            end
        end
    end

    jx .= reduce(+,jxloc)
    jy .= reduce(+,jyloc)

    for i = 1:nx+1
        jx[i, 1] += jx[i, ny+1]
        jx[i, ny+1] = jx[i, 1]
        jy[i, 1] += jy[i, ny+1]
        jy[i, ny+1] = jy[i, 1]
    end
    for j = 1:ny+1
        jx[1, j] += jx[nx+1, j]
        jx[nx+1, j] = jx[1, j]
        jy[1, j] += jy[nx+1, j]
        jy[nx+1, j] = jy[1, j]
    end

end


text:
using AdventOfCode2018
using AdventOfCode2018.Day9

input = get_aoc_input(9)
players, max_points = parse_input(input)

println("Solution part 1: ", part_1(players, max_points))
println("Solution part 1: ", part_1(players, max_points * 100))


text:
s₁ = Scenario(q₁ = -24.0, q₂ = -28.0, d₁ = 500.0, d₂ = 100.0, probability = 0.4)
s₂ = Scenario(q₁ = -28.0, q₂ = -32.0, d₁ = 300.0, d₂ = 300.0, probability = 0.6)
deferred = StochasticProgram([s₁,s₂], solver=GLPKSolverLP())

@first_stage deferred = begin
    @variable(model, x₁ >= 40)
    @variable(model, x₂ >= 20)
    @objective(model, Min, 100*x₁ + 150*x₂)
    @constraint(model, x₁+x₂ <= 120)
end defer

@second_stage deferred = begin
    @decision x₁ x₂
    @uncertain q₁ q₂ d₁ d₂
    @variable(model, 0 <= y₁ <= d₁)
    @variable(model, 0 <= y₂ <= d₂)
    @objective(model, Min, q₁*y₁ + q₂*y₂)
    @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
    @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
end defer


text:
using NearestNeighbors, Random, PyPlot, BenchmarkTools
Random.seed!(12)
# make a function to populate with nearest neighbour values
function getvoronoi!(cellvalues, kdtree, gridpoints, gridvalues)
    idxs, = knn(kdtree, gridpoints, 1)
    for i in 1:length(gridvalues)
        gridvalues[i] = cellvalues[idxs[i][1]]
    end
end
## run one 2D example
nrows, ncols = 100, 100
λ_rows, λ_cols = 1, 10.
ncells = 20
xx, yy = [i for i = 1:nrows, j = 1:ncols], [j for i = 1:nrows, j = 1:ncols]
celllocations = rand(2, ncells).*[nrows/λ_rows, ncols/λ_cols]
cellvalues = rand(ncells, 1)
kdtree = KDTree(celllocations)
gridpoints = [xx[:]/λ_rows yy[:]/λ_cols]'
gridvalues = zeros(size(xx))
getvoronoi!(cellvalues, kdtree, gridpoints, gridvalues)
figure()
imshow(gridvalues, extent=[1,ncols,nrows,1])
scatter(celllocations[2,:]*λ_cols, celllocations[1,:]*λ_rows, color="k")
xlim(extrema(xx))
ylim(extrema(yy))
## now try 3D
nrows, ncols, nnext = 80, 90, 100
λ_rows, λ_cols, λ_next = 50, 10., 2
ncells = 50
xxx, yyy, zzz = [i for i = 1:nrows, j = 1:ncols, k = 1:nnext],
                [j for i = 1:nrows, j = 1:ncols, k = 1:nnext],
                [k for i = 1:nrows, j = 1:ncols, k = 1:nnext]
celllocations = rand(3, ncells).*[nrows/λ_rows, ncols/λ_cols, nnext/λ_next]
cellvalues = rand(ncells, 1)
kdtree = KDTree(celllocations)
gridpoints = [xxx[:]/λ_rows yyy[:]/λ_cols zzz[:]/λ_next]'
gridvalues = zeros(size(xxx))
getvoronoi!(cellvalues, kdtree, gridpoints, gridvalues)
## make slices
ix, iy, iz = 50, 50, 50
fig = figure(figsize=(10,3))
ax = fig.add_subplot(1, 3, 1, projection="3d")
x, y, z = 1:nrows, 1:ncols, 1:nnext
yy, zz = [j for j = 1:ncols, k = 1:nnext],
        [k for j = 1:ncols, k = 1:nnext]
c=(gridvalues[ix,:,:].-minimum(gridvalues))./(maximum(gridvalues)-minimum(gridvalues))
ax.plot_surface(x[ix]*ones(size(yy)), yy, zz, facecolors=plt.cm.jet_r(c), shade=false)

ax = fig.add_subplot(1, 3, 2, projection="3d")
xx, zz = [i for i = 1:nrows, k = 1:nnext],
         [k for i = 1:nrows, k = 1:nnext]
c=(gridvalues[:,iy,:].-minimum(gridvalues))./(maximum(gridvalues)-minimum(gridvalues))
ax.plot_surface(xx, y[iy]*ones(size(zz)), zz, facecolors=plt.cm.jet_r(c), shade=false)

xx, yy = [i for i = 1:nrows, j = 1:ncols],
         [j for i = 1:nrows, j = 1:ncols]
ax = fig.add_subplot(1, 3, 3, projection="3d")
c=(gridvalues[:,:,iz].-minimum(gridvalues))./(maximum(gridvalues)-minimum(gridvalues))
ax.plot_surface(xx, yy, z[iz]*ones(size(yy)), facecolors=plt.cm.jet_r(c), shade=false)

for a in gcf().axes

        # a.scatter3D(celllocations[1,:][:]*λ_rows,
        #             celllocations[2,:][:]*λ_cols,
        #             celllocations[3,:][:]*λ_next,
        #             c=cellvalues[:], vmin=minimum(gridvalues), vmax=maximum(gridvalues),
        #             s=50, cmap="jet_r")
    a.set_xlim(extrema(x))
    a.set_ylim(extrema(y))
    a.set_zlim(extrema(z))
    a.set_xlabel("x")
    a.set_ylabel("y")
    a.set_zlabel("z")
    a.invert_zaxis()
end


text:
using Test
using AsyncFinalizers

include("utils.jl")

function AsyncFinalizers.Internal.chaos_run_finalizers()
    error("!!! CHAOS !!!")
end

FALLBACK_DONE = Channel{Nothing}(1)

function AsyncFinalizers.Internal.queue_fallback_done()
    put!(FALLBACK_DONE, nothing)
end

AsyncFinalizers.Internal.reinit(; nresets = 3)

# Run an async finalizer to hit the fallback code path:
Utils.outlined() do
    AsyncFinalizers.onfinalize(Ref(0)) do _
        function () end
    end
    return
end
GC.gc()
take!(FALLBACK_DONE)

nfinalizers = 5
FINALIZER_CALLED = zeros(Int, nfinalizers)

for i in 1:nfinalizers
    @info "Testing $i-th finalizer"
    function async_finalizer(_)
        FINALIZER_CALLED[i] += 1
    end
    Utils.collect_garbage(; async_finalizer = async_finalizer, check_executor = false)
end

@test FINALIZER_CALLED == ones(Int, nfinalizers)


text:
#TODO make it a function on ExplicitODE and Options

"""

Chooses an initial step-size basing on the equation, initial data,
time span and the order of the method of integration.

"""
function dtinit{T,S}(F, y0::Vector{S}, tspan::Vector{T}, reltol, abstol; order = 1)
    t0 = abs(tspan[1])
    tstop = abs(tspan[end])
    tau = max(reltol*norm(y0, Inf), abstol)
    d0 = norm(y0, Inf)/tau
    f0 = F(t0, y0)
    d1 = norm(f0, Inf)/tau
    if min(d0,d1) < eps(T)^(1/3)
        dt0 = eps(T)^(1/3)/10
    else
        dt0 = (d0/d1)/100
    end
    # perform Euler step
    y1 = similar(y0)
    for d = 1:length(y1)
        y1[d] = y0[d]+dt0*f0[d]
    end
    f1 = F(t0 + dt0, y1)
    # estimate second derivative
    d2 = norm(f1 - f0, Inf)/(tau*dt0)
    if max(d1, d2) <= 10*eps(T)
        dt1 = max(eps(T)^(1/3)/10, dt0/10^3)
    else
        pow = -(2 + log10(max(d1, d2)))/(order+1)
        dt1 = 10^pow
    end
    return min(100*dt0, dt1, abs(tstop-t0))
end

# a scalar version of the above
dtinit(F, y0::Number, args...; kargs...) = dtinit((t,y)->[F(t,y[1])], [y0], args...; kargs...)

"""

A simple bisection algorithm for finding a root of a solution f(x)=0
starting within the range x∈rng, the result is a point x₀ which is
located within the distance eps from the true root of f(x)=0.  For
this algorithm to work we need f(rng[1]) to have a different sign then
f(rng[2]).

"""
function findroot(f,rng,eps)
    xl, xr = rng
    fl, fr = f(xl), f(xr)

    if fl*fr > 0 || xl > xr
        error("Inconsistent bracket")
    end

    while xr-xl > eps
        xm = (xl+xr)/2
        fm = f(xm)

        if fm*fr > 0
            xr = xm
            fr = fm
        else
            xl = xm
            fl = fm
        end
    end

    return (xr+xl)/2
end