text:
include("../src/iPredict.jl")
#include("src/iPredict.jl")
using .iPredict
using Test

#An extremely general systems test - make sure the output is a dictionary!
result = globalPredict("test/data/small_counts.txt",
"test/data/small_gene_isoform_relationships.txt", 1, true)

@test isa(result, Dict) == true

#Check an error is thrown when the gene names don't match
@test_throws ArgumentError globalPredict("test/data/test_counts_matrix.txt",
"test/data/all_isoforms.txt", 1, false)

#Do some general systems checks for oneGenePrediction too
result = oneGenePrediction("test/data/small_counts.txt",
"test/data/small_gene_isoform_relationships.txt", "G1", 10, true)

@test isa(result, Dict) == true
@test length(result) == 1

result = oneGenePrediction("test/data/alt_small_counts.txt",
"test/data/small_gene_isoform_relationships.txt", "G1", 10, true)

@test isa(result, Dict) == true
@test length(result) == 2

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

include(joinpath("helpdb", "Base.jl"))

ext: jl

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
include("BayesianSZNet.jl")
using .BayesianSZNet

p = parse(Float64, ARGS[1])
λ = parse(Float64, ARGS[2])
patience = parse(Int, ARGS[3])

X_tr, z_tr, X_va, z_va = getdr12data()
multitrain(FCNN, X_tr, z_tr, X_va, z_va, constructorname="fcnn"; p, patience, λ)

max_stars_repo_name: podondra/bdl-redshift

max_issues_repo_name: podondra/bdl-redshift

max_forks_repo_name: podondra/bdl-redshift


text:
################################################################################
#
#   FakeAbsOrdQuoRing
#
################################################################################

# The following is a fake implementation for finite O/I to speed up
# multiplication in the quotient ring. This is good for applications where you
# only do computations in O/I and there is not much going back and forth
# between O and O/I.
#
# At the moment, the additive structure is not implemented and there is some
# other restriction, but it would not be hard to implement the full
# ring interface.
#
# We use the SNF basis SI of I, meaning that if a = (x_1,...,x_n) is an element
# with respect to this basis, then reduction modulo I is very cheap.  It is
# just the reduction of the x_i mod n_i.
#
# We represent [a] in O/I by the vector (x_1,...,x_n) *and* an nxn matrix M
# over Z/mZ, where m is the exponent of O/I. The matrix M is the transpose of
# the representation matrix of z -> z * a with respect to SI. Note that this
# matrix M does not identify a uniquely.
#
# Two elements (x, M), (y, N) are multiplied by doing (N * x, N * M). Note that
# N * x has to be reduced modulo the n_i.
#
# The reason the transposing and change of order appears is not fully written
# down yet, but it seems to work.

mutable struct FakeAbsOrdQuoRing{S, T, U, V}
  O::S                      # the order
  id::T                     # the ideal
  B::U                      # the SNF basis of the ideal
  eldiv::V                  # the elementary divisors, Vector{Int}
  adjusted_bmat::fmpz_mat   # the basis matrix of B, from basis(O) to B
  n::Int                    # the exponent of O/id
  ZnZ::NmodRing             # the residue ring Z/nZ
  onevec::Vector{Int}       # vector of the one
end

mutable struct FakeAbsOrdQuoRingElem{S}
  rep_mat::nmod_mat         # reduction of representation matrix mod n
  v::Vector{Int}            # coordinate vector with respect to the SNF basis
                            # this is always reduced
  par::S                    # parent
end

function FakeAbsOrdQuoRing(O::S, id::T) where {S, T}
  B = integral_basis_matrix_wrt(id, O)
  Snf, U, V = snf_with_transform(B)
  # the new basis is given by V^-1, since Snf = U * B * V
  Vinv = inv(V)
  newbas = [ O(_eltseq(Vinv[i, :])) for i in 1:degree(O)]
  @assert O == Order(algebra(O), elem_in_algebra.(newbas))
  eld = diagonal(Snf)
  # I am assuming everything fits Int :)
  n = Int(eld[end])
  ZnZ = ResidueRing(ZZ, n)

  coord = coordinates(one(O)) * V
  coords = Int[Int(mod(coord[i], eld[i])) for i in 1:length(eld)]

  return FakeAbsOrdQuoRing{typeof(O), typeof(id), typeof(newbas), typeof(eld)}(O, id, newbas, eld, Vinv, n, ZnZ, coords)
end

function (R::FakeAbsOrdQuoRing)(x::AlgAssAbsOrdElem)
  repmat = transpose(change_base_ring(R.ZnZ,
                     representation_matrix_wrt(x, R.B, :right)))
  coord = coordinates(x) * inv(R.adjusted_bmat)
  coords = Int[Int(mod(coord[i], R.eldiv[i])) for i in 1:length(R.eldiv)]
  return FakeAbsOrdQuoRingElem{typeof(R)}(repmat, coords, R)
end

# equality
function Base.:(==)(x::FakeAbsOrdQuoRingElem, y::FakeAbsOrdQuoRingElem)
  return x.v == y.v
end

# multiplication
function Base.:(*)(x::FakeAbsOrdQuoRingElem, y::FakeAbsOrdQuoRingElem)
  rmat = y.rep_mat * x.rep_mat
  vv = map(n -> Int(n.data), y.rep_mat * (x.par.ZnZ.(x.v)))
  for i in 1:length(vv)
    vv[i] = mod(vv[i], x.par.eldiv[i])
  end
  return FakeAbsOrdQuoRingElem{typeof(x.par)}(rmat, vv, x.par)
end

# isone
function isone(x::FakeAbsOrdQuoRingElem)
  return x.v == x.par.onevec
end

function one(R::FakeAbsOrdQuoRing)
  return FakeAbsOrdQuoRingElem{typeof(R)}(identity_matrix(R.ZnZ, degree(R.O)), R.onevec, R)
end

# lift, mainly for debugging purposes
function lift(x::FakeAbsOrdQuoRingElem)
  c = fmpz.(x.v)
  return dot(c, x.par.B)
end


text:
module TiledArrays

# Ceiling division
cdiv(a::Integer, b::Integer) = cdiv(promote(a, b)...)
cdiv(a::T, b::T) where {T<:Integer} = one(T) + div(a - one(T), b)

# Opposite of `Base.tail`.
# Keep all but the last element in a tuple.
head(x, y...) = (x, head(y...)...)
head(x) = ()

# # Type traits for describing if there's something scary going on with the strides.
# # For now, assume memory strides are monotonically increasing such that the underlying
# # array describes a contiguous block of memory.
# abstract type AbstractStrides end
# struct DenseStrides <: AbstractStrides end

struct Tile
    dim::Int
    size::Int
end
Base.show(io::IO, tile::Tile) = print(io, "Tile($(tile.dim), $(tile.size))")

# Note that `NonBlockedLayout <: MemoryLayout`, but without any tiling, allowing
# for specific handling of the easy case.
const MemoryLayout{N} = Tuple{Vararg{Union{Int,Tile},N}} where {N}
const NonBlockedLayout{N} = Tuple{Vararg{Int,N}} where {N}

# Utility function for removing all "Tile"s from the front of an array.
function strip_leading_tiles(x::Tuple{Vararg{Union{Int,Tile},N}}) where {N}
    return strip_leading_tiles(x...)
end

function strip_leading_tiles(x::Tile, y::Vararg{Union{Int,Tile},N}) where {N}
    return strip_leading_tiles(y...)
end

strip_leading_tiles(x::Int, y::Vararg{Union{Int,Tile},N}) where {N} = (x, y...)

#####
##### Layout Validation
#####

# Validation
# Invariants to maintain:
#
# 1. We must have exactly one occurance of an Integer index.
# 2. All `Tile` indexes must occur BEFORE the corresponding Integer index.
# 3. Each `Tile` index must have a corresponding Integer index
function validate(layout::MemoryLayout{N}) where {N}
    @nospecialize

    # Pass 1 - Make sure there are no nonpositive indices.
    # Also, find the maximum index.
    maximum_index = 0
    for desc in layout
        if isa(desc, Int)
            desc <= 0 && error("Found non-positive index: $desc")
            maximum_index = max(maximum_index, desc)
        elseif isa(desc, Tile)
            desc.dim <= 0 && error("Found non-positive tile index: $desc")
            maximum_index = max(maximum_index, desc.dim)
        else
            # Should never happen due to type restrictions.
            error("Unreachable reached")
        end
    end

    # Pass 2 - Check invariants.
    for index in Base.OneTo(maximum_index)
        saw_int = false
        saw_tile = false
        for desc in layout
            if isa(desc, Int) && desc == index
                saw_int && error("Integer index $index occurs more than once!")
                saw_int = true
            elseif isa(desc, Tile) && desc.dim == index
                saw_int && error(
                    "Integer index $index occurs before its corresponding tile index!"
                )
                saw_tile = true
            end
        end
        # We saw a tile but not an int, throw error
        saw_tile &&
            !saw_int &&
            error("Found tile with index $index with no corresponding Int")

        # We didn't find anything
        !saw_tile && !saw_int && error("Found no entry with index $index")
    end
    return true
end

#####
##### Adjust indexes for padding.
#####

unchecked_divrem(x::Int, y::Int) = (Base.sdiv_int(x, y), Base.srem_int(x, y))

# Note: not necessarily safe since we are using the unchecked integer division functions.
# Also, `index` should already be converted from index-1 to index-0.
function adjust_for_padding(
    size::NTuple{N,Int}, padded_size::NTuple{N,Int}, index::NTuple{N,Int}
) where {N}
    return ntuple(i -> adjust_for_padding(size[i], padded_size[i], index[i]), Val(N))
end

function adjust_for_padding(size::Int, padded_size::Int, index::Int)
    size == padded_size && return index
    a, b = unchecked_divrem(index, size)
    return (a * padded_size) + b
end

#####
##### Index Conversion Logic
#####

# Expand the index according to the given layout.
# NOTE: `I` should be post padding adjusting if necessary.
function splitindex(::Val{T}, I::Tuple{Vararg{Int,N}}) where {T,N}
    return _splitindex(T, I)
end

# Work from the first index to the last.
# At each step, strip off the index we processed from `layout` and recurse.
#
# When `layout[1]` is an Integer, we just return the corresponding index from `I`.
#
# When `layout[1]` is a Tile, then we grab the corresponding index from `I`,
# `divrem` that value, then reconstruct `I` with the reduced index.
#
# NB: This depends on Julia's super agressive constant propagation heuristics when dealing
# with simplifying tuple recursive functions to generate great assembly code.
@inline function _splitindex(layout::MemoryLayout, I::Tuple{Vararg{Int,N}}) where {N}
    v, _I = process(layout[1], I)
    return (v, _splitindex(Base.tail(layout), _I)...)
end

### Base case - last item should never be a `Tile`.
@inline function _splitindex(layout::Tuple{Int}, I::Tuple{Vararg{Int,N}}) where {N}
    return (I[layout[1]],)
end

function _splitindex(::Tuple{Tile}, I::Tuple{Vararg{Int,N}}) where {N}
    return error("Malformed memory layout with a `Tile` in the tail position!")
end

@inline process(x::Int, I::Tuple{Vararg{Int,N}}) where {N} = (I[x], I)
@inline function process(x::Tile, I::Tuple{Vararg{Int,N}}) where {N}
    outer, inner = unchecked_divrem(I[x.dim], x.size)
    return (inner, ntuple(i -> i == x.dim ? outer : I[i], Val(N)))
end

#####
##### Obtaining full dimensions
#####

# Move layout from type domain to value domain, providing "const" information to the
# compiler.
@inline dims(valT::Val{T}, size::Tuple{Vararg{Int,N}}) where {T,N} = _dims(T, size)

@inline function _dims(layout::MemoryLayout, size::Tuple{Vararg{Int,N}}) where {N}
    v, _size = _dimcheck(layout[1], size)
    return (v, _dims(Base.tail(layout), _size)...)
end

# Bottom of recursion.
@inline _dims(layout::Tuple{Int}, size::Tuple{Vararg{Int,N}}) where {N} = (size[layout[1]],)

@inline _dimcheck(x::Int, size::Tuple{Vararg{Int,N}}) where {N} = (size[x], size)
@inline function _dimcheck(x::Tile, size::Tuple{Vararg{Int,N}}) where {N}
    newdim = cdiv(size[x.dim], x.size)
    return (x.size, ntuple(i -> i == x.dim ? newdim : size[i], Val(N)))
end

#####
##### Users of "dims"
#####

# NB: Cannot extend `Base.strides` because a `MemoryLayout` can potentially just be a tuple
# of Ints, so we'd be committing type piracy.
function _strides(valT::Val{T}, size::Tuple{Vararg{Int,N}}) where {T,N}
    # Drop the last value from `dims`.
    return cumprod((one(Int), head(dims(valT, size)...)...))
end

function fullsize(valT::Val{T}, size::Tuple{Vararg{Int,N}}) where {T,N}
    # Constant propagation to the rescue.
    # If the layout is non-blocked, just standard or permuted, then there is no
    # over-allocation necessary.
    if isa(T, NonBlockedLayout)
        return prod(size)
        # Otherwise, we need to compute the full dimensions, which could include zero padding
        # to keep the inner blocksizes constant.
    else
        return prod(dims(valT, size))
    end
end

#####
##### Linear Offset
#####

function getoffset(
    valT::Val{T},
    size::Tuple{Vararg{Int,N}},     # Padded Size
    I::Tuple{Vararg{Int,N}},        # Post-padded index
) where {T,N}
    Base.@_inline_meta
    strides = _strides(valT, size)
    index = splitindex(valT, I)
    offset = zero(Int)
    for i in eachindex(strides)
        @inbounds offset += strides[i] * index[i]
    end
    return offset
end

# Generates indices into a tiled array.
# Parameters:
#
# `L`: Layout type parameter.
# `N`: Number of Dims.
struct TiledIndexer{L,N}
    size::NTuple{N,Int}
    padded_size::NTuple{N,Int}
end

function TiledIndexer{L}(size::NTuple{N,Int}, padded_size::NTuple{N,Int}) where {L,N}
    return TiledIndexer{L,N}(size, padded_size)
end

function genindex(x::TiledIndexer{L,N}, I::NTuple{N,Int}) where {L,N}
    return 1 + getoffset(
        Val(L), x.padded_size, adjust_for_padding(x.size, x.padded_size, I .- 1)
    )
end

# function genindex(x::TiledIndexer{L,N,Nothing}, I::NTuple{N,Int}) where {L,N}
#     return getoffset(Val(N), x.size, I)
# end

# #####
# ##### Tiled Array
# #####
#
# struct TiledArray{T,L,N} <: AbstractArray{T,N}
#     parent::Vector{T}
#     size::NTuple{N,Int}
# end
# Base.parent(A::TiledArray) = A.parent
#
# function TiledArray{T,L}(::UndefInitializer, size::Tuple{Vararg{Int,N}}) where {T,L,N}
#     # Allocate enough room for the whole layout.
#     parent = Vector{T}(undef, fullsize(Val(L), size))
#     return TiledArray{T,L,N}(parent, size)
# end
#
# layout(::Type{<:TiledArray{<:Any,L}}) where {L} = L
# layout(x::T) where {T<:TiledArray} = layout(T)
# vlayout(x) = Val(layout(x))
#
# # Array interface
# @inline Base.size(x::TiledArray) = x.size
# Base.IndexStyle(::Type{<:TiledArray}) = Base.IndexCartesian()
# Base.@propagate_inbounds function Base.getindex(
#     A::TiledArray{T,L,N}, I::Vararg{Int,N}; strides = nothing
# ) where {T,L,N}
#     @boundscheck checkbounds(A, I...)
#     offset = getoffset(vlayout(A), size(A), I, strides) + one(Int)
#     return @inbounds(getindex(parent(A), offset))
# end
#
# Base.@propagate_inbounds function Base.setindex!(
#     A::TiledArray{T,L,N}, v, I::Vararg{Int,N}; strides = nothing
# ) where {T,L,N}
#     @boundscheck checkbounds(A, I...)
#     offset = getoffset(vlayout(A), size(A), I, strides) + one(Int)
#     return @inbounds setindex!(parent(A), v, offset)
# end

end # module


text:
# =========================================================================== #
# Error Functions
# =========================================================================== #

abstract type AbstractErrorFunction
end

struct LeastSquaresErrorFunction <: AbstractErrorFunction
end

function (error_fcn::LeastSquaresErrorFunction)(meas::T, x::T) where {T <: Number}
    dx = -meas + x
    return dx^2
end
function (error_fcn::LeastSquaresErrorFunction)(meas::Vector{T}, x::Vector{T}) where {T}
    dx = sum( (-meas .+ x).^2 )
    return dx
end
function (error_fcn::LeastSquaresErrorFunction)(meas, x::Vector{T}) where {T}
    dx = sum( (-meas .+ x).^2 )
    return dx
end

function differentiate(error_fcn::LeastSquaresErrorFunction, meas, x)
    return -meas + x
end

function eval(error_fcn::LeastSquaresErrorFunction, meas, x)
    return error_fcn(meas, x)
end

function gradient!(error_fcn::LeastSquaresErrorFunction, y_meas, y_out, delta_e)
    for kk = 1:length(delta_e)
        delta_e[kk] = differentiate(error_fcn, y_meas[kk], y_out[kk])
    end
    return delta_e
end

# function error_gradient!(nn::NeuralNetwork, y_meas, y_out, delta_e)
#     return gradient!(nn.error_function, y_meas, y_out, delta_e)
# end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
p1 = plot(title="Title 1")
@test p1.layout["title"] == "Title 1"

p2 = title!(plot(), "Title 2")
@test p2.layout["title"] == "Title 2"

p3 = plot()
title!("Title 3")
@test p3.layout["title"] == "Title 3"

p4 = plot(xlabel="var x")
@test p4.layout["xaxis"]["title"] == "var x"

p5 = plot()
p6 = ylabel!("var y")

test_is_same_object(p5, p6)
@test p5.layout["yaxis"]["title"] == "var y"


text:
# Generate Julia modules based on Google API Discovery Service.
#    See https://developers.google.com/discovery/
#
module Bootstrap

module Rest
  using HTTP
  using LazyJSON
  using Base.Iterators
  const apisUrl = "https://www.googleapis.com/discovery/v1/apis"
  # draft-zyp-json-schema-03
  # https://tools.ietf.org/html/draft-zyp-json-schema-03
  export JsonSchema
  mutable struct JsonSchema
    name::String
    type::String
    properties
    patternProperties
    additionalProperties
    items::JsonSchema
    additionalItems
    required::Bool
    dependencies
    minimum
    maximum
    exclusiveMinimum
    exclusiveMaximum
    minItems
    maxItems
    uniqueItems
    pattern::String
    minLength
    maxLength
    enum
    default::String
    title::String
    description::String
    format::String
    divisibleBy
    disallow
    extends
    id
    ref # $ref
    schema # $schema
    JsonSchema() = new()
  end
  function Base.convert(::Type{JsonSchema}, x::LazyJSON.Object{Nothing, String})
    # @info "*** convert(::Type{JsonSchema}", x
    if x == nothing return nothing end
    r = JsonSchema()
    for fld in drop(fieldnames(JsonSchema), 1)
      jkey = String(fld)
      if jkey == "ref" || jkey == "schema"
        jkey = "\$"*jkey
      end
      if haskey(x, jkey)
        # @info "$(fld)", fieldtype(JsonSchema, fld)
        setfield!(r, fld, convert(fieldtype(JsonSchema, fld), x[jkey]))
      end
    end
    r
  end
  export Parameter
  mutable struct Parameter
    name::String
    type::String
    format::String
    default::String
    location::String
    required::Bool
    description::String
    ord::Int32 # manually set
    Parameter() = new()
  end
  function Base.isless(x::Parameter, y::Parameter)
    if x.ord > -1 || y.ord > -1
      return x.ord < y.ord
    end
    x.name < y.name
  end
  function Base.convert(::Type{Parameter}, x::LazyJSON.Object{Nothing, String})
    if x == nothing return nothing end
    r = Parameter()
    for fld in drop(fieldnames(Parameter), 1)
      if haskey(x, String(fld))
        setfield!(r, fld, convert(fieldtype(Parameter, fld), x[String(fld)]))
      end
    end
    r
  end
  function Base.convert(::Type{Vector{Parameter}}, x::LazyJSON.Object{Nothing, String})
    if x == nothing return nothing end
    r = Vector{Parameter}()
    for kvp in x
      push!(r, convert(Parameter, kvp.second))
      r[end].name = kvp.first
    end
    r
  end
  export Property
  mutable struct Property
    name::String
    type::String
    format::String
    items::JsonSchema # TODO $ref, object (see Dataset) array of object could be made Dataframe
    # additionalProperties
    required::Bool
    description::String
    default::String
    pattern::String
    ref::String
    Property() = new()
  end
  function Base.convert(::Type{Property}, x::LazyJSON.Object{Nothing, String})
    if x == nothing return nothing end
    r = Property()
    for fld in drop(fieldnames(Property), 1)
      jkey = String(fld)
      if jkey == "ref"
        jkey = "\$"*jkey
      end
      if haskey(x, jkey)
        setfield!(r, fld, convert(fieldtype(Property, fld), x[jkey]))
      end
    end
    r
  end
  function Base.convert(::Type{Vector{Property}}, x::LazyJSON.Object{Nothing, String})
    # @info "convert(Type{Vector{Property}}", x
    if x == nothing return nothing end
    r = Vector{Property}()
    for kvp in x
      push!(r, convert(Property, kvp.second))
      r[end].name = kvp.first
    end
    r
  end
  export Schema
  mutable struct Schema
    name::String
    properties::Vector{Property}
    refs::Vector{Schema}
    expr::Expr
    Schema() = new()
  end
  Base.isless(x::Schema, y::Schema) = x.name < y.name
  function Base.convert(::Type{Schema}, x::LazyJSON.Object{Nothing, String})
    if x == nothing return nothing end
    r = Schema()
    r.refs = Vector{Schema}()
    for fld in drop(fieldnames(Schema), 1)
      if haskey(x, String(fld))
        setfield!(r, fld, convert(fieldtype(Schema, fld), x[String(fld)]))
      end
    end
    r
  end
  function Base.convert(::Type{Dict{String,Schema}}, x::LazyJSON.Object{Nothing, String})
    if x == nothing return nothing end
    r = Dict{String,Schema}()
    for kvp in x
      r[kvp.first] = convert(Schema, kvp.second)
      r[kvp.first].name = kvp.first
    end
    r
  end
  export Method
  mutable struct Method
    name::String
    path::String
    httpMethod::String
    description::String
    parameters::Vector{Parameter}
    parameterOrder::Vector{String}
    request::JsonSchema
    response::JsonSchema
    expr::Expr
    Method() = new()
  end
  Base.isless(x::Method, y::Method) = x.name < y.name
  function Base.convert(::Type{Method}, x::LazyJSON.Object{Nothing,String})
    if x == nothing return nothing end
    r = Method()
    for fld in drop(fieldnames(Method), 1)
      if haskey(x, String(fld))
        setfield!(r, fld, convert(fieldtype(Method, fld), x[String(fld)]))
      end
    end
    r
  end
  function Base.convert(::Type{Vector{Method}}, x::LazyJSON.Object{Nothing,String})
    if x == nothing return nothing end
    r = Vector{Method}()
    for kvp in x
      push!(r, convert(Method, kvp.second))
      r[end].name = kvp.first
    end
    r
  end
  export Resource
  mutable struct Resource
    name::String
    methods::Vector{Method}
    expr::Expr
    Resource() = new()
  end
  function Base.convert(::Type{Resource}, x::LazyJSON.Object{Nothing,String})
    if x == nothing return nothing end
    r = Resource()
    if haskey(x, "methods")
      r.methods = convert(Vector{Method}, x["methods"])
    end
    r
  end
  function Base.convert(::Type{Vector{Resource}}, x::LazyJSON.Object{Nothing,String})
    if x == nothing return nothing end
    r = Vector{Resource}()
    for kvp in x
      push!(r, convert(Resource, kvp.second))
      r[end].name = kvp.first
    end
    r
  end
  export Service
  mutable struct Service
    name::String
    title::String
    description::String
    baseUrl::String
    parameters::Vector{Parameter}
    schemas::Dict{String,Schema}
    resources::Vector{Resource}
    expr::Expr
    Service() = new()
  end
  function Base.convert(::Type{Service}, x::LazyJSON.Object{Nothing,String})
    if x == nothing return nothing end
    r = Service()
    for fld in fieldnames(Service)
      if haskey(x, String(fld))
        setfield!(r, fld, convert(fieldtype(Service, fld), x[String(fld)]))
      end
    end
    r
  end

  function apis(
    ;name::Union{Symbol,Nothing}=nothing
    ,preferred::Union{Bool,Nothing}=nothing
    ,kws...
  )
    query = nothing
    if name != nothing || preferred != nothing
      query = []
      if name != nothing
        push!(query, :name => name)
      end
      if preferred != nothing
        push!(query, :preferred => preferred)
      end
    end
    res = HTTP.get(apisUrl; query=query, kws...)
    LazyJSON.value(String(res.body))
  end
  function get(url::String; kws...)
    res = HTTP.get(url; kws...)
    convert(Service, LazyJSON.value(String(res.body)))
  end
  function getService(name::Symbol; kws...)
    jDirList = apis(;name=name, preferred=true, kws...)
    apiUrl = String(jDirList["items"][1]["discoveryRestUrl"])
    get(apiUrl; kws...)
  end

end # module Rest

module Emit
  using HTTP
  using LazyJSON
  using LazyJSON: convert, value
  using MacroTools: @q, unblock
  using Dates
  using DataFrames
  using ..Rest

  function exprToString(expr::Expr)
    io = IOBuffer()
    print(io, expr)
    String(take!(io))
  end

  # from: https://github.com/JuliaCloud/LazyJSON.jl/blob/035be5ce9bd435157f9291548eba640f507255e6/src/AbstractDict.jl#L64
  isreserved(s::String) = Symbol(s) in (:while, :if, :for, :try, :return, :break, :continue,
  :function, :macro, :quote, :let, :local, :global, :const, :do,
  :struct,
  :type, :immutable, :importall,  # to be deprecated
  :module, :baremodule, :using, :import, :export,
  :end, :else, :catch, :finally, :true, :false)

  # Structs are written to file based dependency chain
  # so that most-dependent Structs are last
  # https://www.geeksforgeeks.org/topological-sorting/
  function topologicalSort(cur::Schema, stack::Vector{Schema})
    if length(cur.refs) != 0
      for nxt in cur.refs
        if cur != nxt
          topologicalSort(nxt, stack)
        end
      end
    end
    if !in(cur, stack)
      push!(stack, cur)
    end
  end

  function getEmit(names::Symbol...; kws...)
    services = Vector{Service}()
    for name in names
      push!(services, Rest.getService(name; kws...))
    end
    emitGcp(services...)
  end

  function emitGcp(services::Service...)
    @info "emitGcp: start ---"
    path = "$(pwd())/src/"
    mkpath(path) # ensure path

    # emit service files
    for service in services
      emitService(service, path)
    end

    eGcpModule = unblock(@q module GCP end)
    # add GCP file includes
    push!(eGcpModule.args[3].args, @q include("util.jl") )
    push!(eGcpModule.args[3].args, @q include("credential.jl") )
    push!(eGcpModule.args[3].args, @q include("session.jl") )
    for service in services
      push!(eGcpModule.args[3].args, @q include($("$(service.name).jl")))
    end # for service file include
    # add usings for generated Services
    for service in services
      push!(eGcpModule.args[3].args, unblock(@q using .$(Symbol(uppercasefirst(service.name)))) )
    end # for: service.name: usings
    # add Project struct
    eProject = unblock(@q mutable struct Project end)
    push!(eGcpModule.args[3].args, eProject)
    # add fields for Project
    push!(eProject.args[3].args, unblock(@q id::String) )
    push!(eProject.args[3].args, unblock(@q credential::Credential) )
    push!(eProject.args[3].args, unblock(@q session::Session) )
    push!(eProject.args[3].args, unblock(@q authorize) )
    for service in services
      fieldName = Symbol(service.name)
      fieldType = Symbol(uppercasefirst(service.name))
      eField = unblock(@q $fieldName::$fieldType.Service)
      push!(eProject.args[3].args, eField)
    end # for: apiName: Project fields
    # add constructor for Project struct
    eProjectCtor = unblock(@q function Project(credential, scopes::AbstractVector{<: AbstractString}) p = new() end)
    push!(eProject.args[3].args, eProjectCtor)
    # initialize resource fields in constructor

    push!(eProjectCtor.args[2].args, unblock(@q begin
      if isa(credential, String) && isfile(credential)
          p.credential = Credential(credential)
      else
          p.credential = credential
      end
      p.session = Session(p.credential, scopes)
      p.id = p.credential.project_id
      p.authorize = function authorize(; kws...)
        auth = GCP.authorize(p.session; kws...)
        headers = Dict{String, String}("Authorization" => "$(auth[:token_type]) $(auth[:access_token])")
        headers
      end
    end ) )
    # push!(eProjectCtor.args[2].args, unblock(@q p.id = credential.project_id ) )
    # push!(eProjectCtor.args[2].args, unblock(@q p.credential = credential ) )
    for service in services
      varName = Symbol(service.name)
      moduleName = Symbol(uppercasefirst(service.name))
      push!(eProjectCtor.args[2].args, unblock(@q p.$varName = $moduleName.Service(p) ) )
    end # for: eProjectCtor: init fields
    push!(eProjectCtor.args[2].args, unblock(@q p))
    # --- End: Generate Project struct ---

    # --- Write GCP module file to disk ---
    filename = "$(path)GCP.jl" # GCP must be uppercase for Pkg test
    @info "emitGcp: write $(filename)"
    open(filename, "w") do f # write file
      write(f, exprToString(eGcpModule))
    end

    @info "emitGcp: done ---"
  end

  # obj input: Property, Parameter, or JsonSchema
  function typeSymbol(obj; schema=nothing, service=nothing)::Union{Symbol, Nothing}
    # type formats: https://developers.google.com/discovery/v1/type-format
    if isdefined(obj, :type)
      if obj.type == "array"
        # @info "TODO $(obj)"
        if isdefined(obj, :items)
          # @info "obj.items", typeof(obj.items).name, typeSymbol(obj.items), obj.items
          elmSym = typeSymbol(obj.items, schema=schema, service=service)
          if elmSym != nothing
            if schema != nothing && service != nothing && haskey(service.schemas, String(elmSym))
              push!(schema.refs, service.schemas[String(elmSym)])
            end
            return Symbol("Vector{$(elmSym)}")
          end
        end
      elseif obj.type == "boolean"
        return Symbol(Bool.name)
      elseif obj.type == "object"
        # @info "TODO $(obj)"
        if isdefined(obj, :properties)
          properties = obj.properties
          if haskey(properties, "kind") && haskey(properties["kind"], "default")
            name = split(properties["kind"]["default"], "#") # "bigquery#dataset"
            name = uppercasefirst(name[2])
            # ensure type is defined; "bigquery#project" does not have a definition
            if haskey(service.schemas, name)
              return Symbol(name)
            end
            return nothing
          end
        end
      elseif obj.type == "string"
        if !isdefined(obj, :format)
          return Symbol(String.name)
        else
          if obj.format == "byte"
            return Symbol(UInt8.name)
          elseif obj.format == "date"
            return Symbol(Date.name)
          elseif obj.format == "date-time"
            return Symbol(DateTime.name)
          elseif obj.format == "int64"
            return Symbol(Int64.name)
          elseif obj.format == "uint64"
            return Symbol(UInt64.name)
          end
        end
      elseif obj.type == "integer"
        if isdefined(obj, :format)
          if obj.format == "int32"
            return Symbol(Int32.name)
          elseif obj.format == "uint32"
            return Symbol(UInt32.name)
          end
        end
      elseif obj.type == "number"
        if isdefined(obj, :format)
          if obj.format == "float"
            return Symbol(Float32.name)
          elseif obj.format == "double"
            return Symbol(Float64.name)
          end
        end
      else
        # @info "TODO $(obj)"
      end
    elseif isdefined(obj, :ref)
      # @info ">>> obj.ref: $(obj.ref)", service == nothing, obj
      if service != nothing && haskey(service.schemas, obj.ref)
        push!(schema.refs, service.schemas[obj.ref])
      end
      return Symbol(obj.ref)
    else # jProperty: no: type, $ref
      # @info "TODO $(obj): no: type, \$ref: TODO"
    end
    nothing
  end

  # function indexof(arr::Vector{Parameter}, itm::Parameter)
  function indexof(arr, itm)
    for n in 1:length(arr)
      if arr[n] == itm
        return n
      end
    end
    -1
  end

  function symbolisprimitivetype(symbol::Symbol)
    try
      return @eval isprimitivetype($symbol)
    catch e
      return false
    end
  end

  function Base.convert(::Type{Expr}, schema::Schema)
    unblock(@q mutable struct $(Symbol(schema.name)) <: GCP.Abstract end )
  end
  function Base.convert(::Type{Expr}, property::Property, schema::Schema, service::Service)
    # @info ">>> $(service.name).$(schema.name).$(property.name)"
    symbol = typeSymbol(property; schema=schema, service=service)
    if isnothing(symbol) return nothing end
    if symbolisprimitivetype(symbol)
      symbol = Symbol("Union{Nothing, $symbol}")
    end
    # uppercasefirst to allow Julia reserve words to be field names
    fldname = property.name
    if isreserved(fldname)
      fldname = fldname * '_'
    end
    unblock(@q $(Symbol(fldname))::$symbol )
  end

  function methodOverride(eFn::Expr, method::Rest.Method, resource::Resource, service::Service)
    if service.name == "bigquery"
      if resource.name == "tabledata"
        if method.name == "insertAll" # bigquery.tabledata.insertAll
          # method parameters
          push!(eFn.args[1].args, unblock(@q val ))
          push!(eFn.args[1].args, unblock(@q datasetId::String ))
          push!(eFn.args[1].args, unblock(@q tableId::String ))
          push!(eFn.args[1].args[2].args, unblock(@q kws... ))
          # method body
          push!(eFn.args[2].args, unblock(@q begin
            io = IOBuffer()
            write(io, "{\"kind\":\"bigquery#tableDataInsertAllRequest\",rows:[")
            if isa(val, Array)
                for n in 1:length(val)
                    if n != 1 write(io, ",") end
                    write(io, "{json:")
                    write(io, GCP.json(val[n]))
                    write(io, "}")
                end
            else
                write(io, "{json:")
                write(io, GCP.json(val))
                write(io, "}")
            end
            write(io, "]}")
          end))
          push!(eFn.args[2].args, unblock(@q headers = project.authorize(; kws...) ))
          push!(eFn.args[2].args, unblock(@q headers["Content-Type"] = "application/json" ))
          push!(eFn.args[2].args, unblock(@q res = HTTP.request("POST", "$(Bigquery.baseUrl)projects/$(project.id)/datasets/$(datasetId)/tables/$(tableId)/insertAll", headers, String(take!(io)); kws...) ))
          push!(eFn.args[2].args, unblock(@q convert(Bigquery.TableDataInsertAllResponse, LazyJSON.value(String(res.body))) ))
          return true
        elseif method.name == "list" # bigquery.tabledata.list
          # method parameters
          push!(eFn.args[1].args, unblock(@q datasetId::String ))
          push!(eFn.args[1].args, unblock(@q tableId::String ))
          push!(eFn.args[1].args[2].args, unblock(@q kws... ))
          # method body
          push!(eFn.args[2].args, unblock(@q df = DataFrame() ))
          push!(eFn.args[2].args, unblock(@q table = project.bigquery.tables.get(datasetId, tableId) ))
          push!(eFn.args[2].args, unblock(@q if table != nothing && isdefined(table, :schema) && isdefined(table.schema, :fields)
            converts = Vector{Function}()
            for fld in table.schema.fields
              if fld.type_ == "TIMESTAMP"
                df[Symbol(fld.name)] = Vector{DateTime}()
                push!(converts, (str) -> unix2datetime( parse(Float64, str ) ) )
              elseif fld.type_ == "FLOAT"
                df[Symbol(fld.name)] = Vector{Float32}()
                push!(converts, (str) -> parse(Float32, str ) )
              elseif fld.type_ == "INTEGER"
                df[Symbol(fld.name)] = Vector{Int32}()
                push!(converts, (str) -> parse(Int32, str ) )
              elseif fld.type_ == "BOOLEAN"
                df[Symbol(fld.name)] = Vector{Bool}()
                push!(converts, (str) -> parse(Bool, String(str) ) )
              elseif fld.type_ == "STRING"
                df[Symbol(fld.name)] = Vector{String}()
                push!(converts, (str) -> String(str) )
              end
            end
          end))
          push!(eFn.args[2].args, unblock(@q headers = project.authorize(; kws...) ))
          push!(eFn.args[2].args, unblock(@q res = HTTP.request("GET", "$(Bigquery.baseUrl)projects/$(project.id)/datasets/$(datasetId)/tables/$(tableId)/data", headers; kws...) ))
          push!(eFn.args[2].args, unblock(@q jres = LazyJSON.value(String(res.body)) ))
          # push!(eFn.args[2].args, unblock(@q  ))
          push!(eFn.args[2].args, unblock(@q dfFill = function dfFill(jres)
            if haskey(jres, "rows")
              dfRow = Vector{Any}()
              for row in jres["rows"]
                empty!(dfRow)
                cols = row["f"]
                for n in 1:length(cols)
                  push!(dfRow, converts[n]( cols[n]["v"] ) )
                end
                push!(df, dfRow)
              end
            end
            @info "filled", size(df)
          end))
          push!(eFn.args[2].args, unblock(@q dfFill(jres) ))
          push!(eFn.args[2].args, unblock(@q while haskey(jres, "pageToken")
            query = [:pageToken => String(jres["pageToken"])]
            headers = project.authorize(; kws...)
            res = HTTP.request("GET", "$(Bigquery.baseUrl)projects/$(project.id)/datasets/$(datasetId)/tables/$(tableId)/data", headers; query=query, kws...)
            jres = LazyJSON.value(String(res.body))
            dfFill(jres)
          end))
          push!(eFn.args[2].args, unblock(@q df  ))
          return true
        end
      end
    end
    false
  end

  function Base.convert(::Type{Expr}, method::Rest.Method, resource::Resource, service::Service)
    eFn = unblock(@q function $(Symbol(method.name))(;) end)

    # check for manual overried implementation
    if methodOverride(eFn, method, resource, service)
      return eFn
    end

     # add request body object parameter (if any)
     if isdefined(method, :request)
      reqSym = typeSymbol(method.request)
      if reqSym != nothing
        push!(eFn.args[1].args, unblock(@q $(Symbol("body"))::$(Symbol(uppercasefirst(service.name))).$reqSym  ))
      else
        push!(eFn.args[1].args, unblock(@q $(Symbol("body"))))
      end
    end

    # add parameters to method (if any)
    if isdefined(method, :parameters)
      requiredParameters = Vector{Parameter}()
      optionalParameters = Vector{Parameter}()
      for parameter in method.parameters
        if parameter.name == "projectId" continue end # set by project.id
        # @info "$(service.name).$(resource.name).$(method.name) $(parameter.name) required:$(parameter.required)"
        if parameter.required
          if isdefined(method, :parameterOrder)
            parameter.ord = indexof(method.parameterOrder, parameter.name)
            # @info "parameter.ord", parameter.ord, parameter.name
          end
          push!(requiredParameters, parameter)
        else
          push!(optionalParameters, parameter)
        end
      end
      # for n in 1:length(requiredParameters)
      #   @info "$(service.name).$(resource.name).$(method.name) requiredParameters-$n: $(requiredParameters[n].name) parameter.required:$(requiredParameters[n].required) isdefined(method, :parameterOrder):$(isdefined(method, :parameterOrder))"
      # end
      # for n in 1:length(optionalParameters)
      #   @info "$(service.name).$(resource.name).$(method.name) optionalParameters-$n: $(optionalParameters[n].name) parameter.required:$(optionalParameters[n].required) isdefined(method, :parameterOrder):$(isdefined(method, :parameterOrder))"
      # end
      # add required parameters to method
      for rp in sort!(requiredParameters)
        prmSym = typeSymbol(rp)
        if prmSym != nothing
          push!(eFn.args[1].args, unblock(@q $(Symbol(rp.name))::$prmSym  ))
        else
          push!(eFn.args[1].args, unblock(@q $(Symbol(rp.name))))
        end
      end
      # add optional parameters to method
      if length(optionalParameters) != 0
        for op in optionalParameters
          push!(eFn.args[1].args[2].args, unblock(@q $(Symbol(op.name))=nothing ))
        end
      end
    end # isdefined(method, :parameters)
    # add keyword args parameter
    push!(eFn.args[1].args[2].args, unblock(@q kws... ))

    # --- Generate method block
    path = method.path
    path = replace(path, "+" => "") # bigquery/models may have '+'
    path = replace(path, "{projectId}" => "\$(project.id)")
    path = replace(path, "{" => "\$(")
    path = replace(path, "}" => ")")

    url = Symbol("\"\$($(uppercasefirst(service.name)).baseUrl)$(path)\"")
    push!(eFn.args[2].args, unblock(@q headers = project.authorize(; kws...) ))
    push!(eFn.args[2].args, unblock(@q res = HTTP.request($(method.httpMethod), $url, headers; kws...) ))

    # add body serialization
    if isdefined(method, :request)
      insert!(eFn.args[2].args, 2, :(headers["Content-Type"] = "application/json"))
      # @info "eFn.args[2].args[3].args[2].args", eFn.args[2].args[3].args[2].args
      push!(eFn.args[2].args[3].args[2].args, :(GCP.json(body)) )
    end

    # add response object deserialization (if any)
    if !isdefined(method, :response)
      push!(eFn.args[2].args, unblock(@q nothing ))
    else
      resSym = typeSymbol(method.response)
      if resSym == nothing
        push!(eFn.args[2].args, unblock(@q nothing ))
      else
        if haskey(service.schemas, String(resSym)) # GCP.Abstract
          resSym = Symbol("$(uppercasefirst(service.name)).$resSym")
        end
        push!(eFn.args[2].args, unblock(@q convert($resSym , LazyJSON.value(String(res.body)) ) ))

        # @info "$fieldName::$resSym hasKey", haskey(service.schemas, String(resSym))
        # if occursin("Union", String(resSym)) # primitive union
        #   m = match(r"^Union{Nothing,([\w\d\s]+)}$", String(resSym))
        #   resSym = Symbol(strip(m.captures[1], [' ']))
        #   push!(eFn.args[2].args, unblock(@q parse($resSym , LazyJSON.value(String(res.body)) ) ))
        # else
        #   if haskey(service.schemas, String(resSym)) # GCP.Abstract
        #     resSym = Symbol("$(uppercasefirst(service.name)).$resSym")
        #   end
        #   push!(eFn.args[2].args, unblock(@q convert($resSym , LazyJSON.value(String(res.body)) ) ))
        # end

        # if resSym == :String || haskey(service.schemas, String(resSym))
        #   if haskey(service.schemas, String(resSym)) # GCP.Abstract
        #     resSym = Symbol("$(uppercasefirst(service.name)).$resSym")
        #   end
        #   push!(eFn.args[2].args, unblock(@q convert($resSym , LazyJSON.value(String(res.body)) ) ))
        # else # primitive: Union{Nothing, Int64} etc...
        #   push!(eFn.args[2].args, unblock(@q parse($resSym , LazyJSON.value(String(res.body)) ) ))
        # end

      end
    end

    eFn
  end

  function Base.convert(::Type{Expr}, resource::Resource, service::Service)
    # resource module
    resource.expr = unblock(@q module $(Symbol(uppercasefirst(resource.name))) end)
    # resource module usings
    push!(resource.expr.args[3].args, unblock(@q using Dates ))
    push!(resource.expr.args[3].args, unblock(@q using LazyJSON ))
    push!(resource.expr.args[3].args, unblock(@q using JSON ))
    push!(resource.expr.args[3].args, unblock(@q using HTTP ))
    push!(resource.expr.args[3].args, unblock(@q using DataFrames )) # for Tabledata.list
    push!(resource.expr.args[3].args, unblock(@q using $(Symbol(".."*uppercasefirst(service.name))) ))
    push!(resource.expr.args[3].args, unblock(@q using ...GCP ))
    # add resource struct with method
    eResourceStruct = unblock(@q mutable struct Resource end)
    push!(resource.expr.args[3].args, eResourceStruct)
    if isdefined(resource, :methods)
      # add struct fields for each method
      for method in sort(resource.methods)
        push!(eResourceStruct.args[3].args, unblock(@q $(Symbol(method.name)) ) )
      end
    end
    # add resource constructor
    eResourceCtor = unblock(@q function Resource(project) r = new() end )
    push!(eResourceStruct.args[3].args, eResourceCtor )
    # add method definitions
    for method in sort(resource.methods)
      method.expr = convert(Expr, method, resource, service)
      eMethod = unblock(@q begin
        r.$(Symbol(method.name)) = $(method.expr)
      end)
      push!(eResourceCtor.args[2].args, eMethod)
    end
    push!(eResourceCtor.args[2].args, unblock(@q r))

    resource.expr
  end

  function Base.convert(::Type{Expr}, service::Service)
    # service module
    service.expr = unblock(@q module $(Symbol(uppercasefirst(service.name))) end)
    # service module usings
    push!(service.expr.args[3].args, unblock(@q using Dates ))
    push!(service.expr.args[3].args, unblock(@q using LazyJSON ))
    push!(service.expr.args[3].args, unblock(@q using ..GCP ))
    # service consts
    push!(service.expr.args[3].args, unblock(@q const baseUrl = $(service.baseUrl) ))
    # --- Generate empty struct expressions for each schema
    for schema in values(service.schemas)
      schema.expr = convert(Expr, schema)
    end
    # --- Generate struct fields
    for schema in values(service.schemas)
      # @info "schema.name", schema.name
      if isdefined(schema, :properties)
        for property in schema.properties
          # @info "property", property
          # if isreserved(property.name)
          #   @info "$(schema.name).$(property.name) SKIPPING"
          #   continue
          # end
          propertySchema = convert(Expr, property, schema, service)
          if isnothing(propertySchema)
            @info "$(schema.name).$(property.name) NOTHING"
            continue
          end
          push!(schema.expr.args[3].args, propertySchema)
        end
      end
      # --- Generate struct constructor (for mutable)
      push!(schema.expr.args[3].args, unblock(@q function $(Symbol(schema.name))() new() end ))
    end

    # --- Add schema structs to module

     # n-ref structs last; sort by references
    # topological sort
    stack = Vector{Schema}()
    for schema in filter(v -> length(v.refs) != 0, collect(values(service.schemas)))
      # @info "try sort", schema.name
      topologicalSort(schema, stack)
      # @info "===> sorting:$(schema.name) refs:$(length(schema.refs))"
    end

    # 0-ref structs first; maybe be referenced by other structs
    for schema in sort(filter(v -> length(v.refs) == 0, collect(values(service.schemas))))
      if !in(schema, stack)
        push!(service.expr.args[3].args, schema.expr)
        # @info "===> 0refs:$(schema.name)"
      end
    end

    for schema in stack
      push!(service.expr.args[3].args, schema.expr)
      # @info "===> stack:$(schema.name)"
    end

    # --- Generate convert methods for each struct. Variants include:
    #   Base.convert(::Type{T}, x::LazyJSON.Object{T, Nothing})
    #   TODO Base.convert(::Type{Union{Nothing, T}}, x::LazyJSON.Object{Nothing,T})
    for key in sort(collect(keys(service.schemas)))
      schema = service.schemas[key]
      type = Symbol(schema.name)
      # Generate method stub
      eConvert = unblock(@q begin
        function Base.convert(::Type{$type}, x::LazyJSON.Object{Nothing,String})
          if x == nothing return nothing end
          r = $type()
        end
      end)
      # Generate set field statements
      # @info "$(schema.name)", schema.expr.args[3].args
      for eField in schema.expr.args[3].args[1:end-1] # skip last arg which is constructor
        # @info "$(schema.name).eField", eField
        fieldName = eField.args[1]
        fieldType = eField.args[2]
        fieldNameStr = String(fieldName)
        # reserved word suffix has '_' at field name end
        fieldNameStr = replace(fieldNameStr, r"_$" => "")
        eField = unblock(@q begin
          if haskey(x, $fieldNameStr)
            r.$fieldName = convert($fieldType, x[$fieldNameStr])
          end
        end)
        # add field set to method
        push!(eConvert.args[2].args, eField)
      end # end convert
      # add return statement
      push!(eConvert.args[2].args, unblock(@q r ))
      # add convert method to service module
      push!(service.expr.args[3].args, eConvert)
    end

    # --- Generate convert methods for each struct. Variants include:
    #   Base.convert(::Type{Vector{T}}, x::LazyJSON.Array{Nothing,T})
    for key in sort(collect(keys(service.schemas)))
      schema = service.schemas[key]
      type = Symbol(schema.name)
      # Generate method stub
      eConvert = unblock(@q begin
        function Base.convert(::Type{Vector{$type}}, x::LazyJSON.Array{Nothing,String})
          if x == nothing return nothing end
          r = Vector{$type}()
          for v in x
            push!(r, convert($type, v))
          end
        end
      end)
      # add return statement
      push!(eConvert.args[2].args, unblock(@q r ))
      # add convert method to service module
      push!(service.expr.args[3].args, eConvert)
    end

    # --- Generate Resource modules
    for resource in service.resources
      # @info "resource.name", resource.name
      push!(service.expr.args[3].args, convert(Expr, resource, service))
    end

    # --- Generate Service struct
    # add usings for Service
    for resource in service.resources
      push!(service.expr.args[3].args, unblock(@q using .$(Symbol(uppercasefirst(resource.name)))) )
    end
    # create Service struct
    eService = unblock(@q mutable struct Service end)
    push!(service.expr.args[3].args, eService)
    # add resource fields to Service struct
    for resource in service.resources
      fieldName = Symbol(resource.name)
      fieldType = Symbol(uppercasefirst(resource.name))
      push!(eService.args[3].args, unblock(@q $fieldName::$fieldType.Resource) )
    end # for: eService: fields
    # add constructor for Service struct
    eResourceCtor = unblock(@q function Service(project) s = new() end)
    push!(eService.args[3].args, eResourceCtor)
    # initialize resource fields in constructor
    for resource in service.resources
      varName = Symbol(resource.name)
      moduleName = Symbol(uppercasefirst(resource.name))
      eVar = unblock(@q s.$varName = $moduleName.Resource(project) )
      push!(eResourceCtor.args[2].args, eVar)
    end # for: eResourceCtor: init fields
    push!(eResourceCtor.args[2].args, unblock(@q s))

    service.expr
  end

  function emitService(service::Service, path::String)
    @info "emitService: $(service.name): start ---"

    eServiceModule = convert(Expr, service)

    # --- Write Service module file to disk ---
    filename = "$(path)$(service.name).jl"
    @info "emitService: $(service.name): write $(filename)"
    open(filename, "w") do f # write file
      write(f, exprToString(eServiceModule))
    end

    @info "emitService: $(service.name): done ---"
  end

end # module Emit

end # module Bootstrap

max_stars_repo_name: ProfLeao/mais

max_issues_repo_name: ProfLeao/mais

max_forks_repo_name: ProfLeao/mais


text:
ds = cf.MovieLens()
@test ds.users == 943
@test ds.items == 1682
@test length(ds.preferences.possibles) == 5
@test ds.preferences.min == 1
@test ds.preferences.max == 5

ds = cf.FilmTrust()
@test ds.users == 1508
@test ds.items == 2071
@test length(ds.preferences.possibles) == 8
@test ds.preferences.min == 0.5
@test ds.preferences.max == 4.0

ds = cf.CiaoDVD()
@test ds.users == 17615
@test ds.items == 16121
@test length(ds.preferences.possibles) == 5
@test ds.preferences.min == 1
@test ds.preferences.max == 5

ds = cf.MovieTweeting()
@test ds.users == 25011
@test ds.items == 14732
@test length(ds.preferences.possibles) == 11
@test ds.preferences.min == 0
@test ds.preferences.max == 10


text:
function Base.:*(A::AbstractArray, s::Space) where {Space}

    @assert ndims(A) == 2
    @assert size(A,2) == numfunctions(s)

    F = eltype(s.fns)
    fns = [similar(F,0) for i in axes(A,1)]
    geo = s.geo

    for i in axes(A,1)
        for j in axes(A,2)
            A[i,j] == 0 && continue
            shapes = [Shape(sh.cellid, sh.refid, A[i,j]*sh.coeff) for sh in s.fns[j]]
            append!(fns[i], shapes)
        end
        fns[i] = collapse_shapes(fns[i])
    end

    pos = fill(zero(eltype(s.pos)), length(fns))
    return Space(geo, fns, pos)
end

function collapse_shapes(shapevec)
    ids = unique([(sh.cellid,sh.refid) for sh in shapevec])

    S = eltype(shapevec)
    collapsed = [S(id..., 0) for id in ids]

    for sh in shapevec
        pos = findfirst(isequal((sh.cellid, sh.refid)), ids)
        @assert pos != nothing
        old = collapsed[pos]
        @assert old.cellid == sh.cellid
        @assert old.refid == sh.refid
        collapsed[pos] = S(old.cellid, old.refid, old.coeff+sh.coeff)
    end

    collapsed
end


text:
# Inlining calls to zchop for reals into zchop for complex
# improves efficiency in this case, by about .5 %
# Profile shows that Julia does not inline the calls.

# Julia v0.4 was not inlining the call to zchop with complex args
# forcing inline increased speed in loop by 30%

using ZChop

const zeps = 1e-14

zchop2{T<:Real}(x::T, eps=zeps) = abs(x) > convert(T,eps) ? x : zero(T)

function zchop2{T<:Real}(x::Complex{T}, eps=zeps)
    rex = real(x)
    iex = imag(x)
    rex =  abs(rex) > convert(T,eps) ? rex : zero(T)
    iex =  abs(iex) > convert(T,eps) ? iex : zero(T)
    complex(rex,iex)
end

zchop2(a::AbstractArray, eps=zeps) = (b = similar(a); for i in 1:length(a)  b[i] = zchop2(a[i],eps) end ; b)
zchop2(x,eps=zeps) = map((x)->zchop2(x,eps),x)

const n = 10^6
const z = Array(Complex{Float64},n)
for i in 1:n z[i] = complex(.1,.1) end

function testzchop()
    m = 500
    sum = 0
    for i in 1:m
        sum += @elapsed zchop(z)
    end
    println("call ", sum/m)
end

function testzchop2()
    m = 500
    sum = 0
    for i in 1:m
        sum += @elapsed zchop2(z)
    end
    println("expl ", sum/m)
end

function runalt()
    for j in 1:100
        testzchop()
        testzchop2()
    end
end

# Return a tuple of n copies of x. copy of bits that is.
# n limited to about 10^5
function mktuple(n,x)
    ex = Expr(:tuple,)
    a = ex.args
    for i in 1:n
        push!(a,x)
    end
    eval(ex)
end

true


text:
module Generate

# using PyCall

# const map   = pyimport("mpl_toolkits.basemap")
# const Globe = map.Basemap()

"""
    sphere(N; R=1)

Generate a spherical point cloud of `N` points with extent of radius `R`.
"""
function sphere(N; R=1)
    θ = π  .* rand(Float64, N)
    ϕ = 2π .* rand(Float64, N)

    return R.*hcat(sin.(θ).*cos.(ϕ), sin.(θ).*sin.(ϕ), cos.(θ))', hcat(θ, ϕ)'
end

function spherical_distance(x; R=1)
    n̂ = x / R
    D = zeros(size(x,2),size(x,2))

    norm² = sum(n̂[i,:].^2 for i ∈ 1:size(x,1))
    cosΔψ = sum(n̂[i,:]' .* n̂[i,:] for i ∈ 1:size(x,1))
    sinΔψ = [norm(cross(n̂[:,i], n̂[:,i])) for i ∈ 1:size(x,2), j ∈ 1:size(x,2)]

    return R*atan.(sinΔψ, cosΔψ)
end

#=
"""
    globe(N)

Generate a spherical point cloud of `N` points where all points are guaranteed to be sampled from land-masses from Earth.
Uses mpl.basemap Python library.
"""
function globe(N)
    Θ = Array{Float64}(undef, N)
    Φ = Array{Float64}(undef, N)
    for n in 1:N
    @label sample
        θ = π*rand(1)[1]
        ϕ = 2π*rand(1)[1]

        x, y = Globe(180*(ϕ.-π)/π, 180*(π/2 .- θ)/π)

        !Globe.is_land(x,y) && @goto sample

        Θ[n] = θ
        Φ[n] = ϕ
    end

    return hcat(sin.(Θ).*cos.(Φ), sin.(Θ).*sin.(Φ), cos.(Θ))', hcat(Θ, Φ)'
end
=#

"""
    swissroll(N; z₀=10, R=1/20)

Generate a point cloud of `N` distributed on a swiss roll manifold with unit radius and length ``\frac{z₀}{R}``.
"""
function swissroll(N; z₀=10, R=1/20)
    z = (z₀/R)*rand(Float64, N)
    ϕ = 1.5π .+ 3π .* rand(Float64, N)

    return hcat(ϕ .* cos.(ϕ), ϕ .* sin.(ϕ), z)' .* R, hcat(ϕ, z)'
end

"""
    torus(N; R=2, r=1)

Generate a point cloud of `N` distributed on a torus, sized inner `r` and outer radius `R` respectively.
"""
function torus(N; R=2, r=1)
    θ = 2π  .* rand(Float64, N)
    ϕ = 2π .* rand(Float64, N)

    return hcat((R .+ r*cos.(θ)) .* cos.(ϕ), (R .+ r*cos.(θ)).* sin.(ϕ), r*sin.(θ))', hcat(ϕ, θ)'
end

end

max_stars_repo_name: nnoll/seqspace

max_issues_repo_name: nnoll/seqspace

max_forks_repo_name: nnoll/seqspace


text:
## Functions for abstract domains from AbstractDomains.jl for inference
## ====================================================================

for finame in ["interval.jl",
               "boxes.jl"]
    include(joinpath("domains", finame))
end


text:
mutable struct FASTQRecord
   sequence::BioSequence{DNAAlphabet{2}}
   quality::Vector{UInt8}
   raw::FASTQ.Record

   FASTQRecord() = new( BioSequence{DNAAlphabet{2}}(),
                        Vector{UInt8}(), FASTQ.Record() )
end

Base.read!( parser::FASTQ.Reader, rec::FASTQRecord ) = read!( parser, rec.raw )

function Base.fill!( rec::FASTQRecord, offset=33 )
   rec.sequence = BioSequences.BioSequence{BioSequences.DNAAlphabet{2}}( rec.raw.data,
                                                                         first(rec.raw.sequence),
                                                                         last(rec.raw.sequence) )
   rec.quality  = rec.raw.data[rec.raw.quality] .- convert(UInt8, offset)
   rec
end

Base.string( rec::FASTQRecord ) = String(rec.raw.data[rec.raw.identifier])

function reverse_complement!(rec::FASTQRecord)
   BioSequences.reverse_complement!(rec.sequence)
   reverse!(rec.quality)
end


text:
# This needs actual performance testing.
@generated function Chunk(::Val{N}) where N
    if 4N <= VectorizationBase.REGISTER_SIZE
        P = N
    else
        p1 = round(Int, 4N / VectorizationBase.REGISTER_SIZE, RoundUp)
        P = round(Int, N / p1, RoundUp)
    end
    ForwardDiff.Chunk{P}()
end
@generated ValP1(::Val{N}) where N = Val{N+1}()
@generated ValM1(::Val{N}) where N = Val{N-1}()

mutable struct GradientDiffResult{V,P,R} <: DiffResults.DiffResult{1,V,Tuple{ConstantFixedSizeVector{P,V,R,R}}}
    grad::ConstantFixedSizeVector{P,V,R,R}
    value::V
    function GradientDiffResult{V,P,R}(::UndefInitializer) where {V, P, R}
        new{V,P,R}()
    end
    @generated function GradientDiffResult{V,P}(::UndefInitializer) where {V, P}
        R = PaddedMatrices.calc_padding(P, V)
        :(GradientDiffResult{$V,$P,$R}(undef))
    end
    function GradientDiffResult(g::ConstantFixedSizeVector{P,V,R,R}, v::V = zero(V)) where {P,V,R}
        GradientDiffResult{P,V,R}(g, v)
    end
end
mutable struct HessianDiffResult{V,P,R,L} <: DiffResults.DiffResult{2,V,Tuple{ConstantFixedSizeVector{P,V,R,R},ConstantFixedSizeMatrix{P,P,V,R,L}}}
    grad::ConstantFixedSizeVector{P,V,R,R}
    hess::ConstantFixedSizeMatrix{P,P,V,R,L}
    value::V
    function HessianDiffResult{V,P,R,L}(::UndefInitializer) where {V,P,R,L}
        new{V,P,R,L}()
    end
    @generated function HessianDiffResult{V,P}(::UndefInitializer) where {V, P}
        N,R,L = PaddedMatrices.calc_NPL((P,P), V)
        :(HessianDiffResult{$V,$P,$R,$L}(undef))
    end
    function GradientDiffResult(g::ConstantFixedSizeVector{P,V,R,R}, H::ConstantFixedSizeMatrix{P,P,V,R,L}, v::V = zero(V)) where {P,V,R,L}
        GradientDiffResult{P,V,R,R}(g, H, v)
    end
end
const GradResult{V,P,R} = Union{GradientDiffResult{V,P,R},HessianDiffResult{V,P,R}}
@inline Base.pointer(result::GradResult{V}) where {V} = Base.unsafe_convert(Ptr{V}, pointer_from_objref(result))
@inline DiffResults.derivative(result::GradResult{V,P,R}) where {V,P,R} = PtrVector{P,V,R}(pointer(result))
@inline DiffResults.derivative(result::GradResult{V,P,R}, ::Type{Val{1}}) where {V,P,R} = PtrVector{P,V,R}(pointer(result))
@inline DiffResults.derivative(result::HessianDiffResult{V,P,R}, ::Type{Val{2}}) where {V,P,R} = PtrVector{P,V,R}(pointer(result) + sizeof(V)*R )
@inline DiffResults.gradient(result::GradResult{V,P,R}) where {V,P,R} = PtrVector{P,V,R}(pointer(result))
@inline DiffResults.hessian(result::HessianDiffResult{V,P,R})  where {V,P,R}= PtrMatrix{P,P,V,R}(pointer(result) + sizeof(V)*R )

@inline DiffResults.value!(r::GradResult, x::Number) = (r.value = x; return r)
# @inline DiffResults.value!(r::GradResult, x::AbstractArray) = (copyto!(value(r), x); return r)
@inline DiffResults.value!(f::F, r::GradResult, x::Number) where {F} = (r.value = f(x); return r)
# @inline DiffResults.value!(f, r::GradResult, x::AbstractArray) = (map!(f, value(r), x); return r)

function DiffResults.derivative!(r::GradResult, x::AbstractArray)
    copyto!(DiffResults.gradient(r), x)
    return r
end
function DiffResults.derivative!(r::GradResult, x::AbstractArray, ::Type{Val{1}})
    copyto!(DiffResults.gradient(r), x)
    return r
end
function DiffResults.derivative!(r::HessianDiffResult, x::AbstractArray, ::Type{Val{2}})
    copyto!(DiffResults.hessian(r), x)
    return r
end
function DiffResults.derivative!(r::GradResult{P,V,R}, x::ConstantFixedSizeVector{P,V,R}) where {P,V,R}
    r.grad = x
    return r
end
function DiffResults.derivative!(r::GradResult{P,V,R}, x::ConstantFixedSizeVector{P,V,R}, ::Type{Val{1}}) where {P,V,R}
    r.grad = x
    return r
end
function DiffResults.derivative!(r::HessianDiffResult{P,V,R,L}, x::ConstantFixedSizeMatrix{P,P,V,R,L}, ::Type{Val{2}}) where {P,V,R,L}
    r.hess = x
    return r
end
# function DiffResults.derivative!(f, r::GradResult, x::Number)
#     r.grad = f(x)
#     return r
# end
# function DiffResults.derivative!(f, r::GradResult, x::AbstractArray)
#     map!(f, r.grad, x)
#     return r
# end

function GradientDiffResult(x::SizedVector{P,T}) where {T,P}
    GradientDiffResult{T,P}(undef)
end
function HessianDiffResult(x::SizedVector{P,T}) where {T,P}
    HessianDiffResult{T,P}(undef)
    # HessianDiffResult(zero(T), similar(x), FixedSizeMatrix{P,P,T}(undef))
end
# function HessianDiffResult(x::MVector{P,T}) where {T,P}
#     HessianDiffResult{T,P}(undef)
#     # HessianDiffResult(zero(T), similar(x), MMatrix{P,P,T}(undef))
# end

abstract type AbstractAutoDiffDifferentiable{P,T} <: AbstractDifferentiableObject{P,T} end
abstract type AbstractOnceDifferentiableObject{P,T} <: AbstractAutoDiffDifferentiable{P,T} end
abstract type AbstractTwiceDifferentiableObject{P,T} <: AbstractAutoDiffDifferentiable{P,T} end
abstract type Configuration{P,V,F} end

struct GradientConfiguration{P,V,F,T,ND,DG,R} <: Configuration{P,V,F}
    f::F
    result::GradientDiffResult{V,P,R}
    gconfig::ForwardDiff.GradientConfig{T,V,ND,DG}
end

struct HessianConfiguration{P,V,F,T,T2,ND,DJ,DG,DG2, R_V, R_D, LSQM} <: Configuration{P,V,F}
    f::F
    result::HessianDiffResult{V,P,R_V,LSQM}
    inner_result::GradientDiffResult{ForwardDiff.Dual{T,V,ND},P,R_D}
    jacobian_config::ForwardDiff.JacobianConfig{T,V,ND,DJ}
    gradient_config::ForwardDiff.GradientConfig{T,ForwardDiff.Dual{T,V,ND},ND,DG}
    gconfig::ForwardDiff.GradientConfig{T2,V,ND,DG2}
end

struct TwiceDifferentiable{P,T,F,C<:Configuration{P,T,F}, SV_V <: SizedVector{P,T}} <: AbstractTwiceDifferentiableObject{P,T}
    x_f::SV_V # x used to evaluate f (stored in F)
    x_df::SV_V # x used to evaluate df (stored in DF)
    x_h::SV_V #??
    config::C
end
struct OnceDifferentiable{P,T,F,C<:GradientConfiguration{P,T,F}, SV_V <: SizedVector{P,T}} <: AbstractOnceDifferentiableObject{P,T}
    x_f::SV_V# x used to evaluate f (stored in F)
    x_df::SV_V # x used to evaluate df (stored in DF)
    config::C
end

function GradientConfiguration(f::F, x::SizedVector{P,T}) where {F,T,P}
    result = GradientDiffResult(x)
    chunk = Chunk(Val{P}())
    tag = ForwardDiff.Tag(f, T)
    gconfig = ForwardDiff.GradientConfig(f, x, chunk, tag)

    GradientConfiguration(f, result, gconfig)
end
extract_tag(::T) where {T} = T
function HessianConfiguration(f::F, x::SizedVector{P,T}) where {F,T,P}
    result = HessianDiffResult(x)
    chunk = Chunk(Val{P}())
    # tag = ForwardDiff.Tag(f, T)

    # here we construct a JacobianConfig so that yduals is a PtrVector to DiffResults.gradient(inner_result)
    tag_instance = ForwardDiff.Tag(f, T)
    tag_type = extract_tag(tag_instance)
    # @show tag
    seeds = ForwardDiff.construct_seeds(ForwardDiff.Partials{P,T})
    # inner_result = GradientDiffResult(jacobian_config.duals[2])
    inner_result = GradientDiffResult{ForwardDiff.Dual{tag_type,T,P},P}(undef)
    duals = (DiffResults.gradient(inner_result), similar(x, ForwardDiff.Dual{tag_type,T,P}))
    jacobian_config = ForwardDiff.JacobianConfig{tag_type,T,P,typeof(duals)}(seeds, duals)

    # jacobian_config = ForwardDiff.JacobianConfig((f,ForwardDiff.gradient), DiffResults.gradient(result), x, chunk, tag)
    # jacobian_config = ForwardDiff.JacobianConfig((f,ForwardDiff.gradient), DiffResults.gradient(result), x, chunk, tag)
    gradient_config = ForwardDiff.GradientConfig(f, jacobian_config.duals[2], chunk, tag_instance)

    gconfig = ForwardDiff.GradientConfig(f, x, chunk, tag_instance)

    HessianConfiguration(f, result, inner_result, jacobian_config, gradient_config, gconfig)
end

TwiceDifferentiable(f::F, ::Val{P}) where {F,P} = TwiceDifferentiable(f, FixedSizeVector{P,Float64}(undef))
function TwiceDifferentiable(f::F, x::SV_V) where {F,T,P,SV_V <: SizedVector{P,T}}
    TwiceDifferentiable(x, similar(x), similar(x), HessianConfiguration(f, x))
end
function TwiceDifferentiable(x_f::SV_V,x_df::SV_V,x_h::SV_V,config::C) where {T,F,C<:Configuration{T,F},P,SV_V <: SizedVector{P,T}}
    TwiceDifferentiable{P,T,F,C}(x_f, x_df, x_h, config)
end

function OnceDifferentiable(f, x_f::SizedVector{P,T}) where {T,P}
    OnceDifferentiable(x_f, similar(x_f), GradientConfiguration(f, x_f))
end
function OnceDifferentiable(x_f::SV_V, config::C) where {T,F,C<:GradientConfiguration{T,F},P,SV_V <: SizedVector{P,T}}
    OnceDifferentiable{P,T,F,C}(x_f, similar(x_f), config)
end

# ProfileDifferentiable(f::F, ::Val{N}) where {F,N} = ProfileDifferentiable(f, Vector{Float64}(undef, N), Val{N}())
# function ProfileDifferentiable(f::F, x::AbstractArray{T}, ::Val{N}) where {F,T,N}
#     x_f = Vector{T}(undef, N-1)
#     ProfileDifferentiable(x_f, similar(x_f), similar(x_f), Configuration(f, x_f, ValM1(Val{N}())), x, Ref(N), Ref{T}(),::Val{N})
# end
# function ProfileDifferentiable(x_f::Vector{T},x_df::Vector{T},x_h::Vector{T}, config::C, x,::A i::RefValue{Int}, v::RefValue{T},::Val{N}) where {T,A<:AbstractArray{T},C,N}
#     ProfileDifferentiable{N,T,A,C}(x_f, x_df, x_h, config, x, i, v)
# end

@inline DiffResults.value!(obj::AbstractAutoDiffDifferentiable, x::Real) = DiffResults.value!(obj.config.result, x)

@inline value(obj::AbstractAutoDiffDifferentiable) = obj.config.result.value
@inline gradient(obj::AbstractAutoDiffDifferentiable) = DiffResults.gradient(obj.config.result)
# @inline gradient(obj::AbstractAutoDiffDifferentiable, i::Integer) = DiffResults.gradient(obj.config.result)[i]
@inline hessian(obj::TwiceDifferentiable) = DiffResults.hessian(obj.config.result)

@inline f(obj::AbstractAutoDiffDifferentiable, x) = obj.config.f(x)

# """
# For DynamicHMC support. Copies data.
# """
# @inline function (obj::OnceDifferentiable)(x)
#     fdf(obj, x)
#     DiffResults.ImmutableDiffResult(-obj.config.result.value, (-1 .* obj.config.result.grad,))
# end
# @inline function (obj::TwiceDifferentiable)(x)
#     fdf(obj, x)
#     DiffResults.ImmutableDiffResult(-obj.config.result.value, (-1 .* obj.config.result.grad,))
# end

@inline function df(obj::AbstractAutoDiffDifferentiable, x)
    ForwardDiff.gradient!(gradient(obj), obj.config.f, x, obj.config.gconfig, Val{false}())
end

@inline function fdf(obj::AbstractAutoDiffDifferentiable, x)
    # obj.config.result.grad = gradient(obj)
    ForwardDiff.gradient!(obj.config.result, obj.config.f, x, obj.config.gconfig, Val{false}())
    # DiffResults.value(obj.config.result)
    obj.config.result.value
end
@inline function scale_fdf(obj::AbstractAutoDiffDifferentiable, x, scale_target)
    # obj.config.result.grad = gradient(obj)
    scale = scale_gradient!(obj.config.result, obj.config.f, x, scale_target, obj.config.gconfig, Val{false}())
    # DiffResults.value(obj.config.result)
    obj.config.result.value, scale
end
@inline function scaled_fdf(obj::AbstractAutoDiffDifferentiable, x, scale)
    # obj.config.result.grad = gradient(obj)
    scaled_gradient!(obj.config.result, obj.config.f, x, scale, obj.config.gconfig, Val{false}())
    # DiffResults.value(obj.config.result)
    obj.config.result.value
end

# function Configuration(f::F, x::AbstractArray{V}, chunk::Chunk = Chunk(x), tag = Tag(f, V)) where {F,V}
#     result = DiffResults.HessianResult(x)
#     jacobian_config = ForwardDiff.JacobianConfig((f,gradient), DiffResults.gradient(result), x, chunk, tag)
#     gradient_config = ForwardDiff.GradientConfig(f, jacobian_config.duals[2], chunk, tag)
#     inner_result = DiffResults.DiffResult(zero(eltype(jacobian_config.duals[2])), jacobian_config.duals[2])
#     gconfig = ForwardDiff.GradientConfig(f, x, chunk, tag)
#     Configuration(f, result, inner_result, jacobian_config, gradient_config, gconfig)
# end

@inline function (c::HessianConfiguration)(y, z)
    # c.inner_result.grad = y #Already true?
    # @show typeof(y), typeof(z)
    ForwardDiff.gradient!(c.inner_result, c.f, z, c.gradient_config, Val{false}())
    DiffResults.value!(c.result, ForwardDiff.value(DiffResults.value(c.inner_result)))
    # copyto!(y, DiffResults.gradient(c.inner_result))
    y
end
@inline function hessian!(c::HessianConfiguration, x::AbstractArray)
    ForwardDiff.jacobian!(DiffResults.hessian(c.result), c, DiffResults.gradient(c.result), x, c.jacobian_config, Val{false}())
    # DiffResults.hessian(c.result)
    DiffResults.hessian(c.result)
end

"""
Force (re-)evaluation of the objective value at `x`.

Returns `f(x)` and stores the value in `obj.F`
"""
@inline function value!!(obj::AbstractAutoDiffDifferentiable, x)
    # obj.f_calls .+= 1
    # copyto!(obj.x_f, x)
    DiffResults.value!(obj, f(obj, x) )
    value(obj)
end
# """
# Evaluates the objective value at `x`.
#
# Returns `f(x)`, but does *not* store the value in `obj.F`
# """
# @inline function value(obj::AbstractAutoDiffDifferentiable, x)
#     if x != obj.x_f
#         # obj.f_calls .+= 1
#         value!!(obj, x)
#     end
#     value(obj)
# end
"""
Evaluates the objective value at `x`.

Returns `f(x)` and stores the value in `obj.F`
"""
@inline function value!(obj::AbstractAutoDiffDifferentiable, x)
    # if x != obj.x_f
        value!!(obj, x)
    # end
    value(obj)
end

"""
Evaluates the gradient value at `x`

This does *not* update `obj.DF`.
"""
@inline function gradient(obj::AbstractAutoDiffDifferentiable, x)
    DF = gradient(obj)
    gradient!!(obj, x)
    DF
end
"""
Evaluates the gradient value at `x`.

Stores the value in `obj.DF`.
"""
@inline function gradient!(obj::AbstractAutoDiffDifferentiable, x)
    # if x != obj.x_df
        gradient!!(obj, x)
    # end
    gradient(obj)
end
"""
Force (re-)evaluation of the gradient value at `x`.

Stores the value in `obj.DF`.
"""
@inline function gradient!!(obj::AbstractAutoDiffDifferentiable, x)
    # obj.df_calls .+= 1
    copyto!(obj.x_df, x)
    df(obj, x)
end

@inline function value_gradient!(obj::AbstractAutoDiffDifferentiable, x)
    # if x != obj.x_f && x != obj.x_df
        value_gradient!!(obj, x)
    # elseif x != obj.x_f
    #     value!!(obj, x)
    # elseif x != obj.x_df
    #     gradient!!(obj, x)
    # end
    value(obj)
end
@inline function value_gradient!!(obj::AbstractAutoDiffDifferentiable, x)
    # obj.f_calls .+= 1
    # obj.df_calls .+= 1
    # copyto!(obj.x_f, x)
    # copyto!(obj.x_df, x)
    DiffResults.value!(obj, fdf(obj, x))
end

@inline function hessian!(obj::AbstractAutoDiffDifferentiable, x)
    # if x != obj.x_h
        hessian!!(obj, x)
    # end
    nothing
end
@inline function hessian!!(obj::AbstractAutoDiffDifferentiable, x)
    # obj.h_calls .+= 1
    # copyto!(obj.x_h, x)
    hessian!(obj.config, x)
end

# function ForwardDiff.vector_mode_dual_eval(f::F, x::FixedSizeVector, cfg::Union{JacobianConfig,GradientConfig}) where F
#     xdual = cfg.duals
#     seed!(xdual, x, cfg.seeds)
#     return f(xdual)
# end

@inline function ForwardDiff.seed!(duals::AbstractMutableFixedSizeVector{P,ForwardDiff.Dual{T,V,N}}, x,
                                    seed::ForwardDiff.Partials{N,V}) where {T,V,N,P}
    @inbounds for i in 1:N
        duals[i] = ForwardDiff.Dual{T,V,N}(x[i], seed)
    end
    return duals
end

@inline function ForwardDiff.seed!(duals::AbstractMutableFixedSizeVector{P,ForwardDiff.Dual{T,V,N}}, x,
                                    seeds::NTuple{N,ForwardDiff.Partials{N,V}}) where {T,V,N,P}
    @inbounds for i in 1:N
        duals[i] = ForwardDiff.Dual{T,V,N}(x[i], seeds[i])
    end
    return duals
end

@inline function ForwardDiff.seed!(duals::AbstractMutableFixedSizeVector{P,ForwardDiff.Dual{T,V,N}}, x, index,
                                    seed::ForwardDiff.Partials{N,V} = zero(ForwardDiff.Partials{N,V})) where {T,V,N,P}
    offset = index - 1
    @inbounds for i in 1:N
        j = i + offset
        duals[j] = ForwardDiff.Dual{T,V,N}(x[j], seed)
    end
    return duals
end

@inline function ForwardDiff.seed!(duals::AbstractMutableFixedSizeVector{P,ForwardDiff.Dual{T,V,N}}, x, index,
                                    seeds::NTuple{N,ForwardDiff.Partials{N,V}}, chunksize = N) where {T,V,N,P}
    offset = index - 1
    @inbounds for i in 1:chunksize
        j = i + offset
        duals[j] = ForwardDiff.Dual{T,V,N}(x[j], seeds[i])
    end
    return duals
end
function ForwardDiff.extract_jacobian!(::Type{T}, result::AbstractMutableFixedSizeMatrix{M,N}, ydual::AbstractArray, n) where {M,N,T}
    # out_reshaped = reshape(result, length(ydual), n)
    # @inbounds for col in 1:size(out_reshaped, 2), row in 1:size(out_reshaped, 1)
    #     out_reshaped[row, col] = partials(T, ydual[row], col)
    # end
    @inbounds for col in 1:N, row in 1:M
        # result[row, col] = ForwardDiff.partials(T, ydual[row], col)
        result[row, col] = ForwardDiff.partials(T, ydual[col], row)
    end
    return result
end


text:
#!/usr/bin/env julia

include("HHChannels.jl")

using JSON
using SIUnits.ShortUnits
using HHChannels

radius = 18.8µm/2
area = 4*pi*radius^2

stim = Stim(10ms, 100ms, 0.1nA/area)
ts, vs = run_hh(100ms, stim=stim, sample_dt=0.025ms)

trace = Dict(
    :name => "membrane voltage",
    :sim => "numeric",
    :model => "soma",
    :units => "mV",
    :data => Dict(
        :time => map(t->t/ms, ts),
        Symbol("soma.mid") => map(v->v/mV, vs)
    )
)

println(JSON.json([trace]))

max_stars_repo_name: kabicm/arbor

max_issues_repo_name: kabicm/arbor

max_forks_repo_name: kabicm/arbor


text:
# This file is part of Kpax3. License is MIT.

function test_mcmc_gibbs_init()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=3, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  # remove unit 5 from cluster 2 and put it into another cluster
  R = [1; 1; 1; 2; 2; 3]

  i = 5
  hi = 2

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  Kpax3.initsupportgibbs!(hi, priorR, state, support)

  @test support.k == 2
  @test support.cl[1:support.k] == [1; 3]
  @test support.t[1:support.k] == [Kpax3.clusterweight(state.v[1], priorR);
                                   Kpax3.clusterweight(state.v[3], priorR)]

  nothing
end

test_mcmc_gibbs_init()

function test_mcmc_gibbs_move()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=4, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  # remove unit 5 from cluster 2 and put it into another cluster
  R = [1; 1; 1; 2; 2; 3]

  i = 5
  hi = 2
  li = 3
  lj = 4

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  support.ni = vec(float(data[:, 4]))
  support.vi = 1

  Kpax3.initsupportgibbs!(hi, priorR, state, support)

  support.t[li] = Kpax3.clusterweight(support.vi, priorR)
  support.t[lj] = Kpax3.clusterweight(1, state.k, priorR)

  g = 0
  l = 0
  y = 0.0
  v = 0
  z = zeros(Float64, 2)
  for b in 1:m
    support.lpi[1, b] = priorC.logγ[1]
    support.lpi[2, b] = priorC.logγ[2]
    support.lpi[3, b] = priorC.logγ[3]
    support.lpi[4, b] = priorC.logγ[3]

    g = 1
    l = 1
    y = state.n1s[g, b] + float(data[b, i])
    v = state.v[g] + 1

    Kpax3.gibbsupdateclustergmove!(y, v, b, l, g, state.k, priorC, support)

    @test support.lq[1, l, b] == Kpax3.logmarglik(y, v, priorC.A[1, b], priorC.B[1, b])
    @test support.lq[2, l, b] == Kpax3.logmarglik(y, v, priorC.A[2, b], priorC.B[2, b])
    @test support.lq[3, l, b] == Kpax3.logmarglik(y, v, priorC.A[3, b], priorC.B[3, b])
    @test support.lq[4, l, b] == Kpax3.logmarglik(y, v, priorC.A[4, b], priorC.B[4, b])

    g = 3
    l = 2
    y = state.n1s[g, b] + float(data[b, i])
    v = state.v[g] + 1

    Kpax3.gibbsupdateclustergmove!(y, v, b, l, g, state.k, priorC, support)

    @test support.lq[1, l, b] == Kpax3.logmarglik(y, v, priorC.A[1, b], priorC.B[1, b])
    @test support.lq[2, l, b] == Kpax3.logmarglik(y, v, priorC.A[2, b], priorC.B[2, b])
    @test support.lq[3, l, b] == Kpax3.logmarglik(y, v, priorC.A[3, b], priorC.B[3, b])
    @test support.lq[4, l, b] == Kpax3.logmarglik(y, v, priorC.A[4, b], priorC.B[4, b])

    @test isapprox(support.lpi[1, b], priorC.logγ[1] + support.lp[1, 1, b] + support.lp[1, 3, b], atol=ε)
    @test isapprox(support.lpi[2, b], priorC.logγ[2] + support.lp[2, 1, b] + support.lp[2, 3, b], atol=ε)

    z[1] = log(exp(priorC.logω[state.k][1] + support.lp[3, 1, b]) + exp(priorC.logω[state.k][2] + support.lp[4, 1, b]))
    z[2] = log(exp(priorC.logω[state.k][1] + support.lp[3, 3, b]) + exp(priorC.logω[state.k][2] + support.lp[4, 3, b]))

    @test isapprox(support.lpi[3, b], priorC.logγ[3] + z[1] + z[2], atol=ε)

    z[1] = support.lq[1, 1, b] + support.lq[1, li, b] - support.lp[1, 1, b]
    z[2] = support.lq[1, 2, b] + support.lq[1, li, b] - support.lp[1, 3, b]
    @test isapprox(support.lr[1, 1, b], z[1], atol=ε)
    @test isapprox(support.lr[1, 2, b], z[2], atol=ε)

    z[1] = support.lq[2, 1, b] + support.lq[2, li, b] - support.lp[2, 1, b]
    z[2] = support.lq[2, 2, b] + support.lq[2, li, b] - support.lp[2, 3, b]
    @test isapprox(support.lr[2, 1, b], z[1], atol=ε)
    @test isapprox(support.lr[2, 2, b], z[2], atol=ε)

    z[1] = log(exp(priorC.logω[state.k][1] + support.lq[3, 1, b]) + exp(priorC.logω[state.k][2] + support.lq[4, 1, b])) - log(exp(priorC.logω[state.k][1] + support.lp[3, 1, b]) + exp(priorC.logω[state.k][2] + support.lp[4, 1, b]))
    z[2] = log(exp(priorC.logω[state.k][1] + support.lq[3, 2, b]) + exp(priorC.logω[state.k][2] + support.lq[4, 2, b])) - log(exp(priorC.logω[state.k][1] + support.lp[3, 3, b]) + exp(priorC.logω[state.k][2] + support.lp[4, 3, b]))

    @test isapprox(support.lr[3, 1, b], z[1], atol=ε)
    @test isapprox(support.lr[3, 2, b], z[2], atol=ε)

    Kpax3.gibbsupdateclusteri!(state.k, b, li, priorC, support)

    y = support.ni[b]
    v = support.vi

    @test support.lq[1, li, b] == Kpax3.logmarglik(y, v, priorC.A[1, b], priorC.B[1, b])
    @test support.lq[2, li, b] == Kpax3.logmarglik(y, v, priorC.A[2, b], priorC.B[2, b])
    @test support.lq[3, li, b] == Kpax3.logmarglik(y, v, priorC.A[3, b], priorC.B[3, b])
    @test support.lq[4, li, b] == Kpax3.logmarglik(y, v, priorC.A[4, b], priorC.B[4, b])

    z[1] = priorC.logω[state.k][1] + support.lq[3, li, b]
    z[2] = priorC.logω[state.k][2] + support.lq[4, li, b]

    @test isapprox(support.lpj[1, b], log(exp(z[1]) + exp(z[2])),  atol=ε)

    Kpax3.gibbsupdateclusterj!(data[b, i], state.k + 1, b, li, lj, priorC, support)

    y = support.ni[b]
    v = 1

    z[1] = priorC.logω[state.k + 1][1] + Kpax3.logmarglik(y, v, priorC.A[3, b], priorC.B[3, b])
    z[2] = priorC.logω[state.k + 1][2] + Kpax3.logmarglik(y, v, priorC.A[4, b], priorC.B[4, b])

    @test isapprox(support.lpj[2, b], log(exp(z[1]) + exp(z[2])), atol=ε)

    @test support.lq[1, lj, b] == Kpax3.logmarglik(data[b, i], 1, priorC.A[1, b], priorC.B[1, b])
    @test support.lq[2, lj, b] == Kpax3.logmarglik(data[b, i], 1, priorC.A[2, b], priorC.B[2, b])
    @test support.lq[3, lj, b] == Kpax3.logmarglik(data[b, i], 1, priorC.A[3, b], priorC.B[3, b])
    @test support.lq[4, lj, b] == Kpax3.logmarglik(data[b, i], 1, priorC.A[4, b], priorC.B[4, b])

    z[1] = priorC.logω[state.k + 1][1] + support.lq[3, lj, b]
    z[2] = priorC.logω[state.k + 1][2] + support.lq[4, lj, b]

    @test isapprox(support.lpj[3, b], log(exp(z[1]) + exp(z[2])), atol=ε)
  end

  Kpax3.gibbscomputeprobi!(li, support)

  for b in 1:support.m
    Kpax3.gibbscomputeprobg!(b, 1, li, support)
    Kpax3.gibbscomputeprobg!(b, 2, li, support)
    Kpax3.gibbscomputeprobj!(b, li, lj, support)
  end

  st = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 1; 3], priorR, priorC, settings)
  su = Kpax3.MCMCSupport(st, priorC)
  t = Kpax3.clusterweight(state.v[1], priorR) + su.logmlik

  @test isapprox(support.t[1], t, atol=ε)

  st = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 3; 3], priorR, priorC, settings)
  su = Kpax3.MCMCSupport(st, priorC)
  t = Kpax3.clusterweight(state.v[3], priorR) + su.logmlik

  @test isapprox(support.t[2], t, atol=ε)

  st = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 2; 3], priorR, priorC, settings)
  su = Kpax3.MCMCSupport(st, priorC)
  t = Kpax3.clusterweight(state.v[2] - 1, priorR) + su.logmlik

  @test isapprox(support.t[li], t, atol=ε)

  st = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 4; 3], priorR, priorC, settings)
  su = Kpax3.MCMCSupport(st, priorC)
  t = Kpax3.clusterweight(1, state.k, priorR) + su.logmlik

  @test isapprox(support.t[lj], t, atol=ε)

  logc = Kpax3.gibbscomputenormconst(lj, support)

  @test isapprox(logc, log(sum(exp.(support.t[1:4]))), atol=ε)

  nothing
end

test_mcmc_gibbs_move()

function test_mcmc_gibbs_merge()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=3, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  # remove unit 6 from cluster 3 and put it into another cluster
  R = [1; 1; 1; 2; 2; 3]

  i = 6
  hi = 3
  li = 3
  lj = 3

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  Kpax3.initsupportgibbs!(hi, priorR, state, support)

  support.t[lj] = Kpax3.clusterweight(1, state.k - 1, priorR)

  g = 0
  l = 0
  y = 0.0
  v = 0
  z = zeros(Float64, 2)
  for b in 1:m
    support.lpi[1, b] = priorC.logγ[1]
    support.lpi[2, b] = priorC.logγ[2]
    support.lpi[3, b] = priorC.logγ[3]

    g = 1
    l = 1
    y = state.n1s[g, b] + float(data[b, i])
    v = state.v[g] + 1

    Kpax3.gibbsupdateclustergmerge!(y, v, b, l, g, support.k, priorC, support)

    @test support.lq[1, l, b] == Kpax3.logmarglik(y, v, priorC.A[1, b], priorC.B[1, b])
    @test support.lq[2, l, b] == Kpax3.logmarglik(y, v, priorC.A[2, b], priorC.B[2, b])
    @test support.lq[3, l, b] == Kpax3.logmarglik(y, v, priorC.A[3, b], priorC.B[3, b])
    @test support.lq[4, l, b] == Kpax3.logmarglik(y, v, priorC.A[4, b], priorC.B[4, b])

    g = 2
    l = 2
    y = state.n1s[g, b] + float(data[b, i])
    v = state.v[g] + 1

    Kpax3.gibbsupdateclustergmerge!(y, v, b, l, g, support.k, priorC, support)

    @test support.lq[1, l, b] == Kpax3.logmarglik(y, v, priorC.A[1, b], priorC.B[1, b])
    @test support.lq[2, l, b] == Kpax3.logmarglik(y, v, priorC.A[2, b], priorC.B[2, b])
    @test support.lq[3, l, b] == Kpax3.logmarglik(y, v, priorC.A[3, b], priorC.B[3, b])
    @test support.lq[4, l, b] == Kpax3.logmarglik(y, v, priorC.A[4, b], priorC.B[4, b])

    @test isapprox(support.lpi[1, b], priorC.logγ[1] + support.lp[1, 1, b] + support.lp[1, 2, b], atol=ε)
    @test isapprox(support.lpi[2, b], priorC.logγ[2] + support.lp[2, 1, b] + support.lp[2, 2, b], atol=ε)

    z[1] = log(exp(priorC.logω[support.k][1] + support.lp[3, 1, b]) + exp(priorC.logω[support.k][2] + support.lp[4, 1, b]))
    z[2] = log(exp(priorC.logω[support.k][1] + support.lp[3, 2, b]) + exp(priorC.logω[support.k][2] + support.lp[4, 2, b]))

    @test isapprox(support.lpi[3, b], priorC.logγ[3] + z[1] + z[2], atol=ε)

    z[1] = support.lq[1, 1, b] + support.lq[1, li, b] - support.lp[1, 1, b]
    z[2] = support.lq[1, 2, b] + support.lq[1, li, b] - support.lp[1, 2, b]
    @test isapprox(support.lr[1, 1, b], z[1], atol=ε)
    @test isapprox(support.lr[1, 2, b], z[2], atol=ε)

    z[1] = support.lq[2, 1, b] + support.lq[2, li, b] - support.lp[2, 1, b]
    z[2] = support.lq[2, 2, b] + support.lq[2, li, b] - support.lp[2, 2, b]
    @test isapprox(support.lr[2, 1, b], z[1], atol=ε)
    @test isapprox(support.lr[2, 2, b], z[2], atol=ε)

    z[1] = log(exp(priorC.logω[support.k][1] + support.lq[3, 1, b]) + exp(priorC.logω[support.k][2] + support.lq[4, 1, b])) + support.lpj[1, b] - log(exp(priorC.logω[support.k][1] + support.lp[3, 1, b]) + exp(priorC.logω[support.k][2] + support.lp[4, 1, b]))
    z[2] = log(exp(priorC.logω[support.k][1] + support.lq[3, 2, b]) + exp(priorC.logω[support.k][2] + support.lq[4, 2, b])) + support.lpj[1, b] - log(exp(priorC.logω[support.k][1] + support.lp[3, 2, b]) + exp(priorC.logω[support.k][2] + support.lp[4, 2, b]))

    @test isapprox(support.lr[3, 1, b], z[1], atol=ε)
    @test isapprox(support.lr[3, 2, b], z[2], atol=ε)
  end

  Kpax3.gibbscomputeprobi!(lj, support)

  for b in 1:support.m
    Kpax3.gibbscomputeprobg!(b, 1, support)
    Kpax3.gibbscomputeprobg!(b, 2, support)
  end

  st = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 2; 1], priorR, priorC, settings)
  su = Kpax3.MCMCSupport(st, priorC)
  t = Kpax3.clusterweight(state.v[1], priorR) + su.logmlik

  @test isapprox(support.t[1], t, atol=ε)

  st = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 2; 2], priorR, priorC, settings)
  su = Kpax3.MCMCSupport(st, priorC)
  t = Kpax3.clusterweight(state.v[2], priorR) + su.logmlik

  @test isapprox(support.t[2], t, atol=ε)

  st = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 2; 3], priorR, priorC, settings)
  su = Kpax3.MCMCSupport(st, priorC)
  t = Kpax3.clusterweight(1, state.k - 1, priorR) + su.logmlik

  @test isapprox(support.t[lj], t, atol=ε)

  logc = Kpax3.gibbscomputenormconst(lj, support)

  @test isapprox(logc, log(sum(exp.(support.t[1:3]))), atol=ε)

  nothing
end

test_mcmc_gibbs_merge()

function test_mcmc_gibbs_sample_cluster()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=4, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  # remove unit 5 from cluster 2 and put it into another cluster
  R = [1; 1; 1; 2; 2; 3]

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  support.t[1] = log(0.001) + log(0.4)
  support.t[2] = log(0.001) + log(0.3)
  support.t[3] = log(0.001) + log(0.2)
  support.t[4] = log(0.001) + log(0.1)

  logc = log(0.001)

  N = 1000000
  g = 0
  counter = zeros(Float64, 6)
  for t in 1:N
    g = Kpax3.gibbssamplecluster(logc, 4, support)
    counter[g] += 1
  end

  counter /= N

  @test maximum(abs.(counter - [0.4; 0.3; 0.2; 0.1; 0.0; 0.0])) <= 0.005

  nothing
end

test_mcmc_gibbs_sample_cluster()

function test_mcmc_gibbs_perform_move()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=4, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  # remove unit 5 from cluster 2 and put it into another cluster
  R = [1; 1; 1; 2; 2; 3]

  i = 5
  hi = 2
  li = 3
  lj = 4

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  support.ni = vec(float(data[:, 4]))
  support.vi = 1

  Kpax3.initsupportgibbs!(hi, priorR, state, support)

  support.t[li] = Kpax3.clusterweight(support.vi, priorR)
  support.t[lj] = Kpax3.clusterweight(1, state.k, priorR)

  g = 0
  l = 0
  y = 0.0
  v = 0
  z = zeros(Float64, 2)
  for b in 1:m
    support.lpi[1, b] = priorC.logγ[1]
    support.lpi[2, b] = priorC.logγ[2]
    support.lpi[3, b] = priorC.logγ[3]
    support.lpi[4, b] = priorC.logγ[3]

    g = 1
    l = 1
    y = state.n1s[g, b] + float(data[b, i])
    v = state.v[g] + 1

    Kpax3.gibbsupdateclustergmove!(y, v, b, l, g, state.k, priorC, support)

    g = 3
    l = 2
    y = state.n1s[g, b] + float(data[b, i])
    v = state.v[g] + 1

    Kpax3.gibbsupdateclustergmove!(y, v, b, l, g, state.k, priorC, support)

    Kpax3.gibbsupdateclusteri!(state.k, b, li, priorC, support)
    Kpax3.gibbsupdateclusterj!(data[b, i], state.k + 1, b, li, lj, priorC, support)
  end

  Kpax3.gibbscomputeprobi!(li, support)

  for b in 1:support.m
    Kpax3.gibbscomputeprobg!(b, 1, li, support)
    Kpax3.gibbscomputeprobg!(b, 2, li, support)
    Kpax3.gibbscomputeprobj!(b, li, lj, support)
  end

  # R = [1; 1; 1; 2; 3; 3]
  l = 2

  Kpax3.gibbsmove!(i, hi, support.cl[l], li, l, data, priorR, support, state)

  st = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 3; 3], priorR, priorC, settings)
  su = Kpax3.MCMCSupport(st, priorC)

  @test state.R == st.R
  @test state.k == st.k

  @test !state.emptycluster[1]
  @test state.cl[1] == st.cl[1]
  @test state.v[1] == st.v[1]
  @test state.n1s[1, :] == st.n1s[1, :]
  @test state.unit[1][1:state.v[1]] == [1; 2; 3]

  @test !state.emptycluster[2]
  @test state.cl[2] == st.cl[2]
  @test state.v[2] == st.v[2]
  @test state.n1s[2, :] == st.n1s[2, :]
  @test state.unit[2][1:state.v[2]] == [4]

  @test !state.emptycluster[3]
  @test state.cl[3] == st.cl[3]
  @test state.v[3] == st.v[3]
  @test state.n1s[3, :] == st.n1s[3, :]
  @test state.unit[3][1:state.v[3]] == [6; 5]

  @test isapprox(state.logpR, st.logpR, atol=ε)

  @test isapprox(support.logmlik, su.logmlik, atol=ε)

  @test support.lp == su.lp

  nothing
end

test_mcmc_gibbs_perform_move()

function test_mcmc_gibbs_perform_merge()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  settings = Kpax3.KSettings(ifile, ofile, maxclust=3, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  # remove unit 6 from cluster 3 and put it into another cluster
  R = [1; 1; 1; 2; 2; 3]

  i = 6
  hi = 3
  li = 3
  lj = 3

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  Kpax3.initsupportgibbs!(hi, priorR, state, support)

  support.t[lj] = Kpax3.clusterweight(1, state.k - 1, priorR)

  g = 0
  l = 0
  y = 0.0
  v = 0
  z = zeros(Float64, 2)
  for b in 1:m
    support.lpi[1, b] = priorC.logγ[1]
    support.lpi[2, b] = priorC.logγ[2]
    support.lpi[3, b] = priorC.logγ[3]

    g = 1
    l = 1
    y = state.n1s[g, b] + float(data[b, i])
    v = state.v[g] + 1

    Kpax3.gibbsupdateclustergmerge!(y, v, b, l, g, support.k, priorC, support)

    g = 2
    l = 2
    y = state.n1s[g, b] + float(data[b, i])
    v = state.v[g] + 1

    Kpax3.gibbsupdateclustergmerge!(y, v, b, l, g, support.k, priorC, support)
  end

  Kpax3.gibbscomputeprobi!(lj, support)

  for b in 1:support.m
    Kpax3.gibbscomputeprobg!(b, 1, support)
    Kpax3.gibbscomputeprobg!(b, 2, support)
  end

  # R = [1; 1; 1; 2; 2; 1]
  l = 1

  Kpax3.gibbsmerge!(i, hi, support.cl[l], l, data, priorR, support, state)

  st = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 2; 1], priorR, priorC, settings)
  su = Kpax3.MCMCSupport(st, priorC)

  @test state.R == st.R
  @test state.k == st.k

  @test !state.emptycluster[1]
  @test state.cl[1] == st.cl[1]
  @test state.v[1] == st.v[1]
  @test state.n1s[1, :] == st.n1s[1, :]
  @test state.unit[1][1:state.v[1]] == [1; 2; 3; 6]

  @test !state.emptycluster[2]
  @test state.cl[2] == st.cl[2]
  @test state.v[2] == st.v[2]
  @test state.n1s[2, :] == st.n1s[2, :]
  @test state.unit[2][1:state.v[2]] == [4; 5]

  @test state.emptycluster[3]

  @test isapprox(state.logpR, st.logpR, atol=ε)

  @test isapprox(support.logmlik, su.logmlik, atol=ε)

  @test support.lp[:, 1:2, :] == su.lp[:, 1:2, :]

  nothing
end

test_mcmc_gibbs_perform_merge()

function test_mcmc_gibbs_perform_split()
  ifile = "data/read_proper_aa.fasta"
  ofile = "../build/test"

  # state should be resized
  settings = Kpax3.KSettings(ifile, ofile, maxclust=3, maxunit=1)

  data = UInt8[0 0 0 0 0 1;
               1 1 1 1 1 0;
               0 0 1 0 1 1;
               1 1 0 1 0 0;
               1 1 0 0 0 0;
               0 0 0 1 1 0;
               1 1 1 0 0 0;
               0 0 0 1 1 1;
               0 0 1 0 0 0;
               1 0 0 1 0 1;
               0 1 0 0 1 0;
               0 0 0 0 0 1;
               1 1 1 0 0 0;
               0 0 0 1 1 0;
               1 1 0 0 1 1;
               0 0 1 1 0 0;
               1 1 0 1 0 0;
               0 0 1 0 1 1]

  (m, n) = size(data)

  priorR = Kpax3.EwensPitman(settings.α, settings.θ)
  priorC = Kpax3.AminoAcidPriorCol(data, settings.γ, settings.r)

  # remove unit 5 from cluster 2 and put it into another cluster
  R = [1; 1; 1; 2; 2; 3]

  i = 5
  hi = 2
  li = 3
  lj = 4

  state = Kpax3.AminoAcidState(data, R, priorR, priorC, settings)
  support = Kpax3.MCMCSupport(state, priorC)

  Kpax3.resizesupport!(support, state.k + 1)

  support.ni = vec(float(data[:, 4]))
  support.vi = 1

  Kpax3.initsupportgibbs!(hi, priorR, state, support)

  support.t[li] = Kpax3.clusterweight(support.vi, priorR)
  support.t[lj] = Kpax3.clusterweight(1, state.k, priorR)

  g = 0
  l = 0
  y = 0.0
  v = 0
  z = zeros(Float64, 2)
  for b in 1:m
    support.lpi[1, b] = priorC.logγ[1]
    support.lpi[2, b] = priorC.logγ[2]
    support.lpi[3, b] = priorC.logγ[3]
    support.lpi[4, b] = priorC.logγ[3]

    g = 1
    l = 1
    y = state.n1s[g, b] + float(data[b, i])
    v = state.v[g] + 1

    Kpax3.gibbsupdateclustergmove!(y, v, b, l, g, state.k, priorC, support)

    g = 3
    l = 2
    y = state.n1s[g, b] + float(data[b, i])
    v = state.v[g] + 1

    Kpax3.gibbsupdateclustergmove!(y, v, b, l, g, state.k, priorC, support)

    Kpax3.gibbsupdateclusteri!(state.k, b, li, priorC, support)
    Kpax3.gibbsupdateclusterj!(data[b, i], state.k + 1, b, li, lj, priorC,
                             support)
  end

  Kpax3.gibbscomputeprobi!(li, support)

  for b in 1:support.m
    Kpax3.gibbscomputeprobg!(b, 1, li, support)
    Kpax3.gibbscomputeprobg!(b, 2, li, support)
    Kpax3.gibbscomputeprobj!(b, li, lj, support)
  end

  # R = [1; 1; 1; 2; 4; 3]
  l = 4

  Kpax3.gibbssplit!(i, hi, li, lj, data, priorR, settings, support, state)

  st = Kpax3.AminoAcidState(data, [1; 1; 1; 2; 4; 3], priorR, priorC, settings)
  su = Kpax3.MCMCSupport(st, priorC)

  @test state.R == st.R
  @test state.k == st.k

  @test !state.emptycluster[1]
  @test state.cl[1] == st.cl[1]
  @test state.v[1] == st.v[1]
  @test state.n1s[1, :] == st.n1s[1, :]
  @test state.unit[1][1:state.v[1]] == [1; 2; 3]

  @test !state.emptycluster[2]
  @test state.cl[2] == st.cl[2]
  @test state.v[2] == st.v[2]
  @test state.n1s[2, :] == st.n1s[2, :]
  @test state.unit[2][1:state.v[2]] == [4]

  @test !state.emptycluster[3]
  @test state.cl[3] == st.cl[3]
  @test state.v[3] == st.v[3]
  @test state.n1s[3, :] == st.n1s[3, :]
  @test state.unit[3][1:state.v[3]] == [6]

  @test !state.emptycluster[4]
  @test state.cl[4] == st.cl[4]
  @test state.v[4] == st.v[4]
  @test state.n1s[4, :] == st.n1s[4, :]
  @test state.unit[4][1:state.v[4]] == [i]

  @test isapprox(state.logpR, st.logpR, atol=ε)

  @test isapprox(support.logmlik, su.logmlik, atol=ε)

  @test support.lp[:, 1:4, :] == su.lp

  nothing
end

test_mcmc_gibbs_perform_split()