text:
module HalfIntegerArrays

using HalfIntegers
import HalfIntegers: HalfIntegerOrInteger
const HalfIntOrInt = Union{HalfInt,Int}
using LinearAlgebra
using Base: Indices, @propagate_inbounds
@static if !isdefined(Base, :IdentityUnitRange)
    const IdentityUnitRange = Base.Slice
else
    using Base: IdentityUnitRange
end
using OffsetArrays
import Base: throw_boundserror

export HalfIntArray
export SpinMatrix
export CartesianIndicesHalfInt
export CartesianIndexHalfInt
export LinearIndicesHalfInt

"""
	AbstractHalfIntegerArray{T,N}

Supertype for `N`-dimensional arrays (or array-like types) with elements of type `T`.
Such arrays may be indexed using integers as well as half-integers.
"""
abstract type AbstractHalfIntegerArray{T,N} <: AbstractArray{T,N} end

"""
    AbstractHalfIntegerWrapper{T,N} <: AbstractHalfIntegerArray{T,N}

Supertype for `N`-dimensional arrays (or array-like types) with elements of type `T`
that wrap a parent array that has integer indices. The wrapper acts as a map between the
indices of the array and those of the parent.

Subtypes of `AbstractHalfIntegerWrapper` must define a method for `Base.parent` that returns
the underlying array.

[`HalfIntArray`](@ref), [`SpinMatrix`](@ref) and other types are subtypes of this.
"""
abstract type AbstractHalfIntegerWrapper{T,N} <: AbstractHalfIntegerArray{T,N} end

const AbstractHalfIntegerVector{T} = AbstractHalfIntegerArray{T,1}
const AbstractHalfIntegerMatrix{T} = AbstractHalfIntegerArray{T,2}
const AbstractHalfIntegerVecOrMat{T} = Union{AbstractHalfIntegerVector{T},
	AbstractHalfIntegerMatrix{T}}
const AdjOrTransAbsHalfIntVecOrMat{T} = LinearAlgebra.AdjOrTrans{T,<:AbstractHalfIntegerVecOrMat}
const AdjOrTransAbsHalfIntMatrix{T} = LinearAlgebra.AdjOrTrans{T,<:AbstractHalfIntegerMatrix}
const AdjOrTransAbsHalfIntVector{T} = LinearAlgebra.AdjOrTrans{T,<:AbstractHalfIntegerVector}

#= This type is useful to extend Base functions for wrapper types
around AbstractHalfIntegerArrays that are defined in Base or stdlib.
This is because Base functions often rely on CartesianIndices and LinearIndices,
whereas this package uses CartesianIndicesHalfInt and LinearIndicesHalfInt
=#
const AbstractHalfIntegerArrayOrWrapper = Union{AbstractHalfIntegerArray,
	AdjOrTransAbsHalfIntVecOrMat}

unsafeInt(d) = floor(Int, d)

include("axes.jl")
include("lincartindexing.jl")
include("halfintarray.jl")
include("linearalgebra.jl")
include("broadcast.jl")
include("subarray.jl")
include("arrayshow.jl")

# Define functions common to the supertype

Base.keys(h::AbstractHalfIntegerArrayOrWrapper) = CartesianIndicesHalfInt(axes(h))
Base.keys(h::AbstractHalfIntegerArray{<:Any,1}) = LinearIndicesHalfInt(h)
Base.keys(r::IdOffsetRange) = LinearIndicesHalfInt(r)

Base.axes1(A::AbstractHalfIntegerArray{<:Any,0}) = Base.OneTo(1)

_maybetail(::Tuple{}) = ()
_maybetail(t::Tuple) = Base.tail(t)

throw_lininderr() = throw(ArgumentError("linear indexing requires an integer index"))
ensureInt(i) = isinteger(i) || throw_lininderr()
ensureInt(::Int) = true

parenttype(A::AbstractHalfIntegerArray) = parenttype(typeof(A))

function trimtoN(I::Tuple, ::Val{N}) where {N}
    J,_ = Base.IteratorsMD.split(I,Val(N))
    Base.fill_to_length(J,one(Int),Val(N))
end

# Convert Real indices to HalfInt
function Base.to_indices(A::AbstractHalfIntegerArrayOrWrapper, inds, I::Tuple{Real, Vararg{Any}})
    (HalfInt(first(I)), to_indices(A, _maybetail(inds), Base.tail(I))...)
end
# Conversely convert HalfInt indices to Integers if the axis has integer indices
# This should really only happen if the index is already an integer in value
function Base.to_indices(A::AbstractHalfIntegerArrayOrWrapper, inds::Tuple{AbstractUnitRange{Int},Vararg{Any}}, I::Tuple{Real, Vararg{Any}})
    (Int(first(I)), to_indices(A, Base.tail(inds), Base.tail(I))...)
end

# zero dim arrays may be indexed with ones
function Base.to_indices(A::AbstractHalfIntegerArrayOrWrapper, inds::Tuple{}, I::Tuple{Real, Vararg{Any}})
    (Int(first(I)), to_indices(A, (), Base.tail(I))...)
end

# AbstractRanges may be handled by to_index
Base.to_index(::AbstractHalfIntegerArrayOrWrapper, a::UnitRange{Int}) = a
Base.to_index(::AbstractHalfIntegerArrayOrWrapper, a::UnitRange{<:Real}) = UnitRange{HalfInt}(a)

for DT in [:Int, :Real]
    @eval Base.error_if_canonical_getindex(::IndexLinear, A::AbstractHalfIntegerArray, ::$DT) =
        error("getindex not defined for ", typeof(A))
    @eval Base.error_if_canonical_getindex(::IndexCartesian, A::AbstractHalfIntegerArray{T,N}, ::Vararg{$DT,N}) where {T,N} =
        error("getindex not defined for ", typeof(A))

    @eval Base.error_if_canonical_setindex(::IndexLinear, A::AbstractHalfIntegerArray, ::$DT) =
        error("setindex! not defined for ", typeof(A))
    @eval Base.error_if_canonical_setindex(::IndexCartesian, A::AbstractHalfIntegerArray{T,N}, ::Vararg{$DT,N}) where {T,N} =
        error("setindex! not defined for ", typeof(A))
end

@propagate_inbounds function Base.getindex(A::AbstractHalfIntegerArrayOrWrapper,
    I::Union{Real,CartesianIndex,CartesianIndexHalfInt}...)

    Base.error_if_canonical_getindex(IndexStyle(A), A, I...)
    J = trimtoN(to_indices(A, I), Val(ndims(A)))
    A[J...]
end

@propagate_inbounds function Base.setindex!(A::AbstractHalfIntegerArrayOrWrapper, v,
    I::Union{Real, CartesianIndex, CartesianIndexHalfInt}...)

    Base.error_if_canonical_setindex(IndexStyle(A), A, I...)
    J = trimtoN(to_indices(A, I), Val(ndims(A)))
    A[J...] = v
    A
end

#=
indicescompatible is meant to be called after checkbounds. This ensures that arrays
with half-integer axes are indexed with half-integers, and arrays with integer axes are
indexed with integers.
=#
indicescompatible(::Tuple{}, ::Tuple{Real, Vararg{Any}}) = true
indicescompatible(::Tuple{HalfInt}, ::Tuple{}) = true
indicescompatible(::Tuple{}, ::Tuple{}) = true
function indicescompatible(offsets::Tuple{HalfInt,Vararg{Any}}, I::Tuple{Real, Vararg{Any}})
    isinteger(first(offsets) + first(I)) && indicescompatible(Base.tail(offsets), Base.tail(I))
end
throw_indicesincompatible() = throw(ArgumentError("Indices provided are incompatible with the array's axes"))

# Linear indexing
for DT in [:Integer, :Real]
    @eval function Base.isassigned(a::AbstractHalfIntegerArrayOrWrapper, i::$DT...)
        try
            a[i...]
            true
        catch e
            if isa(e, BoundsError) || isa(e, UndefRefError) || isa(e,InexactError)
                return false
            else
                rethrow()
            end
        end
    end
end

Base.pairs(::IndexLinear, A::AbstractHalfIntegerArrayOrWrapper) = Iterators.Pairs(A, LinearIndicesHalfInt(A))
Base.pairs(::IndexCartesian, A::AbstractHalfIntegerArrayOrWrapper) = Iterators.Pairs(A, CartesianIndicesHalfInt(axes(A)))

Base.eachindex(::IndexCartesian, A::AbstractHalfIntegerArrayOrWrapper) = CartesianIndicesHalfInt(axes(A))

function _all_match_first_shadow(f::F, inds, A, B...) where F<:Function
    (inds == f(A)) & _all_match_first_shadow(f, inds, B...)
end
_all_match_first_shadow(f::F, inds) where F<:Function = true

@inline function Base.eachindex(::IndexCartesian, A::AbstractHalfIntegerArrayOrWrapper, B::AbstractHalfIntegerArrayOrWrapper...)
    axsA = axes(A)
    _all_match_first_shadow(axes, axsA, B...) || Base.throw_eachindex_mismatch(IndexCartesian(), A, B...)
    CartesianIndicesHalfInt(axsA)
end

function Base.collect(h::AbstractHalfIntegerArrayOrWrapper)
    b = Array{eltype(h)}(undef,size(h))
    @inbounds for (bi,hi) in zip(eachindex(b), eachindex(h))
        b[bi] = h[hi]
    end
    b
end
Base.Array{T,N}(h::AbstractHalfIntegerArray{T,N}) where {T,N} = collect(h)

# Similar and reshape

Base.similar(A::AbstractHalfIntegerArrayOrWrapper, ::Type{T}, dims::Dims) where T =
    similar(unwraphalfint(A), T, dims)

# Need to import OffsetAxis to get around the type-piracy in OffsetArrays
# See https://github.com/JuliaArrays/OffsetArrays.jl/issues/87#issuecomment-581391453
import OffsetArrays: OffsetAxis, OffsetAxisKnownLength
# This lets us avoid ambiguities with Base
const ReshapedAxis = Union{Integer,Base.OneTo}
const ReshapedAxisOneTo = Union{ReshapedAxis, OneTo}
# Axis types internal to this package
const HalfIntAxis = Union{IdOffsetRange,OneTo}
const HalfIntOffsetAxisKnownLength = Union{OffsetAxisKnownLength, HalfIntAxis}
const HalfIntOffsetAxis = Union{OffsetAxis, HalfIntAxis}

# Similar with non-offset axes returns an Array
# OneTo is treated as a non-offset axis
for ax in [:ReshapedAxis, :OneTo, :ReshapedAxisOneTo]
	@eval function Base.similar(A::AbstractHalfIntegerArrayOrWrapper, ::Type{T}, inds::Tuple{$ax,Vararg{$ax}}) where T
	    similar(parent(A), T, map(indexlength,inds))
	end
end

# Similar with offset axes returns a HalfIntArray
for ax in [:HalfIntOffsetAxisKnownLength, :OffsetAxisKnownLength]
    @eval function Base.similar(A::AbstractHalfIntegerArrayOrWrapper, ::Type{T}, inds::Tuple{$ax,Vararg{$ax}}) where T
        B = similar(parent(A), T, map(indexlength, inds))
        return HalfIntArray(B, map(offset, axes(B), inds))
    end
end

for DT in [:OneTo, :ReshapedAxis, :HalfIntOffsetAxisKnownLength, :OffsetAxisKnownLength]
	@eval function Base.similar(A::AbstractHalfIntegerArrayOrWrapper, inds::Tuple{$DT,Vararg{$DT}})
	    similar(A, eltype(A), inds)
	end
end

for DT in [:AbstractArray, :AbstractHalfIntegerArrayOrWrapper]
	@eval function Base.similar(A::$DT, ::Type{T}, inds::Tuple{HalfIntAxis,Vararg{HalfIntAxis}}) where {T}
		P = similar(A, T, map(indexlength, inds))
		HalfIntArray(P, map(offset, axes(P), inds))
	end
end

function Base.similar(::Type{T}, inds::Tuple{HalfIntAxis,Vararg{HalfIntOffsetAxisKnownLength}}) where {T<:AbstractArray}
    P = T(undef, map(indexlength, inds))
    HalfIntArray(P, map(offset, axes(P), inds))
end

for ax in [:HalfIntOffsetAxis, :OffsetAxis]
    @eval Base.reshape(A::AbstractHalfIntegerArrayOrWrapper, inds::Tuple{$ax,Vararg{$ax}}) =
        HalfIntArray(reshape(parent(A), map(indexlength, inds)), map(indexoffset, inds))
end

Base.reshape(A::AbstractHalfIntegerArray, inds::Tuple{OneTo,Vararg{OneTo}}) = reshape(parent(A), Base.OneTo.(inds))

# And for non-offset axes, we can just return a reshape of the parent directly
Base.reshape(A::AbstractHalfIntegerArray, inds::Vararg{Int}) = reshape(A, inds)
Base.reshape(A::AbstractHalfIntegerArray, inds::Vararg{Union{Colon,Int}}) = reshape(A, inds)
Base.reshape(A::AbstractHalfIntegerArray, inds::Vararg{Union{Int,AbstractUnitRange}}) = reshape(A, inds)
Base.reshape(A::AbstractHalfIntegerArray, ::Colon) = A

for DT in [:AbstractHalfIntegerArray, :AbstractHalfIntegerArrayOrWrapper]
	@eval Base.reshape(A::$DT, inds::Tuple{ReshapedAxis,Vararg{ReshapedAxis}}) = reshape(parent(A), inds)
	@eval Base.reshape(A::$DT, inds::Tuple{Int,Vararg{Int}}) = reshape(parent(A), inds)
	@eval Base.reshape(A::$DT, inds::Tuple{Union{Int,Colon},Vararg{Union{Int,Colon}}}) = reshape(parent(A), inds)
	@eval Base.reshape(A::$DT, inds::Tuple{Union{Int,Colon},Vararg{Union{Int,Colon,Base.OneTo}}}) = reshape(parent(A), inds)

	@eval Base.reshape(A::$DT, inds::Tuple{Colon,Vararg{Colon}}) =
    	throw(DimensionMismatch("new dimensions $inds may have at most one omitted dimension specified by `Colon()`"))
	@eval Base.reshape(A::$DT, inds::Tuple{}) =
	    throw(DimensionMismatch("new dimensions () must be consistent with array size $(length(A))"))
end

Base.sizeof(a::AbstractHalfIntegerWrapper) = sizeof(parent(a))

@static if isdefined(Base, :copyto_unaliased!)
    function Base.copyto_unaliased!(deststyle::IndexStyle, dest::AbstractHalfIntegerArrayOrWrapper, srcstyle::IndexStyle, src::AbstractHalfIntegerArrayOrWrapper)
        isempty(src) && return dest
        destinds, srcinds = LinearIndicesHalfInt(dest), LinearIndicesHalfInt(src)
        idf, isf = first(destinds), first(srcinds)
        Δi = idf - isf
        (checkbounds(Bool, destinds, isf+Δi) & checkbounds(Bool, destinds, last(srcinds)+Δi)) ||
            throw(BoundsError(dest, srcinds))
        if deststyle isa IndexLinear
            if srcstyle isa IndexLinear
                # Single-index implementation
                @inbounds for i in srcinds
                    dest[i + Δi] = src[i]
                end
            else
                # Dual-index implementation
                i = idf - 1
                @inbounds for a in src
                    dest[i+=1] = a
                end
            end
        else
            iterdest, itersrc = eachindex(dest), eachindex(src)
            if iterdest == itersrc
                # Shared-iterator implementation
                for I in iterdest
                    @inbounds dest[I] = src[I]
                end
            else
                # Dual-iterator implementation
                ret = iterate(iterdest)
                @inbounds for a in src
                    idx, state = ret
                    dest[idx] = a
                    ret = iterate(iterdest, state)
                end
            end
        end
        return dest
    end
end

end


text:
# --------------------------------------------------------------------------
# ACE.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# All rights reserved.
# --------------------------------------------------------------------------

using JuLIP, ACE, LinearAlgebra

deg = 10
trans = ACE.IdTransform()
J = ACE.TransformedJacobi(deg, trans, 0.5, 3.0)
ship4 = SHIPBasis(SparseSHIP(3, 10; wL = 1.5), J)
J = ship4.J

nargs = 3
function testf(Rs)
   f = 1.0
   for R in Rs
      b = ACE.evaluate(J, norm(R))
      f *= b[3]
   end
   return f
end

# testf(Rs) = prod(norm.(Rs).^2)

ACE.Exp.determine_order(testf, nargs, ship4)

ext: jl


text:
abstract type AbstractStimulusParameter{T} <: AbstractParameter{T} end
abstract type AbstractStimulusAction{T,D} <: AbstractSpaceAction{T,D} end

# Naturally map stimulus arrays
function (stim_params::AbstractArray{<:AbstractStimulusParameter{T}})(space::AbstractSpace{T,N}) where {T,N}
    return map(stim_params) do param
        param(space)
    end
end
# Stimulus actions are applied IN ORDER
function (stims::AbstractArray{<:AbstractStimulusAction})(args...)
    for stim in stims
        stim(args...)
    end
end

struct NoStimulusParameter{T} <: AbstractStimulusParameter{T} end
struct NoStimulusAction{T,N} <: AbstractStimulusAction{T,N} end
function (nostim::NoStimulusParameter{T})(space::AbstractSpace{T,N}) where {T,N}
    NoStimulusAction{T,N}()
end
(ns::NoStimulusAction)(args...) = nothing


text:
using Revise, Pkg, BenchmarkTools

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

Pkg.activate(pwd())
using MinVariance2

#need this since we are uing generated functions
if @isdefined(FIRST_RUN)
  revise(MinVariance2)
end

const FIRST_RUN = true

function minvarscript()
  MinVariance2.getparameters!()
end

@time minvarscript()

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
##################################################
# Test zero factor marginalization
##################################################
bnd = MCBN.BayesNetDAI(3)
MCBN.add_edge!(bnd, 1, 2)
MCBN.add_edge!(bnd, 1, 3)
DAI.clearBackups!(bnd)
n = 0.2
d = 0.1
MCBN.set_factor!(bnd, 1, [0.0, 1.0])
MCBN.set_factor!(bnd, 2, [1-n, d, n, 1-d])
MCBN.set_factor!(bnd, 3, [1-n, d, n, 1-d])
DAI.clearBackups!(bnd)

allvars = DAI.VarSet(DAI.Var(1,2), DAI.Var(2,2))
MCBN.setupIA!(bnd)
jointfac = DAI.marginal(bnd.ia, allvars, true)
margfac = DAI.marginal(bnd.ia, DAI.VarSet(DAI.Var(2,2)), true)
belfac = DAI.belief(bnd.ia, DAI.VarSet(DAI.Var(2,2)))

@test_approx_eq jointfac[1] 0.0
@test_approx_eq jointfac[2] 0.1
@test_approx_eq jointfac[3] 0.0
@test_approx_eq jointfac[4] 0.9

##################################################
# Test CoD from Zhao et al 2011
##################################################
bnd = MCBN.BayesNetDAI(3)
MCBN.add_edge!(bnd, 1, 2)
MCBN.add_edge!(bnd, 1, 3)
DAI.clearBackups!(bnd)
c10 = 0.4
c11 = 1. - c10
n13 = n12 = 0.2
d12 = d13 = 0.1
MCBN.set_factor!(bnd, 1, [c10, c11])
MCBN.set_factor!(bnd, 2, [1-n12, d12, n12, 1-d12])
MCBN.set_factor!(bnd, 3, [1-n13, d13, n13, 1-d13])

MCBN.check_bnd(bnd)
@test_approx_eq MCBN.cod(bnd, 1, [2,3]) 0.675
@test_approx_eq MCBN.cod(bnd, 1, [2]) 0.65
@test_approx_eq MCBN.cod(bnd, 2, [1]) 0.631578947368421

codbnd = bnd
bns = MCBN.BayesNetSampler(codbnd, MCBN.draw_data(codbnd, 50))

@show MCBN.ting_cod(bns, 1, [2,3])
@show MCBN.ting_cod(bns, 1, [2])
@show MCBN.ting_cod(bns, 2, [1])

##################################################
# Test Controllability
##################################################
@test_approx_eq MCBN.controllability(bnd, [2], [1], 1) 0.42
@test_approx_eq MCBN.controllability(bnd, [2], [2], 1) 0.28
@test_approx_eq MCBN.controllability(bnd, [2], [2], 3) 0.0
@test_approx_eq MCBN.controllability(bnd, [1], [2], 2) 0.0
@test_approx_eq MCBN.controllability(bnd, [1,2], [2,2], 3)+1 1.0 # Not much accuracy here
@test_approx_eq MCBN.controllability(bnd, [2,3], [2,2], 1) 0.308


text:
Base.:*(x::Expression,y::Unitful.AbstractQuantity) = x * ModelingToolkit.Constant(y)

instantiate(x::ModelingToolkit.Constant) = x.value
instantiate(x::ModelingToolkit.Variable{Number}) = 1.0
instantiate(x::ModelingToolkit.Variable) = oneunit(1*ModelingToolkit.vartype(x))
function instantiate(x::ModelingToolkit.Operation)
    if x.op isa Variable
        return instantiate(x.op)
    elseif x.op isa Differential
        instantiate(x.args[1])/instantiate(x.args[1].args[1])
    else
        x.op(instantiate.(x.args)...)
    end
end

function validate(eq::ModelingToolkit.Equation)
    try
        return typeof(instantiate(eq.lhs)) == typeof(instantiate(eq.rhs))
    catch
        return false
    end
end

function validate(sys::AbstractODESystem)
    all(validate.(equations(sys)))
end


text:
using DataConvenience

using DataFrames

df = DataFrame(A = 1:4, B = 2:2:8)

dfcor(df)

max_stars_repo_name: xiaodaigh/DataConvenience

max_issues_repo_name: xiaodaigh/DataConvenience

max_forks_repo_name: xiaodaigh/DataConvenience


text:
function newtondist!(yout::Vector{Float64},iters::Vector{Int64},Vs::Float64,vs::Float64,ts::Float64,
    V::Float64,C::Float64,Q::Float64,V0::Float64,beta::Float64,c0::Float64,A0::Float64,W1end::Float64,rho::Float64,
    f::Function,J::Function,maxiter::Int16,epsJ::Float64,epsN::Float64,maxval::Float64,h::Float64)

    # initial guess
    x0 = zeros(2);
    xx = zeros(2);
    x0[2] = V;
    Pest = 1/C.*(x0[2] .- V0);
    Aest = (Pest./beta .+ sqrt.(A0)).^2;
    x0[1] = W1end .+ 8.*(c0 .- sqrt.(0.5.*beta./rho).*Aest.^0.25);

    # non-dimensionalize
    x0[1] /= vs;
    x0[2] /= Vs;

    # setup for iterations
    xx .= x0;
    x = zeros(length(x0));
    N = 1;

    # max step size for line searches
    stpmax = 100*max(sqrt(norm(xx)),length(x0));
    # println(stpmax)

    # Newton iterations
    while N <= maxiter
        # determine Jacobian, check invertibility
        JJ = J(xx,Vs,vs,ts,C,W1end,c0,rho,beta,h);
        D = diagm(maximum!(zeros(length(xx)),abs.(JJ)).^-1);
        # println(JJ)
        # println(D)
        # println(D*JJ)
        # println(cond(D*JJ))
        # println(det(D*JJ))
        if abs(det(D*JJ)) < epsJ
            println(JJ)
            println(D*JJ)
            print(det(D*JJ))
            error("Distal Newton Jacobian is singular.");
        end
        # compute gradient of line search objective function
        fvec = D*f(xx,V,Q,Vs,vs,ts,V0,rho,beta,C,W1end,c0,A0,h);
        # println(f(xx,system,n,state))
        # println(fvec)
        fold = 0.5*dot(fvec,fvec);
        g = (D*JJ)'*fvec;
        # compute newton step
        # println(inv(D*JJ))
        s = -inv(D*JJ)*fvec;
        # line search to update state vector
        fn, xn, check = CVModule.linedist(xx,fold,g,s,stpmax,f,J,maxiter,V,Q,Vs,vs,ts,V0,rho,beta,C,W1end,c0,A0,h);
        # println(xn[1]*vs)
        # check if sufficiently close to root
        JJ = J(xn,Vs,vs,ts,C,W1end,c0,rho,beta,h);
        D = diagm(maximum!(zeros(length(xx)),abs.(JJ)).^-1);
        fvec = D*f(xn,V,Q,Vs,vs,ts,V0,rho,beta,C,W1end,c0,A0,h);
        # println(fvec)
        if norm(fvec) <= epsN
            x = xn;
            x[1] = x[1]*vs;
            x[2] = x[2]*Vs;
            # println(x)
            # println(system.branches.W2root)
            # println(inv(D*JJ))
            # println(fvec)
            iters[1] += N;
            break
        end
        if check
            test = 0.
            den = max(fn,0.5*length(xn));
            for i = 1:length(xn)
                temp = abs(g[i])*max(abs(xn[i]),1.)/den;
                if temp > test
                    test = temp;
                end
            end
            if test < 1e-6
                check = true;
                println("Warning: gradient of objective function close to zero.")
            else
                check = false;
            end
        end
        # check for divergence
        if norm(fvec,Inf) >= maxval
            iters[1] += N;
            println(xn)
            println(fvec)
            error("Newton iteration diverged.");
        end
        N+=1;
        xx = xn;
        if N == maxiter
            println(JJ)
            # println(D)
            # println(D*JJ)
            println(xn)
            println(fvec)
            println(norm(fvec))
            error("Newton iteration failed to converge.");
        end
    end

    # update based on converged solution
    yout[1] = x[2];
    yout[2] = 1./C.*(x[2] .- V0);
    yout[3] = (2.*rho/beta).^2.*((W1end.-x[1])./8 .+ c0).^4;
    yout[4] = yout[3].*0.5.*(W1end .+ x[1]);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
# a simplified CUDA module, adapted from github.com/JuliaGPU/CUDA.jl
export CUDA
module CUDA
export CuPtr
using Compat

@windows? (
begin
  const libcuda = Libdl.find_library(["nvcuda.dll"], [""])
end
: # linux or mac
begin
  const libcuda = Libdl.find_library(["libcuda","libcudart"], [""])
  if isempty(libcuda)
    error("Libcuda not found via Libdl.find_library! Please check installation and ENV configuration")
  end
end)

const driver_error_descriptions = @compat(Dict(
  0 => "Success",
  1 => "Invalid value",
  2 => "Out of memory",
  3 => "Driver not initialized",
  4 => "Driver being shutdown",
  5 => "Profiler disabled",
  6 => "Profiler not initialized",
  7 => "Profiler already started",
  8 => "Profiler already stopped",
  100 => "No CUDA-capable device",
  101 => "Invalid device ordinal",
  200 => "Invalid kernel image",
  201 => "Invalid context",
  202 => "Context already current",
  205 => "Map operation failed",
  206 => "Unmap operation failed",
  207 => "Array mapped",
  208 => "Resource already mapped",
  209 => "No kernel image available/suitable for GPU",
  210 => "Resource already acquired",
  211 => "Resource not mapped",
  212 => "Resource not mapped as array",
  213 => "Resource not mapped as pointer",
  214 => "Uncorrectable ECC error detected",
  215 => "Unsupported limit",
  216 => "Context already in use",
  217 => "Peer access not supported",
  218 => "PTX JIT compilation failed",
  300 => "Invalid kernel source",
  301 => "File not found",
  302 => "Shared object symbol not found",
  303 => "Shared object initialization failed",
  304 => "OS call failed",
  400 => "Invalid handle",
  500 => "Named symbol not found",
  600 => "Not ready",
  700 => "Kernel launch failed",
  701 => "Launch out of resources",
  702 => "Launch timeout",
  703 => "Incompatible texturing mode",
  704 => "Peer access already enabled",
  705 => "Peer access not enabled",
  708 => "Primary context already active",
  709 => "Context destroyed",
  710 => "Assertion triggered failure",
  711 => "Too many peers",
  712 => "Host memory already registered",
  713 => "Host memory not registered",
  714 => "Device stack error",
  715 => "Device illegal instruction",
  716 => "Device load or store instruction on unaligned memory address",
  717 => "Device instruction on address not belonging to an allowed space",
  718 => "Device program counter wrapped its address space",
  719 => "Device exception while executing a kernel",
  800 => "Operation not permitted",
  801 => "Operation not supported",
  999 => "Unknown error"
))

immutable CuDriverError <: Exception
  code::Int
end

import Base.show
show(io::IO, error::CuDriverError) = print(io, driver_error_descriptions[error.code])

macro cucall(fv, argtypes, args...)
  f = eval(fv)
  args = map(esc, args)
  quote
    _curet = ccall( ($(Meta.quot(f)), $libcuda), Cint, $argtypes, $(args...) )
    if _curet != 0
      throw(CuDriverError(round(Int, _curet)))
    end
  end
end

function init()
  @cucall(:cuInit, (Cint,), 0)
end

############################################################
# Box a julia variable as a pointer
############################################################
cubox{T}(x::T) = T[x]

############################################################
# Device and Context
############################################################
immutable CuDevice
  ordinal::Cint
  handle::Cint

  function CuDevice(i::Int)
    ordinal = convert(Cint, i)
    a = Cint[0]
    @cucall(:cuDeviceGet, (Ptr{Cint}, Cint), a, ordinal)
    handle = a[1]
    new(ordinal, handle)
  end
end

immutable CuContext
  handle::Ptr{Void}
end

const CTX_SCHED_AUTO  = 0x00
const CTX_SCHED_SPIN  = 0x01
const CTX_SCHED_YIELD = 0x02
const CTX_SCHED_BLOCKING_SYNC = 0x04
const CTX_MAP_HOST = 0x08
const CTX_LMEM_RESIZE_TO_MAX = 0x10

function create_context(dev::CuDevice, flags::Integer)
  a = Array(Ptr{Void}, 1)
  @cucall(:cuCtxCreate_v2, (Ptr{Ptr{Void}}, Cuint, Cint), a, flags, dev.handle)
  return CuContext(a[1])
end

create_context(dev::CuDevice) = create_context(dev, 0)

function destroy(ctx::CuContext)
  @cucall(:cuCtxDestroy_v2, (Ptr{Void},), ctx.handle)
end

############################################################
# Memory allocation
############################################################
typealias CUdeviceptr Ptr{Void}

type CuPtr
  p::CUdeviceptr

  CuPtr() = new(convert(CUdeviceptr, 0))
  CuPtr(p::CUdeviceptr) = new(p)
end

cubox(p::CuPtr) = cubox(p.p)

function cualloc(T::Type, len::Integer)
  a = CUdeviceptr[0]
  nbytes = round(Int, len) * sizeof(T)
  @cucall(:cuMemAlloc_v2, (Ptr{CUdeviceptr}, Csize_t), a, nbytes)
  return CuPtr(a[1])
end

function free(p::CuPtr)
  @cucall(:cuMemFree_v2, (CUdeviceptr,), p.p)
end

############################################################
# CUDA streams
############################################################
immutable CuStream
  handle::Ptr{Void}
  blocking::Bool
  priority::Int
end

function null_stream()
  CuStream(convert(Ptr{Void}, 0), true, 0)
end

function destroy(s::CuStream)
  @cucall(:cuStreamDestroy_v2, (Ptr{Void},), s.handle)
end

function synchronize(s::CuStream)
  @cucall(:cuStreamSynchronize, (Ptr{Void},), s.handle)
end

############################################################
# PTX Module and Function
############################################################
immutable CuModule
  handle::Ptr{Void}

  function CuModule(filename::AbstractString)
    a = Array(Ptr{Void}, 1)
    @cucall(:cuModuleLoad, (Ptr{Ptr{Void}}, Ptr{Cchar}), a, filename)
    new(a[1])
  end
end

function unload(md::CuModule)
  @cucall(:cuModuleUnload, (Ptr{Void},), md.handle)
end

immutable CuFunction
  handle::Ptr{Void}

  function CuFunction(md::CuModule, name::ASCIIString)
    a = Array(Ptr{Void}, 1)
    @cucall(:cuModuleGetFunction, (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Cchar}),
    a, md.handle, name)
    new(a[1])
  end
end

############################################################
# Launching kernels
############################################################
# This value should be manually synced with the one in kernels/kernels.cu
const THREADS_PER_BLOCK_X = 128
const THREADS_PER_BLOCK_Y = 1
const THREADS_PER_BLOCK_Z = 8

using Compat
get_dim_x(g::Int) = g
get_dim_x(g::@compat(Tuple{Int, Int})) = g[1]
get_dim_x(g::@compat(Tuple{Int, Int, Int})) = g[1]

get_dim_y(g::Int) = 1
get_dim_y(g::@compat(Tuple{Int, Int})) = g[2]
get_dim_y(g::@compat(Tuple{Int, Int, Int})) = g[2]

get_dim_z(g::Int) = 1
get_dim_z(g::@compat(Tuple{Int, Int})) = 1
get_dim_z(g::@compat(Tuple{Int, Int, Int})) = g[3]

using Compat
@compat typealias CuDim Union{Int, Tuple{Int, Int}, Tuple{Int, Int, Int}}

# Stream management

function launch(f::CuFunction, grid::CuDim, block::CuDim, args::Tuple; shmem_bytes::Int=4, stream::CuStream=null_stream())
  gx = get_dim_x(grid)
  gy = get_dim_y(grid)
  gz = get_dim_z(grid)

  tx = get_dim_x(block)
  ty = get_dim_y(block)
  tz = get_dim_z(block)

  kernel_args = [cubox(arg) for arg in args]

  @cucall(:cuLaunchKernel, (
      Ptr{Void},       # function
      Cuint,           # grid dim x
      Cuint,           # grid dim y
      Cuint,           # grid dim z
      Cuint,           # block dim x
      Cuint,           # block dim y
      Cuint,           # block dim z
      Cuint,           # shared memory bytes,
      Ptr{Void},       # stream
      Ptr{Ptr{Void}},  # kernel parameters,
      Ptr{Ptr{Void}}), # extra parameters
      f.handle, gx, gy, gz, tx, ty, tz, shmem_bytes, stream.handle, kernel_args, Ptr{Ptr{Void}}(0))
end

end # module CUDA


text:
using Probably
using Test
using Probably: probabilities

@testset "Probably" begin

a = QBool()
@test probabilities(a.bit.state, 1) ≈ [1,0]
hadamard(a)
@test probabilities(a.bit.state, 1) ≈ [0.5,0.5]

b = QBool(false)
cnot(a, b)
@test probabilities(a.bit.state, 1, 2) ≈ [0.5,0,0,0.5]
cnot(a, b)
@test probabilities(a.bit.state, 1, 2) ≈ [0.5,0,0.5,0]

a = QBool(true)
@test probabilities(a.bit.state, 1) ≈ [0,1]
@test measure(a) == true

a = QBool(true)
b = QBool(false)
@test measure.((a,b)) == (true,false)
swap(a, b)
@test measure.((a,b)) == (false,true)

end


text:
export PigEnv

const PIG_TARGET_SCORE = 100
const PIG_N_SIDES = 6

mutable struct PigEnv{N} <: AbstractEnv
    scores::Vector{Int}
    current_player::Int
    is_chance_player_active::Bool
    tmp_score::Int
end

"""
    PigEnv(;n_players=2)

See [wiki](https://en.wikipedia.org/wiki/Pig_(dice_game)) for explanation of this game.

Here we use it to demonstrate how to write a game with more than 2 players.
"""
PigEnv(; n_players = 2) = PigEnv{n_players}(zeros(Int, n_players), 1, false, 0)

function RLBase.reset!(env::PigEnv)
    fill!(env.scores, 0)
    env.current_player = 1
    env.is_chance_player_active = false
    env.tmp_score = 0
end

RLBase.current_player(env::PigEnv) =
    env.is_chance_player_active ? CHANCE_PLAYER : env.current_player
RLBase.players(env::PigEnv) = 1:length(env.scores)
RLBase.action_space(env::PigEnv, ::Int) = Space((:roll, :hold))
RLBase.action_space(env::PigEnv, ::ChancePlayer) = Space(OneTo(PIG_N_SIDES))

RLBase.prob(env::PigEnv, ::ChancePlayer) = fill(1 / 6, 6)  # TODO: uniform distribution, more memory efficient

RLBase.state(env::PigEnv, ::Observation{Vector{Int}}, p) = env.scores
RLBase.state_space(env::PigEnv, ::Observation, p) = Space(
    SVector(ntuple(_ -> 0 .. (PIG_TARGET_SCORE + PIG_N_SIDES - 1), length(env.scores))),
)

RLBase.is_terminated(env::PigEnv) = any(s >= PIG_TARGET_SCORE for s in env.scores)

function RLBase.reward(env::PigEnv, player)
    winner = findfirst(>=(PIG_TARGET_SCORE), env.scores)
    if isnothing(winner)
        0
    elseif winner == player
        1
    else
        -1
    end
end

function (env::PigEnv)(action, player::Int)
    if action == :roll
        env.is_chance_player_active = true
    else
        env.scores[player] += env.tmp_score
        env.tmp_score = 0
        env.current_player += 1
        if env.current_player > length(players(env))
            env.current_player = 1
        end
    end
end

function (env::PigEnv)(action, ::ChancePlayer)
    env.is_chance_player_active = false
    if action == 1
        env.tmp_score = 0
        env.current_player += 1
        if env.current_player > length(players(env))
            env.current_player = 1
        end
    else
        env.tmp_score += action
    end
end

RLBase.NumAgentStyle(::PigEnv{N}) where {N} = MultiAgent(N)
RLBase.DynamicStyle(::PigEnv) = SEQUENTIAL
RLBase.ActionStyle(::PigEnv) = MINIMAL_ACTION_SET
RLBase.InformationStyle(::PigEnv) = PERFECT_INFORMATION
RLBase.StateStyle(::PigEnv) = Observation{Vector{Int}}()
RLBase.RewardStyle(::PigEnv) = TERMINAL_REWARD
RLBase.UtilityStyle(::PigEnv) = CONSTANT_SUM
RLBase.ChanceStyle(::PigEnv) = EXPLICIT_STOCHASTIC


text:
@params mutable struct Displacement{T} <: AbstractFunction{T}
    u::AbstractVector{T} # displacement vector
    dudx_tmp::AbstractVector # directional derivative
    solver::AbstractDisplacementSolver
    global_dofs::AbstractVector{<:Integer}
    fevals::Int
    maxfevals::Int
end

function Base.show(::IO, ::MIME{Symbol("text/plain")}, ::Displacement)
    return println("TopOpt displacement function")
end

"""
    Displacement()

Construct the Displacement function struct.
"""
function Displacement(solver::AbstractFEASolver; maxfevals=10^8)
    T = eltype(solver.u)
    dh = solver.problem.ch.dh
    k = ndofs_per_cell(dh)
    global_dofs = zeros(Int, k)
    total_ndof = ndofs(dh)
    u = zeros(T, total_ndof)
    dudx_tmp = zeros(T, length(solver.vars))
    return Displacement(u, dudx_tmp, solver, global_dofs, 0, maxfevals)
end

"""
# Arguments
`x` = design variables

# Returns
displacement vector `u`
"""
function (dp::Displacement{T})(x) where {T}
    @unpack solver, global_dofs = dp
    @unpack penalty, problem, xmin = solver
    dp.fevals += 1
    @assert length(global_dofs) == ndofs_per_cell(solver.problem.ch.dh)
    solver.vars .= x
    solver()
    return copy(solver.u)
end

"""
rrule for autodiff.

du/dxe = -K^-1 * dK/dxe * u
d(u)/d(x_e) = - K^-1 * d(K)/d(x_e) * u
            = - K^-1 * (Σ_ei d(ρ_ei)/d(x_e) * K_ei) * u
            = - K^-1 * [d(ρ_e)/d(x_e) * K_e * u]
d(u)/d(x_e)' * Δ = -d(ρ_e)/d(x_e) * u' * K_e * (K^-1 * Δ)

where d(u)/d(x) ∈ (nDof x nCell); d(u)/d(x)^T * Δ = (nCell x nDof) * (nDof x 1) -> nCell x 1
"""
function ChainRulesCore.rrule(dp::Displacement, x)
    @unpack dudx_tmp, solver, global_dofs = dp
    @unpack penalty, problem, u, xmin = solver
    dh = getdh(problem)
    @unpack Kes = solver.elementinfo
    # Forward-pass
    # Cholesky factorisation
    u = dp(x)
    return u, Δ -> begin # v
        solver.rhs .= Δ
        solver(; reuse_chol=true, assemble_f=false)
        dudx_tmp .= 0
        for e in 1:length(x)
            _, dρe = get_ρ_dρ(x[e], penalty, xmin)
            celldofs!(global_dofs, dh, e)
            Keu = bcmatrix(Kes[e]) * u[global_dofs]
            dudx_tmp[e] = -dρe * dot(Keu, solver.lhs[global_dofs])
        end
        return nothing, dudx_tmp # J1' * v, J2' * v
    end
end


text:
using Temporal
using Indicators
using PyPlot

X = tsread("/julia/packages/v0.4/temporal/data/corn.csv")
x = X["2015-03/", :Settle]

plot(x.index, x.values, label="Corn", lw=3, color="blue")
grid(true, ls="-", color="black", alpha=0.25)

# First-order trendlines
maxi = maxima(x)
mini = minima(x)
plot(x[maxi].index, x[maxi].values, label="Resistance", color="red", marker="o")
plot(x[mini].index, x[mini].values, label="Support", color="green", marker="o")

# Second-order trendlines
maxi = [x maxima(x[maxi])][:,2]
mini = [x minima(x[mini])][:,2]
maxi.values[isnan(maxi.values)] = 0.0
mini.values[isnan(mini.values)] = 0.0
maxi = ts(map(bool maxi.values), maxi.index, maxi.fields)
mini = ts(map(bool mini.values), mini.index, mini.fields)
plot(x[maxi].index, x[maxi].values, color="red", marker="o")
plot(x[mini].index, x[mini].values, color="green", marker="o")

# Third-order trendlines
maxi = [x maxima(x[maxi])][:,2]
mini = [x minima(x[mini])][:,2]
maxi.values[isnan(maxi.values)] = 0.0
mini.values[isnan(mini.values)] = 0.0
maxi = ts(map(bool, maxi.values), maxi.index, maxi.fields)
mini = ts(map(bool, mini.values), mini.index, mini.fields)
plot(x[maxi].index, x[maxi].values, color="red", marker="o")
plot(x[mini].index, x[mini].values, color="green", marker="o")

legend(loc="best", frameon=false)

tight_layout()


text:
module BNP

    using Distributions
    using Clustering
    using StatsBase
    using StatsFuns

	import Base.convert

    # source files
    include("distributions.jl")
    include("distfunctions.jl")

		include("show.jl")
    include("utils.jl")
    include("common.jl")

    include("datasets.jl")

    include("dpmm.jl") # Dirichlet Process Mixture Model

    include("vcm_utils.jl")
    include("vcm.jl") # Variable Clustering Model

    include("hdp_utils.jl")
    include("hdp.jl") # Hierarchical Dirichlet Proces

    # analysis functions
    include("psm.jl")
    include("vipointestimate.jl")

    # define models
    abstract ModelType

    type DPM <: ModelType
      H::ConjugatePostDistribution
      α::Float64

      DPM(H::ConjugatePostDistribution; α = 1.0) = new(H, α)
    end

    type VCM <: ModelType
      α::Float64

      VCM(;α = 1.0) = new(α)
    end

    type HDP <: ModelType
      H::ConjugatePostDistribution
      α::Float64
      γ::Float64

      HDP(H::ConjugatePostDistribution; α = 0.1, γ = 5.0) = new(H, α, γ)
    end

    abstract SamplingType
    type Gibbs <: SamplingType
        burnin::Int
        thinout::Int
        maxiter::Int

        Gibbs(; burnin = 0, thinout = 1, maxiter = 100) = new(burnin, thinout, maxiter)

    end

    abstract InitialisationType

    type PrecomputedInitialisation <: InitialisationType
      Z::Array{Int}
      PrecomputedInitialisation(Z::Array{Int}) = new(Z)
    end

    type RandomInitialisation <: InitialisationType
      k::Int

      RandomInitialisation(;k = 2) = new(k)
    end

    type KMeansInitialisation <: InitialisationType
      k::Int

      KMeansInitialisation(;k = 2) = new(k)
    end

    type IncrementalInitialisation <: InitialisationType end

    # define helper functions
    function train(model::DPM, sampler::Gibbs, init::RandomInitialisation, X::AbstractArray)

      # init
      (Z, G) = init_random_dpmm(X, model.H, k = init.k)

      # inference
      return train_cgibbs_dpmm(X, model.H, Z, G, DPMHyperparam(), alpha = model.α, burnin = sampler.burnin, thinout = sampler.thinout, maxiter = sampler.maxiter)

    end

    function train(model::DPM, sampler::Gibbs, init::PrecomputedInitialisation, X::AbstractArray)

      # init
      (Z, G) = init_precomputed_dpmm(X, model.H, init.Z)

      # inference
      return train_cgibbs_dpmm(X, model.H, Z, G, DPMHyperparam(), alpha = model.α, burnin = sampler.burnin, thinout = sampler.thinout, maxiter = sampler.maxiter)

    end

    function train(model::DPM, sampler::Gibbs, init::KMeansInitialisation, X::AbstractArray)

      # init
      (Z, G) = init_kmeans_dpmm(X, model.H, k = init.k)

      # inference
      return train_cgibbs_dpmm(X, model.H, Z, G, DPMHyperparam(), alpha = model.α, burnin = sampler.burnin, thinout = sampler.thinout, maxiter = sampler.maxiter)

    end

    function train(model::VCM, sampler::Gibbs, init::IncrementalInitialisation, X::AbstractArray)

      # init
      B = init_incremental_vcm(X, VCMHyperparam(); α = model.α)

      # inference
      return train_gibbs_vcm(B, VCMHyperparam(), α = model.α, burnin = sampler.burnin, thinout = sampler.thinout, maxiter = sampler.maxiter)

    end

    function train{T <: Real}(model::HDP, sampler::Gibbs, init::RandomInitialisation, X::Vector{Vector{T}})

      # init
      (Z, G) = init_random_hdp(X, model.H, k = init.k)

      # inference

      return train_gibbs_hdp(X, model.H, Z, G, HDPHyperparam(), init.k, α = model.α, γ = model.γ, maxiter = sampler.maxiter)

    end

    # export commands
    export
      # Models
      DPM,
      VCM,
      HDP,

      # Inference Algorithms
      Gibbs,

      # Initialisations
      PrecomputedInitialisation,
      IncrementalInitialisation,
      RandomInitialisation,
      KMeansInitialisation,

      # Parameters
      DPMHyperparam,
      VCMHyperparam,
      HDPHyperparam,

      # Resulting Models
      DPMData,
      VCMData,
      HDPData,

      # training method
      train,

      # expose Distribution functions
      add_data!,
      remove_data!,
      logpred,

      # expose utility functions
      rand_indices,

      # dataset generation
      generateBarsDataset,

	  # analysis
	  compute_psm,
      point_estimate

end


text:
using Pkg
Pkg.activate(@__DIR__)
cd(@__DIR__)
Pkg.update()

push!(LOAD_PATH, "../src/")
const CI = get(ENV, "CI", nothing) == "true"

using Documenter, LimnoSES

makedocs(
    modules = [LimnoSES],
    sitename = "LimnoSES.jl",
    authors = "Tim DuBois",
    format = Documenter.HTML(prettyurls = CI),
    pages = [
        "Introduction" => "index.md",
        "Public API" => "api.md",
        "Developer Docs" => "dev.md",
    ],
)

if CI
    deploydocs(
        repo = "github.com/Libbum/LimnoSES.jl.git",
        target = "build",
        push_preview = true,
    )
end


text:
using GameOfLife
using Test

@testset "GameOfLife.jl" begin
    include("rules.jl")
    include("world.jl")
end

max_stars_repo_name: VasanthManiVasi/Conways-Game-of-Life

max_issues_repo_name: VasanthManiVasi/Conways-Game-of-Life

max_forks_repo_name: VasanthManiVasi/Conways-Game-of-Life


text:
using Plots
x = 1:10
ys=reshape(1:10*5,10,5)
p=plot(x,ys)
plot!(x,sin.(collect(x)))
plot!(p,x,cos.(collect(x)))

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
import VsRecorderBase: vs_tryparse_scene, vs_result
import VsRecorderBase.DefaultStrategyModule: feature_image_and_masks
using VsRecorder.Scenes: feature_image_name

function VsRecorderBase.vs_init!(ctx::PokemonContext{DefaultStrategy})
    invoke(vs_init!, Tuple{VsContext{DefaultStrategy}}, ctx)
    tessdata_path = datapath("tessdata")
    ocr_language = ctx.config.ocr_language
    if isfile(tessdata_path, "$ocr_language.traineddata")
        # Load specialized OCR trained data.
        Tesseract.tess_init(
            ctx.ocr_instances[ocr_language],
            ocr_language,
            tessdata_path
        )
    end
    init_multiple_ocr!(ctx, ["chi_sim+chi_tra", "jpn", "kor", "fra+deu+spa+ita+eng"])
    data = ctx.data
    data.context = ParsingContext()
    Data.initialize_icons()
    ctx
end

function feature_image_and_masks(source::PokemonBattle, ctx::VsContext)
    language = lowercase(string(source.language))

    scene_types = Type[]
    # TODO
    append!(scene_types, Scenes.AvailableScenes)

    map(scene_types) do scene_type
        name = feature_image_name(scene_type, source)
        filename = "$name.png"
        img = load_data("scenes", language, filename)
        mask = try
            load_data("scenes", "mask", filename)
        catch
            ones(Float32, size(img))
        end
        (scene_type, img, mask)
    end
end

function vs_result(ctx::PokemonContext)
    data = ctx.data
    battles = data.battles
    parsed_battles = data.parsed_battles
    PokemonBattleResult(
        battles = battles,
        parsed_battles = parsed_battles
    )
end


text:
##############################################################################
##
## References
##
##############################################################################

# - https://github.com/bitcoin/bitcoin/blob/master/src/base58.cpp

##############################################################################
##
## Constants
##
##############################################################################

# TODO: This would be nice as an associative array so we can just look up by
#       base.
const base58 = "123456789ABCDEFGHJKLMNPQRSTUVWXYZabcdefghijkmnopqrstuvwxyz"

##############################################################################
##
## Function definitions
##
##############################################################################

# We really only need these for base 58 because Julia's default base 58 math
# uses a different alphabet.
function encode58(n::Array{Uint8})
  result = encode58(Crypto.oct2int(n))

  # TODO: this isn't correct in the case of general # of 0 bytes
  for byte in n
    if byte == 0x00
      result = string(base58[1], result) # Add the zero element
    end
  end

  return result
end

# - n: integer we want to convert
function encode58(n::Integer)
  # Require base to be 58 for now
  b = 58

  output = ""

  while n > 0
    n, rem = divrem(n, b)
    output = string(base58[rem + 1], output)
  end

  return output
end

# Decode from base 58 to integer
# - n: base 58 number we want to convert
function decode58(n::String)
  # Require base to be 58 for now
  b = 58

  result = BigInt(0)

  for i = 1:length(n)
    result = result * b + search(base58, n[i]) - 1
  end

  return result
end

function decode58_to_array(n::String)
  result = Crypto.hex2oct(hex(decode58(n)))

  # TODO: this isn't correct in the case of general # of 0s
  # Capture 2 leading 0s
  if n[1] == base58[1] # Get the zero element
    result = [0x00, result]
  end

  return result
end


text:
# # Example 1: Optimization of a State-to-State Transfer in a Two-Level-System

#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`simple_state_to_state.ipynb`](@__NBVIEWER_ROOT_URL__/examples/simple_state_to_state.ipynb).
#md #
#md #     Compare this example against the [same example using Krotov's
#md #     method](https://juliaquantumcontrol.github.io/Krotov.jl/stable/examples/simple_state_to_state/).

#md # ``\gdef\op#1{\hat{#1}}``
#md # ``\gdef\init{\text{init}}``
#md # ``\gdef\tgt{\text{tgt}}``

#nb # $
#nb # \newcommand{tr}[0]{\operatorname{tr}}
#nb # \newcommand{diag}[0]{\operatorname{diag}}
#nb # \newcommand{abs}[0]{\operatorname{abs}}
#nb # \newcommand{pop}[0]{\operatorname{pop}}
#nb # \newcommand{aux}[0]{\text{aux}}
#nb # \newcommand{opt}[0]{\text{opt}}
#nb # \newcommand{tgt}[0]{\text{tgt}}
#nb # \newcommand{init}[0]{\text{init}}
#nb # \newcommand{lab}[0]{\text{lab}}
#nb # \newcommand{rwa}[0]{\text{rwa}}
#nb # \newcommand{bra}[1]{\langle#1\vert}
#nb # \newcommand{ket}[1]{\vert#1\rangle}
#nb # \newcommand{Bra}[1]{\left\langle#1\right\vert}
#nb # \newcommand{Ket}[1]{\left\vert#1\right\rangle}
#nb # \newcommand{Braket}[2]{\left\langle #1\vphantom{#2}\mid{#2}\vphantom{#1}\right\rangle}
#nb # \newcommand{op}[1]{\hat{#1}}
#nb # \newcommand{Op}[1]{\hat{#1}}
#nb # \newcommand{dd}[0]{\,\text{d}}
#nb # \newcommand{Liouville}[0]{\mathcal{L}}
#nb # \newcommand{DynMap}[0]{\mathcal{E}}
#nb # \newcommand{identity}[0]{\mathbf{1}}
#nb # \newcommand{Norm}[1]{\lVert#1\rVert}
#nb # \newcommand{Abs}[1]{\left\vert#1\right\vert}
#nb # \newcommand{avg}[1]{\langle#1\rangle}
#nb # \newcommand{Avg}[1]{\left\langle#1\right\rangle}
#nb # \newcommand{AbsSq}[1]{\left\vert#1\right\vert^2}
#nb # \newcommand{Re}[0]{\operatorname{Re}}
#nb # \newcommand{Im}[0]{\operatorname{Im}}
#nb # $

# This first example illustrates the basic use of the `GRAPE.jl` by solving a simple canonical optimization problem: the transfer of population in a two level system.

using DrWatson
@quickactivate "GRAPETests"
#-
using QuantumControl

#jl using Test
#jl println("")

# ## Two-level Hamiltonian

# We consider the Hamiltonian $\op{H}_{0} = - \frac{\omega}{2} \op{\sigma}_{z}$, representing a simple qubit with energy level splitting $\omega$ in the basis $\{\ket{0},\ket{1}\}$. The control field $\epsilon(t)$ is assumed to couple via the Hamiltonian $\op{H}_{1}(t) = \epsilon(t) \op{\sigma}_{x}$ to the qubit, i.e., the control field effectively drives transitions between both qubit states.
#
# We we will use

ϵ(t) = 0.2 * QuantumControl.Shapes.flattop(t, T=5, t_rise=0.3, func=:blackman);

#-
"""Two-level-system Hamiltonian."""
function hamiltonian(Ω=1.0, ϵ=ϵ)
    σ̂_z = ComplexF64[
        1  0
        0 -1
    ]
    σ̂_x = ComplexF64[
        0  1
        1  0
    ]
    Ĥ₀ = -0.5 * Ω * σ̂_z
    Ĥ₁ = σ̂_x
    return (Ĥ₀, (Ĥ₁, ϵ))
end;
#-

H = hamiltonian();
#jl @test length(H) == 2

# The control field here switches on from zero at $t=0$ to it's maximum amplitude
# 0.2 within the time period 0.3 (the switch-on shape is half a [Blackman pulse](https://en.wikipedia.org/wiki/Window_function#Blackman_window)).
# It switches off again in the time period 0.3 before the
# final time $T=5$). We use a time grid with 500 time steps between 0 and $T$:

tlist = collect(range(0, 5, length=500));

#-
using Plots
Plots.default(
    linewidth               = 3,
    size                    = (550, 300),
    legend                  = :right,
    foreground_color_legend = nothing,
    background_color_legend = RGBA(1, 1, 1, 0.8)
)
#-
function plot_control(pulse::Vector, tlist)
    plot(tlist, pulse, xlabel="time", ylabel="amplitude", legend=false)
end

plot_control(ϵ::T, tlist) where {T<:Function} = plot_control([ϵ(t) for t in tlist], tlist);
#-
fig = plot_control(H[2][2], tlist)
#jl display(fig)

# ## Optimization target

# First, we define a convenience function for the eigenstates.

function ket(label)
    result = Dict("0" => Vector{ComplexF64}([1, 0]), "1" => Vector{ComplexF64}([0, 1]))
    return result[string(label)]
end;

#-
#jl using LinearAlgebra
#jl @test dot(ket(0), ket(1)) ≈ 0
#-

# The physical objective of our optimization is to transform the initial state $\ket{0}$ into the target state $\ket{1}$ under the time evolution induced by the Hamiltonian $\op{H}(t)$.

objectives = [Objective(initial_state=ket(0), generator=H, target_state=ket(1))];

#-
#jl @test length(objectives) == 1
#-

# The full control problem includes this objective, information about the time grid for the dynamics, and the functional to be used (the square modulus of the overlap $\tau$ with the target state in this case).

problem = ControlProblem(
    objectives=objectives,
    tlist=tlist,
    pulse_options=Dict(),
    iter_stop=500,
    J_T=QuantumControl.Functionals.J_T_sm,
    check_convergence=res -> begin
        ((res.J_T < 1e-3) && (res.converged = true) && (res.message = "J_T < 10⁻³"))
    end,
);

# ## Simulate dynamics under the guess field

# Before running the optimization procedure, we first simulate the dynamics under the guess field $\epsilon_{0}(t)$. The following solves equation of motion for the defined objective, which contains the initial state $\ket{\Psi_{\init}}$ and the Hamiltonian $\op{H}(t)$ defining its evolution.

guess_dynamics = propagate_objective(
    objectives[1],
    problem.tlist;
    storage=true,
    observables=(Ψ -> abs.(Ψ) .^ 2,)
)

#-
function plot_population(pop0::Vector, pop1::Vector, tlist)
    fig = plot(tlist, pop0, label="0", xlabel="time", ylabel="population")
    plot!(fig, tlist, pop1; label="1")
end;
#-
fig = plot_population(guess_dynamics[1, :], guess_dynamics[2, :], tlist)
#jl display(fig)

# ## Optimization with LBFGSB

# In the following we optimize the guess field $\epsilon_{0}(t)$ such that the intended state-to-state transfer $\ket{\Psi_{\init}} \rightarrow \ket{\Psi_{\tgt}}$ is solved.

# The GRAPE package performs the optimization by calculating the gradient of $J_T$ with respect to the values of the control field at each point in time. This gradient is then fed into a backend solver that calculates an appropriate update based on that gradient.

# By default, this backend is [LBFGSB.jl](https://github.com/Gnimuc/LBFGSB.jl), a wrapper around the true and tested [L-BFGS-B Fortran library](http://users.iems.northwestern.edu/%7Enocedal/lbfgsb.html). L-BFGS-B is a pseudo-Hessian method: it efficiently estimates the second-order Hessian from the gradient information. The search direction determined from that Hessian dramatically improves convergence compared to using the gradient directly as a search direction. The L-BFGS-B method performs its own linesearch to determine how far to go in the search direction.

# It can be quite instructive to see how the improvement in the pseudo-Hessian search direction compares to the gradient, how the linesearch finds an appropriate step width. For this purpose, we have a [GRAPELinesearchAnalysis](https://github.com/JuliaQuantumControl/GRAPELinesearchAnalysis.jl) package that automatically generates plots in every iteration of the optimization showing the linesearch behavior

using GRAPELinesearchAnalysis

# We feed this into the optimization as part of the `info_hook`.

opt_result_LBFGSB, file = @optimize_or_load(
    datadir("TLS"),
    problem,
    method = :grape,
    force = true,
    filename = "opt_result_LBFGSB.jld2",
    info_hook = chain_infohooks(
        GRAPELinesearchAnalysis.plot_linesearch(datadir("TLS", "Linesearch", "LBFGSB")),
        QuantumControl.GRAPE.print_table,
    )
);
#-
#jl @test opt_result_LBFGSB.J_T < 1e-3
#-

# When going through this tutorial locally, the [generated images for the linesearch](https://github.com/JuliaQuantumControl/GRAPE.jl/tree/data-dump/TLS/Linesearch/LBFGSB) can be found in `docs/TLS/Linesearch/LBFGSB`.

datadir("TLS", "Linesearch", "LBFGSB")
#-
opt_result_LBFGSB

# We can plot the optimized field:

#-
fig = plot_control(opt_result_LBFGSB.optimized_controls[1], tlist)
#jl display(fig)
#-

# ## Optimization via χ

opt_result_LBFGSB_via_χ = optimize(problem; method=:grape, gradient_via=:chi);
#-
opt_result_LBFGSB_via_χ

# ## Optimization with Optim.jl

# As an alternative to the default L-BFGS-B backend, we can also use any of the gradient-based optimizers in [Optiml.jl](https://github.com/JuliaNLSolvers/Optim.jl). This also gives full control over the linesearch method.

using Optim
using LineSearches

# Here, we use the LBFGS implementation that is part of Optim (which is not exactly the same as L-BFGS-B; "B" being the variant of LBFGS with optional additional bounds on the control) with a Hager-Zhang linesearch

opt_result_OptimLBFGS, file = @optimize_or_load(
    datadir("TLS"),
    problem,
    method = :grape,
    filename = "opt_result_OptimLBFGS.jld2",
    info_hook = chain_infohooks(
        GRAPELinesearchAnalysis.plot_linesearch(datadir("TLS", "Linesearch", "OptimLBFGS")),
        QuantumControl.GRAPE.print_table,
    ),
    optimizer = Optim.LBFGS(;
        alphaguess=LineSearches.InitialStatic(alpha=0.2),
        linesearch=LineSearches.HagerZhang(alphamax=2.0)
    )
);

#-
#jl # do a single iteration, to ensure coverage
#jl _optimizer = Optim.LBFGS(;
#jl     alphaguess=LineSearches.InitialStatic(alpha=0.2),
#jl     linesearch=LineSearches.HagerZhang(alphamax=2.0)
#jl )
#jl optimize(problem, method = :grape, iter_stop=1, optimizer=_optimizer)
#-
#jl @test opt_result_OptimLBFGS.J_T < 1e-3
#-

opt_result_OptimLBFGS
#-

# We can plot the optimized field:

fig = plot_control(opt_result_OptimLBFGS.optimized_controls[1], tlist)

# We can see that the choice of linesearch parameters in particular strongly influence the convergence and the resulting field. Play around with different methods and parameters, and compare the different [plots generated by `GRAPELinesearchAnalysis`](https://github.com/JuliaQuantumControl/GRAPE.jl/tree/data-dump/TLS/Linesearch/OptimLBFGS)!
#
# Empirically, we find the default L-BFGS-B to have a very well-behaved linesearch.

# ## Simulate the dynamics under the optimized field

# Having obtained the optimized control field, we can simulate the dynamics to verify that the optimized field indeed drives the initial state $\ket{\Psi_{\init}} = \ket{0}$ to the desired target state $\ket{\Psi_{\tgt}} = \ket{1}$.

opt_dynamics = propagate_objective(
    objectives[1],
    problem.tlist;
    controls_map=IdDict(ϵ => opt_result_LBFGSB.optimized_controls[1]),
    storage=true,
    observables=(Ψ -> abs.(Ψ) .^ 2,)
)

#-
fig = plot_population(opt_dynamics[1, :], opt_dynamics[2, :], tlist)
#jl display(fig)
#-

#-
#jl @test opt_dynamics[2,end] > 0.99
#-


text:
# Copyright(c) Rafael Dias
# MIT Licensed
include("propriedades-secoes.jl")

module AnaliseEstruturas
    include("propriedades-secoes.jl")
    mutable struct Elemento
        coordenadaInicial::PropriedadesSecoes.Coordenada
        coordenadaFinal::PropriedadesSecoes.Coordenada
    end

    mutable struct Estrutura
        elementos::Array{Elemento}
    end
end

max_stars_repo_name: rsd-consultoria/scripts-engenharia

max_issues_repo_name: rsd-consultoria/scripts-engenharia

max_forks_repo_name: rsd-consultoria/scripts-engenharia


text:
# Autogenerated wrapper script for qr_mumps_jll for aarch64-linux-musl-libgfortran5
export libcqrm, libdqrm, libqrm_common, libsqrm, libzqrm

using OpenBLAS32_jll
using CompilerSupportLibraries_jll
using SuiteSparse_jll
using METIS_jll
using SCOTCH_jll
JLLWrappers.@generate_wrapper_header("qr_mumps")
JLLWrappers.@declare_library_product(libcqrm, "libcqrm.so")
JLLWrappers.@declare_library_product(libdqrm, "libdqrm.so")
JLLWrappers.@declare_library_product(libqrm_common, "libqrm_common.so")
JLLWrappers.@declare_library_product(libsqrm, "libsqrm.so")
JLLWrappers.@declare_library_product(libzqrm, "libzqrm.so")
function __init__()
    JLLWrappers.@generate_init_header(OpenBLAS32_jll, CompilerSupportLibraries_jll, SuiteSparse_jll, METIS_jll, SCOTCH_jll)
    JLLWrappers.@init_library_product(
        libcqrm,
        "lib/libcqrm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libdqrm,
        "lib/libdqrm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libqrm_common,
        "lib/libqrm_common.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libsqrm,
        "lib/libsqrm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libzqrm,
        "lib/libzqrm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function dstats(x, rnames="";short=false)
    k = size(x,2)
    if rnames==""
        rnames = 1:k
        rnames = rnames'
    end
    m = mean(x,dims=1)
    mm = median(x,dims=1)
    s = std(x,dims=1)
    sk = m-m
    kt = m-m
    mn = minimum(x,dims=1)
    mx = maximum(x,dims=1)
    q05 = fill(0.0,k)
    q95 = fill(0.0,k)
    if short == false
        for i = 1:size(x,2) q05[i] = quantile(x[:,i], 0.05) end
        for i = 1:size(x,2) q95[i] = quantile(x[:,i], 0.95) end
        cnames = ["  mean", " median","  std", "min", "max", "q05", "q95"]
        stats = [m' mm' s' mn' mx' q05 q95]
        prettyprint(stats, cnames, rnames)
    else
        cnames = ["  mean", " median", "  std", "min", "max"]
        stats = [m' mm' s' mn' mx']
        prettyprint(stats, cnames, rnames)
    end
    return stats
end

max_stars_repo_name: manuelmontesinos/Econometrics

max_issues_repo_name: manuelmontesinos/Econometrics

max_forks_repo_name: manuelmontesinos/Econometrics


text:
# AUTO GENERATED FILE - DO NOT EDIT

export dcc_link

"""
    dcc_link(;kwargs...)
    dcc_link(children::Any;kwargs...)
    dcc_link(children_maker::Function;kwargs...)

A Link component.
Link allows you to create a clickable link within a multi-page app.

For links with destinations outside the current app, `html.A` is a better
component to use.
Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): The children of this component
- `id` (String; optional): The ID of this component, used to identify dash components
in callbacks. The ID needs to be unique across all of the
components in an app.
- `href` (String; required): The URL of a linked resource.
- `refresh` (Bool; optional): Controls whether or not the page will refresh when the link is clicked
- `className` (String; optional): Often used with CSS to style elements with common properties.
- `style` (Dict; optional): Defines CSS styles which will override styles previously set.
- `title` (String; optional): Adds the title attribute to your link, which can contain supplementary
information.
- `target` (String; optional): Specifies where to open the link reference.
- `loading_state` (optional): Object that holds the loading state object coming from dash-renderer. loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
"""
function dcc_link(; kwargs...)
        available_props = Symbol[:children, :id, :href, :refresh, :className, :style, :title, :target, :loading_state]
        wild_props = Symbol[]
        return Component("dcc_link", "Link", "dash_core_components", available_props, wild_props; kwargs...)
end

dcc_link(children::Any; kwargs...) = dcc_link(;kwargs..., children = children)
dcc_link(children_maker::Function; kwargs...) = dcc_link(children_maker(); kwargs...)

max_stars_repo_name: asc-csa/dash-core-components

max_issues_repo_name: asc-csa/dash-core-components

max_forks_repo_name: asc-csa/dash-core-components


text:
const base_url = "http://vrp.atd-lab.inf.puc-rio.br/media/com_vrp/instances/"

const dir_names = [
    "A", "B", "E", "F", "M", "P", "CMT", "tai", "Golden", "Li", "X", "XXL"
]
const XXL_names = [
    "Antwerp", "Brussels", "Flanders", "Ghent", "Leuven"
]

function download_cvrp(name::String)
    is_xxl = any([startswith(name, xxl) for xxl in XXL_names])

    if is_xxl
        dir_name = "XXL"
    else
        for d in dir_names
            if startswith(name, d)
                dir_name = d
                break
            end
        end
    end

    data_dir = joinpath(@__DIR__, "..", "data")
    vrp_url = base_url * dir_name * "/" * name
    vrp_file = joinpath(data_dir, name)
    exts = [".vrp", ".sol"]

    if !isdir(data_dir)
        mkdir(data_dir)
    end

    for ext in exts
        file = vrp_file * ext
        if !isfile(file)
            try
                output = Downloads.download(vrp_url * ext, file)
            catch e
                error("No such data name exists in CVRPLIB: $name")
                print(e)
            end
        end
    end

    return vrp_file .* exts
end

if abspath(PROGRAM_FILE) == @__FILE__
    download_cvrp("Li_21")
end


text:
"""
Water and water vapor properties.

* [`water_density`](@ref)
* [`water_dissoc`](@ref)
* [`latent_heat_vap`](@ref)
* [`latent_heat_sub`](@ref)
* [`e_sat`](@ref)
* [`e_sat_deriv`](@ref)
* [`e_sat_ice`](@ref)
* [`e_sat_ice_deriv`](@ref)
* [`vapor_deficit`](@ref)
* [`vapor_mole_frac`](@ref)
* [`vapor_deficit_mole_frac`](@ref)

# References

## Water density

* [WP02] Wagner, W. and Pruß, A. (2002). The IAPWS Formulation 1995 for the
  Thermodynamic Properties of Ordinary Water Substance for General and
  Scientific Use. *J. Phys. Chem. Ref. Data*, 31, 387--535.
  <https://doi.org/10.1063/1.1461829>.

## Water dissociation constant

* [BL06] Bandura, A. V. and Lvov, S. N. (2006). The ionization constant of
  water over wide ranges of temperature and density. *J. Phys. Chem. Ref.
  Data*, 35, 15--30. <https://doi.org/10.1063/1.1928231>.

## Latent heat

* [H84] Henderson-Sellers, B. (1984). A new formula for latent heat of
  vaporization of water as a function of temperature. *Quart. J. Roy. Meteorol.
  Soc.*, 110(466), 1186--1190. <https://doi.org/10.1002/qj.49711046626>.
* [MK05] Murphy, D. M., and Koop, T. (2005). Review of the vapour pressures of
  ice and supercooled water for atmospheric applications. *Quart. J. Roy.
  Meteorol. Soc.*, 131(608), 1539--1565. <https://doi.org/10.1256/qj.04.94>.

## Saturation vapor pressure

* [GG46] Goff, J. A., and Gratch, S. (1946). Low-pressure properties of water
  from -160 to 212 F, in *Transactions of the American Society of Heating and
  Ventilating Engineers*, pp 95--122, presented at the 52nd Annual Meeting of
  the American Society of Heating and Ventilating Engineers, New York, NY, USA.
* [V16] Vömel, Holger. (2016). *Saturation vapor pressure formulations*.
  Earth Observing Laboratory, National Center for Atmospheric Research,
  Boulder, CO, USA. Retrieved from <http://cires1.colorado.edu/~voemel/vp.html>
  on 2017-10-09.
"""
module Water

include("docstring_style.jl")

using Canopy.Constants: M_w

export water_density,
    water_dissoc,
    latent_heat_vap,
    latent_heat_sub,
    e_sat,
    e_sat_deriv,
    e_sat_ice,
    e_sat_ice_deriv,
    vapor_deficit,
    vapor_mole_frac,
    vapor_deficit_mole_frac

"""
Calculate water density [kg m^-3] as a function of temperature [K].

# Examples

```jldoctest
julia> water_density(298.15)
996.9993666156083
```
"""
function water_density(temp)
    T_crit = 647.096  # critical temperature
    rho_crit = 322.  # critical density
    theta = 1. - temp / T_crit
    rho_ratio = 1. + 1.99274064 * theta ^ (1. / 3.) +
        1.09965342 * theta ^ (2. / 3.) - 0.510839303 * theta ^ (5. / 3.) -
        1.7549349 * theta ^ (16. / 3.) - 45.5170352 * theta ^ (43. / 3.) -
        6.74694450e5 * theta ^ (110. / 3.)
    return rho_crit * rho_ratio
end

"""
Calculate water dissociation constant (pK_w) as a function of temperature [K].

# Examples

```jldoctest
julia> water_dissoc(298.15)
13.994884354781636
```
"""
function water_dissoc(temp)
    n = 6.
    rho_w = water_density(temp) * 1e-3  # convert to g cm^-3 here
    temp_sq = temp * temp
    Z = rho_w * exp(-0.864671 + 8659.19 / temp -
                    22786.2 / temp_sq * rho_w ^ (2. / 3.))
    pK_w_G = 0.61415 + 48251.33 / temp - 67707.93 / temp_sq +
        10102100. / (temp * temp_sq)
    pK_w = -2. * n * (log10(1. + Z) - Z / (Z + 1.) * rho_w *
                      (0.642044 - 56.8534 / temp - 0.375754 * rho_w)) +
        pK_w_G + 2. * log10(M_w)
    return pK_w
end

"""
Calculate the molar latent heat of vaporization of water [J mol^-1]. Applicable
range: 1 atm, -130 -- 40°C.

# Examples

```jldoctest
julia> latent_heat_vap(273.15)  # 0 C
45053.505423916824

julia> latent_heat_vap(373.15)  # 100 C
41816.390994192036

julia> latent_heat_vap(243.15)  # -30 C, supercooled water
46398.07458086907
```

# Known issues

!!! note

    Temperature at 273.15 K is a discontinuity due to small differences in the
    fitted parameters for normal and supercooled water. However, the difference
    is only 1.4 J mol^-1 and can be neglected in most applications.
"""
latent_heat_vap(temp) =
    if temp < 273.15
        # supercooled water
        L_v_sc = 56579. - 42.212 * temp + exp(0.1149 * (281.6 - temp))
    else
        L_v = 1.91846e6 * (temp / (temp - 33.91))^2 * M_w
    end

"""
Calculate the molar latent heat of sublimation of water ice [J mol^-1].
Applicable range: 1 atm, -240 -- 0°C.

# Examples

```jldoctest
julia> latent_heat_sub(273.15)  # 0 C
51059.029186873086

julia> latent_heat_vap(258.15)  # -15 C
45696.76891208987
```
"""
latent_heat_sub(temp) = 46782.5 + 35.8925 * temp - 0.07414 * temp * temp +
    541.5 * exp(- (temp / 123.75)^2)

"""
Calculate saturation vapor pressure of water [Pa] using Goff--Gratch equation.

# Examples

```jldoctest
julia> e_sat(298.15)  # 25 C
3165.1956333836806

julia> e_sat(273.15)  # 0 C
610.3360999334138
```
"""
function e_sat(temp)
    u = 373.16 / temp
    v = temp / 373.16
    log10_e_sat = -7.90298 * (u - 1.) +
        5.02808 * log10(u) -
        1.3816e-7 * (10^(11.344 * (1. - v)) - 1.) +
        8.1328e-3 * (10^(-3.49149 * (u - 1.)) - 1.) +
        log10(1013.246) + 2.  # add 2 to convert from hPa to Pa
    return 10^log10_e_sat
end

"""
Calculate the temperature derivative of saturation vapor pressure of water [Pa
K^-1].

# Examples

```jldoctest
julia> e_sat_deriv(298.15)  # 25 C
188.6862159855073

julia> e_sat_deriv(273.15)  # 0 C
44.35218469842299
```
"""
function e_sat_deriv(temp)
    log_e_sat_deriv = (6790.4984743899386 +
                       exp(12.068003566856145 -
                           3000.0022166762069 / temp)) / (temp * temp) -
        5.02808 / temp +
        exp(8.5004184700093912 - 0.069998191914793798 * temp)
    return e_sat(temp) * log_e_sat_deriv
end

"""
Calculate saturation vapor pressure of water ice [Pa] using Goff--Gratch
equation.

# Examples

```jldoctest
julia> e_sat_ice(258.15)  # -15 C
165.0147739235893

julia> e_sat_ice(273.15)  # 0 C
610.2072697609788
```
"""
function e_sat_ice(temp)
    u = 273.16 / temp
    v = temp / 273.16
    log10_e_sat = -9.09718 * (u - 1.) -
        3.56654 * log10(u) + 0.876793 * (1. - v) +
        log10(6.1071) + 2.  # add 2 to convert from hPa to Pa
    return 10^log10_e_sat
end

"""
Calculate the temperature derivative of saturation vapor pressure of water ice
[Pa K^-1].

# Examples

```jldoctest
julia> e_sat_ice_deriv(258.15)  # -15 C
15.228515762185872

julia> e_sat_ice_deriv(273.15)  # 0 C
50.254185101185534
```
"""
function e_sat_ice_deriv(temp)
    log_e_sat_deriv = 5721.891003334422 / (temp * temp) +
        3.56654 / temp - 0.007390871618983484
    return e_sat_ice(temp) * log_e_sat_deriv
end

"""
Calculate water vapor pressure deficit [Pa] from temperature [K] and relative
humidity [0--1].

# Examples

```jldoctest
julia> vapor_deficit(298.15, 0.5)
1582.5978166918403

julia> vapor_deficit(288.15, 0.6)
681.3124028440158
```
"""
vapor_deficit(temp, RH) = e_sat(temp) * (1. - RH)

"""
Calculate water vapor mole fraction [mol mol^-1] from temperature [K], pressure
[Pa], and relative humidity [0--1].

# Examples

```jldoctest
julia> vapor_mole_frac(298.15, 101325.0, 0.7)
0.021866636500059967

julia> vapor_mole_frac(288.15, 101325.0, 1.0)
0.016810076556723803
```
"""
vapor_mole_frac(temp, pressure, RH) = e_sat(temp) * RH / pressure

"""
Calculate the mole fraction vapor deficit [mol mol^-1] from temperature [K],
pressure [Pa], and relative humidity [0--1].

# Examples

```jldoctest
julia> vapor_deficit_mole_frac(298.15, 101325.0, 0.7)
0.009371415642882845

julia> vapor_deficit_mole_frac(288.15, 5e4, 0.8)
0.006813124028440155
```
"""
vapor_deficit_mole_frac(temp, pressure, RH) =
    vapor_deficit(temp, RH) / pressure

end  # module


text:
# TODO: Remove some T's from output type signatures

#' @description
#'
#' A type that extends the normal Julia Array{T} type by allowing
#' the entries to be missing (i.e. NULL or NA) in addition to
#' taking on values of type `T`.
#'
#' @field data::Array{T, N} The elements of the DataArray are stored here,
#'        unless the given element is missing. In this case, the contents
#'        of the corresponding entry in `data` can be any arbitrary value.
#' @field na::BitArray{N} A 1-bit per entry Boolean mask specifying,
#'        for each entry of the DataArray, whether that entry is `NA`
#'        or not.
#'
#' @examples
#'
#' dv = DataArray([1, 2, 3], [false, false, true])
#'
#' dm = DataArray([1 2; 3 4], [false false; true false])
type DataArray{T, N} <: AbstractDataArray{T, N}
    data::Array{T, N}
    na::BitArray{N}

    function DataArray(d::Array{T, N}, m::BitArray{N})
        # Ensure data values and missingness metadata match
        if size(d) != size(m)
            msg = "Data and missingness arrays must be the same size"
            throw(ArgumentError(msg))
        end
        # additionally check that d does not contain NA entries
        if eltype(d) === Any
            for i in eachindex(d)
                if isassigned(d, i) && isna(d, i)
                    m[i] = true
                end
            end
        elseif eltype(d) <: NAtype
            m = trues(m)
        end
        new(d, m)
    end
end

#' @description
#'
#' An DataVector is an DataArray of order 1.
typealias DataVector{T} DataArray{T, 1}

#' @description
#'
#' An DataMatrix is an DataArray of order 2.
typealias DataMatrix{T} DataArray{T, 2}

#' @description
#'
#' Create a DataArray from a set of data values and missingness
#' mask as an `Array{Bool}`.
#'
#' NB: This definition exists because Julia requires that we redefine
#' inner DataArray constructor as an outer constuctor.
#'
#' @param d::Array The non-NA values of the DataArray being constructed.
#' @param m::BitArray A Boolean mask indicating whether each element of
#'        the DataArray is missing or not.
#'
#' @returns out::DataArray A DataArray based on the input data and
#'          missingness mask.
#'
#' @examples
#'
#' da = DataArray([1, 2, 3], falses(3))
function DataArray{T, N}(d::Array{T, N},
                         m::BitArray{N} = falses(size(d))) # -> DataArray{T}
    return DataArray{T, N}(d, m)
end

#' @description
#'
#' Create a DataArray from a set of data values and missingness
#' mask as an `Array{Bool}`.
#'
#' @param d::Array The non-NA values of the DataArray being constructed.
#' @param m::Array{Bool} A Boolean mask indicating whether each element of
#'        the DataArray is missing or not.
#'
#' @returns out::DataArray A DataArray based on the input data and
#'          missingness mask.
#'
#' @examples
#'
#' da = DataArray([1, 2, 3], [false, false, true])
function DataArray(d::Array, m::Array{Bool}) # -> DataArray{T}
    return DataArray(d, BitArray(m))
end

#' @description
#'
#' Create a DataArray of a given type and dimensionality. All
#' entries will be set to `NA` by default.
#'
#' @param T::Type The type of the output DataArray.
#' @param dims::Integer... The size of the output DataArray.
#'
#' @returns out::DataArray A DataArray of the desired type and size.
#'
#' @examples
#'
#' da = DataArray(Int, 2, 2)
function DataArray(T::Type, dims::Integer...) # -> DataArray{T}
    return DataArray(Array(T, dims...), trues(dims...))
end

#' @description
#'
#' Create a DataArray of a given type and dimensionality. All
#' entries will be set to `NA` by default.
#'
#' @param T::Type The type of the output DataArray.
#' @param dims::NTuple{N, Int} The size of the output DataArray.
#'
#' @returns out::DataArray A DataArray of the desired type and size.
#'
#' @examples
#'
#' da = DataArray(Int, (2, 2))
function DataArray{N}(T::Type, dims::NTuple{N, Int}) # -> DataArray{T}
    return DataArray(Array(T, dims...), trues(dims...))
end

#' @description
#'
#' Create a copy of a DataArray.
#'
#' @param da::DataArray The DataArray that will be copied.
#'
#' @returns out::DataArray A copy of `da`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' dv_new = copy(dv)
Base.copy(d::DataArray) = Base.copy!(similar(d), d) # -> DataArray{T}

function Base.copy!(dest::DataArray, src::DataArray) # -> DataArray{T}
    if isbits(eltype(src)) && isbits(eltype(dest))
        copy!(dest.data, src.data)
    else
        # Elements of src_data are not necessarily initialized, so
        # only copy initialized elements
        dest_data = dest.data
        src_data = src.data
        length(dest_data) >= length(src_data) || throw(BoundsError())
        src_chunks = src.na.chunks
        for i = 1:length(src_data)
            @inbounds if !Base.unsafe_bitgetindex(src_chunks, i)
                dest_data[i] = src_data[i]
            end
        end
    end
    copy!(dest.na, src.na)
    dest
end

function Base.copy!(dest::DataArray, doffs::Integer, src::DataArray) # -> DataArray{T}
    copy!(dest, doffs, src, 1, length(src))
end

# redundant on Julia 0.4
function Base.copy!(dest::DataArray, doffs::Integer, src::DataArray, soffs::Integer) # -> DataArray{T}
    soffs <= length(src) || throw(BoundsError())
    copy!(dest, doffs, src, soffs, length(src)-soffs+1)
end

function Base.copy!(dest::DataArray, doffs::Integer, src::DataArray, soffs::Integer, n::Integer) # -> DataArray{T}
    if n == 0
        return dest
    elseif n < 0
        if VERSION >= v"0.5.0-dev+4711"
            throw(ArgumentError("tried to copy n=$n elements, but n should be nonnegative"))
        else
            throw(BoundsError())
        end
    end
    if isbits(eltype(src))
        copy!(dest.data, doffs, src.data, soffs, n)
    else
        # Elements of src_data are not necessarily initialized, so
        # only copy initialized elements
        dest_data = dest.data
        src_data = src.data
        if doffs < 1 || length(dest_data) - doffs + 1 < n ||
           soffs < 1 || length(src_data) - soffs + 1 < n
            throw(BoundsError())
        end
        src_chunks = src.na.chunks
        for i = 0:(n-1)
            di, si = doffs + i, soffs + i

            @inbounds if !Base.unsafe_bitgetindex(src_chunks, si)
                dest_data[di] = src_data[si]
            end
        end
    end
    copy!(dest.na, doffs, src.na, soffs, n)
    dest
end

Base.fill!(A::DataArray, ::NAtype) = (fill!(A.na, true); A)
Base.fill!(A::DataArray, v) = (fill!(A.data, v); fill!(A.na, false); A)

#' @description
#'
#' Create a deep copy of a DataArray.
#'
#' @param da::DataArray The DataArray that will be deep copied.
#'
#' @returns out::DataArray A deep-copy of `da`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' dv_new = deepcopy(dv)
function Base.deepcopy(d::DataArray) # -> DataArray{T}
    return DataArray(deepcopy(d.data), deepcopy(d.na))
end

function Base.resize!{T}(da::DataArray{T,1}, n::Int)
    resize!(da.data, n)
    oldn = length(da.na)
    resize!(da.na, n)
    da.na[oldn+1:n] = true
    da
end

#' @description
#'
#' Create a new DataArray{T} that is similar to an existing DataArray.
#'
#' @param da::DataArray DataArray based on which a new DataArray
#'        will be created.
#' @param T::Type The element type of the output DataArray.
#' @param dims::Dims The dimensionality of the output DataArray.
#'
#' @returns out::DataArray{T} A new DataArray with the desired type and
#'          dimensionality.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' dv_new = similar(dv, Float64, 2, 2, 2)
function Base.similar(da::DataArray, T::Type, dims::Dims) #-> DataArray{T}
    return DataArray(Array(T, dims), trues(dims))
end

#' @description
#'
#' Find the sizes along each dimension of the DataArray.
#'
#' @param da::DataArray{Bool} DataArray whose size is desired.
#'
#' @returns a::Array{Int} Array containing the indices of all `true` values
#'          in `da`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' inds = find(dv)
Base.size(d::DataArray) = size(d.data) # -> (Int...)

#' @description
#'
#' Determine the number of dimensions (i.e. order) of a DataArray.
#'
#' @param da::DataArray{Bool} DataArray whose dimensionality is desired.
#'
#' @returns d::Int The number of dimensions of `da`.
#'
#' @examples
#'
#' dm = @data [false false; true false]
#' inds = ndims(dm)
Base.ndims(da::DataArray) = ndims(da.data) # -> Int

#' @description
#'
#' Determine the length (e.g. number of elements) of a DataArray.
#'
#' @param da::DataArray The DataArray whose length is desired.
#'
#' @returns n::Int The length of `da`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' i = length(dv)
Base.length(d::DataArray) = length(d.data) # -> Int

#' @description
#'
#' Find the index of the last element of `da`.
#'
#' @param da::DataArray The DataArray whose final element's index is desired.
#'
#' @returns i::Int The index of the final element of `da`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' i = endof(dv)
Base.endof(da::DataArray) = endof(da.data) # -> Int

#' @description
#'
#' Find the indices of all `true` values in a `DataArray{Bool}`.
#'
#' @param da::DataArray{Bool} DataArray whose `true` values will be found.
#'
#' @returns a::Array{Int} Array containing the indices of all `true` values
#'          in `da`.
#'
#' @examples
#'
#' dv = @data [false, false, true, false]
#' inds = find(dv)
function Base.find(da::DataArray{Bool}) # -> Array{Int}
    data = da.data
    ntrue = 0
    @inbounds @bitenumerate da.na i na begin
        ntrue += !na && data[i]
    end
    res = Array(Int, ntrue)
    count = 1
    @inbounds @bitenumerate da.na i na begin
        if !na && data[i]
            res[count] = i
            count += 1
        end
    end
    return res
end

#' @description
#'
#' Convert a DataArray{T} to an Array{S}. Throws an `NAException`
#' if the input contains `NA` values that prohibit conversion.
#'
#' @param da::DataArray{T} The DataArray that will be converted.
#'
#' @returns a::Array{S} The (possibly type-converted) elements of
#'         `da` if none were `NA`.
#'
#' @examples
#'
#' da = @data [1 2; 3 NA]
#' a = convert(Array{Float64}, da)
#'
#' da = @data [1 2; 3 4]
#' a = convert(Array{Float64}, da)
function Base.convert{S, T, N}(::Type{Array{S, N}},
                               x::DataArray{T, N}) # -> Array{S, N}
    if anyna(x)
        err = "Cannot convert DataArray with NA's to desired type"
        throw(NAException(err))
    else
        return convert(Array{S, N}, x.data)
    end
end

function Base.convert{S, T, N}(::Type{Array{S}}, da::DataArray{T, N})
    return convert(Array{S, N}, da)
end

function Base.convert{T}(::Type{Vector}, dv::DataVector{T})
    return convert(Array{T, 1}, dv)
end

function Base.convert{T}(::Type{Matrix}, dm::DataMatrix{T})
    return convert(Array{T, 2}, dm)
end

function Base.convert{T, N}(::Type{Array}, da::DataArray{T, N})
    return convert(Array{T, N}, da)
end

function Base.convert{S, T, N}(
    ::Type{Array{S, N}},
    da::DataArray{T, N},
    replacement::Any
) # -> Array{S, N}
    replacementS = convert(S, replacement)
    res = Array(S, size(da))
    for i in 1:length(da)
        if da.na[i]
            res[i] = replacementS
        else
            res[i] = da.data[i]
        end
    end
    return res
end

function Base.convert{T}(::Type{Vector}, dv::DataVector{T}, replacement::Any)
    return convert(Array{T, 1}, dv, replacement)
end

function Base.convert{T}(::Type{Matrix}, dm::DataMatrix{T}, replacement::Any)
    return convert(Array{T, 2}, dm, replacement)
end

function Base.convert{T, N}(::Type{Array}, da::DataArray{T, N}, replacement::Any)
    return convert(Array{T, N}, da, replacement)
end

#' @description
#'
#' Turn a DataVector into a Vector. Drop any NA's.
#'
#' NB: Because NA's are dropped instead of replaced, this function only
#'     works on DataVector's and will not work on DataArray's of higher
#'     order.
#'
#' @param dv::DataVector{T} DataArray that will be converted to an Array.
#'
#' @returns v::Array{T} Array containing only the non-NA values of `dv`.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' v = dropna(dv)
dropna(dv::DataVector) = copy(dv.data[!dv.na]) # -> Vector

#' @description
#'
#' Determine if the entries of an DataArray are `NA`.
#'
#' @param a::DataArray{T, N} The DataArray whose missingness will
#'        be assessed.
#'
#' @returns na::BitArray{N} Elementwise Boolean whether entry is missing.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' isna(da)
isna(da::DataArray) = copy(da.na) # -> BitArray

#' @description
#'
#' Safe and type-stable way to determine if element `i` of an
#' DataArray is `NA`.
#'
#' @param a::DataArray The DataArray whose missingness will
#'        be assessed.
#' @param inds::Any The indices of the elements to be checked for `NA`.
#'
#' @returns na::Any Are the indexed elements `NA` or not?
#'
#' @examples
#'
#' a = @data([1, 2, 3])
#' isna(a, 1)
#' isna(a, 1:2)
isna(da::DataArray, I::Real) = getindex(da.na, I)

@nsplat N function isna(da::DataArray, I::NTuple{N,Real}...)
    getindex(da.na, I...)
end

#' @description
#'
#' Determine if any of the entries of an DataArray are `NA`.
#'
#' @param da::DataArray{T, N} The DataArray whose elements will
#'        be assessed.
#'
#' @returns out::Bool Are any of the elements of `da` an `NA` value?
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' anyna(da)
anyna(da::DataArray) = any(da.na) # -> Bool

#' @description
#'
#' Determine if all of the entries of an DataArray are `NA`.
#'
#' @param da::DataArray{T, N} The DataArray whose elements will
#'        be assessed.
#'
#' @returns out::Bool Are all of the elements of `da` an `NA` value?
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' allna(da)
allna(da::DataArray) = all(da.na) # -> Bool

#' @description
#'
#' Determine if the entries of an DataArray are `NaN`.
#'
#' @param da::DataArray{T, N} The DataArray whose elements will
#'        be assessed.
#'
#' @returns na::DataArray{Bool} Elementwise Boolean whether entry is `NaN`.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' isnan(da)
function Base.isnan(da::DataArray) # -> DataArray{Bool}
    return DataArray(isnan(da.data), copy(da.na))
end

#' @description
#'
#' Determine if the entries of an DataArray are finite, which means
#' neither `+/-NaN` nor `+/-Inf`.
#'
#' @param da::DataArray{T, N} The DataArray whose elements will
#'        be assessed.
#'
#' @returns na::DataArray{Bool} Elementwise Boolean whether entry is finite.
#'
#' @examples
#'
#' da = @data([1, 2, 3])
#' isfinite(da)
function Base.isfinite(da::DataArray) # -> DataArray{Bool}
    n = length(da)
    res = Array(Bool, size(da))
    for i in 1:n
        if !da.na[i]
            res[i] = isfinite(da.data[i])
        end
    end
    return DataArray(res, copy(da.na))
end

# Promotion rules

# promote_rule{T, T}(::Type{AbstractDataArray{T}},
#                    ::Type{T}) = promote_rule(T, T)
# promote_rule{S, T}(::Type{AbstractDataArray{S}},
#                    ::Type{T}) = promote_rule(S, T)
# promote_rule{T}(::Type{AbstractDataArray{T}}, ::Type{T}) = T

#' @description
#'
#' Convert an Array{T} to a DataArray{S}.
#'
#' @param a::Array{T} The Array that will be converted.
#'
#' @returns da::DataArray{S} The converted DataArray with potential
#'          type-conversion of the elements of `a`.
#'
#' @examples
#'
#' a = [1 2; 3 4]
#' da = convert(DataArray{Float64}, a)
function Base.convert{S, T, N}(::Type{DataArray{S, N}},
                               a::AbstractArray{T, N}) # -> DataArray{S, N}
    return DataArray(convert(Array{S, N}, a), falses(size(a)))
end
Base.convert{S,T,N}(::Type{DataArray{S}}, x::AbstractArray{T,N}) =
    convert(DataArray{S,N}, x)
Base.convert{T, N}(::Type{DataArray}, x::AbstractArray{T, N}) =
    convert(DataArray{T,N}, x)

#' @description
#'
#' Convert a DataArray{T} to a DataArray{S}.
#'
#' @param da::DataArray{T} The DataArray that will be converted.
#'
#' @returns out::DataArray{S} The converted DataArray.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' dv_alt = convert(DataVector{Float64}, dv)
function Base.convert{S, T, N}(::Type{DataArray{S, N}},
                               x::DataArray{T, N}) # -> DataArray{S, N}
    v = similar(x.data, S)
    @inbounds for i = 1:length(x)
        if !x.na[i]
            v[i] = convert(S, x.data[i])
        end
    end
    return DataArray(v, x.na)
end

#' @description
#'
#' Convert an AbstractArray to a DataArray.
#'
#' @param a::AbstractArray{T} The AbstractArray that will be converted to
#'        a DataArray.
#'
#' @returns da::DataArray{T} `a` after conversion to a DataArray.
#'
#' @examples
#'
#' dv = data([1, 2, 3, 4])
#' dm = data([1 2; 3 4])
#' dv = data(1:10)
#' dv = data(falses(3))

data(a::AbstractArray) = convert(DataArray, a)

#' @description
#'
#' Convert a DataArray to an Array of float type.
#'
#' @param da::DataArray{T} The DataArray that will be converted.
#'
#' @returns a::Array{Float64} An Array containing the
#'          type-converted values of `da`.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' v = float(dv)
#
# TODO: Make sure these handle copying correctly
# TODO: Remove these? They have odd behavior, because they convert to Array's.
# TODO: Rethink multi-item documentation approach
for f in (:(Base.float),)
    @eval begin
        function ($f)(da::DataArray) # -> DataArray
            if anyna(da)
                err = "Cannot convert DataArray with NA's to desired type"
                throw(NAException(err))
            else
                ($f)(da.data)
            end
        end
    end
end

#' @internal
#' @description
#'
#' Find the unique values in a DataArray, noting if `NA` occurs in the
#' `DataArray`.
#'
#' @param da::DataArray{T} `DataArray` whose unique values are desired.
#'
#' @returns v::Vector{T} `Vector` containing the unique values from `da`.
#' @returns firstna::Int The index of the first `NA` value, or 0 if
#'          none is present.
#'
#' @examples
#'
#' dv = @data [1, 2, NA, 4]
#' distinct_values, firstna = finduniques(dv)
function finduniques{T}(da::DataArray{T}) # -> Vector{T}, Int
    out = Array(T,0)
    seen = Set{T}()
    n = length(da)
    firstna = 0
    for i in 1:n
        if isna(da, i)
            if firstna == 0
                firstna = length(out) + 1
            else
                continue
            end
        elseif !in(da.data[i], seen)
            push!(seen, da.data[i])
            push!(out, da.data[i])
        end
    end
    return out, firstna
end

#' @description
#'
#' Return a DataVector containing the unique values of a `DataArray`,
#' in the order they appear in the data, including `NA` if any missing entries
#' are encountered.
#'
#' @param da::DataArray{T} `DataArray` whose unique values are desired.
#'
#' @returns dv::DataVector{T} `DataVector` containing the unique values
#'          from `da`, in the order they appear, including `NA` if there are
#'          any missing entries in `da`.
#'
#' @examples
#'
#' dv = @data [1, -2, 1, NA, 4]
#' distinct_values = unique(dv)
function Base.unique{T}(da::DataArray{T}) # -> DataVector{T}
    unique_values, firstna = finduniques(da)
    n = length(unique_values)
    if firstna > 0
        res = DataArray(Array(T, n + 1))
        i = 1
        for val in unique_values
            if i == firstna
                res.na[i] = true
                i += 1
            end
            res.data[i] = val
            i += 1
        end

        if firstna == n + 1
            res.na[n + 1] = true
        end

        return res
    else
        return DataArray(unique_values)
    end
end

#' @description
#'
#' Return a DataVector containing the unique values of a `DataArray`,
#' excluding `NA`.
#'
#'
#' @param da::DataArray{T} `DataArray` whose unique values are desired.
#'
#' @returns dv::DataVector{T} `DataVector` containing the unique values
#'          from `da`, excluding `NA`.
#'
#' @examples
#'
#' dv = @data [1, -2, 1, NA, 4]
#' distinct_values = levels(dv)
function levels(da::DataArray) # -> DataVector{T}
    unique_values, firstna = finduniques(da)
    return DataArray(unique_values)
end

#' @description
#'
#' Return a Vector containing the unique values of an `Array`. This
#' function is identical to `unique` and is only defined for consistency
#' with `DataArray`s and `PooledDataArray`s.
#'
#'
#' @param a::Array{T} `Array` whose unique values are desired.
#'
#' @returns v::Vector{T} `Vector` containing the unique values from `da`.
#'
#' @examples
#'
#' v = [1, -2, 1, 4]
#' distinct_values = levels(v)
function levels(a::AbstractArray) # -> Vector{T}
    return unique(a)
end


text:
__precompile__()

module Xpress

    using Libdl

    # Load in `deps.jl`, complaining if it does not exist
    const depsjl_path = joinpath(@__DIR__, "..", "deps", "deps.jl")

    if isfile(depsjl_path)
        include(depsjl_path)
    elseif !haskey(ENV, "XPRESS_JL_NO_DEPS_ERROR")
        error("XPRESS cannot be loaded. Please run Pkg.build(\"Xpress\").")
    else
        const xpressdlpath = ""
    end

    const libxprs = joinpath(xpressdlpath, string(Sys.iswindows() ? "" : "lib", "xprs", ".", Libdl.dlext))

    ### imports

    import Base.show, Base.copy

    # ### exports

    ### include source files
    module Lib
        import ..Xpress
        const libxprs = Xpress.libxprs
        include("ctypes.jl")
        include("common.jl")
        include("lib.jl")
    end

    include("utils.jl")
    include("helper.jl")
    include("api.jl")
    include("xprs_callbacks.jl")
    include("license.jl")

    const XPRS_ATTRIBUTES = Dict{String, Any}(
        replace(string(name), "XPRS_"=>"") => getfield(Lib, name)
            for name in names(Lib; all = true)
            if startswith(string(name), "XPRS_") && all(isuppercase(c) || isdigit(c) for c in string(name) if c != '_')
    )

    function initialize()
        Libdl.dlopen(libxprs)
        userlic()
        init() # Call XPRSinit for initialization
        # free is not strictly necessary since destroyprob is called
        # inthe finalizer.
        # the user can call free if needed.
        # leaving it uncommented results in many print errors becaus
        # it is called prior to finalizers.
        # atexit(free) # Call free when process terminates
    end

    include("MOI/MOI_wrapper.jl")

    function __init__()
        if !haskey(ENV, "XPRESS_JL_NO_AUTO_INIT") && get(ENV, "JULIA_REGISTRYCI_AUTOMERGE", "false") != "true"
            Xpress.initialize()
        end
    end

end


text:
module Regression

using CatBoost

train_data = [[1, 4, 5, 6], [4, 5, 6, 7], [30, 40, 50, 60]]
eval_data = [[2, 4, 6, 8], [1, 4, 50, 60]]
train_labels = [10, 20, 30]

# Initialize CatBoostRegressor
model = CatBoostRegressor(; iterations=2, learning_rate=1, depth=2)

# Fit model
fit!(model, train_data, train_labels)

# Get predictions
preds = predict(model, eval_data)

end # module


text:
using Profile

function develop_nodes!(network::Network, rules)
    if network.fully_developed return end
    if length(rules) == 0 return end

    apply_to_all(fill_temp_links!, network)

    any_change = false
    num_reachable_nodes = 0
    reachable_nodes = Node[]

    apply_to_all(network) do node
        if MAX_NEURONS <= num_reachable_nodes return end

        old_name = deepcopy(node.string)
        current_reachable_nodes = _develop_node!(node, rules)

        if !any_change && (node.string != old_name || current_reachable_nodes != [node])
            any_change = true
        end

        append!(reachable_nodes, current_reachable_nodes)
        num_reachable_nodes += length(current_reachable_nodes)
    end

    if !any_change || MAX_NEURONS <= num_reachable_nodes
        network.fully_developed = true
    end

    apply_to_all(apply_temp_links!, reachable_nodes, to_temp = true)

    apply_to_all(_remove_deleted, reachable_nodes)
    filter!(node -> node.string != DELETED, reachable_nodes)

    update_inputs_outputs!(network, reachable_nodes)
    remove_non_reachable_nodes!(network)

    network.num_development_steps += 1
    if MAX_DEVELOPMENT_STEPS <= network.num_development_steps
        network.fully_developed = true
    end
end

function _develop_node!(node, rules)
    best_matching_rules = find_best_matches(node, rules)

    if 0 < length(best_matching_rules)
        rule_to_apply = first(best_matching_rules)
        return _apply_rule!(rule_to_apply, node)
    else
        return [node]
    end
end

function _apply_rule!(rule, node) # Assumes that temp are filled before calling

    # 1. replace string with sucessor_1 or remove node completely

    keep_node = (rule.successor_1 != "")

    if keep_node
        node.string = rule.successor_1
    else
        node.string = DELETED
    end

    # if we don't create a successor we're finished

    if rule.successor_2 == ""
        return keep_node ? [node] : []
    end

    # 2. create successor

    successor = Node(rule.successor_2)
    successor.type = node.type

    # 3. create ingoing inhibitory links to successor

    if rule.inhibitory_1_2 && keep_node
        push!(successor.temp_in_inhibitory, node)
    elseif rule.reverse_input
        append!(successor.temp_in_inhibitory, node.in_excitatory)
    else
        append!(successor.temp_in_inhibitory, node.in_inhibitory)
    end

    # 4. create ingoing excitatory links to successor

    if rule.excitatory_1_2 && keep_node
        push!(successor.temp_in_excitatory, node)
    elseif rule.reverse_input
        append!(successor.temp_in_excitatory, node.in_inhibitory)
    else
        append!(successor.temp_in_excitatory, node.in_excitatory)
    end

    # 5. create outgoing inhibitory links from successor

    if rule.inhibitory_2_1 && keep_node
        push!(successor.temp_out_inhibitory, node)
    elseif rule.reverse_output
        append!(successor.temp_out_inhibitory, node.out_excitatory)
    else
        append!(successor.temp_out_inhibitory, node.out_inhibitory)
    end

    # 6. create outgoing excitatory links from successor

    if rule.excitatory_2_1 && keep_node
        push!(successor.temp_out_excitatory, node)
    elseif rule.reverse_output
        append!(successor.temp_out_excitatory, node.out_inhibitory)
    else
        append!(successor.temp_out_excitatory, node.out_excitatory)
    end

    # 7. create reverse connections

    for in_node in successor.temp_in_inhibitory
        push!(in_node.temp_out_inhibitory, successor)
    end
    for out_node in successor.temp_out_inhibitory
        push!(out_node.temp_in_inhibitory, successor)
    end
    for in_node in successor.temp_in_excitatory
        push!(in_node.temp_out_excitatory, successor)
    end
    for out_node in successor.temp_out_excitatory
        push!(out_node.temp_in_excitatory, successor)
    end

    return keep_node ? [node, successor] : [successor]
end

function _remove_deleted(node)
    filter!(x -> x.string != DELETED, node.in_inhibitory)
    filter!(x -> x.string != DELETED, node.in_excitatory)
    filter!(x -> x.string != DELETED, node.out_inhibitory)
    filter!(x -> x.string != DELETED, node.out_excitatory)
end

max_stars_repo_name: tochsner/ALifeBenchmark

max_issues_repo_name: tochsner/ALifeBenchmark

max_forks_repo_name: tochsner/ALifeBenchmark


text:
export ReservoirTrajectory

using MacroTools: @forward
using Random

mutable struct ReservoirTrajectory{B,R<:AbstractRNG} <: AbstractTrajectory
    buffer::B
    n::Int
    capacity::Int
    rng::R
end

@forward ReservoirTrajectory.buffer Base.keys, Base.haskey, Base.getindex

Base.length(x::ReservoirTrajectory) = length(x.buffer[1])

function ReservoirTrajectory(capacity; n = 0, rng = Random.GLOBAL_RNG, kw...)
    buffer = VectorTrajectory(; kw...)
    ReservoirTrajectory(buffer, n, capacity, rng)
end

function Base.push!(b::ReservoirTrajectory; kw...)
    b.n += 1
    if b.n <= b.capacity
        push!(b.buffer; kw...)
    else
        i = rand(b.rng, 1:b.n)
        if i <= b.capacity
            for (k, v) in kw
                b.buffer[k][i] = v
            end
        end
    end
end

function RLBase.update!(
    trajectory::ReservoirTrajectory,
    policy::AbstractPolicy,
    env::AbstractEnv,
    ::PreActStage,
    action,
)
    s = policy isa NamedPolicy ? state(env, nameof(policy)) : state(env)
    if haskey(trajectory.buffer, :legal_actions_mask)
        lasm =
            policy isa NamedPolicy ? legal_action_space_mask(env, nameof(policy)) :
            legal_action_space_mask(env)
        push!(trajectory; :state => s, :action => action, :legal_actions_mask => lasm)
    else
        push!(trajectory; :state => s, :action => action)
    end
end


text:
# outer approximation of a control system example

# imports
using InvariantSetApproximation
# using Plots

# system model. should define x, u and w as arguments
model(x,u,w) = (x[2] + w[1], x[1] + x[2] + u[1] + w[2])

function main(model::Function, iter::Int)
    # state, input and uncertainty bounds
    Xub = [5., 5.]
    Xlb = [-5., -5.]
    Uub = [2.]
    Ulb = [-2.]
    # diameter of cell at 14 iterations is roughly 0.11
    # enlarge W by at least 0.11 (we use 0.15 here) to disturbance to satisfy theorem 4 at the 14th subdivision
    Wub = [0.3, 0.3] .+ 0.15
    Wlb = [-0.3, -0.3] .+ -0.15

    # system
    S = system(model, Xlb, Xub, ulb=Ulb, uub=Uub, wlb=Wlb, wub=Wub)

    # options
    options = Dict(:max_step=>iter, :nXsamples=>5, :Xsampletype=>:face, :nUsamples=>5, :nWsamples=>2)
    O = Options(options)

    # computation
    sol = computeISet(S,O)

    return sol
end

# run main
sol = main(model, 14);
# find convexhull of set (for faster viewing of solution)
iset_cvxh = ConvexHullArray(sol.iset)
# uncomment to plot set
# plot(sol.iset) # slow if large number of cells are present
# plot(iset_cvxh)

println("done!")