text:
# https://www.juliabox.com/notebook/notebooks/tutorials/intro-to-julia/AutoDiff.ipynb
# Dual Number

using InteractiveUtils
versioninfo()

struct D <: Number
    f::Tuple{Float64,Float64}
end

Base.show(io::IO, x::D) = print(io,x.f[1]," + ", x.f[2],"ϵ")

import Base: +, -, *, /, convert, promote_rule
+(x::D, y::D) = D(x.f .+ y.f)
-(x::D, y::D) = D(x.f .- y.f)
*(x::D, y::D) = D((x.f[1] .* y.f[1], x.f[2] * y.f[1] + x.f[1] * y.f[2]))
/(x::D, y::D) = D((x.f[1] / y.f[1], (y.f[1] * x.f[2] - x.f[1] * y.f[2]) / y.f[1]^2))
convert(::Type{D}, x::Real) = D((x,zero(x)))
promote_rule(::Type{D}, ::Type{<:Number}) = D

ϵ = D((0,1))

@show ϵ * ϵ, ϵ^2
@show 1/(1+ϵ)
@show (1+ϵ)^10

function nthroot(x, n=2; t=1, N = 10)
    for i in 1:N
        t += (x/t^(n-1)-t)/n
    end
    t
end

@show nthroot(2,3), ∛2
@show nthroot(7,12), 7^(1/12)
x=17.0; @show nthroot(x+ϵ, 3), ∛x, 1/x^(2/3)/3

#= output
Julia Version 0.7.0
Commit a4cb80f3ed (2018-08-08 06:46 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin14.5.0)
  CPU: Intel(R) Core(TM) M-5Y51 CPU @ 1.10GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, broadwell)
(ϵ * ϵ, ϵ ^ 2) = (0.0 + 0.0ϵ, 0.0 + 0.0ϵ)
1 / (1 + ϵ) = 1.0 + -1.0ϵ
(1 + ϵ) ^ 10 = 1.0 + 10.0ϵ
(nthroot(2, 3), ∛2) = (1.2599210498948732, 1.2599210498948732)
(nthroot(7, 12), 7 ^ (1 / 12)) = (1.1760474285795146, 1.1760474285795146)
(nthroot(x + ϵ, 3), ∛x, (1 / x ^ (2 / 3)) / 3) = (2.571281590658235 + 0.05041728609133794ϵ, 2.571281590658235, 0.05041728609133795)
=#

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
"""
"""

"""----------------------------------------------------------------------------
    Computes the contraction C(1,x) using the indicies of the tensor
  non-zeros to compute the marginalization step and then the vector contraction
  step after.
----------------------------------------------------------------------------"""
function implicit_kronecker_model_marginalization_contraction!(B_indices::Array{Array{Int,2},1},
                                                               B_vals::Array{Array{Float64,1},1},
                                                               H_indices::Array{Array{Int,2},1},
                                                               x::Array{Float64,1},
                                                               y::Array{Float64,1},
                                                               n::Int,m::Int)

    for i in 1:n
        for j in 1:m
           y[(i-1)*m + j] = 0.0
       end
    end

    for (Hn_indices,Bn_indices,Bn_vals) in zip(H_indices,B_indices,B_vals)

        for i = 1:size(Hn_indices,1)

            for p in permutations([1,2,3])
                _,i2,i3  = Hn_indices[i,p]

                for j = 1:size(Bn_indices,1)
                    _,j2,j3 = Bn_indices[j,p]
                    y[(i2-1)*m + j2] +=  Bn_vals[j]*x[(i3-1)*m + j3]

                end
            end
        end
    end
end

#= TODO: verify this can be retired
function implicit_kronecker_mode1_marginalization(Bn_indices::Array{Int,2},
    Bn_vals::Array{Float64,1},Hn_indices::Array{Int,2},n,m)

    ei = Array{Int,1}(undef,6*size(Hn_indices,1)*length(Bn_vals))
    ej = Array{Int,1}(undef,6*size(Hn_indices,1)*length(Bn_vals))
    vals = Array{Float64,1}(undef,6*size(Hn_indices,1)*length(Bn_vals))

    index = 1
    for i = 1:size(Hn_indices,1)

        i1 = Hn_indices[i,1]
        i2 = Hn_indices[i,2]
        i3 = Hn_indices[i,3]

        for j = 1:size(Bn_indices,1)

            #_,j2,j3 = Bn_indices[j,:]
            j2 = Bn_indices[j,2]
            j3 = Bn_indices[j,3]

            ei[index] = (i2-1)*m + j2
            ej[index] = (i3-1)*m + j3
            vals[index] = Bn_vals[j]
            index += 1

            ei[index] = (i3-1)*m + j2
            ej[index] = (i2-1)*m + j3
            vals[index] = Bn_vals[j]
            index += 1

            ei[index] = (i1-1)*m + j2
            ej[index] = (i2-1)*m + j3
            vals[index] = Bn_vals[j]
            index += 1

            ei[index] = (i2-1)*m + j2
            ej[index] = (i1-1)*m + j3
            vals[index] = Bn_vals[j]
            index += 1

            ei[index] = (i1-1)*m + j2
            ej[index] = (i3-1)*m + j3
            vals[index] = Bn_vals[j]
            index += 1

            ei[index] = (i3-1)*m + j2
            ej[index] = (i1-1)*m + j3
            vals[index] = Bn_vals[j]
            index += 1
        end
    end

    return sparse(ei,ej,vals,n*m,n*m)
    #return ei,ej,vals,n*m
end
=#

function sym_mode1_marginalization(indices::Array{Int,2},vals::Array{Float64,1},m::Int)
    n,d = size(indices)

    ei = Array{Int,1}(undef,6*length(vals))
    ej = Array{Int,1}(undef,6*length(vals))
    new_vals = Array{Float64,1}(undef,6*length(vals))

    index = 1
    for i =1:n
        i1,i2,i3 = indices[i,:]

        ei[index] = i1
        ej[index] = i2
        new_vals[index] = vals[i]
        index += 1

        ei[index] = i2
        ej[index] = i1
        new_vals[index] = vals[i]
        index += 1

        ei[index] = i2
        ej[index] = i3
        new_vals[index] = vals[i]
        index += 1

        ei[index] = i3
        ej[index] = i2
        new_vals[index] = vals[i]
        index += 1

        ei[index] = i1
        ej[index] = i3
        new_vals[index] = vals[i]
        index += 1

        ei[index] = i3
        ej[index] = i1
        new_vals[index] = vals[i]
        index += 1
    end

    return sparse(ei,ej,new_vals,m,m)
end

function mode1_marginalization(indices::Array{Int,2},vals::Array{Float64,1},m::Int)
    n,d = size(indices)

    ei = Array{Int,1}(undef, length(vals))
    ej = Array{Int,1}(undef,length(vals))
    new_vals = Array{Float64,1}(undef, length(vals))

    index = 1
    for i =1:n

        ei[index] = indices[index,2]
        ej[index] = indices[index,3]
        new_vals[index] = vals[index]
        index += 1

    end

    return sparse(ei,ej,new_vals,m,m)
end

"""----------------------------------------------------------------------------
    Computes the contraction C(1,x) using list of pairs of the marginalized Hn
  Bn pairs. Hn is typically quite sparse, so it's more efficient to figure out
  which rows of the folded multiplicand should be multiplied to Bn, rather
  than doing the explicit HnXBnᵀ operation.
----------------------------------------------------------------------------"""
function HOFASM_contraction!(tensor_pairs::Array{NTuple{2,SparseMatrixCSC{Float64,Int64}},1},
                             X_k::Array{Float64,2},X_k_1::Array{Float64,2}) #where {F <: AbstractFloat}

    n,m = size(X_k_1)

    for i=1:n
       for j=1:m
           X_k_1[i,j] = 0.0
       end
    end

    for (Hn,Bn) in tensor_pairs
#         X_k_1 .+= (Bn*(Hn*X_k)' )'
        is,js,vs = findnz(Hn)
        for (i,j,v) in zip(is,js,vs)
            if v > 1.0
                X_k_1[i,:]  .+= v*Bn*X_k[j,:]
            else
                X_k_1[i,:]  .+= Bn*X_k[j,:]
            end
        end
    end

end

"""----------------------------------------------------------------------------
    Computes the contraction C(1,x) using list of pairs of the marginalized Hn
  Bn pairs. Hn is typically quite sparse, so it's more efficient to figure out
  which rows of the folded multiplicand should be multiplied to Bn, rather
  than doing the explicit HnXBnᵀ operation. The input here is a preprocessed
  version of each of the Hn, so that the findnz operation isn't called
  repeatedly.
----------------------------------------------------------------------------"""
function HOFASM_contraction!(tensor_pairs::Array{Tuple{Array{Int,1},Array{Int,1},Array{Float64,1},
                                                       SparseMatrixCSC{Float64,Int64}},1},
                             X_k::Array{Float64,2},X_k_1::Array{Float64,2})

    n,m = size(X_k_1)
    for i=1:n
       for j=1:m
           X_k_1[i,j] = 0.0
       end
    end

    for (is,js,vs,Bn) in tensor_pairs
        #is,js,_ = findnz(Hn)
        for (i,j,v) in zip(is,js,vs)
            if v > 1.0
                X_k_1[i,:]  .+= v*Bn*X_k[j,:]
            else
                X_k_1[i,:]  .+= Bn*X_k[j,:]
            end
        end
    end

end

"""----------------------------------------------------------------------------
----------------------------------------------------------------------------"""
function HOFASM_contraction!(marginalized_tensors::Array{SparseMatrixCSC{Float64,Int64},1},
                             x_k::Array{Float64,1},x_k_1::Array{Float64,1})

    for i=1:length(x_k_1)
       x_k_1[i] = 0.0
    end

    for Tn in marginalized_tensors
        x_k_1 .+= Tn*x_k
    end

end

#TODO: need to correct the output
#TODO: check to see if this version is faster now that we don't explicitly compute HnX
function SIMD_sparse_matrix_multiplication(Hn_matrices::Array{SparseMatrixCSC{Float64,Int64},1},
                                    Bn_matrices::Array{SparseMatrixCSC{Float64,Int64},1},
                                    X::Array{Float64,2},res_loc::Array{Float64,2})

    #pack the vectors X into groups of 8
    n,m= size(X)
    num_segs = Int(ceil(m/8))
    X_segments = Array{Array{Float64,2},1}(undef,num_segs)

    for i=1:num_segs
        if i == num_segs

            X_segments[i] = X'[:,(i-1)*8+1:end]
        else
            X_segments[i] = X'[:,(i-1)*8+1:i*8]
        end

    end

    packed_segments = [pack(x) for x in X_segments]

    temp_segs = Array{Array{Float64,2},1}(undef,num_segs)

    for (Hn,Bn) in zip(Hn_matrices,Bn_matrices)

        #compute Bn*X'
        temp = zeros(m,n)
        for i=1:num_segs

            if i == num_segs
                temp[:,(i-1)*8+1:end] = unpack(Bn*packed_segments[i])
            else
                temp[:,(i-1)*8+1:i*8] = unpack(Bn*packed_segments[i])
            end

        end

        #repack temp segments
        for i=1:num_segs
            if i == num_segs
                temp_segs[i] = temp'[:,(i-1)*8+1:end]
            else
               temp_segs[i] = temp'[:,(i-1)*8+1:i*8]
            end
        end
        packed_temp_segs = [pack(x) for x in temp_segs]

        #compute A*(BX')' and store in result location
        for i = 1:num_segs
            if i == num_segs
                res_loc[:,(i-1)*8+1:end] = unpack(Hn*packed_temp_segs[i])
            else
                res_loc[:,(i-1)*8+1:i*8] = unpack(Hn*packed_temp_segs[i])
            end
        end
    end

end

"""----------------------------------------------------------------------------
    Computes the marginalization of a 3rd order tensor with respect to a
  generalized transpose operation.
----------------------------------------------------------------------------"""
#allows the application of a permutation before contracting
function perm_mode1_marginalization(indices::Array{Int,2},vals::Array{Float64,1},
                                    perm::Array{Int,1},m::Int)
    n,d = size(indices)

    ei = Array{Int,1}(undef, length(vals))
    ej = Array{Int,1}(undef,length(vals))
    new_vals = Array{Float64,1}(undef, length(vals))

    index = 1
    for i =1:n

        ei[index] = indices[index,perm[2]]
        ej[index] = indices[index,perm[3]]
        new_vals[index] = vals[index]
        index += 1

    end
    index -= 1 #counter the last extra index
    return sparse(ei[1:index],ej[1:index],new_vals[1:index],m,m)::SparseMatrixCSC{Float64,Int64}
end

"""----------------------------------------------------------------------------
    Computes the marginalization of ∑ Hn^{<σ>} ⨂ Bn^{<σ>} from the indices of
  their tensor representations.
----------------------------------------------------------------------------"""
function perm_marginalize(H_indices::Array{Int64,2},B_indices::Array{Int64,2},
                          B_vals::Array{Float64,1},n::Int,m::Int)

    B_nnz = size(B_indices,1)
    H_nnz = size(H_indices,1)

    ei = zeros(Int,B_nnz*H_nnz*6)
    ej = zeros(Int,B_nnz*H_nnz*6)
    vals = zeros(Float64,B_nnz*H_nnz*6)
    index = 1

    for i in 1:H_nnz
        for p in permutations([1,2,3])

            i1,i2,i3 = H_indices[i,p]

            for j in 1:B_nnz

                (_,j2,j3) = B_indices[j,p]
                ei[index] = (i2-1)*m + j2
                ej[index] = (i3-1)*m + j3
                vals[index] = B_vals[j]
                index += 1

            end
        end
    end

    return sparse(ei,ej,vals,n*m,n*m)
end

ext: jl


text:
# squared Euclidean distance from a set

export SqrDistL2

"""
**Squared distance from a convex set**

    SqrDistL2(ind_S, λ=1.0)

Given `ind_S` the indicator function of a convex set ``S``, and an optional positive parameter `λ`, returns the (weighted) squared Euclidean distance from ``S``, that is function
```math
g(x) = \\tfrac{λ}{2}\\mathrm{dist}_S^2(x) = \\min \\left\\{ \\tfrac{λ}{2}\\|y - x\\|^2 : y \\in S \\right\\}.
```
"""
struct SqrDistL2{R <: Real, T <: ProximableFunction} <: ProximableFunction
  ind::T
  lambda::R
  function SqrDistL2{R,T}(ind::T, lambda::R) where {R <: Real, T<:ProximableFunction}
    if !is_convex(ind) || !is_set(ind)
      error("`ind` must be the indicator of a convex set")
    end
    if lambda < 0
      error("parameter λ must be nonnegative")
    else
      new(ind, lambda)
    end
  end
end

is_prox_accurate(f::SqrDistL2) = is_prox_accurate(f.ind)
is_convex(f::SqrDistL2) = true
is_smooth(f::SqrDistL2) = true
is_quadratic(f::SqrDistL2) = is_affine(f.ind)
is_strongly_convex(f::SqrDistL2) = is_singleton(f.ind)

SqrDistL2(ind::T, lambda::R=1.0) where {R <: Real, T <: ProximableFunction} = SqrDistL2{R, T}(ind, lambda)

function (f::SqrDistL2)(x::AbstractArray{T}) where T <: RealOrComplex
  p, = prox(f.ind, x)
  return (f.lambda/2)*normdiff2(x,p)
end

function prox!(y::AbstractArray{T}, f::SqrDistL2, x::AbstractArray{T}, gamma::Real=1.0) where T <: RealOrComplex
  p, = prox(f.ind, x)
  sqrd = (f.lambda/2)*normdiff2(x,p)
  c1 = 1/(1+f.lambda*gamma)
  c2 = f.lambda*gamma*c1
  for k in eachindex(p)
    y[k] = c1*x[k] + c2*p[k]
  end
  return sqrd*c1^2
end

function gradient!(y::AbstractArray{T}, f::SqrDistL2, x::AbstractArray{T}) where T <: RealOrComplex
  p, = prox(f.ind, x)
  dist2 = normdiff2(x,p)
  y .= f.lambda.*(x .- p)
  return (f.lambda/2)*dist2
end

fun_name(f::SqrDistL2) = "squared Euclidean distance from a convex set"
fun_dom(f::SqrDistL2) = fun_dom(f.ind)
fun_expr(f::SqrDistL2) = "x ↦ (λ/2) inf { ||x-y||^2 : y ∈ S} "
fun_params(f::SqrDistL2) = string("λ = $(f.lambda), S = ", typeof(f.ind))

function prox_naive(f::SqrDistL2, x::AbstractArray{T}, gamma::Real=1.0) where T <: RealOrComplex
  p, = prox(f.ind, x)
  sqrd = (f.lambda/2)*norm(x-p)^2
  gamlam = f.lambda*gamma
  return 1/(1+gamlam)*x + gamlam/(1+gamlam)*p, sqrd/(1+gamlam)^2
end


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
using Seep
using Base.Test
using CUDArt

CUDArt.init(0)

tocpu(x::Real) = x
togpu(x::Real) = x
tocpu(x::Array) = SeepNode(x)
togpu(x::Array) = SeepNode(CudaArray(x))
Base.vec(x::Real) = Float64[x]

function cmp(ex)
  @assert ex.head == :call
  op = ex.args[1]
  args = map(eval, ex.args[2:end])
  cpuargs = map(tocpu, args)

  jl = eval(:($op($args...)))
  cpu = eval(:($op($(cpuargs...))))
  gpu = eval(:($op($(map(togpu, args)...))))

  # compare cpu forward implementation to julia
  @assert norm(vec(forward!(cpu, 1))-vec(jl)) < 1e-6

  # compare gpu forward implementation to julia
  @assert norm(vec(to_host(forward!(gpu, 1)))-vec(jl)) < 1e-6

  kill!(cpu)
  kill!(gpu)
end

srand(1)

const N = 5

for op in [:+, :-]
  @show op
  cmp(:($op(randn(N,N), randn(N,N))))
end

for op in [:.+, :.-, :.*, :./]
  @show op
  cmp(:($op(randn(N,N), randn(N,N))))
  cmp(:($op(randn(1,N), randn(N,N))))
  cmp(:($op(randn(N,N), randn(1,N))))
  cmp(:($op(randn(N,1), randn(N,N))))
  cmp(:($op(randn(N,N), randn(N,1))))
  cmp(:($op(randn(N,1), randn(1,N))))
  cmp(:($op(randn(1,N), randn(N,1))))
end

for a in ("", "t", "c")
  for b in ("", "t", "c")
    op = a == b == "" ? :* : symbol("A", a, "_mul_B", b)
    if !isdefined(Base, op) continue end
    @show op
    cmp(:($op(randn(N,N), randn(N,N))))
  end
end

for op in [:exp, :log]
  @show op
  cmp(:($op(rand(N,N))))
  cmp(:($op(rand(1,N))))
  cmp(:($op(rand(N,1))))
end

for op in [:exp, :tanh, :sum, :softmax]
  @show op
  cmp(:($op(randn(N,N))))
  cmp(:($op(randn(1,N))))
  cmp(:($op(randn(N,1))))
end

for op in [:sum]
  @show op
  for i in 1:2
    cmp(:($op(randn(N,N), $i)))
    cmp(:($op(randn(1,N), $i)))
    cmp(:($op(randn(N,1), $i)))
  end
end


text:
include("stratified_couette_flow.jl")

simulate_stratified_couette_flow(Nxy=128, Nz=128, Ri=0)
simulate_stratified_couette_flow(Nxy=128, Nz=128, Ri=0.01)
simulate_stratified_couette_flow(Nxy=128, Nz=128, Ri=0.04)


text:
module Generic

import Base: Array, abs, asin, asinh, atan, atanh, base, bin, checkbounds,
             conj, convert, cmp, contains, cos, cosh, dec, deepcopy,
             deepcopy_internal, denominator, deserialize, det, div, divrem,
             exp, eye, gcd, gcdx, getindex, hash, hcat, hex, intersect, inv,
             invmod, isequal, isfinite, isless, isqrt, isreal, iszero, lcm,
             ldexp, length, log, lufact, lufact!, mod, ndigits, nextpow2, norm,
             nullspace, numerator, oct, one, parent, parse, precision,
             prevpow2, rand, rank, Rational, rem, reverse, serialize,
             setindex!, show, similar, sign, sin, sinh, size, sqrt, string,
             tan, tanh, trace, trailing_zeros, transpose, transpose!, truncate,
             typed_hvcat, typed_hcat, var, vcat, zero, zeros, +, -, *, ==, ^,
             &, |, $, <<, >>, ~, <=, >=, <, >, //, /, !=

if VERSION >= v"0.7.0-DEV.1144"
import Base: isone
end

import Base: floor, ceil, hypot, sqrt, log, log1p, exp, expm1, sin, cos, sinpi,
             cospi, tan, cot, sinh, cosh, tanh, coth, atan, asin, acos, atanh,
             asinh, acosh, gamma, lgamma, sinpi, cospi, atan2

import Nemo: Integers, Rationals, Ring, RingElem, RingElement, Field,
             FieldElement, promote_rule

using Nemo

include("generic/GenericTypes.jl")

include("generic/PermGroups.jl")

include("generic/YoungTabs.jl")

include("generic/Residue.jl")

include("generic/Poly.jl")

include("generic/MPoly.jl")

include("generic/SparsePoly.jl")

include("generic/RelSeries.jl")

include("generic/AbsSeries.jl")

include("generic/LaurentSeries.jl")

include("generic/Matrix.jl")

include("generic/Fraction.jl")

end # generic


text:
using Test
using WaveProp
using WaveProp.Geometry
using WaveProp.IO
using WaveProp.Mesh

@testset "Sphere geo" begin
    # Test the simple sphere geometry
    fname = joinpath(@__DIR__,"sphere.geo")
    Geometry.clear!()
    Ω,M = read_geo(fname)
    @test length(Ω) == 1
    @test geometric_dimension(Ω) == 3

    Γi = internal_boundary(Ω)
    Γe = external_boundary(Ω)

    @test length(Γi) == 0
    @test length(Γe) == 1
    @test geometric_dimension(Γe) == 2

    Ce = external_boundary(Γe)
    Ci = internal_boundary(Γe)
    S  = skeleton(Γe)
    @test length(Ce) == 2
    @test geometric_dimension(Ce) == 1
end

@testset "Sphere msh" begin
    # Test the simple sphere geometry
    Geometry.clear!()
    fname = joinpath(@__DIR__,"sphere.msh")
    Ω, M  = read_msh(fname)
    etypes(M)
    T = SVector{3,Float64}
    @test etypes(M) == [LagrangeLine{2,T},LagrangeTriangle{3,T},LagrangeTetrahedron{4,T},T] # mesh composed of gmsh simplices
    # Test internal creation of sphere
    Ω, M = WaveProp.IO.gmsh_sphere()
    @test etypes(M) == [LagrangeLine{2,T},LagrangeTriangle{3,T},LagrangeTetrahedron{4,T},T] # mesh composed of gmsh simplices
end

@testset "Disk" begin
    # Test internal creation of disk
    Geometry.clear!()
    Ω, M = WaveProp.IO.gmsh_disk()
    T = SVector{2,Float64}
    @test etypes(M) == [LagrangeLine{2,T},LagrangeTriangle{3,T},T] # mesh composed of gmsh simplices
end

@testset "Element iterator" begin
    Geometry.clear!()
    (lx,ly,lz) = widths = (1.,1.,2.)
    Ω, M  = WaveProp.IO.gmsh_box(;widths=widths)
    idx  = 2
    E    = etypes(M)[idx]
    iter = ElementIterator{E}(M)
    @test eltype(iter) == E
    @test length(iter) == size(M.elements[E],2)
end

@testset "Sub mesh" begin
    Geometry.clear!()
    (lx,ly,lz) = widths = (1.,1.,2.)
    Ω, M  = WaveProp.IO.gmsh_box(;widths=widths)
    subM  = SubMesh(M,external_boundary(Ω))
    idx  = 2
    E    = etypes(M)[idx]
    iter = ElementIterator{E}(subM)
    @test eltype(iter) == E
    @test length(iter) == size(M.elements[E],2)
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
#=
  MathOptInterface wrapper.

  The wrapper allows to access the solver through JuMP.
  This implementation is greatly inspired by the wrapper for
  the SDPA solver since both use a similar structure.
  The latter can be found at:
  https://github.com/jump-dev/SDPA.jl/blob/master/src/MOI_wrapper.jl
=#

const AFFEQ = MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}}
const SupportedSets = Union{MOI.Nonnegatives, MOI.PositiveSemidefiniteConeTriangle}

#=
  Structure describing the optimization problem.
  It stores the information about the variables as well
  as some solver's parameters. The numerical data
  for the problem are stored in the 'problem' inner_problem structure
  for convenience.
=#
mutable struct Optimizer <: MOI.AbstractOptimizer
    objconstant::Float64
    objsign::Int
    blockdims::Vector{Int}
    varmap::Vector{Tuple{Int, Int, Int}} # Variable Index vi -> blk, i, j
    b::Vector{Float64}
    problem::Union{Nothing, inner_problem}
    solve_time::Float64
    silent::Bool
    options::Dict{Symbol, Any}
    function Optimizer(; kwargs...)
		optimizer = new(
            zero(Float64), 1, Int[], Tuple{Int, Int, Int}[], Float64[],
            nothing, NaN, false, Dict{Symbol, Any}())
		for (key, value) in kwargs
			MOI.set(optimizer, MOI.RawParameter(key), value)
		end
		return optimizer
    end
end

# Return the name of the solver.
MOI.get(::Optimizer, ::MOI.SolverName) = "MySDPSolver"

# Return the inner problem structure.
MOI.get(opt::Optimizer, ::MOI.RawSolver) = opt.problem

#=
  This part deals with the solver's parameters.
  The workflow consists of writing three functions:
    (1): MOI.supports: boolean function telling MOI if the solver
	                   supports the wanted parameter.
	(2): MOI.get: a getter for the wanted param (if supported).
	(3): MOI.set: a setter (note that some parameters cannot be set).
=#

const SUPPORTED_RAW_PARAMS = [:maxIter, :ε_relgap, :ε_φ, :verbose]

# MOI.RawParameter: a solver-specific attribute.
function MOI.supports(optimizer::Optimizer, param::MOI.RawParameter)
	return param.name in SUPPORTED_RAW_PARAMS
end

function MOI.get(optimizer::Optimizer, param::MOI.RawParameter)
	if haskey(optimizer.options, param.name)
		return optimizer.options[param.name]
	else
		println("ERROR: $(param.name) is not a valid parameter name.")
	end
end

function MOI.set(optimizer::Optimizer, param::MOI.RawParameter, value)
	if !MOI.supports(optimizer, param)
		throw(MOI.UnsupportedAttribute(param))
	end
	optimizer.options[param.name] = value
end

# MOI.Silent: verbosity of the solver. If set to 'true', no info displayed.
MOI.supports(::Optimizer, ::MOI.Silent) = true

function MOI.set(optimizer::Optimizer, ::MOI.Silent, value::Bool)
	optimizer.silent = value
end

MOI.get(optimizer::Optimizer, ::MOI.Silent) = optimizer.silent

# Base.hash(p, u::UInt64) = hash(convert(Int32, p), u)

# const RAW_STATUS = Dict(
#     noINFO     => "The iteration has exceeded the maxIteration and stopped with no informationon the primal feasibility and the dual feasibility.",
#     pdOPT      => "The normal termination yielding both primal and dual approximate optimal solutions.",
#     pFEAS      => "The primal problem got feasible but the iteration has exceeded the maxIteration and stopped.",
#     dFEAS      => "The dual problem got feasible but the iteration has exceeded the maxIteration and stopped.",
#     pdFEAS     => "Both primal problem and the dual problem got feasible, but the iterationhas exceeded the maxIteration and stopped.",
#     pdINF      => "At least one of the primal problem and the dual problem is expected to be infeasible.",
#     pFEAS_dINF => "The primal problem has become feasible but the dual problem is expected to be infeasible.",
#     pINF_dFEAS => "The dual problem has become feasible but the primal problem is expected to be infeasible.",
#     pUNBD      => "The primal problem is expected to be unbounded.",
#     dUNBD => "The dual problem is expected to be unbounded.")

function MOI.get(optimizer::Optimizer, ::MOI.RawStatusString)
    if optimizer.problem === nothing
        return "`MOI.optimize!` not called."
    end
	return RAW_STATUS[getPhaseValue(optimizer.problem)]
end

# MOI.SolveTime: timing of the solver.
function MOI.get(optimizer::Optimizer, ::MOI.SolveTime)
	return optimizer.solve_time
end

#=
  is_empty checks if the optimizer is empty. An empty
  optimizer has no variables, constraints or model
  attributes, however it can have optimizer attributes.
  Indeed, we can create a model with none of the above
  and set some parameters.
  Ex: model = Model(MySDPSolver.Optimizer) # empty optimizer
      MOI.set(model, MOI.Silent(), true)   # still an empty optimizer
=#
function MOI.is_empty(optimizer::Optimizer)
    return iszero(optimizer.objconstant) &&
        optimizer.objsign == 1 &&
        isempty(optimizer.blockdims) &&
        isempty(optimizer.varmap) &&
        isempty(optimizer.b)
end

# empty! effectively empties the optimizer structure.
function MOI.empty!(optimizer::Optimizer)
    optimizer.objconstant = zero(Float64)
    optimizer.objsign = 1
    empty!(optimizer.blockdims) # sets a Vector to [] in Julia
    empty!(optimizer.varmap)
    empty!(optimizer.b)
    optimizer.problem = nothing
end

function MOI.supports(
    optimizer::Optimizer,
    ::Union{MOI.ObjectiveSense,
            MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}})
    return true
end

MOI.supports_add_constrained_variables(::Optimizer, ::Type{MOI.Reals}) = false

MOI.supports_add_constrained_variables(::Optimizer, ::Type{<:SupportedSets}) = true
function MOI.supports_constraint(
    ::Optimizer, ::Type{MOI.ScalarAffineFunction{Float64}},
    ::Type{MOI.EqualTo{Float64}})
    return true
end

function MOI.copy_to(dest::Optimizer, src::MOI.ModelLike; kws...)
    return MOIU.automatic_copy_to(dest, src; kws...)
end

MOIU.supports_allocate_load(::Optimizer, copy_names::Bool) = !copy_names

function MOIU.allocate(optimizer::Optimizer, ::MOI.ObjectiveSense, sense::MOI.OptimizationSense)
    # To be sure that it is done before load(optimizer, ::ObjectiveFunction, ...), we do it in allocate
    optimizer.objsign = sense == MOI.MIN_SENSE ? -1 : 1
end

function MOIU.allocate(::Optimizer, ::MOI.ObjectiveFunction, ::MOI.ScalarAffineFunction) end

function MOIU.load(::Optimizer, ::MOI.ObjectiveSense, ::MOI.OptimizationSense) end

#=
  Utility function returning the block as well as the indices
  within that block (blk, i, j) of a variable indexed by
  vi by MOI.
=#
varmap(optimizer::Optimizer, vi::MOI.VariableIndex) = optimizer.varmap[vi.value]

# Loads objective coefficient α * vi
function load_objective_term!(optimizer::Optimizer, α, vi::MOI.VariableIndex)
    blk, i, j = varmap(optimizer, vi)
    coef = optimizer.objsign * α
    if i != j
        coef /= sqrt(2.0)
    end
	# println("blk, i, j, coef: ", blk,", ", i,", ",j,", ", coef)

	# push!(optimizer.problem.Data.I_obj, i)
	# push!(optimizer.problem.Data.J_obj, j)
	# push!(optimizer.problem.Data.V_obj, coef)
	# inputElement(optimizer.problem, 0, blk, i, j, float(coef), false)
    # in SDP format, it is max and in MPB Conic format it is min
    # inputElement(optimizer.problem, 0, blk, i, j, float(coef), false)
end

function MOIU.load(optimizer::Optimizer, ::MOI.ObjectiveFunction, f::MOI.ScalarAffineFunction)
	# f is of type MOI.ScalarAffineFunction. We put it in the canonical
	# form, which is represented as a triplet (x, terms, constant) where
	# x is a vector of MOI.VariableIndex, coeff are the coefficient and constant
	# is the independant term.
    obj = MOIU.canonical(f)
    optimizer.objconstant = f.constant # extract the constant

	# Preallocate the vectors for the temp. variables.
	I_s = [Int64[] for i in 1:length(optimizer.blockdims)]   # row indices
	J_s = [Int64[] for i in 1:length(optimizer.blockdims)]   # col indices
	V_s = [Float64[] for i in 1:length(optimizer.blockdims)] # val indices

    for t in obj.terms
        if !iszero(t.coefficient)
			blk, i, j = varmap(optimizer, t.variable_index)
            coef = t.coefficient
			# The canonical repr. sums the same variables,
			# ex. for symmetric matrix var., X[1,2] == X[2,1],
			# and since only one is kept, the coefficient will be doubled.
            if i != j
                coef /= 2
            end
			# Push the element to the corresponding entry
			# in the corresponding block.
			push!(I_s[blk], i)
			push!(J_s[blk], j)
			push!(V_s[blk], coef)
        end
    end
	# Vectorize the matrices and update the inner problem in the Optimizer.
	C_temp = [sparse(I_s[k], J_s[k], V_s[k], abs(optimizer.blockdims[k]), abs(optimizer.blockdims[k])) for k in 1:length(optimizer.blockdims)]
	row_temp = [(optimizer.blockdims[k] > 0) ? sparsevec(svec_lower(C_temp[k])) : sparsevec(diag(C_temp[k])) for k in 1:length(optimizer.blockdims)]
	optimizer.problem.normC = max(1.0, maximum(norm.(row_temp)))
	optimizer.problem.C = -optimizer.objsign*vcat(row_temp...)
	I_s = nothing # Help garbage-collect these
	J_s = nothing
	V_s = nothing
	C_temp = nothing
	row_temp = nothing
end

# Add new blocks to the varmap. If the block is diagonally-constrained,
# the size is negative.
function new_block(optimizer::Optimizer, set::MOI.Nonnegatives)
    push!(optimizer.blockdims, -MOI.dimension(set)) # dimension(set) = nb of vars.
    blk = length(optimizer.blockdims)
    for i in 1:MOI.dimension(set)
        push!(optimizer.varmap, (blk, i, i))
    end
end

function new_block(optimizer::Optimizer, set::MOI.PositiveSemidefiniteConeTriangle)
    push!(optimizer.blockdims, set.side_dimension) # dimension(set) = nxn
    blk = length(optimizer.blockdims)
    for i in 1:set.side_dimension
        for j in 1:i
            push!(optimizer.varmap, (blk, i, j))
        end
    end
end

function MOIU.allocate_constrained_variables(optimizer::Optimizer,
                                             set::SupportedSets)
    offset = length(optimizer.varmap)
    new_block(optimizer, set)
    ci = MOI.ConstraintIndex{MOI.VectorOfVariables, typeof(set)}(offset + 1)
    return [MOI.VariableIndex(i) for i in offset .+ (1:MOI.dimension(set))], ci
end

function MOIU.load_constrained_variables(
    optimizer::Optimizer, vis::Vector{MOI.VariableIndex},
    ci::MOI.ConstraintIndex{MOI.VectorOfVariables},
    set::SupportedSets)
end

function MOIU.load_variables(optimizer::Optimizer, nvars)
    @assert nvars == length(optimizer.varmap)
	N = sum((optimizer.blockdims[k] > 0) ? (optimizer.blockdims[k])*(optimizer.blockdims[k]+1)÷2 : abs(optimizer.blockdims[k]) for k in 1:length(optimizer.blockdims))
    dummy = isempty(optimizer.b)
    if dummy
        optimizer.b = [one(Float64)]
        optimizer.blockdims = [optimizer.blockdims; -1]
    end
	# Instantiate inner_problem with correct structure and data set to 0.
    optimizer.problem = inner_problem(optimizer.blockdims,
	                                  spzeros(Float64, N, length(optimizer.b)),
									  -optimizer.b,
									  Vector{Float64}(undef, N))
	optimizer.problem.normsA = zeros(Float64, length(optimizer.blockdims))

end

function MOIU.allocate_constraint(optimizer::Optimizer,
                                  func::MOI.ScalarAffineFunction{Float64},
                                  set::MOI.EqualTo{Float64})
    push!(optimizer.b, MOI.constant(set))
    return AFFEQ(length(optimizer.b))
end

function MOIU.load_constraint(m::Optimizer, ci::AFFEQ,
                              f::MOI.ScalarAffineFunction, s::MOI.EqualTo)
	# See MOIU.load for comments.
    f = MOIU.canonical(f)

	I_s = [Int64[] for i in 1:length(m.blockdims)]
	J_s = [Int64[] for i in 1:length(m.blockdims)]
	V_s = [Float64[] for i in 1:length(m.blockdims)]

    for t in f.terms
        if !iszero(t.coefficient)
            blk, i, j = varmap(m, t.variable_index)
            coef = t.coefficient
            if i != j
                coef /= 2
            end
			push!(I_s[blk], i)
			push!(J_s[blk], j)
			push!(V_s[blk], -coef)
        end
    end

	Ablocks_i = [sparse(I_s[k], J_s[k], V_s[k], abs(m.blockdims[k]), abs(m.blockdims[k])) for k in 1:length(m.blockdims)]
	temp_row = [(m.blockdims[k] > 0) ? sparsevec(svec_lower(Ablocks_i[k])) : sparsevec(diag(Ablocks_i[k])) for k in 1:length(m.blockdims)]
	m.problem.normsA += (norm.(temp_row)).^2
	m.problem.A[:, ci.value] = vcat(temp_row...)
	I_s = nothing
	J_s = nothing
	V_s = nothing
	Ablocks_i = nothing
	temp_row = nothing
end

function MOI.optimize!(m::Optimizer)
    # Update the norms for scaling
	map!(x -> max(1.0, sqrt(sqrt(x))), m.problem.normsA, m.problem.normsA) # NormsA was generated incrementally
	m.problem.normb = max(1.0, norm(m.problem.b))

    # Generate initial iterate uses unscaled data. We scale them later.
	generate_init!(m.problem)
	scale_inner_problem!(m.problem)

    # Set some params.
    m.problem.verbose = !m.silent
	for (key, value) in m.options
		# MOI.set(optimizer, MOI.RawParameter(key), value)
		if key == :maxIter
			m.problem.maxIter = value
		elseif key == :ε_relgap
			m.problem.ε_relgap = value
		elseif key == :ε_φ
			m.problem.ε_φ = value
		elseif key == :verbose
			m.problem.verbose = value
		end
	end
    # Print some info
    if m.problem.verbose == true
		println("======================================")
        println("========    MySDPSolver    ===========")
		println("======================================\n")
    end
	start_time = time()

	solve!(m.problem)

    m.solve_time = time() - start_time
end

function MOI.supports(::Optimizer, ::MOI.TerminationStatus)
	return true
end

function MOI.supports(::Optimizer, ::MOI.PrimalStatus)
	return true
end

function MOI.supports(::Optimizer, ::MOI.DualStatus)
	return true
end

function MOI.get(m::Optimizer, ::MOI.BarrierIterations)
	return m.problem.num_iters
end

function MOI.get(m::Optimizer, ::MOI.TerminationStatus)
	return MOI.OPTIMAL
end

function MOI.get(m::Optimizer, attr::MOI.PrimalStatus)
	return MOI.FEASIBLE_POINT
end
function MOI.get(m::Optimizer, attr::MOI.DualStatus)
	return MOI.FEASIBLE_POINT
end

MOI.get(m::Optimizer, ::MOI.ResultCount) = m.problem === nothing ? 0 : 1

function MOI.get(m::Optimizer, attr::MOI.ObjectiveValue)
	return (-m.objsign)*dot(m.problem.C, m.problem.X) + m.objconstant
    # return m.objsign * dot(m.problem.C, m.problem.X) + m.objconstant
end

function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal)

end

function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex)
    blk, i, j = varmap(optimizer, vi)
    scaling = (i==j) ? 1.0 : 1.0/sqrt(2.0)
	index = sum([(elem > 0) ? elem*(elem+1)÷2 : abs(elem) for elem in optimizer.blockdims[1:blk-1]])
	if optimizer.blockdims[blk] > 0
		if j > i
			index += sum(1:j) - (j-i)
		elseif i == j
			index += sum(1:j)
		else
			index += sum(1:i) - (i-j)
		end
	else
		index += abs(i)
	end
	return scaling*optimizer.problem.X[index]
end

function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S<:SupportedSets
	blk, i, j = optimizer.varmap[ci.value]
	scaling = (i==j) ? 1.0 : 1.0/sqrt(2.0)
    start_index = Int(sum([(elem > 0) ? elem*(elem+1)÷2 : abs(elem) for elem in optimizer.blockdims[1:blk-1]])) + 1
	length_block = (optimizer.blockdims[blk] > 0) ? optimizer.blockdims[blk]*(optimizer.blockdims[blk]+1)÷2 : abs(optimizer.blockdims[blk])
	temp = optimizer.problem.S[start_index:start_index+length_block-1]
	idx = 1
	for i = 1:optimizer.blockdims[blk]
		for j = 1:i
			if i != j
				temp[idx] /= sqrt(2.0)
			end
			idx += 1
		end
	end
	return temp
end

function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintPrimal,
	             ci::MOI.ConstraintIndex{MOI.VectorOfVariables, S}) where S<:SupportedSets
	return 1.0
end

function MOI.get(optimizer::Optimizer, attr::MOI.ConstraintDual, ci::AFFEQ)
    # MOI.check_result_index_bounds(optimizer, attr)
    return -optimizer.problem.y[ci.value]
end

max_stars_repo_name: TomaszWasTaken/MySDPSolver

max_issues_repo_name: TomaszWasTaken/MySDPSolver

max_forks_repo_name: TomaszWasTaken/MySDPSolver


text:
using FirstOrderLogic
using Base.Test

include("test_parser.jl")
include("test_renamed_quantifiers_form.jl")
include("test_quantified_variables_form.jl")
include("test_prenex_normal_form.jl")
include("test_skolem_normal_form.jl")
include("test_removed_quantifier_form.jl")
include("test_conjunctive_normal_form.jl")
include("test_resolution.jl")


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""A label selector requirement is a selector that contains values, a key, and an operator that relates the key and values.

    IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement(;
        key=nothing,
        operator=nothing,
        values=nothing,
    )

    - key::String : key is the label key that the selector applies to.
    - operator::String : operator represents a key's relationship to a set of values. Valid operators are In, NotIn, Exists and DoesNotExist.
    - values::Vector{String} : values is an array of string values. If the operator is In or NotIn, the values array must be non-empty. If the operator is Exists or DoesNotExist, the values array must be empty. This array is replaced during a strategic merge patch.
"""
mutable struct IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement <: SwaggerModel
    key::Any # spec type: Union{ Nothing, String } # spec name: key
    operator::Any # spec type: Union{ Nothing, String } # spec name: operator
    values::Any # spec type: Union{ Nothing, Vector{String} } # spec name: values

    function IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement(;key=nothing, operator=nothing, values=nothing)
        o = new()
        validate_property(IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement, Symbol("key"), key)
        setfield!(o, Symbol("key"), key)
        validate_property(IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement, Symbol("operator"), operator)
        setfield!(o, Symbol("operator"), operator)
        validate_property(IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement, Symbol("values"), values)
        setfield!(o, Symbol("values"), values)
        o
    end
end # type IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement

const _property_map_IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement = Dict{Symbol,Symbol}(Symbol("key")=>Symbol("key"), Symbol("operator")=>Symbol("operator"), Symbol("values")=>Symbol("values"))
const _property_types_IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement = Dict{Symbol,String}(Symbol("key")=>"String", Symbol("operator")=>"String", Symbol("values")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement }) = collect(keys(_property_map_IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement))
Swagger.property_type(::Type{ IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement[name]))}
Swagger.field_name(::Type{ IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement }, property_name::Symbol) =  _property_map_IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement[property_name]

function check_required(o::IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement)
    (getproperty(o, Symbol("key")) === nothing) && (return false)
    (getproperty(o, Symbol("operator")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApimachineryPkgApisMetaV1LabelSelectorRequirement }, name::Symbol, val)
end


text:
using MPI
using Oceananigans.Grids: topology, size, halo_size, architecture, pop_flat_elements
using Oceananigans.Grids: validate_rectilinear_domain, validate_topology
using Oceananigans.Grids: validate_rectilinear_grid_args, validate_lat_lon_grid_args
using Oceananigans.Grids: generate_coordinate, with_precomputed_metrics
using Oceananigans.Grids: cpu_face_constructor_x, cpu_face_constructor_y, cpu_face_constructor_z

import Oceananigans.Grids: RectilinearGrid, LatitudeLongitudeGrid, with_halo

const DistributedGrid{FT, TX, TY, TZ} = AbstractGrid{FT, TX, TY, TZ, <:MultiArch}

@inline get_local_coords(c::Tuple         , nc, R, index) = (c[1] + (index-1) * (c[2] - c[1]) / R,    c[1] + index * (c[2] - c[1]) / R)
@inline get_local_coords(c::AbstractVector, nc, R, index) = c[1 + (index-1) * nc : 1 + nc * index]

@inline get_global_coords(c::Tuple        , nc, R,  index, arch) = (c[2] - index * (c[2] - c[1]), c[2] - (index - R) * (c[2] - c[1]))

"""
    get_global_coords(c_local::AbstractVector, nc, R, index, arch)

Build a linear global coordinate vector given a local coordinate vector `c_local`
a local number of elements `nc`, number of ranks `R`, rank `index`,
and `arch`itecture.
"""
function get_global_coords(c_local::AbstractVector, nc, R, index, arch)
    c_global = zeros(eltype(c_local), nc*R+1)
    c_global[1 + (index-1) * nc : nc * index] .= c[1:end-1]
    index == R && (c_global[end] = c[end])

    MPI.Allreduce!(c_global, +, arch.communicator)

    return c_global
end

"""
    RectilinearGrid(arch::MultiArch, FT=Float64; kw...)

Return the rank-local portion of `RectilinearGrid` on `arch`itecture.
"""
function RectilinearGrid(arch::MultiArch, FT = Float64;
                         size,
                         x = nothing,
                         y = nothing,
                         z = nothing,
                         halo = nothing,
                         extent = nothing,
                         topology = (Periodic, Periodic, Bounded))

    TX, TY, TZ, size, halo, x, y, z =
        validate_rectilinear_grid_args(topology, size, halo, FT, extent, x, y, z)

    Nx, Ny, Nz = size
    Hx, Hy, Hz = halo

    ri, rj, rk = arch.local_index
    Rx, Ry, Rz = arch.ranks

    # Make sure we can put an integer number of grid points in each rank.
    # Will generalize in the future.
    @assert isinteger(Nx / Rx)
    @assert isinteger(Ny / Ry)
    @assert isinteger(Nz / Rz)

    # Local sizes are denoted with lowercase `n`
    nx, ny, nz = local_size = Nx÷Rx, Ny÷Ry, Nz÷Rz

    xl = get_local_coords(x, nx, Rx, ri)
    yl = get_local_coords(y, ny, Ry, rj)
    zl = get_local_coords(z, nz, Rz, rk)

    Lx, xᶠᵃᵃ, xᶜᵃᵃ, Δxᶠᵃᵃ, Δxᶜᵃᵃ = generate_coordinate(FT, topology[1], nx, Hx, xl, child_architecture(arch))
    Ly, yᵃᶠᵃ, yᵃᶜᵃ, Δyᵃᶠᵃ, Δyᵃᶜᵃ = generate_coordinate(FT, topology[2], ny, Hy, yl, child_architecture(arch))
    Lz, zᵃᵃᶠ, zᵃᵃᶜ, Δzᵃᵃᶠ, Δzᵃᵃᶜ = generate_coordinate(FT, topology[3], nz, Hz, zl, child_architecture(arch))

    architecture = MultiArch(child_architecture(arch), topology = topology, ranks = arch.ranks, communicator = arch.communicator)

    return RectilinearGrid{TX, TY, TZ}(architecture,
                                       nx, ny, nz,
                                       Hx, Hy, Hz,
                                       Lx, Ly, Lz,
                                       Δxᶠᵃᵃ, Δxᶜᵃᵃ, xᶠᵃᵃ, xᶜᵃᵃ,
                                       Δyᵃᶜᵃ, Δyᵃᶠᵃ, yᵃᶠᵃ, yᵃᶜᵃ,
                                       Δzᵃᵃᶠ, Δzᵃᵃᶜ, zᵃᵃᶠ, zᵃᵃᶜ)
end

function LatitudeLongitudeGrid(arch::MultiArch,
                               FT = Float64;
                               precompute_metrics = false,
                               size,
                               latitude,
                               longitude,
                               z,
                               radius = R_Earth,
                               halo = (1, 1, 1))

    Nλ, Nφ, Nz, Hλ, Hφ, Hz, latitude, longitude, topo =
        validate_lat_lon_grid_args(latitude, longitude, size, halo)

    i, j, k = arch.local_index
    Rx, Ry, Rz = arch.ranks

    # Make sure we can put an integer number of grid points in each rank.
    # Will generalize in the future.
    @assert isinteger(Nλ / Rx)
    @assert isinteger(Nφ / Ry)
    @assert isinteger(Nz / Rz)

    nλ, nφ, nz = local_size = Nλ÷Rx, Nφ÷Ry, Nz÷Rz

    λl = get_local_coords(longitude, nx, Rx, i)
    φl = get_local_coords(latitude , ny, Ry, j)
    zl = get_local_coords(z,         nz, Rz, k)

    # Calculate all direction (which might be stretched)
    # A direction is regular if the domain passed is a Tuple{<:Real, <:Real},
    # it is stretched if being passed is a function or vector (as for the VerticallyStretchedRectilinearGrid)
    Lλ, λᶠᵃᵃ, λᶜᵃᵃ, Δλᶠᵃᵃ, Δλᶜᵃᵃ = generate_coordinate(FT, topo[1], nλ, Hλ, λl, arch.child_architecture)
    Lφ, φᵃᶠᵃ, φᵃᶜᵃ, Δφᵃᶠᵃ, Δφᵃᶜᵃ = generate_coordinate(FT, topo[2], nφ, Hφ, φl, arch.child_architecture)
    Lz, zᵃᵃᶠ, zᵃᵃᶜ, Δzᵃᵃᶠ, Δzᵃᵃᶜ = generate_coordinate(FT, topo[3], nz, Hz, zl, arch.child_architecture)

    architecture = MultiArch(child_architecture(arch), grid = grid, ranks = arch.ranks, communicator = arch.communicator)

    preliminary_grid = LatitudeLongitudeGrid(architecture,
                                             Nλ, Nφ, Nz,
                                             Hλ, Hφ, Hz,
                                             Lλ, Lφ, Lz,
                                             Δλᶠᵃᵃ, Δλᶜᵃᵃ, λᶠᵃᵃ, λᶜᵃᵃ,
                                             Δφᵃᶠᵃ, Δφᵃᶜᵃ, φᵃᶠᵃ, φᵃᶜᵃ,
                                             Δzᵃᵃᶠ, Δzᵃᵃᶜ, zᵃᵃᶠ, zᵃᵃᶜ,
                                             (nothing for i=1:10)..., radius)

    return !precompute_metrics ? preliminary_grid : with_precomputed_metrics(preliminary_grid)
end

function reconstruct_global_grid(grid::RectilinearGrid)

    arch = grid.architecture
    i, j, k = arch.local_index

    Rx, Ry, Rz = R = arch.ranks

    nx, ny, nz = n = size(grid)
    Hx, Hy, Hz = H = halo_size(grid)
    Nx, Ny, Nz = n .* R

    TX, TY, TZ = topology(grid)

    x = cpu_face_constructor_x(grid)
    y = cpu_face_constructor_y(grid)
    z = cpu_face_constructor_z(grid)

    xG = get_global_coords(x, nx, Rx, i, arch)
    yG = get_global_coords(y, ny, Ry, j, arch)
    zG = get_global_coords(z, nz, Rz, k, arch)

    child_arch = child_architecture(arch)

    FT = eltype(grid)

    Lx, xᶠᵃᵃ, xᶜᵃᵃ, Δxᶠᵃᵃ, Δxᶜᵃᵃ = generate_coordinate(FT, TX, Nx, Hx, xG, child_arch)
    Ly, yᵃᶠᵃ, yᵃᶜᵃ, Δyᵃᶠᵃ, Δyᵃᶜᵃ = generate_coordinate(FT, TY, Ny, Hy, yG, child_arch)
    Lz, zᵃᵃᶠ, zᵃᵃᶜ, Δzᵃᵃᶠ, Δzᵃᵃᶜ = generate_coordinate(FT, TZ, Nz, Hz, zG, child_arch)

    return RectilinearGrid{TX, TY, TZ}(child_arch,
                                       Nx, Ny, Nz,
                                       Hx, Hy, Hz,
                                       Lx, Ly, Lz,
                                       Δxᶠᵃᵃ, Δxᶜᵃᵃ, xᶠᵃᵃ, xᶜᵃᵃ,
                                       Δyᵃᶠᵃ, Δyᵃᶜᵃ, yᵃᶠᵃ, yᵃᶜᵃ,
                                       Δzᵃᵃᶠ, Δzᵃᵃᶜ, zᵃᵃᶠ, zᵃᵃᶜ)
end

function with_halo(new_halo, grid::DistributedGrid)
    new_grid = with_halo(new_halo, reconstruct_global_grid(grid))
    return scatter_local_grids(architecture(grid), new_grid)
end

function scatter_grid_properties(global_grid)
    # Pull out face grid constructors
    x = cpu_face_constructor_x(global_grid)
    y = cpu_face_constructor_y(global_grid)
    z = cpu_face_constructor_z(global_grid)

    topo = topology(grid)
    sz   = pop_flat_elements(size(global_grid), topo)
    halo = pop_flat_elements(halo_size(global_grid), topo)

    return x, y, z, topo, sz, halo
end

function scatter_local_grids(arch::MultiArch, global_grid::RectilinearGrid)
    x, y, z, topo, sz, halo = scatter_grid_properties(global_grid)
    return RectilinearGrid(arch, eltype(grid); size=sz, x=x, y=y, z=z, halo=halo, topology=topo)
end

function scatter_local_grids(arch::MultiArch, global_grid::LatitudeLongitudeGrid)
    x, y, z, topo, sz, halo = scatter_grid_properties(global_grid)
    return LatitudeLongitudeGrid(arch, eltype(global_grid); size=sz, longitude=x, latitude=y, z=z, halo=halo)
end


text:
#
# Factors IO
#
# Printing IO stuff

# Base.show(io::IO, ϕ::Factor) = show(io, convert(DataFrame, ϕ))

function Base.show(io::IO, ϕ::Factor)
    print(io, "$(length(ϕ)) instantiations:")
    for (d, s) in zip(ϕ.dimensions, size(ϕ))
        println(io, "")
        print(io, "  $d ($s)")
    end
end

# # mimewriting doesn't use a monospace font, which I don't like ...
Base.mimewritable(::MIME"text/html", ϕ::Factor) = true
Base.show(io::IO, a::MIME"text/html", ϕ::Factor) = show(io, a, convert(DataFrame, ϕ))
#    print(io, replace(replace(repr(ϕ), "\n", "<br>"), "\t",
#            "   "))


text:
using Test

using DiffEqBase, ForwardDiff
using DiffEqBase: prob2dtmin, timedepentdtmin

@testset "tspan2dtmin" begin
  tspan2dtmin(tspan; kwargs...) = prob2dtmin(ODEProblem(identity, 1, tspan); kwargs...)
  @test tspan2dtmin((10, 100.0)) === eps(100.0)
  @test tspan2dtmin((-10000.0, 100.0)) === eps(10000.0)
  @test tspan2dtmin((1, 2)) === 0
  @test tspan2dtmin((1//10, 2//10)) === 1//10^10
  @test tspan2dtmin((2//10, Inf)) === eps(1.0)
  @test tspan2dtmin((2//1, Inf)) === eps(2.0)
  @test tspan2dtmin((0, Inf)) === tspan2dtmin((0.0, Inf)) === eps(1.0)
  @test_throws ArgumentError tspan2dtmin((Inf, 100.0))
  @test tspan2dtmin((0f0, 1f5); use_end_time=false) === eps(1f0)
  @test timedepentdtmin(10f0, eps(1f0)) === eps(10f0)
  @test timedepentdtmin(10, eps(1f0)) === eps(1f0)
end

@testset "prob2dtmin" begin
  @test prob2dtmin((0.0, 10.0), 1.0, false) == eps(Float64)
  @test prob2dtmin((0f0, 10f0), 1f0, false) == eps(Float32)
  @test prob2dtmin((0.0, 10.0), ForwardDiff.Dual(1.0), false) == eps(Float64)
end


text:
# seg_sites.jl
# =============
#
# Identify and count segregating sites with BioJulia data types.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/GeneticVariation.jl/blob/master/LICENSE.md

"""
`Segregating` sites are positions which show differences (polymorphisms)
between related genes in a sequence alignment (are not conserved).
Segregating sites include conservative, semi-conservative and non-conservative
mutations.
"""
struct Segregating <: Position end

function Base.count(::Type{Segregating}, seqs::Vector{T}) where T <: Sequence
    cont = true
    i = 0
    count = 0
    @inbounds refseq = seqs[1]
    while cont
        i += 1
        @inbounds refelem = refseq[i]
        @inbounds for j in 2:lastindex(seqs)
            seq = seqs[j]
            cont = ifelse(lastindex(seq) == i, false, true)
            if refelem != seq[i]
                count += 1
                break
            end
        end
    end
    return count, i
end


text:
############################################################################################
### Types
############################################################################################
"""
Connects to `NLopt.jl` as the optimization backend.

# Constructor

    SemOptimizerNLopt(;
        algorithm = :LD_LBFGS,
        options = Dict{Symbol, Any}(),
        local_algorithm = nothing,
        local_options = Dict{Symbol, Any}(),
        equality_constraints = Vector{NLoptConstraint}(),
        inequality_constraints = Vector{NLoptConstraint}(),
        kwargs...)

# Arguments
- `algorithm`: optimization algorithm.
- `options::Dict{Symbol, Any}`: options for the optimization algorithm
- `local_algorithm`: local optimization algorithm
- `local_options::Dict{Symbol, Any}`: options for the local optimization algorithm
- `equality_constraints::Vector{NLoptConstraint}`: vector of equality constraints
- `inequality_constraints::Vector{NLoptConstraint}`: vector of inequality constraints

# Example
```julia
my_optimizer = SemOptimizerNLopt()

# constrained optimization with augmented lagrangian
my_constrained_optimizer = SemOptimizerNLopt(;
    algorithm = :AUGLAG,
    local_algorithm = :LD_LBFGS,
    local_options = Dict(:ftol_rel => 1e-6),
    inequality_constraints = NLoptConstraint(;f = my_constraint, tol = 0.0),
)
```

# Usage
All algorithms and options from the NLopt library are available, for more information see
the NLopt.jl package and the NLopt online documentation.
For information on how to use inequality and equality constraints, see XXX in our online documentation.

# Extended help

## Interfaces
- `algorithm(::SemOptimizerNLopt)`
- `local_algorithm(::SemOptimizerNLopt)`
- `options(::SemOptimizerNLopt)`
- `local_options(::SemOptimizerNLopt)`
- `equality_constraints(::SemOptimizerNLopt)`
- `inequality_constraints(::SemOptimizerNLopt)`

## Implementation

Subtype of `SemOptimizer`.
"""
struct SemOptimizerNLopt{A, A2, B, B2, C} <: SemOptimizer
    algorithm::A
    local_algorithm::A2
    options::B
    local_options::B2
    equality_constraints::C
    inequality_constraints::C
end

Base.@kwdef mutable struct NLoptConstraint
    f
    tol = 0.0
end

############################################################################################
### Constructor
############################################################################################

function SemOptimizerNLopt(;
        algorithm = :LD_LBFGS,
        local_algorithm = nothing,
        options = Dict{Symbol, Any}(),
        local_options = Dict{Symbol, Any}(),
        equality_constraints = Vector{NLoptConstraint}(),
        inequality_constraints = Vector{NLoptConstraint}(),
        kwargs...)
    applicable(iterate, equality_constraints) ||
        (equality_constraints = [equality_constraints])
    applicable(iterate, inequality_constraints) ||
        (inequality_constraints = [inequality_constraints])
    return SemOptimizerNLopt(
        algorithm,
        local_algorithm,
        options,
        local_options,
        equality_constraints,
        inequality_constraints)
end

############################################################################################
### Recommended methods
############################################################################################

update_observed(optimizer::SemOptimizerNLopt, observed::SemObserved; kwargs...) = optimizer

############################################################################################
### additional methods
############################################################################################

algorithm(optimizer::SemOptimizerNLopt) = optimizer.algorithm
local_algorithm(optimizer::SemOptimizerNLopt) = optimizer.local_algorithm
options(optimizer::SemOptimizerNLopt) = optimizer.options
local_options(optimizer::SemOptimizerNLopt) = optimizer.local_options
equality_constraints(optimizer::SemOptimizerNLopt) = optimizer.equality_constraints
inequality_constraints(optimizer::SemOptimizerNLopt) = optimizer.inequality_constraints


text:
using PyPlot

b_ref = DataFrame(bus = collect(1:15), v_squared = [1.000, 0.942,  0.964,  1.000,
             0.997,  0.994,  0.992,  1.041, 1.021,  1.023,
             1.031,  1.034, 0.959,  0.950,  0.944])

function compare_voltages(res_dict1, res_dict2; label1="", label2="" )

u1 = res_dict1[:v_squared]
u2 = res_dict2[:v_squared]
buses = res_dict1[:bus]

v1 = sqrt.(u1)
v2 = sqrt.(u2)

if length(v1) != length(v1)
    warn("Plotter: length of input vectors not the same")
end

clf()
plot(buses, v1, "ro-", label=label1)
plot(buses, v2, "bo-", label=label2)
legend()

end

function compare_voltages(res_dict1, res_dict2, res_dict3; label1="", label2="", label3="" )

u1 = res_dict1[:v_squared]
u2 = res_dict2[:v_squared]
u3 = res_dict3[:v_squared]
buses = res_dict1[:bus]

v1 = sqrt.(u1)
v2 = sqrt.(u2)
v3 = sqrt.(u3)

if length(v1) != length(v1)
    warn("Plotter: length of input vectors not the same")
end

clf()
plot(buses, v1, "ro-", label=label1)
plot(buses, v2, "bo-", label=label2)
plot(buses, v3, "ko-", label=label3)
legend()

end


text:
#################### Binary Metropolised Gibbs Sampler ####################

#################### Types and Constructors ####################

typealias BMGForm Union{Int, Vector{Vector{Int}}}

type BMGTune{F<:BMGForm} <: SamplerTune
  logf::Nullable{Function}
  k::F

  BMGTune() = new()

  BMGTune(x::Vector, k::F) = new(Nullable{Function}(), k)

  BMGTune(x::Vector, k::F, logf::Function) = new(Nullable{Function}(logf), k)
end

typealias BMGIntVariate SamplerVariate{BMGTune{Int}}
typealias BMGVecVariate SamplerVariate{BMGTune{Vector{Vector{Int}}}}

BMGVariate{F<:BMGForm}(x::Vector, logf::Function; k::F=1) =
  SamplerVariate{BMGTune{F}}(x, k, logf)

function validate(v::BMGIntVariate)
  n = length(v)
  v.tune.k <= n || throw(ArgumentError("k exceeds variate length $n"))
  validatebinary(v)
end

function validate(v::BMGVecVariate)
  n = length(v)
  mapreduce(maximum, max, v.tune.k) <= n ||
    throw(ArgumentError("indices in k exceed variate length $n"))
  validatebinary(v)
end

#################### Sampler Constructor ####################

function BMG{F<:BMGForm}(params::ElementOrVector{Symbol}; k::F=1)
  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block)
    v = SamplerVariate(block, k)
    sample!(v, x -> logpdf!(block, x))
    relist(block, v)
  end
  Sampler(params, samplerfx, BMGTune{F}())
end

#################### Sampling Functions ####################

sample!{F<:BMGForm}(v::SamplerVariate{BMGTune{F}}) = sample!(v, v.tune.logf)

function sample!{F<:BMGForm}(v::SamplerVariate{BMGTune{F}}, logf::Function)
  n = length(v)
  probs = Vector{Float64}(n)
  idx = randind(v)

  pbernoulli! = function(x, probs)
    for i in idx
      x_i = x[i]

      x[i] = 0.0
      logf0 = logf(x)
      x[i] = 1.0
      logf1 = logf(x)

      p = invlogit(logf1 - logf0)
      probs[i] = 0.0 < p < 1.0 ? p : 0.5

      x[i] = x_i
    end
    probs
  end

  x = v[:]
  pbernoulli!(x, probs)
  theta = map(p -> rand() < p, probs[idx])

  if n == 1
    v[idx] = theta
  else
    y = v[:]
    y[idx] = theta
    qy = mapreduce(i -> y[i] == 1.0 ? log(probs[i]) : log(1.0 - probs[i]), +,
                   idx)

    pbernoulli!(y, probs)
    qx = mapreduce(i -> x[i] == 1.0 ? log(probs[i]) : log(1.0 - probs[i]), +,
                   idx)

    if rand() < exp((logf(y) - qy) - (logf(x) - qx))
      v[idx] = theta
    end
  end

  v
end

randind(v::BMGIntVariate) = sample(1:length(v), v.tune.k, replace=false)
randind(v::BMGVecVariate) = sample(v.tune.k)


text:
##################################################
# Raw apply/combine
##################################################

unwrap(x) = x
function unwrap(x::Sym{TAG}; context::Context = @__context__) where {TAG}
    hashead(x, :object) && return firstarg(x)
    if hashead(x, :call)
        op = firstarg(x)
        if hashead(op, :fn)
            xs = tailargs(x)
            f = firstarg(op)
            tags = map(tagof, xs)
            if all(tag -> tag isa Type, tags) && hasmethod(f, Tuple{tags...})
                if all(hashead(:object), xs) && all(isfaithful, xs)
                    return invoke(f, Tuple{tags...}, map(firstarg, xs)...)
                elseif all(ismocking, xs) && all(isabstracttype, tags)
                    return mock(f, xs...; context = context)
                end
            end
        end
    end
    return x
end

_apply(args...) = _apply(map(arg -> convert(Sym, arg), args)...)
function _apply(op::Sym, xs::Sym...)
    tags = map(tagof, xs)
    TAG = promote_tag(:call, op, tags...)
    x = Sym{TAG}(:call, op, xs...)
    return x::Sym{TAG}
end

function _combine(head::Symbol, xs...)
    TAG = promote_tag(head, map(tagof, xs)...)
    return Sym{TAG}(head, xs...)
end

##################################################
# Querying apply/combine
##################################################

function apply(op, xs...; context::Context = @__context__)
    return unwrap(query!(_apply(op, xs...), context = context), context = context)
end

function combine(head::Symbol, xs...; context::Context = @__context__)
    unwrap(query!(_combine(head, xs...), context = context), context = context)
end


text:
# © 2016 Massachusetts Institute of Technology.  See LICENSE file for details.
using Base.Test

using ArgParse
s = ArgParseSettings()
@add_arg_table s begin
    "--gpu"
        help = "Use GPU (cuda)"
        action = :store_true
end
const args = parse_args(ARGS, s)
const gpu = args["gpu"]

global compare_cpu_gpu
const T = Array{Float64}

if gpu
    using CUDArt
    CUDArt.init(0)
    device(0)
    Base.rand!(a::CudaArray) = copy!(a, rand(size(a)))
    Base.randn!(a::CudaArray) = copy!(a, randn(size(a)))
    function compare_cpu_gpu(cpu, gpu, node, input; tol=1e-5)
      set_values!(gpu, input)
      gpu()

      #if vecnorm(cpu[node] - to_host(gpu[node])) > tol
      #  @show cpu[node]
      #  @show to_host(gpu[node])
      #  @show cpu[node] - to_host(gpu[node])
      #  @show vecnorm(cpu[node] - to_host(gpu[node]))
      #  @show Seep.sym(node)
      #  Seep.graphviz(STDOUT, node)
      #  for a in gpu.pool.arrays
      #    @show to_host(a)
      #  end
      #end

      @test vecnorm(cpu[node] - to_host(gpu[node])) < tol
    end
else
    ENV["SEEP_NO_GPU"]="1"
    compare_cpu_gpu(x...; kw...) = nothing
end

using Seep
import Seep: sym, arg

Base.randn!(a::Array{Float32}) = copy!(a, randn(size(a)))
mrandn!(as...) = for a in as randn!(a) end

set_values!(inst::Seep.Instance, vals) = for (n,v) in vals copy!(inst[n], v) end

function runtest(y::ANode, x::Vector; fill!::Function=mrandn!, jop=nothing, v_tol=1e-5, d_tol=1e-5, delta=1e-3)
  srand(2)

  grad = gradients(y)
  cpu = instance(T,[y; x; map(a->grad[a], x)])
  igpu = gpu ? instance(CudaArray{eltype(T)}, [y; x; map(a->grad[a], x)]) : nothing

  for i in 1:25
    fill!(map(a->cpu[a], x)...)
    x0 = [(a,copy(cpu[a])) for a in x]
    cpu()

    y0 = copy(cpu[y])

    if jop != nothing
      yref = @eval $(jop)($(map(a->cpu[a], x))...)
      @test !(vecnorm(y0 - yref) > v_tol)
    end

    compare_cpu_gpu(cpu, igpu, y, x0)

    for a in x
      set_values!(cpu, x0)

      a0 = copy(cpu[a])
      gr = copy(cpu[grad[a]])

      cpu()
      compare_cpu_gpu(cpu, igpu, grad[a], x0)

      for j in 1:length(a0)
        set_values!(cpu, x0)

        dx = zeros(eltype(T), size(a0))
        dx[j] += delta*randn()

        copy!(cpu[a], a0+dx)
        cpu()

        if (sum(y0) + sum(gr.*dx) - sum(cpu[y]) > d_tol)
          @show a0
          @show y0
          @show gr
          @show dx
          @show gr.*dx
          @show cpu[y]
          @show sum(y0) + sum(gr.*dx) - sum(cpu[y])
        end
        @test !(sum(y0) + sum(gr.*dx) - sum(cpu[y]) > d_tol)
      end
    end
  end

  if gpu
    close(igpu._pool)
  end
end

println("5x5")
for op in [:+, :.+, :-, :.-, :.*, :*, :At_mul_B, :A_mul_Bt, :At_mul_Bt]
    @show op
    a = ANode(5, 5)
    b = ANode(5, 5)
    y = @eval $op($a, $b)
    runtest(y, [a, b], jop=op)
end

for op in [:./, :.^]
    @show op
    a = ANode(5, 5)
    b = ANode(5, 5)
    y = @eval $op($a, $b)
    runtest(y, [a, b], fill! =(a,b)->(rand!(a); copy!(b, rand(size(b))+1)), jop=op, delta=1e-5)
end

for op in [:dot]
    @show op
    a = ANode(5, 5)
    b = ANode(5, 5)
    y = @eval $op($a, $b)
    runtest(y, [a, b])
end

println("5x5 ⊙ 5x1")
for op in [:.+, :.-, :.*, :*, :At_mul_B]
    @show op
    a = ANode(5, 5)
    b = ANode(5, 1)
    y = @eval $op($a, $b)
    runtest(y, [a, b], jop=op)
end

println("5x5 ⊙ 1x5")
for op in [:.+, :.-, :.*, :At_mul_Bt, :A_mul_Bt]
    @show op
    a = ANode(5, 5)
    b = ANode(1, 5)
    y = @eval $op($a, $b)
    runtest(y, [a, b], jop=op)
end

println("5x1 ⊙ 5x5")
for op in [:.+, :.-, :.*, :At_mul_B, At_mul_Bt]
    @show op
    a = ANode(5, 1)
    b = ANode(5, 5)
    y = @eval $op($a, $b)
    runtest(y, [a, b], jop=op)
end

println("1x5 ⊙ 5x5")
for op in [:.+, :.-, :.*, :*, :A_mul_Bt]
    @show op
    a = ANode(1, 5)
    b = ANode(5, 5)
    y = @eval $op($a, $b)
    runtest(y, [a, b], jop=op)
end

println("vectors")
for op in [:exp, :log, :tanh, :sqrt, :sum, :sigm, :log1p, :expm1, :abs, :sign]
    @show op
    a = ANode(1, 5)
    y = @eval $op($a)
    runtest(y, [a], fill! = (a)->(rand!(a); a[:] += 0.5), jop=op, d_tol=1e-3, v_tol=1e-3)
end

for op in [:softmax]
    @show op
    a = ANode(1, 5)
    y = @eval $op($a)
    runtest(y, [a])
end

println("x,::Real")
for op in [:+, :-, :.+, :.-, :*, :.*]
    @show op

    a = ANode(5, 5)
    b = randn()
    y = @eval $op($a, $b)
    runtest(y, [a])

    a = ANode(5, 5)
    b = randn()
    y2 = @eval $op($b, $a)
    runtest(y2, [a])
end

println("getindex")
for index in (collect(1:5), reshape(25:-1:1, 5, 5))
    a = ANode(5, 5)
    y = a[index]
    runtest(y, [a])
end

println("reduce")
for f in [:sum, :minimum, :maximum]
    @show f

    a = ANode(5, 5)
    y = @eval $f($a)
    runtest(y, [a], delta=1e-5, d_tol=1e-6)

    y = @eval $f($a, 1)
    runtest(y, [a], delta=1e-5, d_tol=1e-6)

    y = @eval $f($a, 2)
    runtest(y, [a], delta=1e-5, d_tol=1e-6)
end

println("extrema")
for f in [:min!, :max!]
    @show f
    a = ANode(5, 5)

    y = @eval Seep.$f($a, 0)
    runtest(y, [a], delta=1e-5, d_tol=1e-5)
end

for f in [:min, :max]
    @show f

    a = ANode(5, 5)
    b = ANode(5, 5)

    y = @eval Seep.$f($a, 0)
    runtest(y, [a], delta=1e-5, d_tol=1e-6)

    y = @eval Seep.$f($a, $b)
    runtest(y, [a, b], delta=1e-5, d_tol=1e-6)
end

let
    a = ANode(1, 10, 10, 9)
    b = ANode(2, 1, 3, 3)

    @show :conv2
    y = conv2(b, a)
    runtest(y, [a,b], delta=1e-5, d_tol=1e-6)

    y = conv2(b, a, padding=(2,2,2))
    runtest(y, [a,b], delta=1e-5, d_tol=1e-6)

    @show :pool
    y = Seep.pool(a, (1,2,2))
    runtest(y, [a], delta=1e-5, d_tol=1e-6)

    y = Seep.pool(a, (1,2,2), padding=(2,2,2))
    runtest(y, [a], delta=1e-5, d_tol=1e-1)
end

if gpu
  CUDArt.close(0)
end


text:
function walk1D(T)
    x = 0
    for _ in 1:T
        x += rand() < 0.5 ? 1 : -1
    end
    return x
end

function create_random_walk!(positions, T)
    for i in 1:length(positions)
        positions[i] = walk1D(T)
    end
end

# using BenchmarkTools;
# positions = zeros(Int, 64);
# T = 10^7;
# @btime create_random_walk!(positions, T);
# @benchmark create_random_walk!(positions, T)
# Threads.nthreads()

max_stars_repo_name: JamieMair/intro-to-julia-talk

max_issues_repo_name: JamieMair/intro-to-julia-talk

max_forks_repo_name: JamieMair/intro-to-julia-talk


text:
using Gen
using PyPlot

# depends on: scenes.jl, path_planner.jl

@gen function model(scene::Scene, times::Vector{Float64})

    # start point of the agent
    start_x = @trace(uniform(0, 1), :start_x)
    start_y = @trace(uniform(0, 1), :start_y)
    start = Point(start_x, start_y)

    # goal point of the agent
    stop_x = @trace(uniform(0, 1), :stop_x)
    stop_y = @trace(uniform(0, 1), :stop_y)
    stop = Point(stop_x, stop_y)

    # plan a path that avoids obstacles in the scene
    maybe_path = plan_path(start, stop, scene, PlannerParams(300, 3.0, 2000, 1.))

    # speed
    speed = @trace(uniform(0, 1), :speed)

    # walk path at constant speed
    if !isnull(maybe_path)
        path = get(maybe_path)
        locations = walk_path(path, speed, times)
    else
        locations = fill(start, length(times))
    end

    # generate noisy observations
    noise = 0.1 * @trace(uniform(0, 1), :noise)
    for (i, point) in enumerate(locations)
        @trace(normal(point.x, noise), i => :x)
        @trace(normal(point.y, noise), i => :y)
    end

    (start, stop, speed, noise, maybe_path, locations)
end

function render(scene::Scene, trace, ax;
                show_measurements=true,
                show_start=true, show_stop=true,
                show_path=true, show_noise=true,
                start_alpha=1., stop_alpha=1., path_alpha=1.)
    (start, stop, speed, noise, maybe_path, locations) = get_call_record(trace).retval
    assignment = get_choices(trace)
    render(scene, ax)
    sca(ax)
    if show_start
        scatter([start.x], [start.y], color="blue", s=100, alpha=start_alpha)
    end
    if show_stop
        scatter([stop.x], [stop.y], color="red", s=100, alpha=stop_alpha)
    end
    if !isnull(maybe_path) && show_path
        path = get(maybe_path)
        for i=1:length(path.points)-1
            prev = path.points[i]
            next = path.points[i+1]
            plot([prev.x, next.x], [prev.y, next.y], color="black", alpha=0.5, linewidth=5)
        end
    end

    # plot locations with measurement noise around them
    scatter([pt.x for pt in locations], [pt.y for pt in locations],
        color="orange", alpha=path_alpha, s=25)
    if show_noise
        for pt in locations
            circle = patches[:Circle]((pt.x, pt.y), noise, facecolor="purple", alpha=0.2)
            ax[:add_patch](circle)
        end
    end

    # plot measured locations
    if show_measurements
        measured_xs = [assignment[i => :x] for i=1:length(locations)]
        measured_ys = [assignment[i => :y] for i=1:length(locations)]
        scatter(measured_xs, measured_ys, marker="x", color="black", alpha=1., s=25)
    end
end