API Reference

control_dependencies(f, ops::Union{Array{PyObject}, PyObject})

Executes all operations in ops before any operations created inside the block.

op1 = tf.print("print op1")
op3 = tf.print("print op3")
control_dependencies(op1) do
    global op2 = tf.print("print op2")
end
run(sess, [op2,op3])

In this example, op1 must be executed before op2. But there is no guarantee when op3 will be executed. There are several possible outputs of the program such as

print op3
print op1
print op2

or

print op1
print op3
print op2

get_collection(name::Union{String, Missing})

Returns the collection with name name. If name is missing, returns all the trainable variables.

get_mpi()

Returns the MPI include directory and shared library.

get_mpirun()

Returns the default mpirun executable.

has_mpi(verbose::Bool = true)

Determines whether MPI is installed.

if_else(condition::Union{PyObject,Array,Bool}, fn1, fn2, args...;kwargs...)

• If condition is a scalar boolean, it outputs fn1 or fn2 (a function with no input argument or a tensor) based on whether condition is true or false.
• If condition is a boolean array, if returns condition .* fn1 + (1 - condition) .* fn2

tensorflow.python.framework.errors_impl.InvalidArgumentError: Retval[0] does not have value

independent(o::PyObject, args...; kwargs...)

Returns o but when computing the gradients, the top gradients will not be back-propagated into dependent variables of o.

reset_default_graph()

Resets the graph by removing all the operators.

tensor(s::String)

Returns the tensor with name s. See tensorname.

tensorname(o::PyObject)

Returns the name of the tensor. See tensor.

while_loop(condition::Union{PyObject,Function}, body::Function, loop_vars::Union{PyObject, Array{Any}, Array{PyObject}};
    parallel_iterations::Int64=10, kwargs...)

Loops over loop_vars while condition is true. This operator only creates one extra node to mark the loops in the computational graph.

Example

The following script computes

\[\sum_{i=1}^{10} i\]

function condition(i, ta)
    i <= 10
end
function body(i, ta)
    u = read(ta, i-1)
    ta = write(ta, i, u+1)
    i+1, ta
end
ta = TensorArray(10)
ta = write(ta, 1, constant(1.0))
i = constant(2, dtype=Int32)
_, out = while_loop(condition, body, [i, ta])
summation = stack(out)[10]

bind(op::PyObject, ops...)

Adding operations ops to the dependencies of op. ops are guaranteed to be executed before op. The function is useful when we want to execute ops but ops is not in the dependency of the final output. For example, if we want to print i each time i is evaluated

i = constant(1.0)
op = tf.print(i)
i = bind(i, op)

Session(args...; kwargs...)

Create an ADCME session. By default, ADCME will take up all the GPU resources at the start. If you want the GPU usage to grow on a need basis, before starting ADCME, you need to set the environment variable via

ENV["TF_FORCE_GPU_ALLOW_GROWTH"] = "true"

Configuration

Session accepts some runtime optimization configurations

• intra: Number of threads used within an individual op for parallelism
• inter: Number of threads used for parallelism between independent operations.
• CPU: Maximum number of CPUs to use.
• GPU: Maximum number of GPU devices to use
• soft: Set to True/enabled to facilitate operations to be placed on CPU instead of GPU

run_profile(args...;kwargs...)

Runs the session with tracing information.

save_profile(filename::String="default_timeline.json")

Save the timeline information to file filename.

• Open Chrome and navigate to chrome://tracing
• Load the timeline file

TensorArray(size_::Int64=0, args...;kwargs...)

Constructs a tensor array for while_loop.

Variable(initial_value;kwargs...)

Constructs a trainable tensor from value.

cell(arr::Array, args...;kwargs...)

Construct a cell tensor.

Example

julia> r = cell([[1.],[2.,3.]])
julia> run(sess, r[1])
1-element Array{Float32,1}:
 1.0
julia> run(sess, r[2])
2-element Array{Float32,1}:
 2.0
 3.0

constant(value; kwargs...)

Constructs a non-trainable tensor from value.

convert_to_tensor(o::Union{PyObject, Number, Array{T}, Missing, Nothing}; dtype::Union{Type, Missing}=missing) where T<:Number
convert_to_tensor(os::Array, dtypes::Array)

Converts the input o to tensor. If o is already a tensor and dtype (if provided) is the same as that of o, the operator does nothing. Otherwise, convert_to_tensor converts the numerical array to a constant tensor or casts the data type. convert_to_tensor also accepts multiple tensors.

Example

convert_to_tensor([1.0, constant(rand(2)), rand(10)], [Float32, Float64, Float32])

get_variable(o::Union{PyObject, Bool, Array{<:Number}}; 
    name::Union{String, Missing} = missing, 
    scope::String = "")

Creates a new variable with initial value o. If name exists, get_variable returns the variable instead of create a new one.

get_variable(dtype::Type;
shape::Union{Array{<:Integer}, NTuple{N, <:Integer}}, 
name::Union{Missing,String} = missing
scope::String = "")

Creates a new variable with initial value o. If name exists, get_variable returns the variable instead of create a new one.

gradient_checkpointing(type::String="speed")

Uses checkpointing scheme for gradients.

• 'speed': checkpoint all outputs of convolutions and matmuls. these ops are usually the most expensive, so checkpointing them maximizes the running speed (this is a good option if nonlinearities, concats, batchnorms, etc are taking up a lot of memory)
• 'memory': try to minimize the memory usage (currently using a very simple strategy that identifies a number of bottleneck tensors in the graph to checkpoint)
• 'collection': look for a tensorflow collection named 'checkpoints', which holds the tensors to checkpoint

gradient_magnitude(l::PyObject, o::Union{Array, PyObject})

Returns the gradient sum

\[\sqrt{\sum_{i=1}^n \|\frac{\partial l}{\partial o_i}\|^2}\]

This function is useful for debugging the training

gradients(ys::PyObject, xs::PyObject; kwargs...)

Computes the gradients of ys w.r.t xs.

• If ys is a scalar, gradients returns the gradients with the same shape as xs.
• If ys is a vector, gradients returns the Jacobian $\frac{\partial y}{\partial x}$

gradients_colocate(loss::PyObject, xs::Union{PyObject, Array{PyObject}}, args...;use_locking::Bool = true, kwargs...)

Computes the gradients of a scalar loss function loss with respect to xs. The gradients are colocated with respect to the forward pass. This function is usually in distributed computing.

hessian(ys::PyObject, xs::PyObject; kwargs...)

hessian computes the hessian of a scalar function f with respect to vector inputs xs.

Example

x = constant(rand(10))
y = 0.5 * sum(x^2)
o = hessian(y, x)

sess = Session(); init(sess)
run(sess, o) # should be an identity matrix

is_variable(o::PyObject)

Determines whether o is a trainable variable.

jacobian(y::PyObject, x::PyObject)

Computes the Jacobian matrix $J_{ij} = \frac{\partial y_i}{\partial x_j}$

ones_like(o::Union{PyObject,Real, Array{<:Real}}, args...; kwargs...)

Returns a all-one tensor, which has the same size as o.

Example

a = rand(100,10)
b = ones_like(a)
@assert run(sess, b)≈ones(100,10)

placeholder(dtype::Type; kwargs...)

Creates a placeholder of the type dtype.

Example

a = placeholder(Float64, shape=[20,10])
b = placeholder(Float64, shape=[]) # a scalar 
c = placeholder(Float64, shape=[nothing]) # a vector

placeholder(o::Union{Number, Array, PyObject}; kwargs...)

Creates a placeholder of the same type and size as o. o is the default value.

tensor(v::Array{T,2}; dtype=Float64, sparse=false) where T

tensor(v::Array{T,2}; dtype=Float64, sparse=false) where T

Convert a generic array v to a tensor. For example,

v = [0.0 constant(1.0) 2.0
    constant(2.0) 0.0 1.0]
u = tensor(v)

u will be a $2\times 3$ tensor.

zeros_like(o::Union{PyObject,Real, Array{<:Real}}, args...; kwargs...)

Returns a all-zero tensor, which has the same size as o.

Example

a = rand(100,10)
b = zeros_like(a)
@assert run(sess, b)≈zeros(100,10)

copy(o::PyObject)

Creates a tensor that has the same value that is currently stored in a variable.

read(ta::PyObject, i::Union{PyObject,Integer})

Reads data from TensorArray at index i.

write(ta::PyObject, i::Union{PyObject,Integer}, obj)

Writes data obj to TensorArray at index i.

categorical(n::Union{PyObject, Integer}; kwargs...)

kwargs has a keyword argument logits, a 2-D Tensor with shape [batch_size, num_classes]. Each slice [i, :] represents the unnormalized log-probabilities for all classes.

choice(inputs::Union{PyObject, Array}, n_samples::Union{PyObject, Integer};replace::Bool=false)

Choose n_samples samples from inputs with/without replacement.

logpdf(dist::T, x) where T<:ADCMEDistribution

Returns the log(prob) for a distribution dist.

SparseTensor

A sparse matrix object. It has two fields

• o: internal data structure

• _diag: true if the sparse matrix is marked as "diagonal".

SparseTensor(A::SparseMatrixCSC)
SparseTensor(A::Array{Float64, 2})

Creates a SparseTensor from numerical arrays.

SparseTensor(I::Union{PyObject,Array{T,1}}, J::Union{PyObject,Array{T,1}}, V::Union{Array{Float64,1}, PyObject}, m::Union{S, PyObject, Nothing}=nothing, n::Union{S, PyObject, Nothing}=nothing) where {T<:Integer, S<:Integer}

Constructs a sparse tensor. Examples:

ii = [1;2;3;4]
jj = [1;2;3;4]
vv = [1.0;1.0;1.0;1.0]
s = SparseTensor(ii, jj, vv, 4, 4)
s = SparseTensor(sprand(10,10,0.3))

Array(A::SparseTensor)

Converts a sparse tensor A to dense matrix.

RawSparseTensor(indices::Union{PyObject,Array{T,2}}, value::Union{PyObject,Array{Float64,1}},
    m::Union{PyObject,Int64}, n::Union{PyObject,Int64}; is_diag::Bool=false) where T<:Integer

A convenient wrapper for making custom operators. Here indices is 0-based.

SparseAssembler(handle::Union{PyObject, <:Integer}, n::Union{PyObject, <:Integer}, tol::Union{PyObject, <:Real}=0.0)

Creates a SparseAssembler for accumulating row, col, val for sparse matrices.

• handle: an integer handle for creating a sparse matrix. If the handle already exists, SparseAssembler return the existing sparse matrix handle. If you are creating different sparse matrices, the handles should be different.
• n: Number of rows of the sparse matrix.
• tol (optional): Tolerance. SparseAssembler will treats any values less than tol as zero.

Example 1

handle = SparseAssembler(100, 5, 1e-8)
op1 = accumulate(handle, 1, [1;2;3], [1.0;2.0;3.0])
op2 = accumulate(handle, 2, [1;2;3], [1.0;2.0;3.0])
J = assemble(5, 5, [op1;op2])

J will be a SparseTensor object.

Example 2

handle = SparseAssembler(0, 5)
op1 = accumulate(handle, 1, [1;2;3], ones(3))
op2 = accumulate(handle, 1, [3], [1.])
op3 = accumulate(handle, 2, [1;3], ones(2))
J = assemble(5, 5, [op1;op2;op3]) # op1, op2, op3 are parallel
Array(run(sess, J))≈[1.0  1.0  2.0  0.0  0.0
                1.0  0.0  1.0  0.0  0.0
                0.0  0.0  0.0  0.0  0.0
                0.0  0.0  0.0  0.0  0.0
                0.0  0.0  0.0  0.0  0.0]

assemble(m::Union{PyObject, <:Integer}, n::Union{PyObject, <:Integer}, ops::PyObject)

Assembles the sparse matrix from the ops created by accumulate. ops is either a single output from accumulate, or concated from several ops

op1 = accumulate(handle, 1, [1;2;3], [1.0;2.0;3.0])
op2 = accumulate(handle, 2, [1;2;3], [1.0;2.0;3.0])
op = [op1;op2] # equivalent to `vcat([op1, op2]...)`

m and n are rows and columns of the sparse matrix.

See SparseAssembler for an example.

compress(A::SparseTensor)

Compresses the duplicated index in A.

Example

using ADCME
indices = [
    1 1 
    1 1
    2 2
    3 3
]
v = [1.0;1.0;1.0;1.0]
A = SparseTensor(indices[:,1], indices[:,2], v, 3, 3)
Ac = compress(A)
sess = Session(); init(sess)

run(sess, A.o.indices) # expected: [0 0;0 0;1 1;2 2]
run(sess, A.o.values) # expected: [1.0;1.0;1.0;1.0]


run(sess, Ac.o.indices) # expected: [0 0;1 1;2 2]
run(sess, Ac.o.values) # expected: [2.0;1.0;1.0]

find(s::SparseTensor)

Returns the row, column and values for sparse tensor s.

scatter_update(A::Union{SparseTensor, SparseMatrixCSC{Float64,Int64}},
i1::Union{Integer, Colon, UnitRange{T}, PyObject,Array{S,1}},
i2::Union{Integer, Colon, UnitRange{T}, PyObject,Array{T,1}},
B::Union{SparseTensor, SparseMatrixCSC{Float64,Int64}})  where {S<:Real,T<:Real}

Adds B to a subblock of a sparse matrix A. Equivalently,

A[i1, i2] += B

scatter_update(A::Union{SparseTensor, SparseMatrixCSC{Float64,Int64}},
i1::Union{Integer, Colon, UnitRange{T}, PyObject,Array{S,1}},
i2::Union{Integer, Colon, UnitRange{T}, PyObject,Array{T,1}},
B::Union{SparseTensor, SparseMatrixCSC{Float64,Int64}})  where {S<:Real,T<:Real}

Updates a subblock of a sparse matrix by B. Equivalently,

A[i1, i2] = B

solve(A_factorized::Tuple{SparseTensor, PyObject}, rhs::Union{Array{Float64,1}, PyObject})

Solves the equation A_factorized * x = rhs using the factorized sparse matrix. See factorize.

spdiag(n::Int64)

Constructs a sparse identity matrix of size $n\times n$, which is equivalent to spdiag(n, 0=>ones(n))

spdiag(m::Integer, pair::Pair...)

Constructs a square $m\times m$ SparseTensor from pairs of the form

offset => array 

Example

Suppose we want to construct a $10\times 10$ tridiagonal matrix, where the lower off-diagonals are all -2, the diagonals are all 2, and the upper off-diagonals are all 3, the corresponding Julia code is

spdiag(10, -1=>-2*ones(9), 0=>2*ones(10), 1=>3ones(9))

spdiag(o::PyObject)

Constructs a sparse diagonal matrix where the diagonal entries are o, which is equivalent to spdiag(length(o), 0=>o)

spzero(m::Int64, n::Union{Missing, Int64}=missing)

Constructs a empty sparse matrix of size $m\times n$. n=m if n is missing

trisolve(a::Union{PyObject, Array{Float64,1}},b::Union{PyObject, Array{Float64,1}},
    c::Union{PyObject, Array{Float64,1}},d::Union{PyObject, Array{Float64,1}})

Solves a tridiagonal matrix linear system. The equation is as follows

\[a_i x_{i-1} + b_i x_i + c_i x_{i+1} = d_i\]

In the matrix format,

\[\begin{bmatrix} b_1 & c_1 & &0 \\ a_2 & b_2 & c_2 & \\ & a_3 & b_3 & &\\ & & & & c_{n-1}\\ 0 & & &a_n & b_n \end{bmatrix}\begin{bmatrix} x_1\\ x_2\\ \vdots \\ x_n \end{bmatrix} = \begin{bmatrix} d_1\\ d_2\\ \vdots\\ d_n\end{bmatrix}\]

zero_out_row(A::Union{SparseMatrixCSC, SparseTensor}, bd::Array{Int64, 1})

Zeros out rows of a sparse tensor A. The counterpart of numerical array function is

\[A[bd,:] .= 0.0\]

\(A::SparseTensor, o::PyObject, method::String="SparseLU")

Solves the linear equation $A x = o$

Method

For square matrices $A$, one of the following methods is available

• auto: using the solver specified by ADCME.options.sparse.solver
• SparseLU
• SparseQR
• SimplicialLDLT
• SimplicialLLT

Base.:\(A_factorized::Tuple{SparseTensor, PyObject}, rhs::Union{Array{Float64,1}, PyObject})

A convenient overload for solve. See factorize.

accumulate(handle::PyObject, row::Union{PyObject, <:Integer}, cols::Union{PyObject, Array{<:Integer}}, vals::Union{PyObject, Array{<:Real}})

Accumulates row-th row. It adds the value to the sparse matrix

for k = 1:length(cols)
    A[row, cols[k]] += vals[k]
end

handle is the handle created by SparseAssembler.

See SparseAssembler for an example.

factorize(A::Union{SparseTensor, SparseMatrixCSC}, max_cache_size::Int64 = 999999)

Factorizes $A$ for sparse matrix solutions. max_cache_size specifies the maximum cache sizes in the C++ kernels, which determines the maximum number of factorized matrices. The function returns the factorized matrix, which is basically Tuple{SparseTensor, PyObject}.

Example

A = sprand(10,10,0.7)
Afac = factorize(A) # factorizing the matrix
run(sess, Afac\rand(10)) # no factorization, solving the equation
run(sess, Afac\rand(10)) # no factorization, solving the equation

argsort(o::PyObject; 
stable::Bool = false, rev::Bool=false, dims::Integer=-1, name::Union{Nothing,String}=nothing)

Returns the indices of a tensor that give its sorted order along an axis.

batch_matmul(o1::PyObject, o2::PyObject)

Computes o1[i,:,:] * o2[i, :] or o1[i,:,:] * o2[i, :, :] for each index i.

clip(o::Union{Array{Any}, Array{PyObject}}, vmin, vmax, args...;kwargs...)

Clips the values of o to the range [vmin, vmax]

Example

a = constant(3.0)
a = clip(a, 1.0, 2.0)
b = constant(rand(3))
b = clip(b, 0.5, 1.0)

rvec(o::PyObject; kwargs...)

Vectorizes the tensor o to a column vector, assuming column major.

pad(o::PyObject, paddings::Array{Int64, 2}, args...; kwargs...)

Pads o with values on the boundary.

Example

o = rand(3,3)
o = pad(o, [1 4      # first dimension
             2 3])   # second dimension
run(sess, o)

Expected:

8×8 Array{Float64,2}:
 0.0  0.0  0.0       0.0       0.0       0.0  0.0  0.0
 0.0  0.0  0.250457  0.666905  0.823611  0.0  0.0  0.0
 0.0  0.0  0.23456   0.625145  0.646713  0.0  0.0  0.0
 0.0  0.0  0.552415  0.226417  0.67802   0.0  0.0  0.0
 0.0  0.0  0.0       0.0       0.0       0.0  0.0  0.0
 0.0  0.0  0.0       0.0       0.0       0.0  0.0  0.0
 0.0  0.0  0.0       0.0       0.0       0.0  0.0  0.0
 0.0  0.0  0.0       0.0       0.0       0.0  0.0  0.0

pmap(fn::Function, o::Union{Array{PyObject}, PyObject})

Parallel for loop. There should be no data dependency between different iterations.

Example

x = constant(ones(10))
y1 = pmap(x->2.0*x, x)
y2 = pmap(x->x[1]+x[2], [x,x])
y3 = pmap(1:10, x) do z
    i = z[1]
    xi = z[2]
    xi + cast(Float64, i)
end
run(sess, y1)
run(sess, y2)
run(sess, y3)

rollmean(u, window::Int64)

Returns the rolling mean given a window size m

\[o_k = \frac{\sum_{i=k}^{k+m-1} u_i}{m}\]

Rolling functions in ADCME:

• rollmean: rolling mean
• rollsum: rolling sum
• rollvar: rolling variance
• rollstd: rolling standard deviation

rollstd(u, window::Int64)

Returns the rolling standard deviation given a window size m

\[o_k = \sqrt{\frac{\sum_{i=k}^{k+m-1} (u_i - m_i)^2}{m-1}}\]

Here $m_i$ is the rolling mean computed using rollmean

Rolling functions in ADCME:

• rollmean: rolling mean
• rollsum: rolling sum
• rollvar: rolling variance
• rollstd: rolling standard deviation

rollsum(u, window::Int64)

Returns the rolling sum given a window size m

\[o_k = \sum_{i=k}^{k+m-1} u_i\]

Rolling functions in ADCME:

• rollmean: rolling mean
• rollsum: rolling sum
• rollvar: rolling variance
• rollstd: rolling standard deviation

rollvar(u, window::Int64)

Returns the rolling variance given a window size m

\[o_k = \frac{\sum_{i=k}^{k+m-1} (u_i - m_i)^2}{m-1}\]

Here $m_i$ is the rolling mean computed using rollmean

Rolling functions in ADCME:

• rollmean: rolling mean
• rollsum: rolling sum
• rollvar: rolling variance
• rollstd: rolling standard deviation

rvec(o::PyObject; kwargs...)

Vectorizes the tensor o to a row vector, assuming column major.

scatter_add(A::PyObject, 
    xind::Union{Colon, Int64, Array{Int64}, BitArray{1}, Array{Bool,1}, UnitRange{Int64}, StepRange{Int64, Int64}, PyObject},
    yind::Union{Colon, Int64, Array{Int64}, BitArray{1}, Array{Bool,1}, UnitRange{Int64}, StepRange{Int64, Int64}, PyObject},
    updates::Union{Array{<:Real}, Real, PyObject})

A[xind, yind] += updates

scatter_add(a::PyObject, 
    indices::Union{Colon, Int64, Array{Int64}, BitArray{1}, Array{Bool,1}, UnitRange{Int64}, StepRange{Int64, Int64}, PyObject},
    updates::Union{Array{<:Real}, Real, PyObject})

Updates array add

a[indices] += updates

Example

Julia:

A[[1;2;3]] += rand(3)
A[2] += 1.0

ADCME:

A = scatter_add(A, [1;2;3], rand(3))
A = scatter_add(A, 2, 1.0)

scatter_add(A::PyObject, 
    xind::Union{Colon, Int64, Array{Int64}, BitArray{1}, Array{Bool,1}, UnitRange{Int64}, StepRange{Int64, Int64}, PyObject},
    yind::Union{Colon, Int64, Array{Int64}, BitArray{1}, Array{Bool,1}, UnitRange{Int64}, StepRange{Int64, Int64}, PyObject},
    updates::Union{Array{<:Real}, Real, PyObject})

A[xind, yind] -= updates

scatter_sub(a::PyObject, 
    indices::Union{Colon, Int64, Array{Int64}, BitArray{1}, Array{Bool,1}, UnitRange{Int64}, StepRange{Int64, Int64}, PyObject},
    updates::Union{Array{<:Real}, Real, PyObject})

Updates array a

a[indices] -= updates

Example

Julia:

A[[1;2;3]] -= rand(3)
A[2] -= 1.0

ADCME:

A = scatter_sub(A, [1;2;3], rand(3))
A = scatter_sub(A, 2, 1.0)

scatter_update(A::PyObject, 
    xind::Union{Colon, Int64, Array{Int64}, BitArray{1}, Array{Bool,1}, UnitRange{Int64}, StepRange{Int64, Int64}, PyObject},
    yind::Union{Colon, Int64, Array{Int64}, BitArray{1}, Array{Bool,1}, UnitRange{Int64}, StepRange{Int64, Int64}, PyObject},
    updates::Union{Array{<:Real}, Real, PyObject})

A[xind, yind] = updates

scatter_update(a::PyObject, 
    indices::Union{Colon, Int64, Array{Int64}, BitArray{1}, Array{Bool,1}, UnitRange{Int64}, StepRange{Int64, Int64}, PyObject},
    updates::Union{Array{<:Real}, Real, PyObject})

Updates array a

a[indices] = updates

Example

Julia:

A[[1;2;3]] = rand(3)
A[2] = 1.0

ADCME:

A = scatter_update(A, [1;2;3], rand(3))
A = scatter_update(A, 2, 1.0)

set_shape(o::PyObject, s::Union{Array{<:Integer}, Tuple{Vararg{<:Integer, N}}}) where N
set_shape(o::PyObject, s::Integer...)

Sets the shape of o to s. s must be the actual shape of o. This function is used to convert a tensor with unknown dimensions to a tensor with concrete dimensions.

Example

a = placeholder(Float64, shape=[nothing, 10])
b = set_shape(a, 3, 10)
run(sess, b, a=>rand(3,10)) # OK 
run(sess, b, a=>rand(5,10)) # Error
run(sess, b, a=>rand(10,3)) # Error

softmax_cross_entropy_with_logits(logits::Union{Array, PyObject}, labels::Union{Array, PyObject})

Computes softmax cross entropy between logits and labels

logits is typically the output of a linear layer. For example,

logits = [
    0.124575  0.511463   0.945934
    0.538054  0.0749339  0.187802
    0.355604  0.052569   0.177009
    0.896386  0.546113   0.456832
]
labels = [2;1;2;3]

The predicted labels associated with logits is

argmax(softmax(logits), dims = 2)

Labels can also be one hot vectors

labels = [0 1
          1 0
          0 1
          0 1]

solve_batch(A::Union{PyObject, Array{<:Real, 2}}, rhs::Union{PyObject, Array{<:Real,2}})

Solves $Ax = b$ for a batch of right hand sides.

• A: a $m\times n$ matrix, where $m\geq n$
• rhs: a $n_b\times m$ matrix. Each row is a new right hand side to solve.

The returned value is a $n_b\times n$ matrix.

Example

a = rand(10,5)
b = rand(100, 10)
sol = solve_batch(a, b)
@assert run(sess, sol) ≈ (a\b')'

stack(o::PyObject)

Convert a TensorArray o to a normal tensor. The leading dimension is the size of the tensor array.

topk(o::PyObject, k::Union{PyObject,Integer}=1;
    sorted::Bool=true, name::Union{Nothing,String}=nothing)

Finds values and indices of the k largest entries for the last dimension. If sorted=true the resulting k elements will be sorted by the values in descending order.

vector(i::Union{Array{T}, PyObject, UnitRange, StepRange}, v::Union{Array{Float64},PyObject},s::Union{Int64,PyObject})

Returns a vector V with length s such that

V[i] = v

adjoint(o::PyObject; kwargs...)

Returns the conjugate adjoint of o. When the dimension of o is greater than 2, only the last two dimensions are permuted, i.e., permutedims(o, [1,2,...,n,n-1])

map(fn::Function, o::Union{Array{PyObject},PyObject};
kwargs...)

Applies fn to each element of o.

• o∈Array{PyObject} : returns [fn(x) for x in o]
• o∈PyObject : splits o according to the first dimension and then applies fn.

Example

a = constant(rand(10,5))
b = map(x->sum(x), a) # equivalent to `sum(a, dims=2)`

a = constant(ones(10))
b = constant(ones(10))
fn = x->x[1]+x[2]
c = map(fn, [a, b], dtype=Float64)

reshape(o::PyObject, s::Union{Array{<:Integer}, Tuple{Vararg{<:Integer, N}}}) where N 
reshape(o::PyObject, s::Integer; kwargs...)
reshape(o::PyObject, m::Integer, n::Integer; kwargs...)
reshape(o::PyObject, ::Colon, n::Integer)
reshape(o::PyObject, n::Integer, ::Colon)

Reshapes the tensor according to row major if the "TensorFlow style" syntax is used; otherwise reshaping according to column major is assumed.

Example

reshape(a, [10,5]) # row major 
reshape(a, 10, 5) # column major 

reverse(o::PyObject, kwargs...)

Given a tensor o, and an index dims representing the set of dimensions of tensor to reverse.

Example

a = rand(10,2)
A = constant(a)
@assert run(sess, reverse(A, dims=1)) == reverse(a, dims=1)
@assert run(sess, reverse(A, dims=2)) == reverse(a, dims=2)
@assert run(sess, reverse(A, dims=-1)) == reverse(a, dims=2)

Base.:sort(o::PyObject; 
rev::Bool=false, dims::Integer=-1, name::Union{Nothing,String}=nothing)

Sort a multidimensional array o along the given dimension.

• rev: true for DESCENDING and false (default) for ASCENDING
• dims: -1 for last dimension.

split(o::PyObject, 
    num_or_size_splits::Union{Integer, Array{<:Integer}, PyObject}; kwargs...)

Splits o according to num_or_size_splits

Example 1

a = constant(rand(10,8,6))
split(a, 5)

Expected output:

5-element Array{PyCall.PyObject,1}:
 PyObject <tf.Tensor 'split_5:0' shape=(2, 8, 6) dtype=float64>
 PyObject <tf.Tensor 'split_5:1' shape=(2, 8, 6) dtype=float64>
 PyObject <tf.Tensor 'split_5:2' shape=(2, 8, 6) dtype=float64>
 PyObject <tf.Tensor 'split_5:3' shape=(2, 8, 6) dtype=float64>
 PyObject <tf.Tensor 'split_5:4' shape=(2, 8, 6) dtype=float64>

Example 2

a = constant(rand(10,8,6))
split(a, [4,3,1], dims=2)

Expected output:

3-element Array{PyCall.PyObject,1}:
 PyObject <tf.Tensor 'split_6:0' shape=(10, 4, 6) dtype=float64>
 PyObject <tf.Tensor 'split_6:1' shape=(10, 3, 6) dtype=float64>
 PyObject <tf.Tensor 'split_6:2' shape=(10, 1, 6) dtype=float64>

Example 3

a = constant(rand(10,8,6))
split(a, 3, dims=3)

Expected output:

3-element Array{PyCall.PyObject,1}:
 PyObject <tf.Tensor 'split_7:0' shape=(10, 8, 2) dtype=float64>
 PyObject <tf.Tensor 'split_7:1' shape=(10, 8, 2) dtype=float64>
 PyObject <tf.Tensor 'split_7:2' shape=(10, 8, 2) dtype=float64>

vec(o::PyObject;kwargs...)

Vectorizes the tensor o assuming column major.

svd(o::PyObject, args...; kwargs...)

Returns a TFSVD structure which holds the following data structures

S::PyObject
U::PyObject
V::PyObject
Vt::PyObject

We have the equality $o = USV'$

Example

A = rand(10,20)
r = svd(constant(A))
A2 = r.U*diagm(r.S)*r.Vt # The value of `A2` should be equal to `A`

Diary(suffix::Union{String, Nothing}=nothing)

Creates a diary at a temporary directory path. It returns a writer and the corresponding directory path

activate(sw::Diary, port::Int64=6006)

Running Diary at http://localhost:port.

load(sess::PyObject, file::String, vars::Union{PyObject, Nothing, Array{PyObject}}=nothing, args...; kwargs...)

Loads the values of variables to the session sess from the file file. If vars is nothing, it loads values to all the trainable variables. See also save, load

load(sw::Diary, dirp::String)

Loads Diary from dirp.

logging(file::Union{Nothing,String}, o::PyObject...; summarize::Int64 = 3, sep::String = " ")

Logging o to file. This operator must be used with bind.

pload(file::String)

Loads a Python objection from file. See also psave

print_tensor(in::Union{PyObject, Array{Float64,2}})

Prints the tensor in

psave(o::PyObject, file::String)

Saves a Python objection o to file. See also pload

save(sess::PyObject, file::String, vars::Union{PyObject, Nothing, Array{PyObject}}=nothing, args...; kwargs...)

Saves the values of vars in the session sess. The result is written into file as a dictionary. If vars is nothing, it saves all the trainable variables. See also save, load

save(sw::Diary, dirp::String)

Saves Diary to dirp.

scalar(o::PyObject, name::String)

Returns a scalar summary object.

write(sw::Diary, step::Int64, cnt::Union{String, Array{String}})

Writes to Diary.

AdadeltaOptimizer(learning_rate=1e-3;kwargs...)

See AdamOptimizer for descriptions.

AdagradDAOptimizer(learning_rate=1e-3; global_step, kwargs...)

See AdamOptimizer for descriptions.

AdagradOptimizer(learning_rate=1e-3;kwargs...)

See AdamOptimizer for descriptions.

AdamOptimizer(learning_rate=1e-3;kwargs...)

Constructs an ADAM optimizer.

Example

learning_rate = 1e-3
opt = AdamOptimizer(learning_rate).minimize(loss)
sess = Session(); init(sess)
for i = 1:1000
    _, l = run(sess, [opt, loss])
    @info "Iteration $i, loss = $l")
end

Dynamical Learning Rate

We can also use dynamical learning rate. For example, if we want to use a learning rate $l_t = \frac{1}{1+t}$, we have

learning_rate = placeholder(1.0)
opt = AdamOptimizer(learning_rate).minimize(loss)
sess = Session(); init(sess)
for i = 1:1000
    _, l = run(sess, [opt, loss], lr = 1/(1+i))
    @info "Iteration $i, loss = $l")
end

The usage of other optimizers such as GradientDescentOptimizer, AdadeltaOptimizer, and so on is similar: we can just replace AdamOptimizer with the corresponding ones.

BFGS!(sess::PyObject, loss::PyObject, max_iter::Int64=15000; 
vars::Array{PyObject}=PyObject[], callback::Union{Function, Nothing}=nothing, method::String = "L-BFGS-B", kwargs...)

BFGS! is a simplified interface for L-BFGS-B optimizer. See also ScipyOptimizerInterface. callback is a callback function with signature

callback(vs::Array, iter::Int64, loss::Float64)

vars is an array consisting of tensors and its values will be the input to vs.

Example 1

a = Variable(1.0)
loss = (a - 10.0)^2
sess = Session(); init(sess)
BFGS!(sess, loss)

Example 2

θ1 = Variable(1.0)
θ2 = Variable(1.0)
loss = (θ1-1)^2 + (θ2-2)^2
cb = (vs, iter, loss)->begin 
    printstyled("[#iter $iter] θ1=$(vs[1]), θ2=$(vs[2]), loss=$loss\n", color=:green)
end
sess = Session(); init(sess)
cb(run(sess, [θ1, θ2]), 0, run(sess, loss))
BFGS!(sess, loss, 100; vars=[θ1, θ2], callback=cb)

Example 3

Use bounds to specify upper and lower bound of a variable.

x = Variable(2.0)    
loss = x^2
sess = Session(); init(sess)
BFGS!(sess, loss, bounds=Dict(x=>[1.0,3.0]))

BFGS!(sess, loss, method = "BFGS")

BFGS!(value_and_gradients_function::Function, initial_position::Union{PyObject, Array{Float64}}, max_iter::Int64=50, args...;kwargs...)

Applies the BFGS optimizer to value_and_gradients_function

BFGS!(sess::PyObject, loss::PyObject, grads::Union{Array{T},Nothing,PyObject}, 
vars::Union{Array{PyObject},PyObject}; kwargs...) where T<:Union{Nothing, PyObject}

Running BFGS algorithm $\min_{\texttt{vars}} \texttt{loss}(\texttt{vars})$ The gradients grads must be provided. Typically, grads[i] = gradients(loss, vars[i]). grads[i] can exist on different devices (GPU or CPU).

Example 1

import Optim # required
a = Variable(0.0)
loss = (a-1)^2
g = gradients(loss, a)
sess = Session(); init(sess)
BFGS!(sess, loss, g, a)

Example 2

import Optim # required
a = Variable(0.0)
loss = (a^2+a-1)^2
g = gradients(loss, a)
sess = Session(); init(sess)
cb = (vs, iter, loss)->begin 
    printstyled("[#iter $iter] a = $vs, loss=$loss\n", color=:green)
end
BFGS!(sess, loss, g, a; callback = cb)

CustomOptimizer(opt::Function, name::String)

creates a custom optimizer with struct name name. For example, we can integrate Optim.jl with ADCME by constructing a new optimizer

CustomOptimizer("Con") do f, df, c, dc, x0, x_L, x_U
    opt = Opt(:LD_MMA, length(x0))
    bd = zeros(length(x0)); bd[end-1:end] = [-Inf, 0.0]
    opt.lower_bounds = bd
    opt.xtol_rel = 1e-4
    opt.min_objective = (x,g)->(g[:]= df(x); return f(x)[1])
    inequality_constraint!(opt, (x,g)->( g[:]= dc(x);c(x)[1]), 1e-8)
    (minf,minx,ret) = NLopt.optimize(opt, x0)
    minx
end

Here

∘ f: a function that returns $f(x)$

∘ df: a function that returns $\nabla f(x)$

∘ c: a function that returns the constraints $c(x)$

∘ dc: a function that returns $\nabla c(x)$

∘ x0: initial guess

∘ nineq: number of inequality constraints

∘ neq: number of equality constraints

∘ x_L: lower bounds of optimizable variables

∘ x_U: upper bounds of optimizable variables

Then we can create an optimizer with

opt = Con(loss, inequalities=[c1], equalities=[c2])

To trigger the optimization, use

minimize(opt, sess)

Note thanks to the global variable scope of Julia, step_callback, optimizer_kwargs can actually be passed from Julia environment directly.

GradientDescentOptimizer(learning_rate=1e-3;kwargs...)

See AdamOptimizer for descriptions.

NonlinearConstrainedProblem(f::Function, L::Function, θ::PyObject, u0::Union{PyObject, Array{Float64}}; options::Union{Dict{String, T}, Missing}=missing) where T<:Integer

Computes the gradients $\frac{\partial L}{\partial \theta}$

\[\min \ L(u) \quad \mathrm{s.t.} \ F(\theta, u) = 0\]

u0 is the initial guess for the numerical solution u, see newton_raphson.

Caveats: Assume r, A = f(θ, u) and θ are the unknown parameters, gradients(r, θ) must be defined (backprop works properly)

Returns: It returns a tuple (L: loss, C: constraints, and Graidents)

\[\left(L(u), u, \frac{\partial L}{\partial θ}\right)\]

Example

We want to solve the following constrained optimization problem $\begin{aligned}\min_\theta &\; L(u) = (u-1)^3\\ \text{s.t.} &\; u^3 + u = \theta\end{aligned}$ The solution is $\theta = 2$. The Julia code is

function f(θ, u)
    u^3 + u - θ, spdiag(3u^2+1) 
end
function L(u) 
    sum((u-1)^2)
end
pl = Variable(ones(1))
l, θ, dldθ = NonlinearConstrainedProblem(f, L, pl, ones(1))

We can coupled it with a mathematical optimizer

using Optim 
sess = Session(); init(sess)
BFGS!(sess, l, dldθ, pl) 

Optimize!(sess::PyObject, loss::PyObject, max_iter::Int64 = 15000;
vars::Union{Array{PyObject},PyObject, Missing} = missing, 
grads::Union{Array{T},Nothing,PyObject, Missing} = missing, 
optimizer = missing,
callback::Union{Function, Missing}=missing,
x_tol::Union{Missing, Float64} = missing,
f_tol::Union{Missing, Float64} = missing,
g_tol::Union{Missing, Float64} = missing, kwargs...) where T<:Union{Nothing, PyObject}

An interface for using optimizers in the Optim package or custom optimizers.

• sess: a session;

• loss: a loss function;

• max_iter: maximum number of max_iterations;

• vars, grads: optimizable variables and gradients

• optimizer: Optim optimizers (default: LBFGS)

• callback: callback after each linesearch completion (NOT one step in the linesearch)

Other arguments are passed to Options in Optim optimizers.

We can also construct a custom optimizer. For example, to construct an optimizer out of Ipopt:

import Ipopt
x = Variable(rand(2))
loss = (1-x[1])^2 + 100(x[2]-x[1]^2)^2

function opt(f, g, fg, x0, kwargs...)
    prob = createProblem(2, -100ones(2), 100ones(2), 0, Float64[], Float64[], 0, 0,
                     f, (x,g)->nothing, (x,G)->g(G, x), (x, mode, rows, cols, values)->nothing, nothing)
    prob.x = x0 
    Ipopt.addOption(prob, "hessian_approximation", "limited-memory")
    status = Ipopt.solveProblem(prob)
    println(Ipopt.ApplicationReturnStatus[status])
    println(prob.x)
    Ipopt.freeProblem(prob)
    nothing
end

sess = Session(); init(sess)
Optimize!(sess, loss, optimizer = opt)

RMSPropOptimizer(learning_rate=1e-3;kwargs...)

See AdamOptimizer for descriptions.

ScipyOptimizerInterface(loss; method="L-BFGS-B", options=Dict("maxiter"=> 15000, "ftol"=>1e-12, "gtol"=>1e-12), kwargs...)

A simple interface for Scipy Optimizer. See also ScipyOptimizerMinimize and BFGS!.

ScipyOptimizerMinimize(sess::PyObject, opt::PyObject; kwargs...)

Minimizes a scalar Tensor. Variables subject to optimization are updated in-place at the end of optimization.

Note that this method does not just return a minimization Op, unlike minimize; instead it actually performs minimization by executing commands to control a Session https://www.tensorflow.org/api_docs/python/tf/contrib/opt/ScipyOptimizerInterface. See also ScipyOptimizerInterface and BFGS!.

• feed_dict: A feed dict to be passed to calls to session.run.
• fetches: A list of Tensors to fetch and supply to loss_callback as positional arguments.
• step_callback: A function to be called at each optimization step; arguments are the current values of all optimization variables packed into a single vector.
• loss_callback: A function to be called every time the loss and gradients are computed, with evaluated fetches supplied as positional arguments.
• run_kwargs: kwargs to pass to session.run.

newton_raphson(func::Function, 
    u0::Union{Array,PyObject}, 
    θ::Union{Missing,PyObject, Array{<:Real}}=missing,
    args::PyObject...) where T<:Real

Newton Raphson solver for solving a nonlinear equation. ∘ func has the signature

• func(θ::Union{Missing,PyObject}, u::PyObject)->(r::PyObject, A::Union{PyObject,SparseTensor}) (if linesearch is off)
• func(θ::Union{Missing,PyObject}, u::PyObject)->(fval::PyObject, r::PyObject, A::Union{PyObject,SparseTensor}) (if linesearch is on)

where r is the residual and A is the Jacobian matrix; in the case where linesearch is on, the function value fval must also be supplied. ∘ θ are external parameters. ∘ u0 is the initial guess for u ∘ args: additional inputs to the func function ∘ kwargs: keyword arguments to func

The solution can be configured via ADCME.options.newton_raphson

• max_iter: maximum number of iterations (default=100)
• rtol: relative tolerance for termination (default=1e-12)
• tol: absolute tolerance for termination (default=1e-12)
• LM: a float number, Levenberg-Marquardt modification $x^{k+1} = x^k - (J^k + \mu^k)^{-1}g^k$ (default=0.0)
• linesearch: whether linesearch is used (default=false)

Currently, the backtracing algorithm is implemented. The parameters for linesearch are supplied via options.newton_raphson.linesearch_options

• c1: stop criterion, $f(x^k) < f(0) + \alpha c_1 f'(0)$
• ρ_hi: the new step size $\alpha_1\leq \rho_{hi}\alpha_0$
• ρ_lo: the new step size $\alpha_1\geq \rho_{lo}\alpha_0$
• iterations: maximum number of iterations for linesearch
• maxstep: maximum allowable steps
• αinitial: initial guess for the step size $\alpha$

newton_raphson_with_grad(f::Function, 
u0::Union{Array,PyObject}, 
θ::Union{Missing,PyObject, Array{<:Real}}=missing,
args::PyObject...) where T<:Real

Differentiable Newton-Raphson algorithm. See newton_raphson.

Use ADCME.options.newton_raphson to supply options.

Example

function f(θ, x)
    x^3 - θ, 3spdiag(x^2)
end

θ = constant([2. .^3;3. ^3; 4. ^3])
x = newton_raphson_with_grad(f, constant(ones(3)), θ)
run(sess, x)≈[2.;3.;4.]
run(sess, gradients(sum(x), θ))

BatchNormalization(dims::Int64=2; kwargs...)

Creates a batch normalization layer.

Example

b = BatchNormalization(2)
x = rand(10,2)
training = placeholder(true)
y = b(x, training)
run(sess, y)

Conv1D(filters, kernel_size, strides, activation, args...;kwargs...)

c = Conv1D(32, 3, 1, "relu")
x = rand(100, 6, 128) # 128-length vectors with 6 timesteps ("channels")
y = c(x) # shape=(100, 4, 32)

Conv2D(filters, kernel_size, strides, activation, args...;kwargs...)

The arrangement is (samples, rows, cols, channels) (dataformat='channelslast')

Conv2D(32, 3, 1, "relu")

Conv3D(filters, kernel_size, strides, activation, args...;kwargs...)

The arrangement is (samples, rows, cols, channels) (dataformat='channelslast')

c = Conv3D(32, 3, 1, "relu")
x = constant(rand(100, 10, 10, 10, 16))
y = c(x)
# shape=(100, 8, 8, 8, 32)

Dense(units::Int64, activation::Union{String, Function, Nothing} = nothing,
    args...;kwargs...)

Creates a callable dense neural network.

Resnet1D(out_features::Int64, hidden_features::Int64;
    num_blocks::Int64=2, activation::Union{String, Function, Nothing} = "relu", 
    dropout_probability::Float64 = 0.0, use_batch_norm::Bool = false, name::Union{String, Missing} = missing)

Creates a 1D residual network. If name is not missing, Resnet1D does not create a new entity.

Example

resnet = Resnet1D(20)
x = rand(1000,10)
y = resnet(x)

Example: Digit recognition

using MLDatasets
using ADCME

# load data 
train_x, train_y = MNIST.traindata()
train_x = reshape(Float64.(train_x), :, size(train_x,3))'|>Array
test_x, test_y = MNIST.testdata()
test_x = reshape(Float64.(test_x), :, size(test_x,3))'|>Array

# construct loss function 
ADCME.options.training.training = placeholder(true)
x = placeholder(rand(64, 784))
l = placeholder(rand(Int64, 64))
resnet = Resnet1D(10, num_blocks=10)
y = resnet(x)
loss = mean(sparse_softmax_cross_entropy_with_logits(labels=l, logits=y))

# train the neural network 
opt = AdamOptimizer().minimize(loss)
sess = Session(); init(sess)
for i = 1:10000
    idx = rand(1:60000, 64)
    _, loss_ = run(sess, [opt, loss], feed_dict=Dict(l=>train_y[idx], x=>train_x[idx,:]))
    @info i, loss_
end

# test 
for i = 1:10
    idx = rand(1:10000,100)
    y0 = resnet(test_x[idx,:])
    y0 = run(sess, y0, ADCME.options.training.training=>false)
    pred = [x[2]-1 for x in argmax(y0, dims=2)]
    @info "Accuracy = ", sum(pred .== test_y[idx])/100
end

ae(x::PyObject, output_dims::Array{Int64}, scope::String = "default";
    activation::Union{Function,String} = "tanh")

Alias: fc, ae

Creates a neural network with intermediate numbers of neurons output_dims.

ae(x::Union{Array{Float64}, PyObject}, 
    output_dims::Array{Int64}, 
    θ::Union{Array{Array{Float64}}, Array{PyObject}};
    activation::Union{Function,String} = "tanh")

Alias: fc, ae

Constructs a neural network with given weights and biases θ

Example

x = constant(rand(10,30))
θ = ae_init([30, 20, 20, 5])
y = ae(x, [20, 20, 5], θ) # 10×5

ae(x::Union{Array{Float64}, PyObject}, output_dims::Array{Int64}, θ::Union{Array{Float64}, PyObject};
activation::Union{Function,String, Nothing} = "tanh")

Alias: fc, ae

Creates a neural network with intermediate numbers of neurons output_dims. The weights are given by θ

Example 1: Explicitly construct weights and biases

x = constant(rand(10,2))
n = ae_num([2,20,20,20,2])
θ = Variable(randn(n)*0.001)
y = ae(x, [20,20,20,2], θ)

Example 2: Implicitly construct weights and biases

θ = ae_init([10,20,20,20,2]) 
x = constant(rand(10,10))
y = ae(x, [20,20,20,2], θ)

See also ae_num, ae_init.

ae_init(output_dims::Array{Int64}; T::Type=Float64, method::String="xavier")
fc_init(output_dims::Array{Int64})

Return the initial weights and bias values by TensorFlow as a vector. The neural network architecture is

\[o_1 (\text{Input layer}) \rightarrow o_2 \rightarrow \ldots \rightarrow o_n (\text{Output layer})\]

Three types of random initializers are provided

• xavier (default). It is useful for tanh fully connected neural network.

W^l_i \sim \mathcal{N}\left(0, \sqrt{\frac{1}{n_{l-1}}}\right)

• xavier_avg. A variant of xavier

\[W^l_i \sim \mathcal{N}\left(0, \sqrt{\frac{2}{n_l + n_{l-1}}}\right)\]

• he. This is the activation aware initialization of weights and helps mitigate the problem

of vanishing/exploding gradients.

\[W^l_i \sim \mathcal{N}\left(0, \sqrt{\frac{2}{n_{l-1}}}\right)\]

Example

x = constant(rand(10,30))
θ = fc_init([30, 20, 20, 5])
y = fc(x, [20, 20, 5], θ) # 10×5

ae_num(output_dims::Array{Int64})
fc_num(output_dims::Array{Int64})

Estimates the number of weights and biases for the neural network. Note the first dimension should be the feature dimension (this is different from ae since in ae the feature dimension can be inferred), and the last dimension should be the output dimension.

Example

x = constant(rand(10,30))
θ = ae_init([30, 20, 20, 5])
@assert ae_num([30, 20, 20, 5])==length(θ)
y = ae(x, [20, 20, 5], θ)

ae_to_code(file::String, scope::String; activation::String = "tanh")

Return the code string from the feed-forward neural network data in file. Usually we can immediately evaluate the code string into Julia session by

eval(Meta.parse(s))

If activation is not specified, tanh is the default.

bn(args...;center = true, scale=true, kwargs...)

bn accepts a keyword parameter is_training.

Example

bn(inputs, name="batch_norm", is_training=true)

update_ops = get_collection(UPDATE_OPS)
control_dependencies(update_ops) do 
    global train_step = AdamOptimizer().minimize(loss)
end 

dense(inputs::Union{PyObject, Array{<:Real}}, units::Int64, args...; 
    activation::Union{String, Function} = nothing, kwargs...)

Creates a fully connected layer with the activation function specified by activation

dropout(x::Union{PyObject, Real, Array{<:Real}}, 
rate::Union{Real, PyObject}, training::Union{PyObject,Bool} = true; kwargs...)

Randomly drops out entries in x with a rate of rate.

ae(x::PyObject, output_dims::Array{Int64}, scope::String = "default";
    activation::Union{Function,String} = "tanh")

Alias: fc, ae

Creates a neural network with intermediate numbers of neurons output_dims.

ae(x::Union{Array{Float64}, PyObject}, output_dims::Array{Int64}, θ::Union{Array{Float64}, PyObject};
activation::Union{Function,String, Nothing} = "tanh")

Alias: fc, ae

Creates a neural network with intermediate numbers of neurons output_dims. The weights are given by θ

Example 1: Explicitly construct weights and biases

x = constant(rand(10,2))
n = ae_num([2,20,20,20,2])
θ = Variable(randn(n)*0.001)
y = ae(x, [20,20,20,2], θ)

Example 2: Implicitly construct weights and biases

θ = ae_init([10,20,20,20,2]) 
x = constant(rand(10,10))
y = ae(x, [20,20,20,2], θ)

See also ae_num, ae_init.

ae(x::Union{Array{Float64}, PyObject}, 
    output_dims::Array{Int64}, 
    θ::Union{Array{Array{Float64}}, Array{PyObject}};
    activation::Union{Function,String} = "tanh")

Alias: fc, ae

Constructs a neural network with given weights and biases θ

Example

x = constant(rand(10,30))
θ = ae_init([30, 20, 20, 5])
y = ae(x, [20, 20, 5], θ) # 10×5

ae_init(output_dims::Array{Int64}; T::Type=Float64, method::String="xavier")
fc_init(output_dims::Array{Int64})

Return the initial weights and bias values by TensorFlow as a vector. The neural network architecture is

o1 (\text{Input layer}) \rightarrow o2 \rightarrow \ldots \rightarrow o_n (\text{Output layer}) $

Three types of random initializers are provided

• xavier (default). It is useful for tanh fully connected neural network.

W^l_i \sim \mathcal{N}\left(0, \sqrt{\frac{1}{n_{l-1}}}\right)

• xavier_avg. A variant of xavier

W^li \sim \mathcal{N}\left(0, \sqrt{\frac{2}{nl + n_{l-1}}}\right) $

• he. This is the activation aware initialization of weights and helps mitigate the problem

of vanishing/exploding gradients.

W^li \sim \mathcal{N}\left(0, \sqrt{\frac{2}{n{l-1}}}\right) $

Example

x = constant(rand(10,30))
θ = fc_init([30, 20, 20, 5])
y = fc(x, [20, 20, 5], θ) # 10×5

ae_num(output_dims::Array{Int64})
fc_num(output_dims::Array{Int64})

Estimates the number of weights and biases for the neural network. Note the first dimension should be the feature dimension (this is different from ae since in ae the feature dimension can be inferred), and the last dimension should be the output dimension.

Example

x = constant(rand(10,30))
θ = ae_init([30, 20, 20, 5])
@assert ae_num([30, 20, 20, 5])==length(θ)
y = ae(x, [20, 20, 5], θ)

fcx(x::Union{Array{Float64,2},PyObject}, output_dims::Array{Int64,1}, 
θ::Union{Array{Float64,1}, PyObject};
activation::String = "tanh")

Creates a fully connected neural network with output dimension o and inputs $x\in \mathbb{R}^{m\times n}$.

\[x \rightarrow o_1 \rightarrow o_2 \rightarrow \ldots \rightarrow o_k\]

θ is the weights and biases of the neural network, e.g., θ = ae_init(output_dims).

fcx outputs two tensors:

• the output of the neural network: $u\in \mathbb{R}^{m\times o_k}$.

• the sensitivity of the neural network per sample: $\frac{\partial u}{\partial x}\in \mathbb{R}^{m \times o_k \times n}$

GAN(dat::Union{Array,PyObject}, generator::Function, discriminator::Function,
loss::Union{Missing, Function}=missing; latent_dim::Union{Missing, Int64}=missing,
    batch_size::Int64=32)

Creates a GAN instance.

• dat $\in \mathbb{R}^{n\times d}$ is the training data for the GAN, where $n$ is the number of training data, and $d$ is the dimension per training data.
• generator$:\mathbb{R}^{d'} \rightarrow \mathbb{R}^d$ is the generator function, $d'$ is the hidden dimension.
• discriminator$:\mathbb{R}^{d} \rightarrow \mathbb{R}$ is the discriminator function.
• loss is the loss function. See klgan, rklgan, wgan, lsgan for examples.
• latent_dim (default=$d$, the same as output dimension) is the latent dimension.
• batch_size (default=32) is the batch size in training.

Example: Constructing a GAN

dat = rand(10000,10)
generator = (z, gan)->10*z
discriminator = (x, gan)->sum(x)
gan = GAN(dat, generator, discriminator, "wgan_stable")

Example: Learning a Gaussian random variable

using ADCME 
using PyPlot
using Distributions
dat = randn(10000, 1) * 0.5 .+ 3.0
function gen(z, gan)
    ae(z, [20,20,20,1], "generator_$(gan.ganid)", activation = "relu")
end
function disc(x, gan)
    squeeze(ae(x, [20,20,20,1], "discriminator_$(gan.ganid)", activation = "relu"))
end
gan = GAN(dat, gen, disc, g->wgan_stable(g, 0.001); latent_dim = 10)

dopt = AdamOptimizer(0.0002, beta1=0.5, beta2=0.9).minimize(gan.d_loss, var_list=gan.d_vars)
gopt = AdamOptimizer(0.0002, beta1=0.5, beta2=0.9).minimize(gan.g_loss, var_list=gan.g_vars)
sess = Session(); init(sess)
for i = 1:5000
    batch_x = rand(1:10000, 32)
    batch_z = randn(32, 10)
    for n_critic = 1:1
        global _, dl = run(sess, [dopt, gan.d_loss], 
                feed_dict=Dict(gan.ids=>batch_x, gan.noise=>batch_z))
    end
    _, gl, gm, dm, gp = run(sess, [gopt, gan.g_loss, 
        gan.STORAGE["g_grad_magnitude"], gan.STORAGE["d_grad_magnitude"], 
        gan.STORAGE["gradient_penalty"]],
        feed_dict=Dict(gan.ids=>batch_x, gan.noise=>batch_z))
    mod(i, 100)==0 && (@info i, dl, gl, gm, dm, gp)
end

hist(run(sess, squeeze(rand(gan,10000))), bins=50, density = true)
nm = Normal(3.0,0.5)
x0 = 1.0:0.01:5.0
y0 = pdf.(nm, x0)
plot(x0, y0, "g")

build!(gan::GAN)

Builds the GAN instances. This function returns gan for convenience.

jsgan(gan::GAN)

Computes the vanilla GAN loss function.

klgan(gan::GAN)

Computes the KL-divergence GAN loss function.

lsgan(gan::GAN)

Computes the least square GAN loss function.

predict(gan::GAN, input::Union{PyObject, Array})

Predicts the GAN gan output given input input.

rklgan(gan::GAN)

Computes the reverse KL-divergence GAN loss function.

sample(gan::GAN, n::Int64)
rand(gan::GAN, n::Int64)

Samples n instances from gan.

wgan(gan::GAN)

Computes the Wasserstein GAN loss function.

wgan_stable(gan::GAN, λ::Float64)

Returns the discriminator and generator loss for the Wasserstein GAN loss with penalty parameter $\lambda$

The objective function is

\[L = E_{\tilde x\sim P_g} [D(\tilde x)] - E_{x\sim P_r} [D(x)] + \lambda E_{\hat x\sim P_{\hat x}}[(||\nabla_{\hat x}D(\hat x)||^2-1)^2]\]

MCMCSimple(obs::Array{Float64, 1}, h::Function, 
σ::Float64, θ0::Array{Float64,1}, lb::Float64, ub::Float64)

A very simple yet useful interface for MCMC simulation in many scientific computing problems.

• obs: Observations
• h: Forward computation function
• σ: Noise standard deviation for the observed data
• ub, lb: upper and lower bound
• θ0: Initial guess

The mathematical model is

\[y_{obs} = h(\theta)\]

and we have a hard constraint lb\leq \theta \leq ub.

cmake(DIR::String=".."; CMAKE_ARGS::Union{Array{String}, String} = "")

The built-in Cmake command for building C/C++ libraries. If extra Cmake arguments are needed, please specify it through CMAKE_ARGS.

Example

ADCME.cmake(CMAKE_ARGS=["SHARED=YES", "STAITC=NO"])

The executed command might be:

/home/darve/kailaix/.julia/adcme/bin/cmake -G Ninja -DCMAKE_MAKE_PROGRAM=/home/darve/kailaix/.julia/adcme/bin/ninja -DJULIA=/home/darve/kailaix/julia-1.3.1/bin/julia -DCMAKE_C_COMPILER=/home/darve/kailaix/.julia/adcme/bin/x86_64-conda_cos6-linux-gnu-gcc -DCMAKE_CXX_COMPILER=/home/darve/kailaix/.julia/adcme/bin/x86_64-conda_cos6-linux-gnu-g++ SHARED=YES STATIC=NO ..

compile(s::String; force::Bool=false)

Compiles the library given by path deps/s. If force is false, compile first check whether the binary product exists. If the binary product exists, return 2. Otherwise, compile tries to compile the binary product, and returns 0 if successful; it return 1 otherwise.

compile()

Compile a custom operator in the current directory. A CMakeLists.txt must be present.

customop(;with_mpi::Bool = false)

Create a new custom operator. Typically users call customop twice: the first call generates a customop.txt, users edit the content in the file; the second all generates C++ source code, CMakeLists.txt, and gradtest.jl from customop.txt.

Example

julia> customop() # create an editable `customop.txt` file
[ Info: Edit custom_op.txt for custom operators
julia> customop() # after editing `customop.txt`, call it again to generate interface files.

Options

• with_mpi: Whether the custom operator uses MPI

debug(libfile::String = "")

Loading custom operator shared library. If the loading fails, detailed error message is printed.

debug(sess::PyObject, o::PyObject)

In the case a session run yields an error from the TensorFlow backend, this function can help print the exact error. For example, you might encounter InvalidArgumentError() with no detailed error information, and this function can be useful for debugging.

doctor()

Reports health of the current installed ADCME package. If some components are broken, possible fix is proposed.

get_library_symbols(file::Union{String, PyObject})

Returns the symbols in the custom op library file.

get_placement()

Returns the operation placements.

install(s::String; force::Bool = false, islocal::Bool = false)

Install a custom operator from a URL, a directory (when islocal is true), or a string. In any of the three case, install copy the folder to /home/runner/work/ADCME.jl/ADCME.jl/deps/CustomOps/Plugin. When s is a string, s is converted to

https://github.com/ADCMEMarket/<s>

load_library(filename::String)

Load custom operator libraries. If used with

load_op(oplibpath::Union{PyObject, String}, opname::String; verbose::Union{Missing, Bool} = missing)

Loads the operator opname from library oplibpath.

load_op_and_grad(oplibpath::Union{PyObject, String}, opname::String; multiple::Bool=false)

Loads the operator opname from library oplibpath; gradients are also imported. If multiple is true, the operator is assumed to have multiple outputs.

load_system_op(opname::String, grad::Bool=true; multiple::Bool=false)

Loads custom operator from CustomOps directory (shipped with ADCME instead of TensorFlow) For example

s = "SparseOperator"
oplib = "libSO"
grad = true

this will direct Julia to find library CustomOps/SparseOperator/libSO.dylib on MACOSX

make_library(Libdir::String)

Make shared library in Libdir. The structure of the source codes files are

- Libdir 
  - *.cpp 
  - *.h 
  - CMakeLists
  - build (Optional)

nnuq(H::Array{Float64,2}, invR::Union{Float64, Array{Float64,2}}, invQ::Union{Float64, Array{Float64,2}})

Returns the variance matrix for the Baysian inversion.

The negative log likelihood function is

\[l(s) =\frac{1}{2} (y-h(s))^T R^{-1} (y-h(s)) + \frac{1}{2} s^T Q^{-1} s\]

The covariance matrix is computed by first linearizing $h(s)$

\[h(s)\approx h(s_0) + \nabla h(s_0) (s-s_0)\]

and then computing the second order derivative

\[V = \left(\frac{\partial^2 l}{\partial s^T\partial s}\right)^{-1} = (H^T R^{-1} H + Q^{-1})^{-1}\]

Note the result is independent of $s_0$, $y_0$, and only depends on $\nabla h(s_0)$

register(forward::Function, backward::Function; multiple::Bool=false)

Register a function forward with back-propagated gradients rule backward to the backward. ∘ forward: it takes $n$ inputs and outputs $m$ tensors. When $m>1$, the keyword multiple must be true. ∘ backward: it takes $\tilde m$ top gradients from float/double output tensors of forward, $m$ outputs of the forward, and $n$ inputs of the forward. backward outputs $n$ gradients for each input of forward. When input $i$ of forward is not float/double, backward should return nothing for the corresponding gradients.

Example

forward = x->log(1+exp(x))
backward = (dy, y, x)->dy*(1-1/(1+y))
f = register(forward, backward)

sleep_for(t::Union{PyObject, <:Real})

Sleeps for t seconds.

test_gpu()

Tests the GPU ultilities

timestamp(deps::Union{PyObject, <:Real, Missing}=missing)

These functions are usually used with bind for profiling. Note the timing is not very accurate in a multithreaded environment.

• deps: deps is always executed before returning the timestamp.

Example

a = constant(3.0)
t0 = timestamp(a)
sleep_time = sleep_for(a)
t1 = timestamp(sleep_time)
sess = Session(); init(sess)
t0_, t1_ = run(sess, [t0, t1])
time = t1_ - t0_

xavier_init(size, dtype=Float64)

Returns a matrix of size size and its values are from Xavier initialization.

precompile(force::Bool=false)

Precompile the built-in custom operators.

animate(update::Function, frames; kwargs...)

Creates an animation using update function update.

Example

θ = LinRange(0, 2π, 100)
x = cos.(θ)
y = sin.(θ)
pl, = plot([], [], "o-")
t = title("0")
xlim(-1.2,1.2)
ylim(-1.2,1.2)
function update(i)
    t.set_text("$i")
    pl.set_data([x[1:i] y[1:i]]'|>Array)
end
animate(update, 1:100)

gradview(sess::PyObject, pl::PyObject, loss::PyObject, u0; scale::Float64 = 1.0)

Visualizes the automatic differentiation and finite difference convergence converge. For correctly implemented differentiable codes, the convergence rate for AD should be 2 and for FD should be 1 (if not evaluated at stationary point).

• scale: you can control the step size for perturbation.

jacview(sess::PyObject, f::Function, θ::Union{Array{Float64}, PyObject, Missing}, 
u0::Array{Float64}, args...)

Performs gradient test for a vector function. f has the signature

f(θ, u) -> r, J

Here θ is a nuisance parameter, u is the state variables (w.r.t. which the Jacobian is computed), r is the residual vector, and J is the Jacobian matrix (a dense matrix or a SparseTensor).

Example 1

u0 = rand(10)
function verify_jacobian_f(θ, u)
    r = u^3+u - u0
    r, spdiag(3u^2+1.0)
end
jacview(sess, verify_jacobian_f, missing, u0)

Example 2

u0 = rand(10)
rs = rand(10)
function verify_jacobian_f(θ, u)
    r = [u^2;u] - [rs;rs]
    r, [spdiag(2*u); spdiag(10)]
end
jacview(sess, verify_jacobian_f, missing, u0); close("all")

lineview(sess::PyObject, pl::PyObject, loss::PyObject, θ1, θ2=nothing; n::Integer = 10)

Plots the function

\[h(α) = f((1-α)θ_1 + αθ_2)\]

Example

pl = placeholder(Float64, shape=[2])
l = sum(pl^2-pl*0.1)
sess = Session(); init(sess)
lineview(sess, pl, l, rand(2))

saveanim(anim::PyObject, filename::String; kwargs...)

Saves the animation produced by animate

test_gradients(f::Function, x0::Array{Float64, 1}; scale::Float64 = 1.0, showfig::Bool = true)

Testing the gradients of a vector function f: y, J = f(x) where y is a scalar output and J is the vector gradient.

test_hessian(f::Function, x0::Array{Float64, 1}; scale::Float64 = 1.0)

Testing the Hessian of a scalar function f: g, H = f(x) where y is a scalar output, g is a vector gradient output, and H is the Hessian.

test_jacobian(f::Function, x0::Array{Float64, 1}; scale::Float64 = 1.0, showfig::Bool = true)

Testing the gradients of a vector function f: y, J = f(x) where y is a vector output and J is the Jacobian.

Database(filename::Union{Missing, String} = missing; 
    commit_after_execute::Bool = true)

Creates a database from filename. If filename is not provided, an in-memory database is created. If commit_after_execute is false, no commit operation is performed after each execute.

• do block syntax:

Database() do db
    execute(db, "create table mytable (a real, b real)")
end

The database is automatically closed after execution. Therefore, if execute is a query operation, users need to store the results in a global variable.

• Query meta information

keys(db) # all tables 
keys(db, "mytable") # all column names in `db.mytable` 

commit(db::Database)

Commits changes to db.

execute(db::Database, sql::String, args...)

Executes the SQL statement sql in db. Users can also use the do block syntax.

execute(db) do 
    "create table mytable (a real, b real)"
end

execute can also be used to insert a batch of records

t1 = rand(10)
t2 = rand(10)
param = collect(zip(t1, t2))
execute(db, "INSERT TO mytable VALUES (?,?)", param)

or

execute(db, "INSERT TO mytable VALUES (?,?)", t1, t2)

ExplicitNewmark(M::Union{SparseTensor, SparseMatrixCSC}, Z1::Union{Missing, SparseTensor, SparseMatrixCSC}, Z2::Union{Missing, SparseTensor, SparseMatrixCSC}, Δt::Float64)

An explicit Newmark integrator for

\[M \ddot{\mathbf{d}} + Z_1 \dot{\mathbf{d}} + Z_2 \mathbf{d} + f = 0\]

The numerical scheme is

\[\left(\frac{1}{\Delta t^2} M + \frac{1}{2\Delta t}Z_1\right)d^{n+1} = \left(\frac{2}{\Delta t^2} M - Z_2\right)d^n - \left(\frac{1}{\Delta t^2} M - \frac{1}{2\Delta t}Z_1\right) d^{n-1} - f\]

To use this integrator,

en = ExplicitNewmark(M, Z1, Z2, Δt)
d2 = step(en, d0, d1, f)

Newmark(M::Union{SparseTensor, PyObject, Array{Float64, 2}, SparseMatrixCSC}, 
                    C::Union{SparseTensor, PyObject, Array{Float64, 2}, SparseMatrixCSC}, 
                    K::Union{SparseTensor, PyObject, Array{Float64, 2}, SparseMatrixCSC}, 
                    Δt::Float64; β::Float64 = 1/4, γ::Float64 = 1/2, factorize_S::Bool = true)

Creates a Newmark algorithm object. This function uses "a-form" implementation for

\[M\ddot \mathbf{d} + C \dot\mathbf{d} + K \mathbf{d} + f = 0\]

The numerical scheme is

\[\begin{aligned} \tilde d_{n+1} &= d_n + \Delta t v_n + \Delta^2/2(1-2\beta)a_n \\ \tilde v_{n+1} &= v_n + (1-\gamma)\Delta a_n \\ (M+\gamma\Delta C+\beta \Delta^2 K)a_{n+1} &= F_{n+1} - C\tilde v_{n+1} - K\tilde d_{n+1}\\ d_{n+1} &= \tilde d_{n+1} + \beta \Delta^2 a_{n+1}\\ v_{n+1} &= \tilde v_{n+1} + \gamma \Delta ta_n \end{aligned}\]

For unconditional stability $2\beta \geq \gamma \geq \frac12$

For an explicit form, see ExplicitNewmark, which corresponds to $\beta = 0, \gamma = \frac12$.

A common choice for Newmark algorithm is $\beta = \frac14, \gamma = \frac12$ (a.k.a., average acceleration, or Trapezoidal rule), which is unconditionally stable, has second order of accuracy.

To perform time marching, use

step(nm::Newmark, d::PyObject, v::PyObject, a::PyObject, 
        F::Union{Array{Float64, 1}, PyObject})

TR_BDF2(D0::Union{SparseTensor, SparseMatrixCSC}, 
    D1::Union{SparseTensor, SparseMatrixCSC}, 
    Δt::Float64)

Constructs a TR-BDF2 (the Trapezoidal Rule with Second Order Backward Difference Formula) handler for the DAE

\[D_1 \dot y + D_0 y = f\]

The struct is a functor, which performs one step simulation

(tr::TR_BDF2)(y::Union{PyObject, Array{Float64, 1}}, 
    f1::Union{PyObject, Array{Float64, 1}}, 
    f2::Union{PyObject, Array{Float64, 1}}, 
    f3::Union{PyObject, Array{Float64, 1}})

Here f1, f2, and f3 correspond to the right hand side at time step $n$, $n+\frac12$, and $n+1$.

Or we can pass a batched F defined as a (2NT+1) × DOF array

(tr::TR_BDF2)(y0::Union{PyObject, Array{Float64, 1}}, 
    F::Union{PyObject, Array{Float64, 2}})

The output will be the entire solution of size (NT+1) × DOF.

constant(tr::TR_BDF2)

Converts tr to a symbolic solver.

ode45(y::Union{PyObject, Float64, Array{Float64}}, T::Union{PyObject, Float64}, 
            NT::Union{PyObject,Int64}, f::Function, θ::Union{PyObject, Missing}=missing)

Solves

\[\frac{dy}{dt} = f(y, t, \theta)\]

with six-stage, fifth-order, Runge-Kutta method.

rk4(y::Union{PyObject, Float64, Array{Float64}}, T::Union{PyObject, Float64}, 
            NT::Union{PyObject,Int64}, f::Function, θ::Union{PyObject, Missing}=missing)

Solves

\[\frac{dy}{dt} = f(y, t, \theta)\]

with Runge-Kutta (order 4) method.

runge_kutta(f::Function, T::Union{PyObject, Float64}, 
            NT::Union{PyObject,Int64}, y::Union{PyObject, Float64, Array{Float64}}, θ::Union{PyObject, Missing}=missing; method::String="rk4")

Solves

\[\frac{dy}{dt} = f(y, t, \theta)\]

with Runge-Kutta method.

For example, the default solver, RK4, has the following numerical scheme per time step

\[\begin{aligned} k_1 &= \Delta t f(t_n, y_n, \theta)\\ k_2 &= \Delta t f(t_n+\Delta t/2, y_n + k_1/2, \theta)\\ k_3 &= \Delta t f(t_n+\Delta t/2, y_n + k_2/2, \theta)\\ k_4 &= \Delta t f(t_n+\Delta t, y_n + k_3, \theta)\\ y_{n+1} &= y_n + \frac{k_1}{6} +\frac{k_2}{3} +\frac{k_3}{3} +\frac{k_4}{6} \end{aligned}\]

αscheme(M::Union{SparseTensor, SparseMatrixCSC}, 
    C::Union{SparseTensor, SparseMatrixCSC}, 
    K::Union{SparseTensor, SparseMatrixCSC}, 
    Force::Union{Array{Float64}, PyObject}, 
    d0::Union{Array{Float64, 1}, PyObject}, 
    v0::Union{Array{Float64, 1}, PyObject}, 
    a0::Union{Array{Float64, 1}, PyObject}, 
    Δt::Array{Float64}; 
    solve::Union{Missing, Function} = missing,
    extsolve::Union{Missing, Function} = missing, 
    ρ::Float64 = 1.0)

Generalized α-scheme. $M u_{tt} + C u_{t} + K u = F$

Force must be an array of size n×p, where d0, v0, and a0 have a size p Δt is an array (variable time step).

The generalized α scheme solves the equation by the time stepping

\[\begin{aligned} \bf d_{n+1} &= \bf d_n + h\bf v_n + h^2 \left(\left(\frac{1}{2}-\beta_2 \right)\bf a_n + \beta_2 \bf a_{n+1} \right)\\ \bf v_{n+1} &= \bf v_n + h((1-\gamma_2)\bf a_n + \gamma_2 \bf a_{n+1})\\ \bf F(t_{n+1-\alpha_{f_2}}) &= M \bf a _{n+1-\alpha_{m_2}} + C \bf v_{n+1-\alpha_{f_2}} + K \bf{d}_{n+1-\alpha_{f_2}} \end{aligned}\]

where

\[\begin{aligned} \bf d_{n+1-\alpha_{f_2}} &= (1-\alpha_{f_2})\bf d_{n+1} + \alpha_{f_2} \bf d_n\\ \bf v_{n+1-\alpha_{f_2}} &= (1-\alpha_{f_2}) \bf v_{n+1} + \alpha_{f_2} \bf v_n \\ \bf a_{n+1-\alpha_{m_2} } &= (1-\alpha_{m_2}) \bf a_{n+1} + \alpha_{m_2} \bf a_n\\ t_{n+1-\alpha_{f_2}} & = (1-\alpha_{f_2}) t_{n+1 + \alpha_{f_2}} + \alpha_{f_2}t_n \end{aligned}\]

Here the parameters are computed using

\[\begin{aligned} \gamma_2 &= \frac{1}{2} - \alpha_{m_2} + \alpha_{f_2}\\ \beta_2 &= \frac{1}{4} (1-\alpha_{m_2}+\alpha_{f_2})^2 \\ \alpha_{m_2} &= \frac{2\rho_\infty-1}{\rho_\infty+1}\\ \alpha_{f_2} &= \frac{\rho_\infty}{\rho_\infty+1} \end{aligned}\]

∘ solve: users can provide a solver function, solve(A, rhs) for solving Ax = rhs ∘ extsolve: similar to solve, but the signature has the form

extsolve(A, rhs, i)

This provides the users with more control, e.g., (time-dependent) Dirichlet boundary conditions. See Generalized α Scheme for details.

αscheme_time(Δt::Array{Float64}; ρ::Float64 = 1.0)

Returns the integration time $t_{i+1-\alpha_{f_2}}$ between $[t_i, t_{i+1}]$ using the alpha scheme. If $\Delta t$ has length $n$, the output will also have length $n$.

function RBF2D(xc::Union{PyObject, Array{Float64, 1}}, yc::Union{PyObject, Array{Float64, 1}}; 
    c::Union{PyObject, Array{Float64, 1}, Missing} = missing, 
    eps::Union{PyObject, Array{Float64, 1}, Real, Missing} = missing,
    d::Union{PyObject, Array{Float64, 1}} = zeros(0), 
    kind::Int64 = 0)

Constructs a radial basis function representation on a 2D domain

\[f(x, y) = \sum_{i=1}^N c_i \phi(r; \epsilon_i) + d_0 + d_1 x + d_2 y\]

Here d can be either 0, 1 (only $d_0$ is present), or 3 ($d_0$, $d_1$, and $d_2$ are all present).

kind determines the type of radial basis functions

• 0:Gaussian

\[\phi(r; \epsilon) = e^{-(\epsilon r)^2}\]

• 1:Multiquadric

\[\phi(r; \epsilon) = \sqrt{1+(\epsilon r)^2}\]

• 2:Inverse quadratic

\[\phi(r; \epsilon) = \frac{1}{1+(\epsilon r)^2}\]

• 3:Inverse multiquadric

\[\phi(r; \epsilon) = \frac{1}{\sqrt{1+(\epsilon r)^2}}\]

Returns a callable struct, i.e. to evaluates the function at locations $(x, y)$ (x and y are both vectors), run

rbf(x, y)

RBF3D(xc::Union{PyObject, Array{Float64, 1}}, yc::Union{PyObject, Array{Float64, 1}},
    zc::Union{PyObject, Array{Float64, 1}}; 
    c::Union{PyObject, Array{Float64, 1}, Missing} = missing, 
    eps::Union{PyObject, Array{Float64, 1}, Real, Missing} = missing,
    d::Union{PyObject, Array{Float64, 1}} = zeros(0), 
    kind::Int64 = 0)

Constructs a radial basis function representation on a 3D domain

\[f(x, y, z) = \sum_{i=1}^N c_i \phi(r; \epsilon_i) + d_0 + d_1 x + d_2 y + d_3 z\]

Here d can be either 0, 1 (only $d_0$ is present), or 4 ($d_0$, $d_1$, $d_2$, and $d_3$ are all present).

kind determines the type of radial basis functions

• 0:Gaussian

\[\phi(r; \epsilon) = e^{-(\epsilon r)^2}\]

• 1:Multiquadric

\[\phi(r; \epsilon) = \sqrt{1+(\epsilon r)^2}\]

• 2:Inverse quadratic

\[\phi(r; \epsilon) = \frac{1}{1+(\epsilon r)^2}\]

• 3:Inverse multiquadric

\[\phi(r; \epsilon) = \frac{1}{\sqrt{1+(\epsilon r)^2}}\]

Returns a callable struct, i.e. to evaluates the function at locations $(x, y, z)$ (x, y, and z are all vectors), run

rbf(x, y, z)

interp1(x::Union{Array{Float64, 1}, PyObject},v::Union{Array{Float64, 1}, PyObject},xq::Union{Array{Float64, 1}, PyObject})

returns interpolated values of a 1-D function at specific query points using linear interpolation. Vector x contains the sample points, and v contains the corresponding values, v(x). Vector xq contains the coordinates of the query points.

Example

x = sort(rand(10))
y = constant(@. x^2 + 1.0)
z = [x[1]; x[2]; rand(5) * (x[end]-x[1]) .+ x[1]; x[end]]
u = interp1(x,y,z)

dtw(s::Union{PyObject, Array{Float64}}, t::Union{PyObject, Array{Float64}}, 
    use_fast::Bool = false)

Computes the dynamic time wrapping (DTW) distance between two time series s and t. Returns the distance and path. use_fast specifies whether fast algorithm is used. Note fast algorithm may not be accurate.

emd(a::Union{PyObject, Array{Float64}}, b::Union{PyObject, Array{Float64}}, M::Union{PyObject, Array{Float64}};
iter::Int64 = 1000, tol::Float64 = 1e-9, returnall::Bool=false)

Computes the Earth Mover's Distance, which is defined as

\[D(M) = \sum_{i=1}^m \sum_{j=1}^n M_{ij} d_{ij}\]

Here $M \in \mathbb{R}^{m\times n}$ is the ground distance matrix. The algorithm solves the following optimization problem

\[\begin{aligned}\min_{M} &\ D(M)\\\text{s.t.} & \ \sum_{i=1}^m M_{ij} = b_j\\ &\ \sum_{j=1}^n M_{ij} = a_i \end{aligned}\]

The internal solver for the optimization problem is a netflow solver. The algorithm requires $\sum_i a_i = \sum_j b_j = 1$.

empirical_emd(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}};
    iter::Int64 = 1000, tol::Float64 = 1e-9, dist::Union{Integer,Function}=2, returnall::Bool=false)

Same as empirical_sinkhorn, except that the Earth Mover Distance is computed.

empirical_sinkhorn(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}};
    reg::Union{PyObject,Float64} = 1.0, iter::Int64 = 1000, tol::Float64 = 1e-9, method::String="sinkhorn", dist::Function=dist, returnall::Bool=false)

Computes the empirical Sinkhorn distance with sinkhorn algorithm. Here $x$ and $y$ are samples from two distributions.

• reg (default = 1.0): entropy regularization parameter
• tol (default = 1e-9), iter (default = 1000): tolerance and max iterations for the Sinkhorn algorithm
• dist (default = 2): Integer or Function, the distance function between two samples; if dist is integer, $L-dist$ norm is used.
• returnall: returns (TransportMatrix, Loss) if true; otherwise, only Loss is returned.

The implementation are adapted from https://github.com/rflamary/POT.

ot_dist(x::Union{PyObject, Array{Float64}}, y::Union{PyObject, Array{Float64}}, order::Union{Int64, PyObject}=2)

Computes the distance function with norm order. dist returns a $n\times m$ matrix, where $x\in \mathbb{R}^{n\times d}$ and $y\in \mathbb{R}^{m\times d}$, and the return $M\in \mathbb{R}^{n\times m}$

\[M_{ij} = ||x_i - y_j||_{o}\]

ot_plot1D(a::Array{Float64, 1}, b::Array{Float64, 1}, M::Array{Float64, 2})

Plots the optimal transport matrix for 1D distributions.

sinkhorn(a::Union{PyObject, Array{Float64}}, b::Union{PyObject, Array{Float64}}, M::Union{PyObject, Array{Float64}};
reg::Float64 = 1.0, iter::Int64 = 1000, tol::Float64 = 1e-9, method::String="sinkhorn")

Computes the optimal transport with Sinkhorn algorithm. The mathematical formulation is

\[\begin{aligned} \arg\min_P &\ \left(P, M\right) + \lambda \Omega(\Gamma)\\ \text{s.t.} &\ \Gamma 1 = a\\ &\ \Gamma^T 1 = b\\ & \Gamma \geq 0 \end{aligned}\]

Here $\Omega$ is the entropic regularization. Note if $\lambda$ is very small, the algorithm may encounter numerical instabilities.

The implementation are adapted from https://github.com/rflamary/POT.

mutable struct mpi_SparseTensor
    rows::PyObject 
    ncols::PyObject
    cols::PyObject 
    values::PyObject 
    ilower::Int64 
    iupper::Int64 
    N::Int64
    oplibpath::String
end

A structure to hold local data of a sparse matrix. The global matrix is assumed to be a $M\times N$ square matrix. The current processor owns rows from ilower to iupper (inclusive). The data is specified by

• rows: an array indicating the rows that contain nonzero values. Note rows ≥ ilower.
• ncols: an array indicating the number of nonzero values for each row in rows.
• cols: the column indices for nonzero values. Its length is $\sum_{i=1}^{\mathrm{ncols}} \mathrm{ncols}_i$
• vals: the nonzero values corresponding to each column index in cols
• oplibpath: the backend library (returned by ADCME.load_plugin_MPITensor)

All data structure are 0-based. Note if we work with a linear solver, $M=N$.

For example, consider the sparse matrix

[  1 0 0 1  ]
[  0 1 2 1  ]

We have

rows = Int32[0;1]
ncols = Int32[2;3]
cols = Int32[0;3;1,2,3]
values = [1.;1.;1.;2.;1.]
iupper = ilower + 2

mpi_SparseTensor(sp::Union{SparseTensor, SparseMatrixCSC{Float64,Int64}}, 
    ilower::Union{Int64, Missing} = missing,
    iupper::Union{Int64, Missing} = missing)

Constructing mpi_SparseTensor from a SparseTensor or a sparse Array.

mpi_SparseTensor(rows::Union{Array{Int32,1}, PyObject}, ncols::Union{Array{Int32,1}, PyObject}, cols::Union{Array{Int32,1}, PyObject},
    vals::Union{Array{Float64,1}, PyObject}, ilower::Int64, iupper::Int64, N::Int64)

Create a $N\times N$ distributed sparse tensor A for the current MPI processor. The current MPI processor owns rows with indices [ilower, iupper]. The submatrix is specified using the CSR format.

• rows: an array indicating the rows that contain nonzero values. Note rows ≥ ilower.
• ncols: an array indicating the number of nonzero values for each row in rows.
• cols: the column indices for nonzero values. Its length is $\sum_{i=1}^{\mathrm{ncols}} \mathrm{ncols}_i$
• vals: the nonzero values corresponding to each column index in cols

Note that by default the indices are zero-based.

mpi_barrier(o::PyObject)

Returns 0. Impose a barrier for object o.

using ADCME 

mpi_init()
s = constant(rand(10))
t = mpi_bcast(s)
t = t + 1.0 + mpi_barrier(t)
q = mpi_bcast(t)
sess = Session(); init(sess)
run(sess, q)

mpi_bcast(a::Union{Array{Float64}, Float64, PyObject}, root::Int64 = 0)

Broadcast a from processor root to all other processors.

mpi_finalize()

Finalize the MPI call.

mpi_finalized()

Returns a boolean indicating whether the current MPI session is finalized.

mpi_gather(u::Union{Array{Float64, 1}, PyObject}, deps::Union{Missing, PyObject} = missing)

Gathers all the vectors from different processes to the root process. The function returns a long vector which concatenates of local vectors in the order of process IDs.

mpi_halo_exchange(u::Union{Array{Float64, 2}, PyObject},m::Int64,n::Int64; deps::Union{Missing, PyObject} = missing,
fill_value::Float64 = 0.0, tag::Union{PyObject, Int64} = 0)

Perform Halo exchnage on u (a $k \times k$ matrix). The output has a shape $(k+2)\times (k+2)$

• fill_value: value used for the boundaries
• tag: message tag
• deps: a scalar tensor; it can be used to serialize the MPI calls

mpi_halo_exchange2(u::Union{Array{Float64, 2}, PyObject},m::Int64,n::Int64; deps::Union{Missing, PyObject} = missing,
fill_value::Float64 = 0.0, tag::Union{PyObject, Int64} = 0)

Similar to mpi_halo_exchange, but the reach is 2, i.e., for a $N\times N$ matrix $u$, the output will be a $(N+4)\times (N+4)$ matrix.

mpi_init()

Initialized the MPI session. mpi_init must be called before any run(sess, ...).

mpi_initialized()

Returns a boolean indicating whether the current MPI session is initialized.

mpi_optimize(_f::Function, _g!::Function, x0::Array{Float64};
    method::String = "LBFGS", options = missing)

Performs the MPI-enabled optimization.

Example

mpirun -n 4 julia adcme_inverse.jl

You can find the file adcme_inverse.jl here.

mpi_rank()

Returns the rank of current MPI process (rank 0 based).

mpi_recv(a::Union{Array{Float64}, Float64, PyObject}, src::Int64, tag::Int64 = 0)

Receives an array from processor src. mpi_recv requires an input for gradient backpropagation. Typically we can write

r = mpi_rank()
a = constant(Float64(r))
if r==1
    a = mpi_send(a, 0)
end
if r==0
    a = mpi_recv(a, 1)
end

Then a=1 on both processor 0 and processor 1.

mpi_send(a::Union{Array{Float64}, Float64, PyObject}, dest::Int64,root::Int64 = 0)

Sends a to processor dest. a itself is returned so that the send action can be added to the computational graph.

mpi_sendrecv(a::Union{Array{Float64}, Float64, PyObject}, dest::Int64, src::Int64, tag::Int64=0)

A convenient wrapper for mpi_send followed by mpi_recv.

mpi_size()

Returns the size of MPI world.

mpi_sum(a::Union{Array{Float64}, Float64, PyObject}, root::Int64 = 0)

Sum a on the MPI processor root.

mpi_sync!(message::Array{Int64,1}, root::Int64 = 0)
mpi_sync!(message::Array{Float64,1}, root::Int64 = 0)

Sync message across all MPI processors.

require_mpi()

Throws an error if mpi_init() has not been called.

adjoint(A::mpi_SparseTensor)

Returns the adjoint of A, i.e., A'. Each MPI rank owns the same number of rows.

change_directory(directory::Union{Missing, AbstractString})

Change the current working directory to directory. If directory does not exist, it is made.

If directory is missing, the default is ADCME.PREFIXDIR.

copy_file(src::String, dest::String)

Copy file src to dest

get_conda()

Returns the conda executable location.

get_gfortran()

Install a gfortran compiler if it does not exist.

get_library(filename::AbstractString)

Returns a valid library file. For example, for filename = "adcme", we have

• On MacOS, the function returns libadcme.dylib
• On Linux, the function returns libadcme.so
• On Windows, the function returns adcme.dll

get_library_name(filename::AbstractString)

Returns the OS-dependent library name

Example

get_library_name("mylibrary")

• Windows: mylibrary.dll
• MacOS: libmylibrary.dylib
• Linux: libmylibrary.so

get_pip()

Returns the location for pip

git_repository(url::AbstractString, file::AbstractString)

Clone a repository url and rename it to file.

http_file(url::AbstractString, file::AbstractString)

Download a file from url and rename it to file.

link_file(target::AbstractString, link::AbstractString)

Make a symbolic link link -> target

make_directory(directory::AbstractString)

Make a directory if it does not exist.

read_with_env(cmd::Cmd, env::Union{Missing, Dict} = missing)

Similar to run_with_env, but returns a string containing the output.

require_cmakecache(func::Function, DIR::String = ".")

Check if cmake has output something. If not, func is executed.

require_file(f::Function, file::Union{String, Array{String}})

If any of the files/links/directories in file does not exist, execute f.

require_import(s::Symbol)

Checks whether the package s is imported in the Main namespace. Returns the package handle.

require_library(func::Function, filename::AbstractString)

If the library file filename does not exist, func is executed.

run_with_env(cmd::Cmd, env::Union{Missing, Dict} = missing)

Running the command with the default environment and an extra environment variables env

uncompress(zipfile::AbstractString, file::AbstractString)

Uncompress a zip file zipfile to file (a directory). Note this function does not check that the uncompressed content has the name file. It is used as a hint to skip uncompress action.

Users may use mv uncompress_file file to enforce the consistency.

get_gpu()

Returns the compiler information for GPUs. An examplary output is

(NVCC = "/usr/local/cuda/bin/nvcc", LIB = "/usr/local/cuda/lib64", INC = "/usr/local/cuda/include", TOOLKIT = "/home/kailaix/.julia/adcme/pkgs/cudatoolkit-10.0.130-0/lib", CUDNN = "/home/kailaix/.julia/adcme/pkgs/cudnn-7.6.5-cuda10.0_0/lib")

gpu_info()

Returns the CUDA and GPU information. An examplary output is

- NVCC: /usr/local/cuda/bin/nvcc
- CUDA library directories: /usr/local/cuda/lib64
- Latest supported version of CUDA: 11000
- CUDA runtime version: 10010
- CUDA include_directories: /usr/local/cuda/include
- CUDA toolkit directories: /home/kailaix/.julia/adcme/pkgs/cudatoolkit-10.0.130-0/lib:/home/kailaix/.julia/adcme/pkgs/cudnn-7.6.5-cuda10.0_0/lib
- Number of GPUs: 7

has_gpu()

Check if the TensorFlow backend is using CUDA GPUs. Operators that have GPU implementations will be executed on GPU devices. See also get_gpu