text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)
module SQLite

using DataStreams, DataFrames, WeakRefStrings, LegacyStrings
import LegacyStrings: UTF16String

export Data, DataFrame

if !isdefined(Core, :String)
    typealias String UTF8String
end

if Base.VERSION < v"0.5.0-dev+4631"
    unsafe_string = bytestring
    unsafe_wrap{A<:Array}(::Type{A}, ptr, len) = pointer_to_array(ptr, len)
    unsafe_wrap(::Type{String}, ptr, len) = unsafe_string(ptr, len)
end

type SQLiteException <: Exception
    msg::AbstractString
end

include("consts.jl")
include("api.jl")

immutable NullType end
# custom NULL value for interacting with the SQLite database
const NULL = NullType()
show(io::IO,::NullType) = print(io,"#NULL")

#TODO: Support sqlite3_open_v2
# Normal constructor from filename
sqliteopen(file,handle) = sqlite3_open(file,handle)
sqliteopen(file::UTF16String,handle) = sqlite3_open16(file,handle)
sqliteerror() = throw(SQLiteException(unsafe_string(sqlite3_errmsg())))
sqliteerror(db) = throw(SQLiteException(unsafe_string(sqlite3_errmsg(db.handle))))

"""
represents an SQLite database, either backed by an on-disk file or in-memory

Constructors:

* `SQLite.DB()` => in-memory SQLite database
* `SQLite.DB(file)` => file-based SQLite database
"""
type DB
    file::String
    handle::Ptr{Void}
    changes::Int

    function DB(f::AbstractString)
        handle = Ref{Ptr{Void}}()
        f = isempty(f) ? f : expanduser(f)
        if @OK sqliteopen(f,handle)
            db = new(f,handle[],0)
            register(db, regexp, nargs=2, name="regexp")
            finalizer(db, _close)
            return db
        else # error
            sqlite3_close(handle[])
            sqliteerror()
        end
    end
end
DB() = DB(":memory:")

function _close(db::DB)
    sqlite3_close_v2(db.handle)
    db.handle = C_NULL
    return
end

Base.show(io::IO, db::SQLite.DB) = print(io, string("SQLite.DB(",db.file == ":memory:" ? "in-memory" : "\"$(db.file)\"",")"))

"""
`SQLite.Stmt(db::DB, sql::AbstractString)` creates and prepares an SQLite statement
"""
type Stmt
    db::DB
    handle::Ptr{Void}

    function Stmt(db::DB,sql::AbstractString)
        handle = Ref{Ptr{Void}}()
        sqliteprepare(db,sql,handle,Ref{Ptr{Void}}())
        stmt = new(db,handle[])
        finalizer(stmt, _close)
        return stmt
    end
end

function _close(stmt::Stmt)
    sqlite3_finalize(stmt.handle)
    stmt.handle = C_NULL
    return
end

sqliteprepare(db,sql,stmt,null) = @CHECK db sqlite3_prepare_v2(db.handle,sql,stmt,null)

include("UDF.jl")
export @sr_str, @register, register

"""
`SQLite.bind!(stmt::SQLite.Stmt, values)`

bind `values` to parameters in a prepared SQL statement. Values can be:

* `Vector`; where each element will be bound to an SQL parameter by index order
* `Dict`; where dict values will be bound to named SQL parameters by the dict key

Additional methods exist for working individual SQL parameters:

* `SQLite.bind!(stmt, name, val)`: bind a single value to a named SQL parameter
* `SQLite.bind!(stmt, index, val)`: bind a single value to a SQL parameter by index number
"""
function bind! end

function bind!(stmt::Stmt, values::Vector)
    nparams = sqlite3_bind_parameter_count(stmt.handle)
    @assert nparams == length(values) "you must provide values for all placeholders"
    for i in 1:nparams
        @inbounds bind!(stmt, i, values[i])
    end
end
function bind!{V}(stmt::Stmt, values::Dict{Symbol, V})
    nparams = sqlite3_bind_parameter_count(stmt.handle)
    @assert nparams == length(values) "you must provide values for all placeholders"
    for i in 1:nparams
        name = unsafe_string(sqlite3_bind_parameter_name(stmt.handle, i))
        @assert !isempty(name) "nameless parameters should be passed as a Vector"
        # name is returned with the ':', '@' or '$' at the start
        name = name[2:end]
        bind!(stmt, i, values[Symbol(name)])
    end
end
# Binding parameters to SQL statements
function bind!(stmt::Stmt,name::AbstractString,val)
    i::Int = sqlite3_bind_parameter_index(stmt.handle,name)
    if i == 0
        throw(SQLiteException("SQL parameter $name not found in $stmt"))
    end
    return bind!(stmt,i,val)
end
bind!(stmt::Stmt,i::Int,val::AbstractFloat)  = (sqlite3_bind_double(stmt.handle,i,Float64(val)); return nothing)
bind!(stmt::Stmt,i::Int,val::Int32)          = (sqlite3_bind_int(stmt.handle,i,val); return nothing)
bind!(stmt::Stmt,i::Int,val::Int64)          = (sqlite3_bind_int64(stmt.handle,i,val); return nothing)
bind!(stmt::Stmt,i::Int,val::NullType)       = (sqlite3_bind_null(stmt.handle,i); return nothing)
bind!(stmt::Stmt,i::Int,val::AbstractString) = (sqlite3_bind_text(stmt.handle,i,val); return nothing)
bind!(stmt::Stmt,i::Int,val::WeakRefString{UInt8})   = (sqlite3_bind_text(stmt.handle,i,val.ptr,val.len); return nothing)
bind!(stmt::Stmt,i::Int,val::WeakRefString{UInt16})  = (sqlite3_bind_text16(stmt.handle,i,val.ptr,val.len*2); return nothing)
bind!(stmt::Stmt,i::Int,val::UTF16String)    = (sqlite3_bind_text16(stmt.handle,i,val); return nothing)
function bind!(stmt::Stmt,i::Int,val::WeakRefString{UInt32})
    A = UTF32String(pointer_to_array(val.ptr, val.len+1, false))
    return bind!(stmt, i, convert(String,A))
end
# We may want to track the new ByteVec type proposed at https://github.com/JuliaLang/julia/pull/8964
# as the "official" bytes type instead of Vector{UInt8}
bind!(stmt::Stmt,i::Int,val::Vector{UInt8})  = (sqlite3_bind_blob(stmt.handle,i,val); return nothing)
# Fallback is BLOB and defaults to serializing the julia value

# internal wrapper type to, in-effect, mark something which has been serialized
immutable Serialization
    object
end

function sqlserialize(x)
    t = IOBuffer()
    # deserialize will sometimes return a random object when called on an array
    # which has not been previously serialized, we can use this type to check
    # that the array has been serialized
    s = Serialization(x)
    serialize(t,s)
    return takebuf_array(t)
end
# fallback method to bind arbitrary julia `val` to the parameter at index `i` (object is serialized)
bind!(stmt::Stmt,i::Int,val) = bind!(stmt,i,sqlserialize(val))

# magic bytes that indicate that a value is in fact a serialized julia value, instead of just a byte vector
const SERIALIZATION = UInt8[0x11,0x01,0x02,0x0d,0x53,0x65,0x72,0x69,0x61,0x6c,0x69,0x7a,0x61,0x74,0x69,0x6f,0x6e,0x23]
function sqldeserialize(r)
    ret = ccall(:memcmp, Int32, (Ptr{UInt8},Ptr{UInt8}, UInt),
            SERIALIZATION, r, min(18,length(r)))
    if ret == 0
        v = deserialize(IOBuffer(r))
        return v.object
    else
        return r
    end
end
#TODO:
 #int sqlite3_bind_zeroblob(sqlite3_stmt*, int, int n);
 #int sqlite3_bind_value(sqlite3_stmt*, int, const sqlite3_value*);

"""
`SQLite.execute!(stmt::SQLite.Stmt)` => `Void`
`SQLite.execute!(db::DB, sql::String)` => `Void`

Execute a prepared SQLite statement, not checking for or returning any results.
"""
function execute! end

function execute!(stmt::Stmt)
    r = sqlite3_step(stmt.handle)
    if r == SQLITE_DONE
        sqlite3_reset(stmt.handle)
    elseif r != SQLITE_ROW
        sqliteerror(stmt.db)
    end
    return r
end
function execute!(db::DB,sql::AbstractString)
    stmt = Stmt(db,sql)
    return execute!(stmt)
end

"""
`SQLite.esc_id(x::Union{String,Vector{String}})`

Escape SQLite identifiers (e.g. column, table or index names). Can be either
a string, or a vector of strings (note does not check for null characters).
A vector of identifiers will be separated by commas.
"""
function esc_id end

esc_id(x::AbstractString) = "\""*replace(x,"\"","\"\"")*"\""
esc_id{S<:AbstractString}(X::AbstractVector{S}) = join(map(esc_id,X),',')

# Transaction-based commands
"""
`SQLite.transaction(db, mode="DEFERRED")`
`SQLite.transaction(func, db)`

Begin a transaction in the specified `mode`, default = "DEFERRED".

If `mode` is one of "", "DEFERRED", "IMMEDIATE" or "EXCLUSIVE" then a
transaction of that (or the default) type is started. Otherwise a savepoint
is created whose name is `mode` converted to AbstractString.

In the second method, `func` is executed within a transaction (the transaction being committed upon successful execution)
"""
function transaction end

function transaction(db, mode="DEFERRED")
    execute!(db,"PRAGMA temp_store=MEMORY;")
    if uppercase(mode) in ["", "DEFERRED", "IMMEDIATE", "EXCLUSIVE"]
        execute!(db, "BEGIN $(mode) TRANSACTION;")
    else
        execute!(db, "SAVEPOINT $(mode);")
    end
end
function transaction(f::Function, db)
    # generate a random name for the savepoint
    name = string("SQLITE", randstring(10))
    execute!(db, "PRAGMA synchronous = OFF;")
    transaction(db, name)
    try
        f()
    catch
        rollback(db, name)
        rethrow()
    finally
        # savepoints are not released on rollback
        commit(db, name)
        execute!(db, "PRAGMA synchronous = ON;")
    end
end

"""
`SQLite.commit(db)`
`SQLite.commit(db, name)`

commit a transaction or named savepoint
"""
function commit end

commit(db) = execute!(db, "COMMIT TRANSACTION;")
commit(db, name) = execute!(db, "RELEASE SAVEPOINT $(name);")

"""
`SQLite.rollback(db)`
`SQLite.rollback(db, name)`

rollback transaction or named savepoint
"""
function rollback end

rollback(db) = execute!(db, "ROLLBACK TRANSACTION;")
rollback(db, name) = execute!(db, "ROLLBACK TRANSACTION TO SAVEPOINT $(name);")

"""
`SQLite.drop!(db, table; ifexists::Bool=true)`

drop the SQLite table `table` from the database `db`; `ifexists=true` will prevent an error being thrown if `table` doesn't exist
"""
function drop!(db::DB, table::AbstractString; ifexists::Bool=false)
    exists = ifexists ? "IF EXISTS" : ""
    transaction(db) do
        execute!(db,"DROP TABLE $exists $(esc_id(table))")
    end
    execute!(db,"VACUUM")
    return
end

"""
`SQLite.dropindex!(db, index; ifexists::Bool=true)`

drop the SQLite index `index` from the database `db`; `ifexists=true` will not return an error if `index` doesn't exist
"""
function dropindex!(db::DB,index::AbstractString;ifexists::Bool=false)
    exists = ifexists ? "IF EXISTS" : ""
    transaction(db) do
        execute!(db,"DROP INDEX $exists $(esc_id(index))")
    end
    return
end

"""
`SQLite.createindex!(db, table, index, cols; unique::Bool=true, ifnotexists::Bool=false)`

create the SQLite index `index` on the table `table` using `cols`, which may be a single column or vector of columns.
`unique` specifies whether the index will be unique or not.
`ifnotexists=true` will not throw an error if the index already exists
"""
function createindex!{S<:AbstractString}(db::DB,table::AbstractString,index::AbstractString,cols::Union{S,AbstractVector{S}}
                      ;unique::Bool=true,ifnotexists::Bool=false)
    u = unique ? "UNIQUE" : ""
    exists = ifnotexists ? "IF NOT EXISTS" : ""
    transaction(db) do
        execute!(db,"CREATE $u INDEX $exists $(esc_id(index)) ON $(esc_id(table)) ($(esc_id(cols)))")
    end
    execute!(db,"ANALYZE $index")
    return
end

"""
`SQLite.removeduplicates!(db, table, cols::Vector)`

removes duplicate rows from `table` based on the values in `cols` which is an array of column names
"""
function removeduplicates!{T <: AbstractString}(db,table::AbstractString,cols::AbstractArray{T})
    colsstr = ""
    for c in cols
       colsstr = colsstr * esc_id(c) * ","
    end
    colsstr = chop(colsstr)
    transaction(db) do
        execute!(db,"DELETE FROM $(esc_id(table)) WHERE _ROWID_ NOT IN (SELECT max(_ROWID_) from $(esc_id(table)) GROUP BY $(colsstr));")
    end
    execute!(db,"ANALYZE $table")
    return
 end

"`SQLite.Source` implements the `Source` interface in the `DataStreams` framework"
type Source <: Data.Source
    schema::Data.Schema
    stmt::Stmt
    status::Cint
end

"SQLite.Sink implements the `Sink` interface in the `DataStreams` framework"
type Sink <: Data.Sink # <: IO
    schema::Data.Schema
    db::DB
    tablename::String
    stmt::Stmt
    transaction::String
end

include("Source.jl")
include("Sink.jl")

end # module


text:
module IGA

# import dependent packages
using FastGaussQuadrature, PyCall, PyPlot

# export types
export Degree, KnotVector, BezierPatch, NURBS, QuadRule

# export functionality from other libs
export ndgrid, allcomb, linearspace, boundary
export buildvector, uniqueknots, globalrefine, findspan
export dimsplinespace, grevillepoints, onebasisfuneval, bsplinebasiseval, dersbsplinebasisfuns, dersbsplinebasisfunsinterpolationmatrix
export knotinsertionoperator!, knotinsertionoperator, degreeelevationoperator, bezierdecompoperator
export dim

# export functionality from nurbs.jl
export nsdim, nbezierpatches, ienarray
export refine!, globalrefine!, degreeelevate!
export boundary, ∂

# export functionality from this module
export dersbezierbasisfuns, assembly, plot, plot_stress_contours, L2residual

include("functions/baseextension.jl")
include("functions/knotvector.jl")
include("functions/bsplinebasisfuns.jl")
include("functions/knotrefinement.jl")
include("functions/nurbs.jl")
include("functions/quadrature.jl")

######################################################################################################
####################################### Shape function routines ######################################
######################################################################################################

# bernstein basis functions on [-1,1] and their derivatives
dersbezierbasisfuns(p::Int, u::Float64, nout::Int) = dersbsplinebasisfuns(p, buildvector([-1.0;1.0],[p+1;p+1]), p+1, u, nout)'

# bernstein basis functions on [-1,1] and their derivatives at a set of points u
function dersbezierbasisfuns(p::Int, u::Vector{Float64}, nout::Int)
    ϕ = zeros(p+1,length(u),nout)
    for k in 1:length(u)
        ϕ[:,k,:] = dersbezierbasisfuns(p, u[k], nout)
    end
    return ϕ
end

# Compute 1D shape function and derivatives in a single point
function dersbezierbasisfuns(S::BezierPatch{1}, ϕ::Vector{Float64}, dϕ::Vector{Float64})

    # initialize
    C = S.dec
    W = S.wpts
    P = S.cpts

    # NURBS are B-splines in projective three-space
    Bʷ  = broadcast(*, W, C[1]'* ϕ)
    ∇Bʷ = broadcast(*, W, C[1]'* dϕ)

    # compute weightfunction and derivative
    w  = sum(Bʷ)
    ∇w = sum(∇Bʷ)

    # compute NURBS basisfunction and derivative
    R  = Bʷ / w
    ∇R = (∇Bʷ * w - Bʷ * ∇w) / w^2

    # return quadtities in physical space
    return R, ∇R
end

# Compute 1D shape function and derivatives in an array of points
function dersbezierbasisfuns(S::BezierPatch{1}, ϕ::Matrix{Float64}, dϕ::Matrix{Float64})
    R, ∇R = zeros(ϕ), zeros(dϕ)
    for k in 1:size(R,2)
        R[:,k], ∇R[:,k] = dersbezierbasisfuns(S, ϕ[:,k], dϕ[:,k])
    end
    return R, ∇R
end

# Compute 2D shape function, derivatives and determinant in a single point
function dersbezierbasisfuns(S::BezierPatch{2}, ϕ::NTuple{2,Vector{Float64}}, dϕ::NTuple{2,Vector{Float64}})

    # initialize
    C = S.dec
    W = S.wpts
    P = S.cpts

    # compute B-spline basis functions
    ϕ  = ntuple(k-> C[k]'* ϕ[k], 2)
    dϕ = ntuple(k-> C[k]'*dϕ[k], 2)

    # NURBS are B-splines in projective 3-space
    Bʷ  = broadcast(*, W,  kron(ϕ[2], ϕ[1]))
    ∇Bʷ = broadcast(*, W, [kron(ϕ[2], dϕ[1]) kron(dϕ[2], ϕ[1])])

    # compute weightfunction
    w  = sum(Bʷ)
    ∇w = sum(∇Bʷ, 1)

    # compute NURBS functions and derivatives
    R  = Bʷ / w
    ∇R = (∇Bʷ * w - Bʷ * ∇w) / w^2

    # return quadtities in physical space
    return R, ∇R
end

# Compute 2D shape function and derivatives in an array of points
function dersbezierbasisfuns(S::BezierPatch{2}, ϕ::NTuple{2,Matrix{Float64}}, dϕ::NTuple{2,Matrix{Float64}})
    n = prod([size(ϕ[k],1) for k in 1:2]); m = ntuple(k -> size(ϕ[k],2), 2)
    R, ∇R = zeros(n,prod(m)), zeros(n,prod(m),2)
    i = 1
    for k in 1:m[2]
        for l in 1:m[1]
            R[:,i], ∇R[:,i,:] = dersbezierbasisfuns(S, (ϕ[1][:,k], ϕ[2][:,l]), (dϕ[1][:,k], dϕ[2][:,l]))
            i+=1
        end
    end
    return R, ∇R
end

######################################################################################################
########################################## Assembly routines #########################################
######################################################################################################

typealias Forcing  Function
typealias Traction Function
typealias BoundaryDisplacement Function
typealias Material Function

function assembly(S::NURBS, D::Material, f::Forcing, g::BoundaryDisplacement, h::Traction)

    #########################  Initialize #############################
    p     = S.deg                            # polynomial degree
    nsd   = nsdim(S)                         # number of space dimensions
    ndofs = dimsplinespace(S)                # dimension of the spline space
    dims  = nbezierpatches(S)                # number of elements in each direction

    IEN   = S.ien                            # set IEN-array
    ID    = S.id                             # set ID-array

    neq   = maximum(ID)                      # number of equations
    K     = spzeros(neq,neq)                 # allocate space for the stiffness matrix
    F     = zeros(neq)                       # allocate space for the righthandside forcing

    ########################  Precompute data #########################

    # precompute Gauss-Legendre nodes and weights
    qr = ntuple(k -> QuadRule(p[k]+1,"legendre"), nsd)

    # pre-compute univariate Bernstein polynomials and their derivatives
    ϕ = ntuple(k -> dersbezierbasisfuns(p[k], qr[k].points, 2), nsd)

    ########################  Assembly loop over elements #########################

    # integration over domain Ω
    for e in 1:size(IEN,2)

        # initialize
        A = IEN[:,e]                       # global node numbers
        P = ID[:,A]                          # global equation numbers
        i,j = ind2sub(dims,e)                # element subscripts

        # construct Bezier Patch
        Sb = BezierPatch(tuple(p...), S.cpts[A,:], S.wpts[A], (S.dec[1][:,:,i], S.dec[2][:,:,j]))

        # compute element contribution of the stiffness matrix and the forcing vector
        Kₑ, Fₑ = element_stiffness_and_forcing(Sb, ϕ, qr, D, f)

        # add contributions of element stiffness matrix and right hand side forcing
        I = P.!=0                  # index internal displacements
        K[P[I],P[I]] += Kₑ[I[:],I[:]]
        F[P[I]]      += Fₑ[I[:]]

        # subtract Dirichlet boundary displacements
        for i in 1:2
            for a in 1:size(P,2)
                if P[i,a]==0
                    F[P[I]] -= Kₑ[I,a] * g(vec(S.cpts[A[a],:]))[i]
                end
            end
        end
    end

    # integration over boundary ∂Ω
    for dir = 1:2
        I = 1:nsd.!=dir
        for k in 1:2
            C = ∂(S,dir,k) # boundary k, direction dir, of S
            for e in 1:size(C.ien,2)

                # initialize
                A = C.ien[:,e]                       # global node numbers
                P = ID[:,A][:]                       # global equation numbers
                J = P.!=0                            # dofs with Neumann data

                # construct Bezier Patch
                Cb = BezierPatch(tuple(C.deg...), S.cpts[A,:], S.wpts[A], (C.dec[1][:,:,e],))

                # compute element contribution of the stiffness matrix and the forcing vector
                F[P[J]] += neumann_boundary_condition(Cb, ϕ[I], qr[I], h)[J]

            end
        end
    end

    return K, F
end

# computation of element stiffness and forcing
function element_stiffness_and_forcing(Sb::BezierPatch{2}, ϕ::NTuple{2,Array{Float64,3}}, qr::NTuple{2,QuadRule}, D::Material, f::Forcing)

    # initialize
    nen = prod(Sb.deg+1)                 # number of element nodes
    ned = 2                              # number of element dofs per node
    Kₑ = zeros(nen*ned, nen*ned)         # allocate space for element stiffness matrix
    Fₑ = zeros(nen*ned)                  # allocate space for element forcing vector

    # loop over quadrature points
    for j in 1:dim(qr[2])
        for i in 1:dim(qr[1])

            # shape function routine
            R, ∇R = dersbezierbasisfuns(Sb,(ϕ[1][:,i,1],ϕ[2][:,j,1]), (ϕ[1][:,i,2],ϕ[2][:,j,2]))

            # compute geometrical properties
            J = ∇R' * Sb.cpts                                # Jacobian
            x = vec(R'  * Sb.cpts)                           # quadrature point in physical space
            α = qr[1].weights[i]*qr[2].weights[j]*det(J)     # quadrature weight in physical space

            # compute strain components
            N = (inv(J) * ∇R')'
            B = Matrix{Float64}[[N[a,1]  0.0;
                                 0.0    N[a,2];
                                 N[a,2] N[a,1]] for a in 1:nen]

            # compute element stiffness contribution of quadrature point i,j
            for a in 1:nen
                p = ned*(a-1)+(1:ned)
                for b in a:nen
                    q = ned*(b-1)+(1:ned)
                    Kₑ[p,q] += (B[a]' * D(x) * B[b]) * α
                end
            end

            # compute element forcing contribution of quadrature point i,j
            Fₑ += ((f(x) * α) * R')[:]
        end
    end

    # Use symmetry to determine the remainder of the element stiffness matrix
    for a in 1:nen
        p = ned*(a-1)+(1:ned)
        for b in a+1:nen
            q = ned*(b-1)+(1:ned)
            Kₑ[q,p] = Kₑ[p,q]'
        end
    end
    return Kₑ, Fₑ
end

# computation of element tractions
function neumann_boundary_condition(Cb::BezierPatch{1}, ϕ::Tuple{Array{Float64,3}}, qr::Tuple{QuadRule}, h::Function)

    # initialize
    nen = prod(Cb.deg+1)                 # number of element nodes
    ned = 2                              # number of element dofs per node
    Fₑ = zeros(nen*ned)                  # allocate space for element forcing vector

    # loop over quadrature points
    for i in 1:dim(qr[1])

        # shape function routine
        R, ∇R = dersbezierbasisfuns(Cb, ϕ[1][:,i,1], ϕ[1][:,i,2])

        # compute geometrical properties
        t = vec(∇R' * Cb.cpts)                           # Tangent vector
        x = vec(R'  * Cb.cpts)                           # quadrature point in physical space
        α = qr[1].weights[i]*norm(t)                     # quadrature weight in physical space

        # compute element forcing contribution of quadrature point i,j
        Fₑ += ((h(x) * α) * R')[:]
    end

    return Fₑ
end

######################################################################################################
########################################## Postprocessing routines ###################################
######################################################################################################

import PyPlot.plot

# plot the NURBS geometry
function plot(S::NURBS, gran=(3,3))
    # initialize
    if isempty(S.ien)==true
        S.ien = ienarray(S)
    end
    if isempty(S.dec[1])==true
        S.dec = bezierdecompoperator(S)
    end
    dims = ntuple(i -> size(S.dec[i],3), 2)
    u    = ntuple(i -> linearspace(-1.0,1.0,gran[i]), 2)
    ϕ    = ntuple(i -> dersbezierbasisfuns(S.deg[i], u[i], 2), 2)

    # loop over elements
    for e in 1:prod(dims)
        i,j = ind2sub(dims,e)          # element subscripts
        A   = S.ien[:,e]                 # global node numbers

        # construct Bezier Patch
        Sb  = BezierPatch(tuple(S.deg...), S.cpts[A,:], S.wpts[A], (S.dec[1][:,:,i], S.dec[2][:,:,j]))

        # compute 2-dimensional rational basis functions
        R, ∇R = dersbezierbasisfuns(Sb, (ϕ[1][:,:,1], ϕ[2][:,:,1]), (ϕ[1][:,:,2], ϕ[2][:,:,2]))

        # plot surface
        X = ntuple(i -> reshape(R' * Sb.cpts[:,i], gran), size(Sb.cpts,2))
        contourf(X[1], X[2], zeros(X[1]))

        # plot boundary
        plot(vec(X[1][:,1]),   vec(X[2][:,1]),"k")
        plot(vec(X[1][:,end]), vec(X[2][:,end]),"k")
        plot(vec(X[1][1,:]),   vec(X[2][1,:]),"k")
        plot(vec(X[1][end,:]), vec(X[2][end,:]),"k")
    end
end

# plot the contourlevels of stress
function plot_stress_contours(G::NURBS{2}, S::NURBS{2}, ii::Int, D::Material, levels::Vector{Float64}, gran=(4,4))

    # initialize
    p   = G.deg
    nen = prod(p+1)
    dims = ntuple(i -> size(S.dec[i],3), 2)

    # sample points
    u = ntuple(k -> linearspace(-1.0,1.0,gran[k]), 2)

    # pre-compute univariate Bernstein polynomials and their derivatives
    ϕ = ntuple(k -> dersbezierbasisfuns(p[k], u[k], 2), 2)

    # loop over elements
    for e in 1:size(S.ien,2)
        i,j = ind2sub(dims,e)          # element subscripts
        A   = G.ien[:,e]               # global node numbers

        # construct Bezier Patch
        Sb = BezierPatch(tuple(p...), G.cpts[A,:], G.wpts[A], (G.dec[1][:,:,i], G.dec[2][:,:,j]))

        # compute the stress at the sample points
        σd = zeros(gran[1],gran[2],3)
        X = zeros(gran); Y = zeros(gran)
        for j in 1:gran[2]
            for i in 1:gran[1]

                # shape function routine
                R, ∇R = dersbezierbasisfuns(Sb,(ϕ[1][:,i,1],ϕ[2][:,j,1]), (ϕ[1][:,i,2],ϕ[2][:,j,2]))

                # compute geometrical properties
                J = ∇R' * Sb.cpts                                # Jacobian
                x = vec(R'  * Sb.cpts)                           # quadrature point in physical space

                if det(J)==0.0
                    J = eye(2)
                end

                # compute strain basis functions
                N = (inv(J) * ∇R')'
                B = Matrix{Float64}[[N[a,1]  0.0;
                                     0.0    N[a,2];
                                     N[a,2] N[a,1]] for a in 1:nen]

                # compute the stress in point x
                for a in 1:nen
                    save = D(x)*B[a]*S.cpts[A[a],:]'
                    for k in 1:3
                        σd[i,j,k] += save[k]
                    end
                end

                # save the coordinates
                X[i,j] = x[1]; Y[i,j] = x[2]
            end
        end

        cp = contourf(X, Y, σd[:,:,ii],levels)
        if e==1
            colorbar(cp)
        end
    end
end

typealias Geometry NURBS{2}
typealias Solution NURBS{2}
typealias ExactSolution Function

# Compute the L2 residial in stress
function L2residual(G::Geometry, S::Solution, σ::ExactSolution, D::Material)

    # initialize
    p     = S.deg                            # polynomial degree
    nsd   = nsdim(S)                         # number of space dimensions
    nen   = prod(p+1)
    ndofs = dimsplinespace(S)                # dimension of the spline space
    dims  = nbezierpatches(S)                # number of elements in each direction
    IEN   = S.ien                            # set IEN-array

    # precompute Gauss-Legendre nodes and weights
    qr = ntuple(k -> QuadRule(p[k]+1,"legendre"), nsd)

    # pre-compute univariate Bernstein polynomials and their derivatives
    ϕ = ntuple(k -> dersbezierbasisfuns(p[k], qr[k].points, 2), nsd)

    # integration over domain Ω
    L²res = zeros(3)
    L21 = zeros(3)
    L22 = zeros(3)
    for e in 1:size(IEN,2)

        # initialize
        A = IEN[:,e]                         # global node numbers
        i,j = ind2sub(dims,e)                # element subscripts

        # construct Bezier Patch
        Sb = BezierPatch(tuple(p...), G.cpts[A,:], G.wpts[A], (G.dec[1][:,:,i], G.dec[2][:,:,j]))

        # loop over quadrature points
        for j in 1:dim(qr[2])
            for i in 1:dim(qr[1])

                # shape function routine
                R, ∇R = dersbezierbasisfuns(Sb,(ϕ[1][:,i,1],ϕ[2][:,j,1]), (ϕ[1][:,i,2],ϕ[2][:,j,2]))

                # compute geometrical properties
                J = ∇R' * Sb.cpts                                # Jacobian
                x = vec(R'  * Sb.cpts)                           # quadrature point in physical space
                α = qr[1].weights[i]*qr[2].weights[j]*det(J)     # quadrature weight in physical space

                # compute strain components
                N = (inv(J) * ∇R')'

                # compute the stress in point x
                σd = zeros(3)
                for a in 1:nen
                    # strain basis function `a'
                    B = [N[a,1]    0.0;
                            0.0 N[a,2];
                         N[a,2] N[a,1]]

                    # compute the stress in quadrature point i,j
                    σd += D(x) * B * S.cpts[A[a],:]'
                end
                # compute element forcing contribution of quadrature point i,j
                L²res += (σd - σ(x)).^2 * α
            end
        end

    end
    return sqrt(L²res)
end

end

ext: jl

max_stars_repo_name: hiemstar/iga-demos

max_issues_repo_name: hiemstar/iga-demos

max_forks_repo_name: hiemstar/iga-demos


text:
let
    @test isapprox(deltaangle(0.0, 0.0),  0.0)
    @test isapprox(deltaangle(0.0, 1.0),  1.0)
    @test isapprox(deltaangle(0.0, 2π),  0.0, atol=1e-10)
    @test isapprox(deltaangle(0.0, π + 1.0), -(π- 1.0))

    @test isapprox(angledist(0.0, 0.0), 0.0)
    @test isapprox(angledist(0.0, 1.0), 1.0)
    @test isapprox(angledist(0.0, 2π), 0.0, atol=1e-10)
    @test isapprox(angledist(0.0, π + 1.0), π- 1.0)

    @test isapprox(lerp_angle(0.0, 0.0, 1.0), 0.0)
    @test isapprox(lerp_angle(0.0, 2π, 0.5), 0.0, atol=1e-10)
    @test isapprox(lerp_angle(0.0, 2.0, 0.75),  1.5)
    @test isapprox(lerp_angle(0.0, 2π-2, 0.75), -1.5)

    @test are_collinear(VecE2(0.0,0.0), VecE2(1.0,1.0), VecE2(2.0,2.0))
    @test are_collinear(VecE2(0.0,0.0), VecE2(1.0,1.0), VecE2(-2.0,-2.0))
    @test are_collinear(VecE2(0.0,0.0), VecE2(0.5,1.0), VecE2(1.0,2.0))
    @test are_collinear(VecE2(1.0,2.0), VecE2(0.0,0.0), VecE2(0.5,1.0))
    @test !are_collinear(VecE2(0.0,0.0), VecE2(1.0,0.0), VecE2(0.0,1.0))
end

let
    seg = LineSegment1D(0,1)
    seg2 = seg + 0.5
    @test isapprox(seg2.a, 0.5)
    @test isapprox(seg2.b, 1.5)

    seg2 = seg - 0.3
    @test isapprox(seg2.a, -0.3)
    @test isapprox(seg2.b,  0.7)

    seg = LineSegment1D(2,5)
    @test isapprox(get_distance(seg, 1.5), 0.5)
    @test isapprox(get_distance(seg, -1.5), 3.5)
    @test isapprox(get_distance(seg,  2.5), 0.0)
    @test isapprox(get_distance(seg,  2.0), 0.0)
    @test isapprox(get_distance(seg,  5.0), 0.0)
    @test isapprox(get_distance(seg,  5.5), 0.5)
    @test isapprox(get_distance(seg, 15.5), 10.5)

    @test !in(-1.0, seg)
    @test !in(1.0, seg)
    @test in(2.0, seg)
    @test in(3.0, seg)
    @test in(5.0, seg)
    @test !in(5.5, seg)

    @test in(LineSegment1D(3.0, 3.5), seg)
    @test in(LineSegment1D(3.0, 3.0), seg)
    @test !in(LineSegment1D(3.0, 13.0), seg)
    @test !in(LineSegment1D(-3.0, 3.0), seg)
    @test !in(LineSegment1D(-3.0, -2.0), seg)
    @test !in(LineSegment1D(13.0, 15.0), seg)

    @test intersects(seg, LineSegment1D(3.0, 3.5))
    @test intersects(seg, LineSegment1D(3.0, 3.0))
    @test intersects(seg, LineSegment1D(3.0, 13.0))
    @test intersects(seg, LineSegment1D(-3.0, 3.0))
    @test !intersects(seg, LineSegment1D(-3.0, -2.0))
    @test !intersects(seg, LineSegment1D(13.0, 15.0))
end

let
    L = Line(VecE2(0,0), π/4)
    L2 = L + VecE2(-1,1)
    @test isapprox(L2.C, VecE2(-1,1))
    @test isapprox(L2.θ, π/4)
    L2 = L - VecE2(-1,1)
    @test isapprox(L2.C, VecE2(1,-1))
    @test isapprox(L2.θ, π/4)
    L2 = rot180(L)
    @test isapprox(L2.θ, π/4 + π)
    @test isapprox(get_polar_angle(L), atan(1,1))
    @test isapprox(get_distance(L, VecE2(0,0)), 0)
    @test isapprox(get_distance(L, VecE2(1,1)), 0, atol=1e-10)
    @test isapprox(get_distance(L, VecE2(1,0)), √2/2)
    @test isapprox(get_distance(L, VecE2(0,1)), √2/2)
    @test get_side(L, VecE2(0,0)) ==  0
    @test get_side(L, VecE2(1,1)) ==  0
    @test get_side(L, VecE2(0,1)) ==  1
    @test get_side(L, VecE2(1,0)) == -1
end

let
    L = LineSegment(VecE2(0,0), VecE2(1,1))
    L2 = L + VecE2(-1,1)
    @test isapprox(L2.A, VecE2(-1,1))
    @test isapprox(L2.B, VecE2( 0,2))
    L2 = L - VecE2(-1,1)
    @test isapprox(L2.A, VecE2(1,-1))
    @test isapprox(L2.B, VecE2(2, 0))
    @test isapprox(get_polar_angle(L), atan(1,1))
    @test isapprox(get_distance(L, VecE2(0,0)), 0)
    @test isapprox(get_distance(L, VecE2(1,1)), 0)
    @test isapprox(get_distance(L, VecE2(1,0)), √2/2)
    @test isapprox(get_distance(L, VecE2(0,1)), √2/2)
    @test get_side(L, VecE2(1,1)) ==  0
    @test get_side(L, VecE2(0,1)) ==  1
    @test get_side(L, VecE2(1,0)) == -1

    @test angledist(L, L2) ≈ 0.0
    @test parallel(L, L2)

    L2 = LineSegment(VecE2(0,0), VecE2(-1,1))
    @test angledist(L, L2) ≈ π/2
    @test !parallel(L, L2)
    @test intersects(L, L2)

    L2 = LineSegment(VecE2(0,0), VecE2(1,-1))
    @test angledist(L, L2) ≈ π/2
    @test !parallel(L, L2)
    @test intersects(L, L2)

    L2 = LineSegment(VecE2(0,0), VecE2(1,-1)) + VecE2(5,-7)
    @test angledist(L, L2) ≈ π/2
    @test !parallel(L, L2)
    @test !intersects(L, L2)

    @test  intersects(L, LineSegment(VecE2(1,1), VecE2(2,-1)))
    @test !intersects(L, LineSegment(VecE2(2,2), VecE2(5,6)))
    @test !intersects(L, LineSegment(VecE2(-5,-5), VecE2(5,-5)))
end

let
    @test isapprox(intersect(Ray(-1,0,0), Ray(0,-1,π/2)), VecE2(0.0,0.0))

    @test  intersects(Ray(0,0,0), Ray(1,-1,π/2))
    @test !intersects(Ray(0,0,0), Ray(1,-1,-π/2))
    # @test  intersects(Ray(0,0,0), Ray(1,0,0)) # TODO: get this to work

    @test  intersects(Ray(0,0,0),   Line(VecE2(0,0), VecE2(1,1)))
    @test !intersects(Ray(0,0,0),   Line(VecE2(0,1), VecE2(1,1)))
    @test  intersects(Ray(0,0,0),   Line(VecE2(1,0), VecE2(2,0)))
    @test  intersects(Ray(0,0,0),   Line(VecE2(1,-1), VecE2(1,1)))
    @test  intersects(Ray(0,0,π/2), Line(VecE2(-1,1), VecE2(1,1)))
    @test !intersects(Ray(0,0,π/2), Line(VecE2(-1,1), VecE2(-1,2)))
    @test  intersects(Ray(0,0,π/2), Line(VecE2(-1,1), VecE2(-1.5,1.5)))

    @test  intersects(Ray(0,0,0), LineSegment(VecE2(0,0), VecE2(1,1)))
    @test !intersects(Ray(0,0,0), LineSegment(VecE2(0,1), VecE2(1,1)))
    @test  intersects(Ray(0,0,0), LineSegment(VecE2(1,0), VecE2(2,0)))
    @test !intersects(Ray(0,0,0), LineSegment(VecE2(-1,0), VecE2(-2,0)))

    @test isapprox(intersect(VecSE2(4,0,3pi/4), LineSegment(VecE2(5.6,0), VecE2(0,3.6))),
                   VecE2(1.12, 2.88), atol=1e-3)
end

let
    proj1 = Projectile(VecSE2(0.0,0.0,0.0), 1.0)
    proj2 = Projectile(VecSE2(1.0,1.0,1.0), 1.0)
    @test isapprox(VecE2(propagate(proj1, 2.0).pos), VecE2(2.0,0.0))
    @test isapprox(VecE2(propagate(proj2, 0.0).pos), VecE2(1.0,1.0))

    t_PCA, d_PCA = closest_time_of_approach_and_distance(proj1, Projectile(VecSE2(0.0,1.0,0.0), 1.0))
    @test isapprox(t_PCA, 0.0)
    @test isapprox(d_PCA, 1.0)

    t_PCA, d_PCA = closest_time_of_approach_and_distance(proj1, Projectile(VecSE2(0.0,2.0,0.0), 1.0))
    @test isapprox(t_PCA, 0.0)
    @test isapprox(d_PCA, 2.0)

    t_PCA, d_PCA = closest_time_of_approach_and_distance(proj1, Projectile(VecSE2(0.0,1.0,0.0), 2.0))
    @test isapprox(t_PCA, 0.0)
    @test isapprox(d_PCA, 1.0)

    t_PCA, d_PCA = closest_time_of_approach_and_distance(Projectile(VecSE2(0.0,-0.0,0.2), 1.0), Projectile(VecSE2(0.0,1.0,-0.2), 1.0))
    @test t_PCA > 0.0
    @test isapprox(d_PCA, 0.0)

    @test isapprox(get_intersection_time(proj1, LineSegment(VecE2(0, 0), VecE2(0,1))), 0.0)
    @test isapprox(get_intersection_time(proj1, LineSegment(VecE2(1,-1), VecE2(1,1))), 1.0)
    @test    isinf(get_intersection_time(proj1, LineSegment(VecE2(-1,0), VecE2(-1,1))))
    @test isapprox(get_intersection_time(proj1, LineSegment(VecE2(3,0), VecE2(2,0))), 2.0)

    t_PCA, d_PCA = closest_time_of_approach_and_distance(proj1, LineSegment(VecE2(0, 0), VecE2(0,1)))
    @test isapprox(t_PCA, 0.0)
    @test isapprox(d_PCA, 0.0)
    t_PCA, d_PCA = closest_time_of_approach_and_distance(proj1, LineSegment(VecE2(1,-1), VecE2(1,1)))
    @test isapprox(t_PCA, 1.0)
    @test isapprox(d_PCA, 0.0)
    t_PCA, d_PCA = closest_time_of_approach_and_distance(proj1, LineSegment(VecE2(1,0.5), VecE2(1,1)))
    @test isapprox(t_PCA, 1.0)
    @test isapprox(d_PCA, 0.5)
    t_PCA, d_PCA = closest_time_of_approach_and_distance(proj1, LineSegment(VecE2(1,-0.5), VecE2(1,-1)))
    @test isapprox(t_PCA, 1.0)
    @test isapprox(d_PCA, 0.5)
    t_PCA, d_PCA = closest_time_of_approach_and_distance(proj1, LineSegment(VecE2(1,-0.5), VecE2(2,-1)))
    @test isapprox(t_PCA, 1.0)
    @test isapprox(d_PCA, 0.5)
    t_PCA, d_PCA = closest_time_of_approach_and_distance(proj1, LineSegment(VecE2(2,-1), VecE2(1,-0.5)))
    @test isapprox(t_PCA, 1.0)
    @test isapprox(d_PCA, 0.5)
    t_PCA, d_PCA = closest_time_of_approach_and_distance(proj1, LineSegment(VecE2(2,-1), VecE2(1,-0.5)), true)
    @test isinf(t_PCA)
    @test isinf(d_PCA)
    # t_PCA, d_PCA = closest_time_of_approach_and_distance(proj1, LineSegment(VecE2(-1,1), VecE2(1,1)))
    # println((t_PCA, d_PCA ))
    # @test isapprox(t_PCA, 0.0)
    # @test isapprox(d_PCA, 1.0)
end

let
    @test  in(VecE2(0,0), Circ(0,0,1.0))
    @test !in(VecE2(0,0), Circ(2,0,1.0))
    @test  in(VecE2(1.5, 0), Circ(2,0,1.0))
    @test  in(VecE3(1.5,0,3), Circ(2,0,3,1.0))

    box = AABB(VecE2(0.0, 0.5), 2.0, 5.0)
    @test  in(VecE2( 0.0,0.0), box)
    @test  in(VecE2(-1.0,0.0), box)
    @test !in(VecE2(-2.0,0.0), box)
    @test !in(VecE2( 1.0,3.1), box)

    box = AABB(VecE2(-1.0, -2.0), VecE2(1.0, 3.0))
    @test  in(VecE2( 0.0,0.0), box)
    @test  in(VecE2(-1.0,0.0), box)
    @test !in(VecE2(-2.0,0.0), box)
    @test !in(VecE2( 1.0,3.1), box)
end

let
    box = OBB(VecSE2(0.0, 0.5, 0.0), 2.0, 5.0)
    @test  in(VecE2( 0.0,0.0), box)
    @test  in(VecE2(-1.0,0.0), box)
    @test !in(VecE2(-2.0,0.0), box)
    @test !in(VecE2( 1.0,3.1), box)

    box = OBB(VecSE2(0.0, 0.0, π/4), 2.0, 2.0)
    @test  in(VecE2( 0.0,0.0), box)
    @test  in(VecE2( 1.0,0.0), box)
    @test  in(VecE2( √2/2-0.1,√2/2-0.1), box)
    @test !in(VecE2( 1.0,1.0), box)
    @test !in(VecE2( √2/2+0.1,√2/2+0.1), box)

    box = OBB(VecSE2(0.0, 0.0, 0.0), 1.0, 1.0)
    box2 = inertial2body(box, VecSE2(1.0,0.0,0.0))
    @test isapprox(box2.aabb.center.x, -1.0)
    @test isapprox(box2.aabb.center.y,  0.0)
    @test isapprox(box2.θ,  0.0)

    box2 = inertial2body(box, VecSE2(1.0,0.0,π))
    @test isapprox(box2.aabb.center.x,  1.0)
    @test isapprox(box2.aabb.center.y,  0.0, atol=1e-10)
    @test isapprox(angledist(box2.θ,  π), 0.0, atol=1e-10)
end

let
    plane = Plane3()
    @test plane.normal == VecE3(1.0,0.0,0.0)
    @test plane.offset == 0.0

    plane = Plane3(VecE3(2.0,0.0,0.0), 0.5)
    @test plane.normal == VecE3(1.0,0.0,0.0)
    @test plane.offset == 0.5

    plane = Plane3(VecE3(2.0,0.0,0.0), VecE3(0.0,0.0,0.0))
    @test plane.normal == VecE3(1.0,0.0,0.0)
    @test plane.offset == 0.0

    plane = Plane3(VecE3(2.0,0.0,0.0), VecE3(1.0,1.0,1.0))
    @test plane.normal == VecE3(1.0,0.0,0.0)
    @test isapprox(plane.offset, -1.0, atol=1e-10)

    plane = Plane3(VecE3(1.0,1.0,0.0), VecE3(1.0,0.0,0.0), VecE3(0.0,0.0,0.0))
    @test norm(plane.normal - VecE3(0.0,0.0,1.0)) < 1e-8
    @test isapprox(plane.offset, 0.0, atol=1e-10)

    plane = Plane3(VecE3(0.0,0.0,0.0), VecE3(1.0,0.0,0.0), VecE3(1.0,1.0,0.0))
    @test norm(plane.normal - VecE3(0.0,0.0,-1.0)) < 1e-8
    @test isapprox(plane.offset, 0.0, atol=1e-10)

    plane = Plane3()
    @test isapprox(get_signed_distance(plane, VecE3( 0.0,0.0,0.0)),  0.0, atol=1e-10)
    @test isapprox(get_signed_distance(plane, VecE3( 1.0,0.0,0.0)),  1.0, atol=1e-10)
    @test isapprox(get_signed_distance(plane, VecE3(-1.0,0.0,0.0)), -1.0, atol=1e-10)
    @test isapprox(get_signed_distance(plane, VecE3( 1.0,0.5,0.7)),  1.0, atol=1e-10)
    @test isapprox(get_signed_distance(plane, VecE3(-1.0,0.7,0.5)), -1.0, atol=1e-10)
    @test isapprox(get_distance(plane, VecE3(-1.0,0.7,0.5)), 1.0, atol=1e-10)
    @test norm(proj(VecE3(-1.0,0.7,0.5), plane) - VecE3(0.0,0.7,0.5)) < 1e-10
    @test get_side(plane, VecE3( 0.0,0.0,0.0)) ==  0
    @test get_side(plane, VecE3( 1.0,0.0,0.0)) ==  1
    @test get_side(plane, VecE3(-1.0,0.0,0.0)) == -1
    @test get_side(plane, VecE3( 1.0,0.5,0.7)) ==  1
    @test get_side(plane, VecE3(-1.0,0.7,0.5)) == -1
end


text:
#!/bin/julia

###################################
# Garamond search server CLI tool #
###################################

module gars  # The Garamond server

using Pkg
project_root_path = "/" * joinpath(split(@__FILE__, "/")[1:end-4]...)
Pkg.activate(project_root_path)
using Garamond
using Sockets
using Logging
using ArgParse

# Function that parses Garamond's server arguments
function get_server_commandline_arguments(args::Vector{String})
	s = ArgParseSettings()
	@add_arg_table! s begin
        "--data-config", "-d"
            help = "data configuration file"
            arg_type = String
        "--env-cache", "-c"
            help = "search environment cache file"
            arg_type = String
        "--log-level"
            help = "logging level"
            default = "info"
        "--log", "-l"
            help = "logging stream"
            default = "stdout"
        "--unix-socket", "-u"
            help = "UNIX socket for data communication"
            arg_type = String
        "--web-socket-port", "-w"
            help = "WEB socket data communication port"
            arg_type = UInt16
        "--web-socket-ip"
            help = "WEB socket data communication IP"
            default = "127.0.0.1"
        "--http-port", "-p"
            help = "HTTP port for REST services"
            arg_type = Int
        "--http-ip"
            help = "HTTP IP for REST services"
            default = "0.0.0.0"
        "--search-server-port", "-i"
            help = "Internal TCP port for the search server"
            arg_type = Int
            default = 9_000
	end
	return parse_args(args,s)
end

# Tests whether an ip is valid
function isvalidip(ip::AbstractString)
    try
        IPv4(ip), IPv6(ip)
        return true
    catch
        return false
    end
end

########################
# Main module function #
########################
function julia_main()::Cint
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function real_main()
    @info "~ GARAMOND ~ $(Garamond.printable_version())\n"

    # Parse command line arguments
    args = get_server_commandline_arguments(ARGS)

    # Get the argument values
    log_level = args["log-level"]
    logging_stream = args["log"]

    # Logging
    logger = Garamond.build_logger(logging_stream, log_level)
    global_logger(logger)

    # Get IP's
    ws_ip = args["web-socket-ip"]
    http_ip = args["http-ip"]
    if !isvalidip(ws_ip)
        @warn "Web-socket IP $ws_ip is not valid. Exiting..."
        return 0
    elseif !isvalidip(http_ip)
        @warn "HTTP IP $http_ip is not valid. Exiting..."
        return 0
    end

    # Get sockets/ports and check
    unixsocket = args["unix-socket"]
    ws_port = args["web-socket-port"]
    http_port = args["http-port"]
    io_port = args["search-server-port"]

    if unixsocket == nothing && ws_port == nothing && http_port == nothing
        @warn """At least a UNIX-socket, WEB-socket port or HTTP port
                 have to be specified. Use the -u, -w or -p options.
                 Exiting..."""
        return 0
    end

    if io_port < 0
        @warn "Search server I/O port value error. Defaulting to 9000."
        io_port = 9000
    end

    # Check data path and cache
    data_config_path = args["data-config"]
    env_cache_path = args["env-cache"]
    if data_config_path === nothing && env_cache_path === nothing
        @warn "Specify a data configuration or cache file using the -d or -c options. Exiting..."
        return 0
    end

    # Start I/O server(s)
    sserver_ready = Condition()  # search server ready
    unixsocket != nothing && @async unix_socket_server(unixsocket, io_port, sserver_ready)
    ws_port != nothing && @async web_socket_server(ws_port, io_port, sserver_ready, ipaddr=ws_ip)
    http_port != nothing && @async rest_server(http_port, io_port, sserver_ready, ipaddr=http_ip)

    # Start Search server
    Garamond.search_server(data_config_path, io_port, sserver_ready; cache_path=env_cache_path)
    return 0
end

################################
# Start main Garamond function #
################################
if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end

end  # module


text:
# common Base functionality

function Base.repeat(a::AbstractGPUVecOrMat, m::Int, n::Int = 1)
    o, p = size(a, 1), size(a, 2)
    b = similar(a, o*m, p*n)
    if length(b) == 0
        return b
    end
    gpu_call(b, a, o, p, m, n; total_threads=n) do ctx, b, a, o, p, m, n
        j = linear_index(ctx)
        j > n && return
        d = (j - 1) * p + 1
        @inbounds for i in 1:m
            c = (i - 1) * o + 1
            for r in 1:p
                for k in 1:o
                    b[k - 1 + c, r - 1 + d] = a[k, r]
                end
            end
        end
        return
    end
    return b
end

function Base.repeat(a::AbstractGPUVector, m::Int)
    o = length(a)
    b = similar(a, o*m)
    if length(b) == 0
        return b
    end
    gpu_call(b, a, o, m; total_threads=m) do ctx, b, a, o, m
        i = linear_index(ctx)
        i > m && return
        c = (i - 1)*o + 1
        @inbounds for i in 1:o
            b[c + i - 1] = a[i]
        end
        return
    end
    return b
end

## PermutedDimsArrays

using Base: PermutedDimsArrays

# PermutedDimsArrays' custom copyto! doesn't know how to deal with GPU arrays
function PermutedDimsArrays._copy!(dest::PermutedDimsArray{T,N,<:Any,<:Any,<:AbstractGPUArray}, src) where {T,N}
    dest .= src
    dest
end


text:
using Test
using WinEncoding

let
    @testset "codepage 1252" begin
        @test decode1252([0xf0])=="ð"
        @test decode1252([0x33,0x34,0x35])=="345"
        @test decode1252([0xa2,0xe7])=="¢ç"
    end
    @testset "codepage 950" begin
        @test WinEncoding.decode950([0xa4,0x48])=="人"
        @test WinEncoding.decode950([0xa4])=="�" #"\ufffd"
    end
    @testset "codepage 936" begin
        @test WinEncoding.decode936([0x81,0x40])=="丂"
        @test WinEncoding.decode936([0x80])=="€"
        @test WinEncoding.decode936([0xff])=="\uf8f5"
    end
    @testset "codepage 932" begin
        @test WinEncoding.decode932([0x80])=="\u80"
        @test WinEncoding.decode932([0xff])=="\uf8f3"
        @test WinEncoding.decode932([0xa0])=="\uf8f0"
        @test WinEncoding.decode932([0x82,0xa0])=="あ"
        @test WinEncoding.decode932([0x89,0xbe])=="伽"
    end
end


text:
include("./artificial_data.jl")

using Distributions
using CellwiseEditingDifferentiation
using Lora

### Targets

hyperpars = Dict{Symbol, Any}(:λ=>fill(10.0, nsites))

prior = Dict{Symbol, Any}(:v=>Function[(m::Float64, v::Float64, a::Float64, b::Float64)->
                          vpcprior(m, v, a, b, hyperpars[:λ][i]) for i in 1:nsites])

target = Dict{Symbol, Any}(:w=>Array(Function, nsites), :p=>Dict{Int, Vector{Function}}())

for i in 1:nsites
  target[:w][i] = function (w::Vector{Float64}, p::Vector{Float64})
    v = data[:m][i]*(1-data[:m][i])*inv_logit(w[1])
    a, b = beta_pars_from_mv(data[:m][i], v)
    sum([logpdf(Beta(a, b), p[m]) for m in 1:ncells[i]])+prior[:v][i](data[:m][i], v, a, b)+logdvdw(w[1], data[:m][i])
  end

  target[:p][i] = Array(Function, ncells[i])
  for j in 1:ncells[i]
    target[:p][i][j] = function (w::Float64)
      a, b = beta_pars_from_mv(data[:m][i], data[:m][i]*(1-data[:m][i])*inv_logit(w))
      Beta(a+data[:edited][i, cells[i][j]], b+data[:coverage][i, cells[i][j]]-data[:edited][i, cells[i][j]])
    end
  end
end

### Initial conditions

init = Dict{Symbol, Any}(:w=>Float64[0. for i in 1:nsites])
init[:p] = Dict{Int, Vector{Float64}}()
for i in 1:nsites
  init[:p][i] = Array(Float64, ncells[i])
  for j in 1:ncells[i]
    init[:p][i][j] = rand(target[:p][i][j](init[:w][i]))
  end
end

### Modellers

modeller = Dict{Symbol, Any}(:w=>Array(Function, nsites))
for i in 1:nsites
  modeller[:w][i] = function (p::Vector{Float64}, init::Vector{Float64})
    model(w->target[:w][i](w, p), init=init)
  end
end

### Samplers

sampler = Dict{Symbol, Any}(:w=>[RAM() for i in 1:nsites])

### Runners

runner = Dict{Symbol, Any}(:w=>[SerialMC(burnin=0, nsteps=1, thinning=1) for i in 1:nsites])

### Jobs

job = Dict{Symbol, Any}(:w=>Array(Function, nsites), :p=>Dict{Int, Vector{Function}}())

for i in 1:nsites
  job[:w][i] = function (p::Vector{Float64}, init::Vector{Float64})
    MCJob(modeller[:w][i](p, init), sampler[:w][i], runner[:w][i])
  end

  job[:p][i] = Array(Function, ncells[i])
  for j in 1:ncells[i]
    job[:p][i][j] = function (w::Float64)
      DistJob(target[:p][i][j](w))
    end
  end
end

job[:gibbs] = SerialMC(burnin=1000, nsteps=11000, thinning=1)

### Run Gibbs sampler

mcchain = gibbs(data, init, job)

### Write output

for i in 1:nsites
  VOUTFILE = joinpath(OUTDIR, @sprintf("vchain_%s_site%02d.txt", string(simulationid), i))
  writedlm(VOUTFILE,  map(w->data[:m][i]*(1-data[:m][i])*inv_logit(w), mcchain[:w][:, i]), ' ')

  POUTFILE = joinpath(OUTDIR, @sprintf("pchain_%s_site%02d.txt", string(simulationid), i))
  writedlm(POUTFILE, mcchain[:p][i], ' ')
end


text:
# Julia wrapper for header: /usr/local/cuda-7.5/include/nvrtc.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function nvrtcGetErrorString(result)
    ccall((:nvrtcGetErrorString,libnvrtc),Ptr{UInt8},(nvrtcResult,),result)
end

function nvrtcVersion(major,minor)
    check_nvrtcresult(ccall((:nvrtcVersion,libnvrtc),nvrtcResult,(Ptr{Cint},Ptr{Cint}),major,minor))
end

function nvrtcCreateProgram(prog,src,name,numHeaders,headers,includeNames)
    check_nvrtcresult(ccall((:nvrtcCreateProgram,libnvrtc),nvrtcResult,(Ptr{nvrtcProgram},Ptr{UInt8},Ptr{UInt8},Cint,Ptr{Ptr{UInt8}},Ptr{Ptr{UInt8}}),prog,src,name,numHeaders,headers,includeNames))
end

function nvrtcDestroyProgram(prog)
    check_nvrtcresult(ccall((:nvrtcDestroyProgram,libnvrtc),nvrtcResult,(Ptr{nvrtcProgram},),prog))
end

function nvrtcCompileProgram(prog,numOptions,options)
    check_nvrtcresult(ccall((:nvrtcCompileProgram,libnvrtc),nvrtcResult,(nvrtcProgram,Cint,Ptr{Ptr{UInt8}}),prog,numOptions,options))
end

function nvrtcGetPTXSize(prog,ptxSizeRet)
    check_nvrtcresult(ccall((:nvrtcGetPTXSize,libnvrtc),nvrtcResult,(nvrtcProgram,Ptr{Csize_t}),prog,ptxSizeRet))
end

function nvrtcGetPTX(prog,ptx)
    check_nvrtcresult(ccall((:nvrtcGetPTX,libnvrtc),nvrtcResult,(nvrtcProgram,Ptr{UInt8}),prog,ptx))
end

function nvrtcGetProgramLogSize(prog,logSizeRet)
    check_nvrtcresult(ccall((:nvrtcGetProgramLogSize,libnvrtc),nvrtcResult,(nvrtcProgram,Ptr{Csize_t}),prog,logSizeRet))
end

function nvrtcGetProgramLog(prog,log)
    check_nvrtcresult(ccall((:nvrtcGetProgramLog,libnvrtc),nvrtcResult,(nvrtcProgram,Ptr{UInt8}),prog,log))
end


text:
function synthesise_coupler(n::Int, _K, _BW_REL, POLYNOMIAL = :Butterworth)
    # POLYNOMIAL = :Butterworth
    # POLYNOMIAL = :Chebyshev

    # Parse inputs to desired precision (input string in ord)
    K = isa(_K, String) ? parse(BigFloat, _K) : convert(BigFloat, _K)
    BW_REL = isa(_BW_REL, String) ? parse(BigFloat, _BW_REL) : convert(BigFloat, _BW_REL)

    # Polynomial for pole extraction
    P_x̄ = -design_polynomial(POLYNOMIAL, n-1, BigFloat, :x̄)

    # Calculate related constants
    tc = tand(45 * (1 - BW_REL/2))      # Always a real number
	xc = real(x(im*tc))                 # Always a real number
    K_ = parse(BigFloat,"0.5") * K * P_x̄.coeffs[end];

	# Design function for calculations
	Γ_T(t) = parse(BigFloat,"0.5") * K * t * P_x̄(x(t)/xc)

    ## REFLECTION COEFFICIENT
    # Polynomial roots
    r̄ = roots(P_x̄)

    # Numerator zeros
    Z_n = begin
        # Preallocate zeros (for both functions and add )
        local Z_n = zeros(Complex{BigFloat}, 2n-1)
        for i in eachindex(r̄)
            Z_n[2i:(2i+1)] = [im, -im] .* sqrt((1-xc*r̄[i])/(1+xc*r̄[i]))
        end
        Z_n
    end

    # Denominator zeros
    LHS = Polynomial(xc .* BigFloat[1,0,-1], :t)^(n-1)
    RHS = begin
        local RHS_ = K_ * Polynomial(BigFloat[0,1], :t);
        for i in eachindex(r̄)
            RHS_ = RHS_ *
                (Polynomial(BigFloat[1,0,1], :t) - xc*r̄[i]*Polynomial(BigFloat[1,0,-1], :t))
        end
        RHS_
    end
    Z_d = begin
        Z_d_a = roots(LHS - RHS);
        Z_d_b = roots(LHS + RHS);
        # Join roots from both parentheses
        Z_d_ = vcat(Z_d_a, Z_d_b);
        # Select only left half-plane
        Z_d = Z_d_[sortperm(real.(Z_d_))[1:2n-1]]
    end

    # Coefficient reconstruction
    Γ_N = -1 * Polynomial(real.(fromroots(Complex{BigFloat}.(Z_n)).coeffs), :t)
    Γ_D = Polynomial(real.(fromroots(Complex{BigFloat}.(Z_d)).coeffs), :t)
    Γ(t) = (Γ_N(t)/Γ_D(t))

    ## TRANSMISSION COEFFICIENT
    _Kt = sqrt(abs(((LHS+RHS)*(LHS-RHS)).coeffs[end]))
	T(t) = -LHS(t)/(_Kt * Γ_D(t))

    ## DARLINGTON SYNTHESIS
    # Preallocate vectors
    a = Vector{BigFloat}(undef, n)
    b = Vector{BigFloat}(undef, n-1)

    # Initial admittance functions
    Y_N = Γ_D - Γ_N;    # Numerator of admittance
    Y_D = Γ_D + Γ_N;    # Denominator of admittance
    dY_in = polynomial_fraction_derivative(Y_N, Y_D) # Derivative of admittance

    # Extract first branch as an open-circuit stub
    a[1] = dY_in(1)
    a[end] = a[1]

    # Update input admittance and derivative
    Y_N = Y_N - a[1] * Polynomial(BigFloat[0,1], :t) * Y_D
    Y_D = Y_D
    dY_in = polynomial_fraction_derivative(Y_N, Y_D)

    # Iterate through all "transmission line followed by shunt open-circuit stub" sections
    for i in 2:Int(floor(n/2)+1)
        # Extract transmission line
        b[i-1] = Y_N(1)/Y_D(1)
        b[end-i+2] = b[i-1]

        # Get new admittance
        # Calculate numerator and denominator naïvely
        Y_N_ = Y_N * Polynomial(BigFloat[1,0,1], :t) * b[i-1] -
            2 * Y_D * Polynomial(BigFloat[0,1], :t) * (b[i-1]^2)
        Y_D_ = Y_D * Polynomial(BigFloat[1,0,1], :t) * b[i-1] -
            2 * Y_N * Polynomial(BigFloat[0,1], :t)
        # Remove common roots
        rn = roots(Y_N_) # Roots of naïve numerator
        rd = roots(Y_D_) # Roots of naïve denominator
        _goodi = fill(true,length(rn)); # Flags for non-common root
        for j in eachindex(rn)
            # Find first root in common within the range of precision
            #   of reconstructed function
            k = findfirst(isapprox.(rn[j], rd, atol = eps(Float64(abs(rn[j])))))
            # If there exists, pop from denominator and flag as common root for numerator
            isnothing(k) || (_goodi[j] = false; popat!(rd,k))
        end
        rn = rn[_goodi]; # Remove numerator common roots
        # Get coefficients with potentially spurious imaginary parts
        c_Y_N = (Y_N_.coeffs[end] * fromroots(Complex{BigFloat}.(rn); var = :t)).coeffs
        c_Y_D = (Y_D_.coeffs[end] * fromroots(Complex{BigFloat}.(rd); var = :t)).coeffs
        # Build admittance functions without spurious imaginary parts
        Y_N = Polynomial(real.(c_Y_N), :t)
        Y_D = Polynomial(real.(c_Y_D), :t)
        dY_in = polynomial_fraction_derivative(Y_N, Y_D)

        # Extract shunt line
        a[i] = dY_in(1)
        a[end-i+1] = a[i]

        # Prepare next admittance
        Y_N = Y_N - a[i] * Polynomial(BigFloat[0,1], :t) * Y_D
        Y_D = Y_D
        dY_in = polynomial_fraction_derivative(Y_N, Y_D)
    end

    return Coupler(a, b)
end

@inline x(t) = (1 + t^2) / (1 - t^2)

@inline t²(x) = (x-1)/(x+1)

@inline t(x) = sqrt(Complex(t²(x)))

function design_polynomial(family, degree, type_precision = Float64, polyvar = :x)
	p = zeros(type_precision, degree+1)
	p[end] = 1
    if family == :Butterworth
        return Polynomial(p, polyvar)
    elseif family == :Chebyshev
        return convert(Polynomial{type_precision}, ChebyshevT(p, polyvar))
    end
end

function polynomial_fraction_derivative(N, D)
	dN = derivative(N)
	dD = derivative(D)
	return (t) -> ((dN(t)*D(t) - dD(t)*N(t))/((D(t))^2))
end


text:
using ThinningAndShift ; global const T=ThinningAndShift
using Test
using Distributions
using Random ; Random.seed!(0)

@testset "Poisson generator" begin
  t_tot = 567.8
  rate = 123.45
  trains_poisson = T.make_poisson_samples(rate,t_tot)
  rate_num = length(trains_poisson)/t_tot
  @test isapprox(rate,rate_num;rtol=0.1)
  isis = diff(trains_poisson)
  # ISI is exponential. mean same as variance same as 1/rate
  @test isapprox(mean(isis),inv(rate);rtol=0.1)
  @test isapprox(std(isis),mean(isis);rtol=0.2)
end

@testset "Indepdentend processes" begin
  markings = [ [1,],[2,],[3,]]
  p1,p2,p3 = 0.6,0.3,0.1
  motherrat  = 234.5
  markings_probs = [p1,p2,p3]
  jitters_d = fill(T.NoJitter(),3)
  gtas_test = T.GTAS(motherrat,markings,markings_probs,jitters_d)
  t_tot = 500.0
  (t1,t2,t3),_ = T.make_samples_with_parent(gtas_test,t_tot)
  r1_num = length(t1)/t_tot
  r2_num = length(t2)/t_tot
  r3_num = length(t3)/t_tot
  @test isapprox(r1_num,p1*motherrat;rtol=0.1)
  @test isapprox(r3_num,p3*motherrat;rtol=0.1)
  isis1 = diff(t1)
  isis3 = diff(t3)
  @test isapprox(std(isis1),mean(isis1);rtol=0.2)
  @test isapprox(std(isis3),mean(isis3);rtol=0.2)
  rat_an = T.get_expected_rates(gtas_test)
  @test all(isapprox.([r1_num,r2_num,r3_num],rat_an,;rtol=0.1))
end

@testset "cross-correlations, 2D, Gaussian jitter" begin
  markings = [ [1,2],]
  markings_probs = [1.0]
  mySigma = let σ1 = 0.17,σ2 = 0.2, ρ = 0.7
      [ σ1^2 ρ*σ1*σ2 ; ρ*σ1*σ2  σ2^2 ]
  end
  jitters = [T.JitterDistribution( MultivariateNormal(mySigma)) ,]
  lambda_parent = 100.0
  gtas_test = T.GTAS(lambda_parent,markings,markings_probs,jitters)
  trains,t_parent,attr = T.make_samples_with_parent(gtas_test,10_000.0)
  timescov,covboth = T.covariance_density_numerical(trains,0.15,6.0)
  std_diff = sqrt(mySigma[1,1]+mySigma[2,2] - 2*mySigma[1,2])
  density_an = Normal(0.0,std_diff)
  density_an_vals = lambda_parent .* pdf.(density_an,timescov)
  density_num =  @. 0.5(covboth[:,1,2]+covboth[:,2,1])
  @test all(isapprox.(density_an_vals,density_num;atol=10.,rtol=0.2))
end

@testset "some internal stuff for negative corr" begin
  t_tot = 10.0
  rate = 25.0
  train = T.make_poisson_samples(rate,t_tot)
  antiker = T.AntiExponential(0.1)
  d = Exponential(0.1)
  t_now = 0.5 + rand()*3.0
  Cmax = 0.9999
  testidx,testcdf= T.compute_forward_cutprob(t_now,antiker,train,
    50,Cmax)
  @test train[testidx] > t_now
  @test train[testidx-1] < t_now
  idx_end = testidx + length(testcdf) - 1
  @test cdf(d, train[idx_end] - t_now ) > Cmax
  @test cdf(d, train[idx_end-1] - t_now ) < Cmax

  antiker = T.AntiExponential(0.1)
  timehor = T.antikernel_horizon(antiker,1E-5)

  @test isapprox(pdf(Exponential(0.1),timehor),1E-5)
  @test pdf(Exponential(0.1),timehor*0.999) > 1E-5
end

function get_ff_num(_train,Ttot)
  dtbin = 8.0
  rates = T.bin_spikes(_train,dtbin,Ttot)
  var(rates)/mean(rates)
end

@testset "negative corr 1D" begin

  myrate = 100.0
  killratio = 0.59
  rates = [myrate,killratio*myrate]

  markings = [T.Marking([1]),T.AntiMarking([1])]
  jitters = [T.NoJitter(),T.AntiJitterExpSequential(1000.000)]

  gen = T.sAGTAS(1,rates,markings,jitters)

  Ttot = 1E3
  trains = T.make_samples(gen,Ttot);

  rate_num = length(trains[1])/Ttot
  rate_expected = myrate * (1-killratio)

  @test isapprox(rate_num,rate_expected;rtol=0.2)

  # test FF
  ff_expected = sum(rates) / (rates[1]-rates[2])
  ff_num = get_ff_num(trains[1],Ttot)
  @test isapprox(ff_expected,ff_num;rtol=0.2)

  # repeat with the non-sequential one
  jitters = [T.NoJitter(),T.AntiJitterExp(30E-3)]
  gen = T.sAGTAS(1,rates,markings,jitters)
  Ttot = 1E3
  trains = T.make_samples(gen,Ttot);
  rate_num = length(trains[1])/Ttot
  rate_expected = myrate * (1-killratio)
  @test isapprox(rate_num,rate_expected;rtol=0.2)

  # test FF
  ff_expected = sum(rates) / (rates[1]-rates[2])
  ff_num = get_ff_num(trains[1],Ttot)
  @test isapprox(ff_expected,ff_num;rtol=0.2)

  # repeat with interval
  jitters = [T.NoJitter(),T.AntiJitterStepSequential(1000.000)]
  gen = T.sAGTAS(1,rates,markings,jitters)
  Ttot = 1E3
  trains = T.make_samples(gen,Ttot);
  rate_num = length(trains[1])/Ttot
  rate_expected = myrate * (1-killratio)
  @test isapprox(rate_num,rate_expected;rtol=0.2)

  # test FF
  ff_expected = sum(rates) / (rates[1]-rates[2])
  ff_num = get_ff_num(trains[1],Ttot)
  @test isapprox(ff_expected,ff_num;rtol=0.2)

end


text:
x = 1

function times2()
  x * 2
end

println(times2())

for i = 1:1
  global a = 3
  println(a)
end

println(a)
println(a += 1)

# Global variables not optimized by Julia compiler - can be a huge performance overhead

max_stars_repo_name: cadamsmith/julia-programming

max_issues_repo_name: cadamsmith/julia-programming

max_forks_repo_name: cadamsmith/julia-programming


text:
using Jennrich, Test, Random, Manifolds

@testset "CPD" begin
	Random.seed!(1)
	r = 3
	𝒯 = reshape(sum([kron(randn(9), randn(8), randn(7)) for i = 1:r]), (7,8,9))
	A, B, C = jennrich(𝒯, r)
	@test vec(𝒯) ≈ sum([kron(C[:,i], B[:,i], A[:,i]) for i = 1:r])
end

@testset "(Lᵣ,Lᵣ,1)-BTD" begin
	Random.seed!(1)
	n = (9,8,7)
	L = [2,3]
	ℳ = map(l -> Tucker(n, (l,l,1)), L)
	𝒯 = sum([embed(ℳ[i], TuckerPoint(randn(n...), (L[i],L[i],1))) for i = 1:length(L)])
	A,B,C = jennrich(𝒯, L)
	@test size(A) == (n[1], sum(L))
	@test size(B) == (n[2], sum(L))
	@test size(C) == (n[3], length(L))
	@test reshape(𝒯, n[1], :)' ≈ [kron(C[:,1], B[:,1:L[1]]) kron(C[:,2], B[:,L[1]+1:end])] * A'
end


text:
"""
    makehomogeneous(mag::NIVolume; σ_mm=7, nbox=15)

Homogeneity correction for NIVolume from NIfTI files.

###  Keyword arguments:

- `σ_mm`: σ size for smoothing to obtain bias field. Takes NIfTI voxel size into account
- `nbox`: Number of boxes in each dimension for the box-segmentation step.

"""
function makehomogeneous(mag::NIVolume, datatype=eltype(mag); σ_mm=7, nbox=15)
    return makehomogeneous!(datatype.(mag); σ=mm_to_vox(σ_mm, mag), nbox)
end

"""
    makehomogeneous(mag; σ, nbox=15)

Homogeneity correction of 3D arrays. 4D volumes are corrected using the first 3D volume to
obtain the bias field.

###  Keyword arguments:

- `σ`: σ size in voxel for smoothing to obtain bias field. (mandatory)
- `nbox`: Number of boxes in each dimension for the box-segmentation step.

Larger σ-values make the bias field smoother, but might not be able to catch the
inhomogeneity. Smaller values can catch fast varying inhomogeneities but new inhomogeneities
might be created. The stronger the bias field, the more boxes are required for segmentation.
With too many boxes, it can happen that big darker structures are captured and appear
overbrightened.

"""
function makehomogeneous(mag, datatype=eltype(mag); σ, nbox=15)
    return makehomogeneous!(datatype.(mag); σ, nbox)
end
function makehomogeneous!(mag; σ, nbox=15)
    lowpass = getsensitivity(mag; σ, nbox)
    if eltype(mag) <: AbstractFloat
        mag ./= lowpass
    else # Integer doesn't support NaN
        lowpass[isnan.(lowpass) .| (lowpass .<= 0)] .= typemax(eltype(lowpass))
        mag .= div.(mag, lowpass ./ 2048) .|> x -> min(x, typemax(eltype(mag)))
    end
    mag
end

function getpixdim(nii::NIVolume)
    pixdim = nii.header.pixdim[2:(1+ndims(nii))]
    if all(pixdim .== 1)
        println("Warning! All voxel dimensions are 1 in NIfTI header, maybe they are wrong.")
    end
    return pixdim
end

mm_to_vox(mm, nii::NIVolume) = mm_to_vox(mm, getpixdim(nii))
mm_to_vox(mm, pixdim) = mm ./ pixdim

function getsensitivity(mag::NIVolume, datatype=eltype(mag); kw...)
    return getsensitivity(datatype.(mag), getpixdim(mag); kw...)
end
function getsensitivity(mag, pixdim; σ_mm=7, nbox=15)
    return getsensitivity(mag; σ=mm_to_vox(σ_mm, pixdim), nbox)
end
function getsensitivity(mag; σ, nbox=15)
    # segmentation
    firstecho = view(mag,:,:,:,1)
    mask = robustmask(firstecho)
    segmentation = boxsegment(firstecho, mask, nbox)
    # smoothing
    σ1, σ2 = getsigma(σ)
    lowpass = gaussiansmooth3d(firstecho, σ1; mask=segmentation, nbox=8)
    fillandsmooth!(lowpass, mean(firstecho[mask]), σ2)

    return lowpass
end

# split sigma in two parts
function getsigma(σ)
    factorfinalsmoothing = 0.7
    σ1 = sqrt(1 - factorfinalsmoothing^2) .* σ
    σ2 = factorfinalsmoothing .* σ
    return σ1, σ2
end

function fillandsmooth!(lowpass, stablemean, σ2)
    stablethresh = stablemean / 4
    lowpassmask = (lowpass .< stablethresh) .| isnan.(lowpass) .| (lowpass .> 10 * stablemean)
    lowpass[lowpassmask] .= 3 * stablemean
    lowpassweight = 1.2 .- lowpassmask
    gaussiansmooth3d!(lowpass, σ2; weight=lowpassweight)
end

#threshold(image) = threshold(image, robustmask(image))
function threshold(image, mask; width=0.1)
    m = try quantile(skipmissing(image[mask]), 0.9) catch; 0 end
    return ((1 - width) * m .< image .< (1 + width) * m) .& mask
end

function boxsegment!(image::AbstractArray{<:AbstractFloat}, mask, nbox)
    image[boxsegment(image, mask, nbox)] .= NaN
    return image
end

function boxsegment(image, mask, nbox)
    N = size(image)
    dim = ndims(image)
    boxshift = ceil.(Int, N ./ nbox)

    segmented = zeros(UInt8, size(mask))
    for center in Iterators.product([1:boxshift[i]:N[i] for i in 1:dim]...)
        boxside(d) = max(1, center[d] - boxshift[d]):min(center[d] + boxshift[d], N[d])
        I = CartesianIndices(ntuple(boxside, dim))
        segmented[I] .+= threshold(image[I], mask[I])
    end
    return segmented .* mask .>= 2
end


text:
using Test
using DataDeps
using ParameterEstimocean

using Statistics

@testset "Transformations tests" begin
    data_path = datadep"two_day_suite_2m/free_convection_instantaneous_statistics.jld2";
    field_names = :u

    raw_observations = SyntheticObservations(data_path; field_names = field_names)
    Nz = raw_observations.grid.Nz
    times = raw_observations.times

    map = ConcatenatedOutputMap()

    @info "  Test space and time slicing..."
    time_indices = 101
    z_indices = [collect(20:2:Int(Nz/2)); Nz-4]

    time_transformation = TimeIndices(t=time_indices)
    space_transformation = SpaceIndices(y=:, z=z_indices)

    transformation = Transformation(space = space_transformation,
                                    time = time_transformation,
                                    normalization = nothing)

    observations = SyntheticObservations(data_path; field_names = field_names, transformation)

    sliced_raw_observations = raw_observations.field_time_serieses.u[:, :, z_indices, time_indices][1, 1, :, :]
    ny, nz = size(sliced_raw_observations)
    sliced_raw_observations = reshape(sliced_raw_observations, (ny * nz, 1))

    @test parent(observation_map(map, observations)) ≈ parent(sliced_raw_observations)

    time_indices = [100, 200]
    z_indices = 20:40

    time_transformation = TimeIndices(t=time_indices)
    space_transformation = SpaceIndices(y=:, z=z_indices)

    transformation = Transformation(space = space_transformation,
                                    time = time_transformation,
                                    normalization = ZScore())

    observations = SyntheticObservations(data_path; field_names = field_names, transformation)

    sliced_raw_observations = raw_observations.field_time_serieses.u[:, :, z_indices, time_indices][1, 1, :, :]
    ny, nz = size(sliced_raw_observations)
    sliced_raw_observations = reshape(sliced_raw_observations, (ny * nz, 1))

    # normalize
    μ, σ =  mean(sliced_raw_observations), sqrt(cov(sliced_raw_observations; corrected=false))
    @. sliced_raw_observations = (sliced_raw_observations - μ) / σ

    @test parent(observation_map(map, observations)) ≈ parent(sliced_raw_observations)

    @info "  Test applying space and time weights..."

    time_scale_factor = 0.55
    time_weights = time_scale_factor * ones(size(times))

    depth_scale = 100
    z = raw_observations.grid.zᵃᵃᶜ[1:Nz]
    z_weights = @. exp(z / depth_scale)
    space_weights = reshape(z_weights, 1, 1, Nz)

    transformation = Transformation(space = space_weights,
                                    time = time_weights,
                                    normalization = nothing)

    observations = SyntheticObservations(data_path; field_names = field_names, transformation)

    weighted_raw_observations = raw_observations.field_time_serieses.u[1, 1, 1:Nz, :]
    @. weighted_raw_observations *= z_weights * time_scale_factor
    ny, nz = size(weighted_raw_observations)
    weighted_raw_observations = reshape(weighted_raw_observations, (ny * nz, 1))

    @test parent(observation_map(map, observations)) ≈ parent(weighted_raw_observations)
end


text:
using ReadOnlyArrays, Test, SparseArrays

@testset "ReadOnlyArrays" begin
    x = rand(5, 10)
    vx = view(x, 2:4, 3:8)
    s = sparse(x)
    vs = view(s, 2:4, 3:8)

    for a in [x, vx, s, vs]
        @testset "tests for $(typeof(a))" begin
            r = ReadOnlyArray(a)

            # Check the ability to construct it with type parameters.
            @test ReadOnlyArray{Float64}(a) == r
            @test ReadOnlyArray{Float64,2}(a) == r
            @test ReadOnlyArray{Float64,2,typeof(a)}(a) == r

            @test Base.IteratorSize(typeof(r)) == Base.IteratorSize(typeof(a))
            @test Base.IteratorEltype(typeof(r)) == Base.IteratorEltype(typeof(a))
            @test eltype(typeof(r)) == eltype(r) == eltype(a)
            @test size(r) == size(a)
            @test size(r, 1) == size(a, 1)
            @test size(r, 2) == size(a, 2)
            for (p, q) in zip(r, a)
                @test p == q
            end
            @test r[end] == a[end]
            @test r[1:2, end] == a[1:2, end]
            @test_throws BoundsError r[0]
            @test_throws BoundsError r[end+1]

            @test firstindex(r) == firstindex(a)
            @test lastindex(r) == lastindex(a)

            @test IndexStyle(typeof(r)) == IndexStyle(typeof(a))
            @test iterate(r) == iterate(a)
            @test iterate(r, firstindex(r)) == iterate(a, firstindex(a))
            @test iterate(r, lastindex(r)) == iterate(a, lastindex(a))
            @test iterate(r, iterate(r, lastindex(r))[2]) ==
                  iterate(a, iterate(a, lastindex(a))[2])
            @test length(r) == length(a)

            @test axes(r) == axes(a)
            @test parent(r) === a

            if !issparse(parent(a))
                @test strides(r) == strides(a)
                @test Base.unsafe_convert(Ptr{Float64}, r) ==
                      Base.unsafe_convert(Ptr{Float64}, a)
                @test stride(r, 1) == stride(a, 1)
            end

            # Test the implicit conversion from mutable array to read-only.
            f()::ReadOnlyArray{Float64,2} = a
            fa = f()
            @test fa isa ReadOnlyArray{Float64,2}
            @test fa == a
        end
    end

    @testset "Aliases" begin
        @test ReadOnlyVector([4, 5]) isa ReadOnlyArray{Int, 1, Vector{Int}}
        @test ReadOnlyVector{Int}([4, 5]) isa ReadOnlyArray{Int, 1, Vector{Int}}
        @test ReadOnlyVector{Int,Vector{Int}}([4, 5]) isa ReadOnlyVector{Int, Vector{Int}}
        @test ReadOnlyMatrix([1 2 ; 3 4]) isa ReadOnlyArray{Int, 2, Matrix{Int}}
        @test ReadOnlyMatrix{Int}([1 2 ; 3 4]) isa ReadOnlyArray{Int, 2, Matrix{Int}}
        @test ReadOnlyMatrix{Int,Matrix{Int}}([1 2 ; 3 4]) isa ReadOnlyArray{Int, 2, Matrix{Int}}
    end
end


text:
include("parser.jl")

## accidentally tokenizer works on it
t1=tokenizer("any.x,z.P(x) or not Q(x) and ( R(x) or Q(z))")

## This is my intension
t2= tokenizer("∀.x.∃.z.P(x) | ¬Q(x) & ( R(x) | Q(z))")

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
# ===== VSB Verifications: RunTime_Function Building =====

println(" === BEGINNING Verifications3.jl === ")

# === Imports ===
println("   IMPORTING PACKAGES...")
using VSB
using SimpleVPM
using LinearAlgebra
using Printf

cd("/Users/Damyn/Documents/BYU/FLOW Lab/VSB/verifications")

# === Boundary Geometry ===
println("   GENERATING GEOMETRY...")
NPTS = 10                               # Number of Points on Boundary
del_theta = (2*pi)/(NPTS    )           # Step in theta around body
R = 1.0                                 # Cylinder Radius
k_hat = [0.0, 0.0, 1.0]                 # Z direction unit vector
X_coor = zeros(NPTS)                    # X-coordinates
Y_coor = zeros(NPTS)                    # Y-coordinates
X_tHat = zeros(NPTS)                    # X Tangent Vector
Y_tHat = zeros(NPTS)                    # X Tangent Vector
X_nHat = zeros(NPTS)                    # Y Normal Vector
Y_nHat = zeros(NPTS)                    # Y Normal Vector
overlap = 2.0                           # Particle overlap
sigma = overlap * del_theta             # Guassian Spreading
phase = 0.0 * del_theta
for i=1:NPTS
    X_coor[i] = R*cos((i-1)*del_theta + phase)
    Y_coor[i] = R*sin((i-1)*del_theta + phase)

    X_tHat[i] = sin((i-1)*del_theta + phase)
    Y_tHat[i] = -cos((i-1)*del_theta + phase)

    X_nHat[i] = cos((i-1)*del_theta + phase)
    Y_nHat[i] = sin((i-1)*del_theta + phase)
end
body_pts = [[X_coor[i], Y_coor[i], 0.0] for i in 1:NPTS]
t_hats = [[X_tHat[i], Y_tHat[i], 0.0] for i in 1:NPTS]
n_hats = [[X_nHat[i], Y_nHat[i], 0.0] for i in 1:NPTS]

body = VSB.Boundary(body_pts,t_hats,n_hats)

# === Parameters ===
file_path = splitdir(@__FILE__)[1]
save_path = joinpath(file_path, "sims/verification3_08/")
run_name = "verifications3_08"
nu = 1e0                                   # Kinematic Viscosity
magU_inf = 1.0                              # Magnitude of Free-Stream Velocity
del_t = 0.1                                 # Time step
NSTEPS = 100                                # Number of time steps
Re = magU_inf * R / nu                      # Reynolds Number

U_inf(t) = magU_inf * [1.0, 0.0, 0.0]       # Free-Stream velocity

# === Set up SimpleVPM ===
println("   CREATING PARTICLE FIELD...")
pfield = SimpleVPM.ParticleField(nu,
                                 zeta_gauserf,                # Basis function
                                 g_gauserf;                   # Regularizing function of zeta
                                 transposed=true,             # Transposed scheme
                                 relax=!true,                 # Pedrizzeti relaxation scheme
                                 rlxf=0.3,                    # Relaxation factor
                                 integration="euler",         # Time integration scheme
                                 Uinf=U_inf)                  # Free-Stream velocity function

# === Velocity function ===
function U_slip(X)
    if length(pfield.particles) == 0
        return U_inf(0)
    else
        return SimpleVPM.U(pfield, X) + U_inf(0)
    end
end

# === Initial vorticity field ===
params = VSB.Parameters(body, U_slip, verbose = false, sigma=sigma)
alphas = params.alphas
betas = VSB.CalcVSDiffusionCoefs(body, alphas, del_t, pfield.nu, sigma=sigma)
betas[1] = 0.0
betas[NPTS] = 0.0

offset = 0.01             # Normal direction offset for new particle locations
for i in 1:NPTS
    loc = body_pts[i] + offset .* n_hats[i]
    Gamma = (1/sigma^3) * betas[i] .* k_hat
    vol = 1.0
    particle = SimpleVPM.Particle(loc, Gamma, sigma, vol)
    SimpleVPM.addparticle(pfield, particle)
end

# === RunTimeFunction ===
function RunTimeFunction(pfield, t, del_t)

    if (t != 0)

        for p in pfield.particles
            p.Gamma[1] = 0.0
            p.Gamma[2] = 0.0
        end

        # === Calculate parameters ===
        params = VSB.Parameters(body, U_slip, verbose = false)
        alphas = params.alphas
        betas = VSB.CalcVSDiffusionCoefs(body, alphas, del_t, pfield.nu)

        # === Add particles on surface ===
        offset = 0             # Normal direction offset for new particle locations
        for i in 1:NPTS
            loc = body_pts[i] + offset .* n_hats[i]
            Gamma = (1/sigma^3) * betas[i] .* k_hat
            vol = 1.0
            particle = SimpleVPM.Particle(loc, Gamma, sigma, vol)
            SimpleVPM.addparticle(pfield, particle)
        end

        pfield_gamma = [p.Gamma[3] for p in pfield.particles]
        println("  MAX GAMMA = ",maximum(pfield_gamma))

    end

    # === Returns ===
    return false

end

# === Simulation Parameters ===
println("   RUNNING SIMULATION...")

println("***********************************")
println(" ===== SIMULATION PARAMETERS ===== ")
println(" NU:      ",nu)
println(" U_INF:   ",magU_inf)
println(" NSTEPS:  ",NSTEPS)
println(" NPTS:    ",NPTS)
println(" DEL_T:   ",del_t)
println(" RE:      ",Re)
println("***********************************")

NSTEPS_relax = -1
sigma0 = sigma

paraview = true

SimpleVPM.run_vpm!(pfield, del_t, NSTEPS,
                   runtime_function = RunTimeFunction,
                   nsteps_relax = NSTEPS_relax,
                   rbf_itmax = 1000,
                   sgm0 = sigma0,
                   save_path = save_path,
                   run_name = run_name)

if paraview
    println("   CALLING PARAVIEW...")
    strn = run_name * "...xmf;"
    run(`paraview --data="$save_path$strn"`)
end

println(" === END OF Verifications3.jl === ")


text:
using ADCME
using LinearAlgebra
using AdFem
using DelimitedFiles

A = readdlm("fenics/A.txt")
mesh = Mesh(8,8,1/8)
ρ = constant(ones(get_ngauss(mesh)))
B1 = compute_fem_laplace_matrix1(ones(get_ngauss(mesh)), mesh)
B = compute_fem_laplace_matrix1(ρ, mesh)
sess = Session(); init(sess)
B0 = run(sess, B)
B0 = Array(B0)
@show norm(A - B0)
@info norm(B1-B0)

A = readdlm("fenics/A2.txt")
mesh = Mesh(8, 8, 1.0 / 8,degree=2)
ρ = constant(ones(get_ngauss(mesh)))
B1 = compute_fem_laplace_matrix1(ones(get_ngauss(mesh)), mesh)

B = compute_fem_laplace_matrix1(ρ, mesh)
sess = Session(); init(sess)
B0 = run(sess, B)
B0 = Array(B0)
@info norm(B1-B0)

E = Int64.(readdlm("fenics/edges.txt"))
Edof = get_edge_dof(E, mesh)
DOF = [1:mesh.nnode; Edof .+ mesh.nnode]
B0 = B0[DOF, DOF]

@show norm(A - B0)


text:
const dylibpath="builddir/example.dylib"
@show ccall((:empty,dylibpath),Cint,(Ptr{Cvoid},),C_NULL)

@info "call init_jl_runtime"
@show ccall((:init_jl_runtime,dylibpath),Cvoid,(Ptr{Cvoid},),C_NULL)
@info "try to call bar_segfault"

@show ccall((:bar_segfault,dylibpath),Cint,(Cint,),Cint(100))

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
function gen_sys_ltc_orig_mtk(wiring; T=Float32, name)
  @variables t
  D = Differential(t)

  ϵ = T(1e-5)

  n_neurons = wiring.n_total
  n_in = s_in = wiring.s_in
  w_sens = wiring.matrices[:w_sens]
  w_syns = wiring.matrices[:w_syns]

  sts = @variables begin
    (v[1:n_neurons](t) = rand_uniform(T, -0.2, 0.2, n_neurons))#, [lower=T(-1), upper=T(1)]
    I_sens[1:n_neurons](t)
    I_syns[1:n_neurons](t)
    I_L[1:n_neurons](t)
  end

  ps = @parameters begin
    stim[1:n_in] = fill(T(1), n_in)
    C[1:n_neurons] = rand_uniform(T, 0.4, 0.6, n_neurons), [lower=T(0), upper=T(Inf)]
    G_L[1:n_neurons] = rand_uniform(T, 0.001, 1, n_neurons), [lower=T(0), upper=T(Inf)]
    E_L[1:n_neurons] = rand_uniform(T, -0.2, 0.2, n_neurons)#, [lower=T(-1), upper=T(1)]
  end

  sts = reduce(vcat, collect.(sts))
  ps = reduce(vcat, collect.(ps))

  eqs = Equation[]
  for i in 1:n_neurons
    push!(eqs, I_L[i] ~ G_L[i] * (v[i] - E_L[i]))
    push!(eqs, D(v[i]) ~ -(I_L[i] + I_sens[i] + I_syns[i]) / (C[i] + ϵ))
  end

  systems = ODESystem[]

  I_senss = reshape(Num[0 for _ in 1:s_in*n_neurons], s_in, n_neurons)
  I_synss = reshape(Num[0 for _ in 1:n_neurons*n_neurons], n_neurons, n_neurons)

  for d in 1:n_neurons
    for s in 1:n_in
      w_sens[s,d] == 0 && continue

      ss = Symbol("G_sens$(d)")
      G = @parameters $ss = rand_uniform(T, 0.001, 1), [lower=T(0), upper=T(Inf)]
      ss = Symbol("E_sens$(d)")
      E = @parameters $ss = wiring.matrices[:E_sens][s,d]#, [lower=T(-1), upper=T(1)]
      ss = Symbol("vh_sens$(s)_$(d)")
      vh = @parameters $ss = rand_uniform(T, 0.3, 0.8)#, [lower=T(0), upper=T(1)]
      ss = Symbol("k_sens$(s)_$(d)")
      k = @parameters $ss = rand_uniform(T, 3, 8)#, [lower=T(1), upper=T(10)]

      I_senss[s,d] = G[1] * mysigm(k[1] * (fff(stim[s]) - vh[1])) * (v[d] - E[1])

      push!(ps, G[1], E[1], vh[1], k[1])
    end

    for s in 1:n_neurons
      w_syns[s,d] == 0 && continue

      ss = Symbol("G_syns$(d)")
      G = @parameters $ss = rand_uniform(T, 0.001, 1), [lower=T(0), upper=T(Inf)]
      ss = Symbol("E_syns$(d)")
      E = @parameters $ss = wiring.matrices[:E_sens][s,d]#, [lower=T(-1), upper=T(1)]
      ss = Symbol("vh_syns$(s)_$(d)")
      vh = @parameters $ss = rand_uniform(T, 0.3, 0.8)#, [lower=T(0), upper=T(1)]
      ss = Symbol("k_syns$(s)_$(d)")
      k = @parameters $ss = rand_uniform(T, 3, 8)#, [lower=T(1), upper=T(10)]

      I_senss[s,d] = G[1] * mysigm(k[1] * (v[s] - vh[1])) * (v[d] - E[1])

      push!(ps, G[1], E[1], vh[1], k[1])
    end
  end

  for i in 1:n_neurons
    push!(eqs, I_sens[i] ~ sum(I_senss[:,i]))
    push!(eqs, I_syns[i] ~ sum(I_synss[:,i]))
  end

  ODESystem(eqs, t, sts, ps; systems, name)
end

#Zygote.@nograd gen_sys_ltc_syn_state

function LTCOrigMTK(wiring, solver, sensealg;
  T=Float32, tspan=T.((0, 1)),
  #wiring=nothing,
  tgrad=true, gen_jac=false, kwargs...)

  @named sys = gen_sys_ltc_mtk(wiring)
  ssys = structural_simplify(sys)
  defs = ModelingToolkit.get_defaults(ssys)

  prob = ODEProblem(ssys, defs, tspan, tgrad=tgrad, jac=gen_jac)

  lb, ub = BCTRNN.get_bounds(ssys, wiring.s_in; T)

  rnncell = BCTRNNCell(wiring, solver, sensealg, prob, lb, ub; kwargs...)
  MyRecur(rnncell), defs, ssys
end


text:
rocblas_size(t::Char, M::ROCVecOrMat) = (size(M, t=='N' ? 1 : 2), size(M, t=='N' ? 2 : 1))

const ROCBLASArray{T<:ROCBLASFloat} = ROCArray{T}

###########
#
# BLAS 1
#
###########

LinearAlgebra.rmul!(x::ROCArray{<:ROCBLASFloat}, k::Number) =
  scal!(length(x), convert(eltype(x), k), x, 1)

LinearAlgebra.rmul!(x::ROCArray{<:ROCBLASComplex}, k::Real) =
  scal!(length(x), convert(real(eltype(x)), k), x, 1)

LinearAlgebra.rmul!(x::ROCArray{<:ROCBLASComplex}, k::Number) =
  scal!(length(x), convert(eltype(x), k), x, 1)

# Work around ambiguity with GPUArrays wrapper
LinearAlgebra.rmul!(x::ROCArray{<:ROCBLASFloat}, k::Real) =
  invoke(rmul!, Tuple{typeof(x), Number}, x, k)

function LinearAlgebra.BLAS.dot(DX::ROCArray{T}, DY::ROCArray{T}) where T<:Union{Float32,Float64}
    n = length(DX)
    n==length(DY) || throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))"))
    dot(n, DX, 1, DY, 1)
end

function LinearAlgebra.BLAS.dotc(DX::ROCArray{T}, DY::ROCArray{T}) where T<:Union{ComplexF32,ComplexF64}
    n = length(DX)
    n==length(DY) || throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))"))
    dotc(n, DX, 1, DY, 1)
end

function LinearAlgebra.BLAS.dot(DX::ROCArray{T}, DY::ROCArray{T}) where T<:Union{ComplexF32,ComplexF64}
    dotc(DX, DY)
end

function LinearAlgebra.BLAS.dotu(DX::ROCArray{T}, DY::ROCArray{T}) where T<:Union{ComplexF32,ComplexF64}
    n = length(DX)
    n==length(DY) || throw(DimensionMismatch("dot product arguments have lengths $(length(DX)) and $(length(DY))"))
    dotu(n, DX, 1, DY, 1)
end

#LinearAlgebra.norm(x::ROCBLASArray) = nrm2(x)
LinearAlgebra.BLAS.asum(x::ROCBLASArray) = asum(length(x), x, 1)

function LinearAlgebra.axpy!(alpha::Number, x::ROCArray{T}, y::ROCArray{T}) where T<:ROCBLASFloat
    length(x)==length(y) || throw(DimensionMismatch(""))
    axpy!(length(x), convert(T,alpha), x, 1, y, 1)
end

Base.argmin(xs::ROCBLASArray{<:ROCBLASReal}) = iamin(xs)
Base.argmax(xs::ROCBLASArray{<:ROCBLASReal}) = iamax(xs)

############
#
# BLAS 2
#
############

#########
# GEMV
##########
function gemv_wrapper!(y::ROCVector{T}, tA::Char, A::ROCMatrix{T}, x::ROCVector{T},
                       alpha = one(T), beta = zero(T)) where T<:ROCBLASFloat
    mA, nA = rocblas_size(tA, A)
    if nA != length(x)
        throw(DimensionMismatch("second dimension of A, $nA, does not match length of x, $(length(x))"))
    end
    if mA != length(y)
        throw(DimensionMismatch("first dimension of A, $mA, does not match length of y, $(length(y))"))
    end
    if mA == 0
        return y
    end
    if nA == 0
        return rmul!(y, 0)
    end
    gemv!(tA, alpha, A, x, beta, y)
end

LinearAlgebra.mul!(Y::ROCVector{T}, A::ROCMatrix{T}, B::ROCVector{T}) where T<:ROCBLASFloat = gemv_wrapper!(Y, 'N', A,  B)
LinearAlgebra.lmul!(Y::ROCVector{T}, A::LinearAlgebra.Transpose{<:Any, ROCMatrix{T}}, B::ROCVector{T}) where T<:ROCBLASFloat = gemv_wrapper!(Y, 'T', A.parent, B)
LinearAlgebra.lmul!(Y::ROCVector{T}, A::LinearAlgebra.Adjoint{<:Any, ROCMatrix{T}}, B::ROCVector{T}) where T<:ROCBLASFloat = gemv_wrapper!(Y, 'T', A.parent, B)
LinearAlgebra.lmul!(Y::ROCVector{T}, A::LinearAlgebra.Adjoint{<:Any, ROCMatrix{T}}, B::ROCVector{T}) where T<:ROCBLASComplex = gemv_wrapper!(Y, 'C', A.parent, B)

############
#
# BLAS 3
#
############

########
# GEMM
########
function gemm_wrapper!(C::ROCVecOrMat{T}, tA::Char, tB::Char,
                   A::ROCVecOrMat{T},
                   B::ROCVecOrMat{T},
                   alpha = one(T),
                   beta = zero(T)) where T <: ROCBLASFloat
    mA, nA = rocblas_size(tA, A)
    mB, nB = rocblas_size(tB, B)

    if nA != mB
        throw(DimensionMismatch("A has dimensions ($mA,$nA) but B has dimensions ($mB,$nB)"))
    end

    if C === A || B === C
        throw(ArgumentError("output matrix must not be aliased with input matrix"))
    end

    if mA == 0 || nA == 0 || nB == 0
        if size(C) != (mA, nB)
            throw(DimensionMismatch("C has dimensions $(size(C)), should have ($mA,$nB)"))
        end
        return LinearAlgebra.rmul!(C, 0)
    end

    gemm!(tA, tB, alpha, A, B, beta, C)
end

# Mutating
LinearAlgebra.mul!(C::ROCMatrix{T}, A::ROCVecOrMat{T}, B::ROCVecOrMat{T}) where T<:ROCBLASFloat = gemm_wrapper!(C, 'N', 'N', A, B)
LinearAlgebra.mul!(C::ROCMatrix{T}, trA::LinearAlgebra.Transpose{<:Any, <:ROCMatrix{T}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    gemm_wrapper!(C, 'T', 'N', parent(trA), B)
LinearAlgebra.mul!(C::ROCMatrix{T}, A::ROCMatrix{T}, trB::LinearAlgebra.Transpose{<:Any, <:ROCMatrix{T}}) where T<:ROCBLASFloat =
    gemm_wrapper!(C, 'N', 'T', A, parent(trB))
LinearAlgebra.mul!(C::ROCMatrix{T}, trA::LinearAlgebra.Transpose{<:Any, <:ROCMatrix{T}}, trB::LinearAlgebra.Transpose{<:Any, <:ROCMatrix{T}}) where T<:ROCBLASFloat =
    gemm_wrapper!(C, 'T', 'T', parent(trA), parent(trB))
LinearAlgebra.mul!(C::ROCMatrix{T}, adjA::LinearAlgebra.Adjoint{<:Any, ROCMatrix{T}}, B::ROCMatrix{T}) where T<:ROCBLASReal =
    gemm_wrapper!(C, 'T', 'N', parent(adjA), B)
LinearAlgebra.mul!(C::ROCMatrix{T}, adjA::LinearAlgebra.Adjoint{<:Any, <:ROCMatrix{T}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    gemm_wrapper!(C, 'C', 'N', parent(adjA), B)
LinearAlgebra.mul!(C::ROCMatrix{T}, A::ROCMatrix{T}, adjB::LinearAlgebra.Adjoint{<:Any, ROCMatrix{T}}) where T<:ROCBLASReal =
    gemm_wrapper!(C, 'N', 'T', A, parent(adjB))
LinearAlgebra.mul!(C::ROCMatrix{T}, A::ROCMatrix{T}, adjB::LinearAlgebra.Adjoint{<:Any, <:ROCMatrix{T}}) where T<:ROCBLASFloat =
    gemm_wrapper!(C, 'N', 'C', A, parent(adjB))
LinearAlgebra.mul!(C::ROCMatrix{T}, adjA::LinearAlgebra.Adjoint{<:Any, ROCMatrix{T}}, adjB::LinearAlgebra.Adjoint{<:Any, ROCMatrix{T}}) where T<:ROCBLASReal =
    gemm_wrapper!(C, 'T', 'T', parent(adjA), parent(adjB))
LinearAlgebra.mul!(C::ROCMatrix{T}, adjA::LinearAlgebra.Adjoint{<:Any, <:ROCMatrix{T}}, adjB::LinearAlgebra.Adjoint{<:Any, <:ROCMatrix{T}}) where T<:ROCBLASFloat =
    gemm_wrapper!(C, 'C', 'C', parent(adjA), parent(adjB))
LinearAlgebra.mul!(C::ROCMatrix{T}, trA::LinearAlgebra.Transpose{<:Any, <:ROCMatrix{T}}, adjB::LinearAlgebra.Adjoint{T, <:ROCMatrix{T}}) where T<:ROCBLASReal =
    gemm_wrapper!(C, 'T', 'T', parent(trA), parent(adjB))
LinearAlgebra.mul!(C::ROCMatrix{T}, trA::LinearAlgebra.Transpose{<:Any, <:ROCMatrix{T}}, adjB::LinearAlgebra.Adjoint{<:Any, <:ROCMatrix{T}}) where T<:ROCBLASFloat =
    gemm_wrapper!(C, 'T', 'C', parent(trA), parent(adjB))
LinearAlgebra.mul!(C::ROCMatrix{T}, adjA::LinearAlgebra.Adjoint{T, <:ROCMatrix{T}}, trB::LinearAlgebra.Transpose{<:Any, <:ROCMatrix{T}}) where T<:ROCBLASReal =
    gemm_wrapper!(C, 'T', 'T', parent(adjA), parent(trB))
LinearAlgebra.mul!(C::ROCMatrix{T}, adjA::LinearAlgebra.Adjoint{<:Any, <:ROCMatrix{T}}, trB::LinearAlgebra.Transpose{<:Any, <:ROCMatrix{T}}) where T <: ROCBLASFloat =
    gemm_wrapper!(C, 'C', 'T', parent(adjA), parent(trB))

########
# TRSM
########

# ldiv!
## No transpose/adjoint
LinearAlgebra.ldiv!(A::UpperTriangular{T,ROCMatrix{T}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('L', 'U', 'N', 'N', one(T), parent(A), B)
LinearAlgebra.ldiv!(A::UnitUpperTriangular{T,ROCMatrix{T}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('L', 'U', 'N', 'U', one(T), parent(A), B)
LinearAlgebra.ldiv!(A::LowerTriangular{T,ROCMatrix{T}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('L', 'L', 'N', 'N', one(T), parent(A), B)
LinearAlgebra.ldiv!(A::UnitLowerTriangular{T,ROCMatrix{T}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('L', 'L', 'N', 'U', one(T), parent(A), B)
## Adjoint
LinearAlgebra.ldiv!(A::Adjoint{T,UpperTriangular{T,ROCMatrix{T}}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('L', 'U', 'C', 'N', one(T), parent(parent(A)), B)
LinearAlgebra.ldiv!(A::Adjoint{T,UnitUpperTriangular{T,ROCMatrix{T}}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('L', 'U', 'C', 'U', one(T), parent(parent(A)), B)
LinearAlgebra.ldiv!(A::Adjoint{T,LowerTriangular{T,ROCMatrix{T}}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('L', 'L', 'C', 'N', one(T), parent(parent(A)), B)
LinearAlgebra.ldiv!(A::Adjoint{T,UnitLowerTriangular{T,ROCMatrix{T}}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('L', 'L', 'C', 'U', one(T), parent(parent(A)), B)
## Transpose
LinearAlgebra.ldiv!(A::Transpose{T,UpperTriangular{T,ROCMatrix{T}}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('L', 'U', 'T', 'N', one(T), parent(parent(A)), B)
LinearAlgebra.ldiv!(A::Transpose{T,UnitUpperTriangular{T,ROCMatrix{T}}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('L', 'U', 'T', 'U', one(T), parent(parent(A)), B)
LinearAlgebra.ldiv!(A::Transpose{T,LowerTriangular{T,ROCMatrix{T}}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('L', 'L', 'T', 'N', one(T), parent(parent(A)), B)
LinearAlgebra.ldiv!(A::Transpose{T,UnitLowerTriangular{T,ROCMatrix{T}}}, B::ROCMatrix{T}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('L', 'L', 'T', 'U', one(T), parent(parent(A)), B)

# rdiv!
## No transpose/adjoint
LinearAlgebra.rdiv!(A::ROCMatrix{T}, B::UpperTriangular{T,ROCMatrix{T}}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('R', 'U', 'N', 'N', one(T), parent(B), A)
LinearAlgebra.rdiv!(A::ROCMatrix{T}, B::UnitUpperTriangular{T,ROCMatrix{T}}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('R', 'U', 'N', 'U', one(T), parent(B), A)
LinearAlgebra.rdiv!(A::ROCMatrix{T}, B::LowerTriangular{T,ROCMatrix{T}}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('R', 'L', 'N', 'N', one(T), parent(B), A)
LinearAlgebra.rdiv!(A::ROCMatrix{T}, B::UnitLowerTriangular{T,ROCMatrix{T}}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('R', 'L', 'N', 'U', one(T), parent(B), A)
## Adjoint
LinearAlgebra.rdiv!(A::ROCMatrix{T}, B::Adjoint{T,UpperTriangular{T,ROCMatrix{T}}}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('R', 'U', 'C', 'N', one(T), parent(parent(B)), A)
LinearAlgebra.rdiv!(A::ROCMatrix{T}, B::Adjoint{T,UnitUpperTriangular{T,ROCMatrix{T}}}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('R', 'U', 'C', 'U', one(T), parent(parent(B)), A)
LinearAlgebra.rdiv!(A::ROCMatrix{T}, B::Adjoint{T,LowerTriangular{T,ROCMatrix{T}}}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('R', 'L', 'C', 'N', one(T), parent(parent(B)), A)
LinearAlgebra.rdiv!(A::ROCMatrix{T}, B::Adjoint{T,UnitLowerTriangular{T,ROCMatrix{T}}}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('R', 'L', 'C', 'U', one(T), parent(parent(B)), A)
## Transpose
LinearAlgebra.rdiv!(A::ROCMatrix{T}, B::Transpose{T,UpperTriangular{T,ROCMatrix{T}}}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('R', 'U', 'T', 'N', one(T), parent(parent(B)), A)
LinearAlgebra.rdiv!(A::ROCMatrix{T}, B::Transpose{T,UnitUpperTriangular{T,ROCMatrix{T}}}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('R', 'U', 'T', 'U', one(T), parent(parent(B)), A)
LinearAlgebra.rdiv!(A::ROCMatrix{T}, B::Transpose{T,LowerTriangular{T,ROCMatrix{T}}}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('R', 'L', 'T', 'N', one(T), parent(parent(B)), A)
LinearAlgebra.rdiv!(A::ROCMatrix{T}, B::Transpose{T,UnitLowerTriangular{T,ROCMatrix{T}}}) where T<:ROCBLASFloat =
    rocBLAS.trsm!('R', 'L', 'T', 'U', one(T), parent(parent(B)), A)


text:
__precompile__()

module Animated_CMA_ES

using Plots, LinearAlgebra, StaticArrays

export cma_es, cma_es_animated , new_cma_es

include("CMA_ES.jl")
include("cma_es_animated.jl")
include("new_cma_es.jl")

end # module


text:
module Models

function colpitts(t, u, p)
   M = 2

   c = [p["c21"], p["c32"], p["c13"]]

   du = zeros(3*M)

   for i=0:M-1
      x1, x2, x3 = u[i*3 + 1:i*3 + 3]
      du[3*i + 1] = p["p1"]*x2 + (c[i%3 + 1])*(u[(3*(i + 1) + 1) % (3*M)] - x1)
      du[3*i + 2] = -p["p2"]*(x1 + x3) - p["p4"]*x2
      du[3*i + 3] = p["p3"][i%3 + 1]*(x2 + 1 - exp(-x1))
   end
   return du
end

colpitts_true = (t, u)->colpitts(t, u, Dict("p1" => 5.0, "p2" => 0.0797,
                                            "p3" => 3*[3.0, 3.5, 4.0],
                                            "p4" => 0.6898, "c21" => 0.05,
                                            "c32" => 0.1, "c13" => 0.15))

colpitts_err = (t, u)->colpitts(t, u, Dict("p1" => 5.0 + 0.1,
                                           "p2" => 0.0797 + 0.01,
                                           "p3" => 3*[3.0, 3.5, 4.0],
                                           "p4" => 0.6898, "c21" => 0.05,
                                           "c32" => 0.1, "c13" => 0.15))

function chua(t, u, p)
   x, y, z = u

   du = zeros(3)

   f = p["m_1"]*x + 0.5*(p["m_0"] - p["m_1"])*(abs(x + 1) - abs(x - 1))

   du[1] = p["α"]*(y - x - f)
   du[2] = x - y + z
   du[3] = -p["β"]*y

   return du
end

chua_true = (t, u)->chua(t, u, Dict("α" => 15.6, "β" => 25.58, "m_1" => -5/7,
                                    "m_0" => -8/7))

chua_err = (t, u)->chua(t, u, Dict("α" => 15.7, "β" => 24.58, "m_1" => -5/7,
                                    "m_0" => -8/7))

function osc(t, u, p)
   x, y, z, uu, v = u

   du = zeros(5)
   du[1] = p["σ"]*(u[2]-u[1]) + p["c"]*uu
   du[2] = u[1]*(p["ρ"]-u[3]) - u[2]
   du[3] = u[1]*u[2] - p["β"]*u[3]
   du[4] = v
   du[5] = -p["Ω"]^2*uu

   return du
end

osc_true = (t, u)->osc(t, u, Dict("Ω" => 0.3, "σ" => 10, "β" => 8/3,
                                  "ρ" => 28, "c" => 5))

osc_err = (t, u)->osc(t, u, Dict("Ω" => 0.32, "σ" => 10, "β" => 8/3,
                                 "ρ" => 28, "c" => 5.1))

end


text:
push!(LOAD_PATH,joinpath(dirname(@__FILE__),"../src/"))
using Documenter
using ModelicaScriptingTools
using OMJulia

makedocs(
    sitename="ModelicaScriptingTools.jl",
    pages = [
        "index.md",
        "api.md",
        "documentation.md"
    ]
)
deploydocs(
    repo = "github.com/THM-MoTE/ModelicaScriptingTools.jl.git",
    devbranch = "main",
    versions = ["v^", "v#.#", "stable" => "v^"]
)

max_stars_repo_name: THM-MoTE/MoST

max_issues_repo_name: THM-MoTE/MoST

max_forks_repo_name: THM-MoTE/MoST


text:
using Test
using BinaryWrappers
import lib4ti2_jll

function run_4ti2(dir::String, tool::String, example::String; wrapper=false)
    mktempdir() do tmp
        cp(dir, joinpath(tmp, example))
        tmpexample = joinpath(tmp, example, example)
        # at least one output file should be named in the outputname file
        outfile = readchomp(joinpath(tmp, example, "outputname"))
        chomp(outfile)

        if wrapper
            @test success(pipeline(`$tool $tmpexample`,devnull, devnull))
        else
            jlltool = getfield(lib4ti2_jll,Symbol(tool))
            if typeof(jlltool) == String
               # shell script wont work like this
               @test_broken success(pipeline(`$(jlltool()) $tmpexample`, devnull, devnull))
               return nothing
            else
               @test success(pipeline(`$(jlltool()) $tmpexample`, devnull, devnull))
            end
        end
        return readchomp(joinpath(tmp, example, outfile))
    end
end

@testset "lib4ti2 examples" begin
    binpath = @generate_wrappers(lib4ti2_jll)
    ENV["PATH"] = binpath * ":" * ENV["PATH"]

    @test isfile(joinpath(binpath, "circuits"))
    @test isfile(joinpath(binpath, "graver"))
    @test isfile(joinpath(binpath, "groebner"))
    @test isfile(joinpath(binpath, "hilbert"))
    @test isfile(joinpath(binpath, "zsolve"))

    path = joinpath(@__DIR__, "examples", "4ti2")

    for tool in readdir(path)
        for example in readdir(joinpath(path, tool))
            dir = joinpath(path, tool, example)

            @testset "$tool: $example" begin
                jllout = run_4ti2(dir, tool, example)
                wrapperout = run_4ti2(dir, tool, example; wrapper=true)
                if jllout != nothing && wrapperout != nothing
                    @test jllout == wrapperout
                end
            end
        end
    end
end


text:
module BitTwiddle

include("tools.jl")
include("adt.jl")

end # module


text:
"""
    readvcf("vcffile.vcf.gz")

Read a vcf file and return a vcf data array of all lines.
The vcf file may be compressed (bgzip) or uncompressed.
"""
function readbigvcf(vcf)

    isfile(vcf) || error("$vcf does not exist")
    bn = splitext(vcf)
    n = bn[1]

    if occursin(r".gz$", vcf)
        f = bgzip_file(vcf)
        d = readlines(f);
        bgzip_file("$n")
    else
        d = readlines(vcf);
    end

    return d;

end

max_stars_repo_name: ScottWatkins/STRdiff

max_issues_repo_name: ScottWatkins/STRdiff

max_forks_repo_name: ScottWatkins/STRdiff