text:
function gillespie_single(mutfunc, # virulence mutation function
                   T::Real; # max time
                   # variable parameters
                   S::Integer, # susceptible host population size
                   I::AbstractVector, # infected host population sizes
                   R::AbstractVector, # reoverd host population sizes
                   v::AbstractVector, # virulence
                   # fixed parameters
                   β::Real, # infect coefficient
                   b::Real, # birth rate
                   d::Real, # death rate
                   r::Real, # recovery rate
                   μ::Real, # immunity loss rate
                   m::Real, # matation rate
                   maxepoch::Integer = 10_000_000)
    length(R) == length(I) == length(v) || error("R, I and v must have same length!")
    t = 0.0 # initial time
    epoch = 0x0000000000000000 # init epoch

    dynamics_S = Dynamics(t, S)

    @apply copy I R v

    sum_I = sum(I)
    sum_R = sum(R)
    p_I = I ./ sum_I
    E_v = transpose(p_I) * v
    dynamics_S = Dynamics(t, S)
    dynamics_I = Dynamics(t, sum_I)
    dynamics_R = Dynamics(t, sum_R)
    dynamics_v = Dynamics(t, E_v)

    maxepoch = UInt(maxepoch)
    while t <= T && epoch <= maxepoch && length(I) >= 0
        epoch += 0x0000000000000001

        # calculate reactions

        reproduction_S = b * S
        reproduction_I = b * I
        reproduction_R = b * R

        death_S = d * S
        death_I = d * I
        death_R = d * R

        infect_S = β * S * I # I_i infect S
        infect_R = β * I * transpose(R) # I_i infect R_j
        for i in 1:size(infect_R, 1)
            infect_R[i, i] = 0
        end
        recovery = r * I
        immuneloss = μ * R

        kill = v .* I

        pathogenmutation = m * I

        τ, r_idx, idx = choosereaction([reproduction_S, death_S], reproduction_I, reproduction_R,
                                       death_I, death_R, infect_S,
                                       infect_R, recovery, immuneloss,
                                       kill, pathogenmutation)

        t += τ # time goes

        if r_idx == 1 # S reproduction and death
            if idx::Int == 1
                S += 1
            elseif idx == 2
                S -= 1
            else
                error("Unexpected idx at reaction $r_idx.")
            end
        elseif r_idx == 2 || r_idx == 3 # I and R reproduction
            S += 1
        elseif r_idx == 4 # I death
            n = I[idx] -= 1
            if n <= 0 && R[idx] <= 0
                deleteat!(I, idx)
                deleteat!(R, idx)
                deleteat!(v, idx)
            end
        elseif r_idx == 5 # R death
            n = R[idx] -= 1
            if n <= 0 && I[idx] <= 0
                deleteat!(I, idx)
                deleteat!(R, idx)
                deleteat!(v, idx)
            end
        elseif r_idx == 6 # I_i infect S
            S -= 1
            n = I[idx] += 1
        elseif r_idx == 7 # I_i infect R_j
            i, j = idx
            n = I[i] += 1
            n = R[j] -= 1
            if n <= 0 && I[j] <= 0
                deleteat!(I, j)
                deleteat!(R, j)
                deleteat!(v, j)
            end
        elseif r_idx == 8 # recovery
            n = I[idx] -= 1
            n = R[idx] += 1
        elseif r_idx == 9 # immunity loss
            S += 1

            n = R[idx] -= 1
            if n <= 0 && I[idx] <= 0
                deleteat!(I, idx)
                deleteat!(R, idx)
                deleteat!(v, idx)
            end
        elseif r_idx == 10 # pathogen kill host
            n = I[idx] -= 1
            if n <= 0 && R[idx] <= 0
                deleteat!(I, idx)
                deleteat!(R, idx)
                deleteat!(v, idx)
            end
        elseif r_idx == 11 # pathogen mutation
            n = I[idx] -= 1
            new_dynamics_I = Dynamics(t, 1)
            new_dynamics_R = Dynamics(t, 0)
            push!(I, 1)
            push!(R, 0)
            new_v = mutfunc(v[idx])
            push!(v, new_v)
            if n <= 0 && R[idx] <= 0
                deleteat!(I, idx)
                deleteat!(R, idx)
                deleteat!(v, idx)
            end
        end
        sum_I = sum(I)
        sum_R = sum(R)
        p_I = I ./ sum_I
        E_v = transpose(p_I) * v
        record!(dynamics_S, t, S)
        record!(dynamics_I, t, sum_I)
        record!(dynamics_R, t, sum_R)
        record!(dynamics_v, t, E_v)
    end
    println("Simulation ends at epoch ", Int(epoch), "!")
    return dynamics_S, dynamics_I, dynamics_R, dynamics_v
end


text:
using FinancialMonteCarlo
path1 = joinpath(dirname(pathof(FinancialMonteCarlo)), "..", "test")
test_listTmp = readdir(path1);

BlackList = ["REQUIRE", "runtests.jl", "Project.toml", "Manifest.toml", "cuda", "af", "wip", "test_black_rev_diff.jl", "test_black_pre.jl"];
func_scope(x::String) = include(x);
return nothing;
test_list = [test_element for test_element in test_listTmp if !any(x -> x == test_element, BlackList)]
println("Running tests:\n")
for (current_test, i) in zip(test_list, 1:length(test_list))
    println("------------------------------------------------------------")
    println("  * $(current_test) *")
    func_scope(joinpath(path1, current_test))
    println("------------------------------------------------------------")
    if (i < length(test_list))
        println("")
    end
end

ext: jl


text:
using ChainRulesCore
using ChainRulesTestUtils
using QuadGK
using StableRNGs
using Transits
using Test
using Unitful, UnitfulAstro

Unitful.preferunits(u"Rsun,Msun,d"...)
ENV["UNITFUL_FANCY_EXPONENTS"] = false

# Numpy version of `isapprox`
# https://stackoverflow.com/questions/27098844/allclose-how-to-check-if-two-arrays-are-close-in-julia/27100515#27100515
allclose(a, b; rtol=1e-5, atol=1e-8) = mapreduce((a,b) -> abs(a - b) ≤ (atol + rtol*abs(b)), &, a, b)

const PLOT = get(ENV, "TEST_PLOTS", "false") == "true"
PLOT && include("plots.jl")

rng = StableRNG(2752)

@testset "Transits" begin
    include("Mn_integral.jl")
    include("distributions.jl")
    include("distributions.jl")
    include("elliptic.jl")
    include("grads.jl")
    include("orbits/keplerian.jl")
    include("orbits/simple.jl")
    include("orbits/solvers.jl")
    include("poly.jl")
    include("show.jl")
end


text:
export iauPb06
"""
This function forms three Euler angles which implement general
precession from epoch J2000.0, using the IAU 2006 model.  Frame
bias (the offset between ICRS and mean J2000.0) is included.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2  double   TT as a 2-part Julian Date (Note 1)

Returned:
   bzeta        double   1st rotation: radians cw around z
   bz           double   3rd rotation: radians cw around z
   btheta       double   2nd rotation: radians ccw around y

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The traditional accumulated precession angles zeta_A, z_A,
   theta_A cannot be obtained in the usual way, namely through
   polynomial expressions, because of the frame bias.  The latter
   means that two of the angles undergo rapid changes near this
   date.  They are instead the results of decomposing the
   precession-bias matrix obtained by using the Fukushima-Williams
   method, which does not suffer from the problem.  The
   decomposition returns values which can be used in the
   conventional formulation and which include frame bias.

   3. The three angles are returned in the conventional order, which
   is not the same as the order of the corresponding Euler
   rotations.  The precession-bias matrix is
   R_3(-z) x R_2(+theta) x R_3(-zeta).

   4. Should zeta_A, z_A, theta_A angles be required that do not
   contain frame bias, they are available by calling the SOFA
   function iauP06e.

Called:
   iauPmat06    PB matrix, IAU 2006
   iauRz        rotate around Z-axis

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPb06(double date1, double date2,
#              double *bzeta, double *bz, double *btheta)

function iauPb06(date1::Real, date2::Real)
   ref_bzeta  = Ref{Float64}(0.0)
   ref_bz     = Ref{Float64}(0.0)
   ref_btheta = Ref{Float64}(0.0)

   ccall((:iauPb06, libsofa_c), Cvoid,
       (Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
       convert(Float64, date1),
       convert(Float64, date2),
       ref_bzeta, ref_bz, ref_btheta)

   return ref_bzeta[], ref_bz[], ref_btheta[]
end


text:
# QCQP

#using Everest

include("src/loadcode.jl");
test_problem_nlo("BRANDY")

##############################################################
#############################################################
testmat = rand(2,2);
testmat[1,1] = 1e8;
testmat[1,2] = 1e8;
testmat[2,1] = 1e8 + 1;
testmat[2,2] = 1e8
testrhs = rand(2)
sol1 = testmat \ testrhs
norm(testmat*sol1-testrhs);

sol2 = lufact(testmat) \ testrhs
norm(testmat*sol2-testrhs);
norm(testmat*inv(testmat)*testrhs-testrhs);

function test(R::SparseMatrixCSC{Float64,Int64},v::Array{Float64,1})
  println("cholfac")
  @time F = cholfact(R + 1e-7 * speye(size(R,1)));

  println("solve")

  @time begin
    for i = 1:100
        sol = F \ v
    end
  end
end

function test2(A::SparseMatrixCSC{Float64,Int64},v::Array{Float64,1})
  println("lufac")
  @time F = lufact([speye(size(A,2))]);

  println("solve")

  @time begin
    for i = 1:50
        F \ v # compare with mumps
    end
  end
end

function test3(R::SparseMatrixCSC{Float64,Int64},V::AbstractMatrix)
  println("cholfac")
  @time F = cholfact(R + 1e-7 * speye(size(R,1)));

  println("solve")

  @time sol = F \ V
end

A,b,c = get_netlib_problem("QAP15")
@time R = A*A';

v = ones(size(R,1));
test(R,v)

siz = 100;
rhs = zeros(size(R,1),siz);
rhs[1,:] = ones(siz)
test3(R,rhs)

using MUMPS

max_forks_repo_name: ohinder/Everest


text:
module Smile

##############################################################################
##
## Dependencies
##
##############################################################################

##############################################################################
##
## Ensure libsmilejl can be found
##
##############################################################################

pathtoadd = joinpath(dirname(@__FILE__), "..", "deps", "downloads")

if isempty(Libdl.find_library(collect(["libsmilejl"]), collect([pathtoadd])))
	warn("Could not communicate with libsmilejl. Make sure it exists and that it is on the path.")
else

	global const LIB_SMILE = Libdl.find_library(collect(["libsmilejl"]), collect([pathtoadd]))

	import Base: getindex, setindex!, print

	include("low_level/types.jl")
	include("low_level/network.jl")
	include("low_level/pattern.jl")
	include("low_level/node.jl")
	include("low_level/nodedefinition.jl")
	include("low_level/intarray.jl")
	include("low_level/idarray.jl")
	include("low_level/doublearray.jl")
	include("low_level/syscoordinates.jl")
	include("low_level/dataset.jl")
	include("low_level/nodevalue.jl")
	include("low_level/dmatrix.jl")
	include("low_level/learning.jl")

	include("high_level/learning.jl")
	include("high_level/sampling.jl")
end

end # end module


text:
"""
$(TYPEDEF)
$(TYPEDFIELDS)

One or more Petri Nets and an ID Registry shared by all nets.
"""
struct PnmlModel
    nets::Vector{PnmlNet}
    namespace::String
    reg::IDRegistry # Shared by all nets.
    xml::XMLNode
end

"""
$(TYPEDSIGNATURES)
"""

PnmlModel(net::PnmlNet) = PnmlModel([net])
PnmlModel(nets::Vector{PnmlNet}) = PnmlModel(nets, pnml_ns, IDRegistry(), nothing)
PnmlModel(nets::Vector{PnmlNet}, ns, reg::IDRegistry) = PnmlModel(nets, ns, reg, nothing)

"""
$(TYPEDSIGNATURES)

Return all `nets` of `model`.
"""
nets(model::PnmlModel) = model.nets
namespace(model::PnmlModel) = model.namespace
idregistry(model::PnmlModel) = model.reg
xmlnode(model::PnmlModel) = model.xml

"""
$(TYPEDSIGNATURES)

Build a PnmlModel from a string containing XML.
"""
function parse_str(str::AbstractString)
    reg = IDRegistry()
    # Good place for debugging.
    parse_pnml(root(EzXML.parsexml(str)); reg)
end

"""
$(TYPEDSIGNATURES)

Build a PnmlModel from a file containing XML.
"""
function parse_file(fname::AbstractString)
    reg = IDRegistry()
    parse_pnml(root(EzXML.readxml(fname)); reg)
end

"""
$(TYPEDSIGNATURES)
Return nets matching pntd `type` given as string or symbol.
See [`PnmlTypes.pntd_symbol`](@ref), [`PnmlTypes.pnmltype`](@ref).
"""
function find_nets end
find_nets(model, type::AbstractString) = find_nets(model, pntd_symbol(type))
find_nets(model, type::Symbol) = find_nets(model, pnmltype(type))
find_nets(model, ::PNTD) where {PNTD <: PnmlType} =
    filter(n->typeof(n.type) <: PNTD, nets(model))

"""
$(TYPEDSIGNATURES)

Return `PnmlNet` having `id` or `nothing``.
"""
function find_net end

function find_net(model, id::Symbol)
    getfirst(net->pid(net) === id, nets(model))
end

"""
$(TYPEDSIGNATURES)

Return first net contained by `doc`.
"""
first_net(model) = first(nets(model))


text:
using Test
using TightlyBound

function test1()

    @testset "testing sym1" begin
        abc = rand(3)
        lmn = rand(3)
        lmn = lmn / sqrt(sum(lmn.^2))

        lmnx = [1, 0, 0]
        lmny = [0, 1, 0]

        tot1 = 0.0
        tot2 = 0.0
        for orb1 = [:s, :pz, :px, :py, :dz2,:dxz,:dyz,:dx2_y2,:dxy]
            for orb2 = [:s, :pz, :px, :py, :dz2,:dxz,:dyz,:dx2_y2,:dxy]

                t = TightlyBound.CalcTB.symmetry_factor(orb1, orb2, lmn, abc)
                t2 = TightlyBound.CalcTB.symmetry_factor_fit(orb1, orb2, lmn)

#                println(orb1," ", orb2, "   ", t, "   " , sum(t2.*abc) )
#                if isapprox(t, sum(t2.*abc), rtol=1e-9) == false
#                    println("ERROR $orb1 $orb2")
#                end
                @test isapprox(t, sum(t2.*abc), rtol=1e-9)

            end
        end
        for orb1 = [:s, :pz, :px, :py, :dxy, :dz2, :dxz, :dyz]
            for orb2 = [:s, :pz, :px, :py, :dxy, :dz2, :dxz, :dyz]
                tot1 += TightlyBound.CalcTB.symmetry_factor(orb1, orb2, lmnx, abc)
                tot2 += TightlyBound.CalcTB.symmetry_factor(orb1, orb2, lmny, abc)
            end
        end

#        println("tot1 $tot1 tot2 $tot2 ")
        @test isapprox(tot1, tot2, rtol=1e-5)

    end
end

test1()
#test2()


text:
function fourier_quadrature(w::Matrix{CT}, v::Union{Array{CT,2},SharedArray{CT,2}}, u::Union{Array{CT,2},SharedArray{CT,2}}, grid::Grid2d{M,N,RT}) where {M,N,RT,CT}
    @assert size(u) == size(v) == size(w)

    local result::CT = 0

    @inbounds for j in 1:size(u,2)
        for i in 1:size(u,1)
            result += w[i,j] * conj(v[i,j]) * u[i,j]
        end
    end

    result * grid.normalisation
end

function fourier_quadrature(w::Matrix{Matrix{CT}}, v::Union{Array{CT,2},SharedArray{CT,2}}, u::Union{Array{CT,2},SharedArray{CT,2}}, grid::Grid2d{M,N,RT}) where {M,N,RT,CT}
    local result::Matrix{CT} = zeros(CT, size(w,1), size(w,2))

    for l in 1:size(w,2)
        for k in 1:size(w,1)
            result[k,l] = fourier_quarature(w[k,l], v, u, grid)
        end
    end

    result
end

function fourier_quadrature(w::Matrix{Matrix{CT}}, v::Union{Vector{Array{CT,2}},Vector{SharedArray{CT,2}}}, u::Union{Array{CT,2},SharedArray{CT,2}}, grid::Grid2d{M,N,RT}) where {M,N,RT,CT}
    @assert size(w,2) == length(v)

    local result::Vector{CT} = zeros(CT, size(w,1))

    for k in 1:size(w,1)
        for l in 1:size(w,2)
            result[k] += fourier_quarature(w[k,l], v[l], u, grid)
        end
    end

    result
end

# function fourier_quadrature(w::Matrix{CT}, v::Union{Array{CT,3},SharedArray{CT,3}}, l, u::Union{Array{CT,2},SharedArray{CT,2}}, grid::Grid2d{RT}) where {RT,CT}
#     @assert size(u) == (size(v,1), size(v,2)) == size(w)
#     @assert l ≤ size(v,3)
#
#     local result::CT = 0
#
#     @inbounds for j in 1:size(u,2)
#         for i in 1:size(u,1)
#             result += w[i,j] * conj(v[i,j,l]) * u[i,j]
#         end
#     end
#
#     result * grid.normalisation
# end


text:
################################################################################
#
#       GrpAb/GrpAbFinGen.jl : Finitely generated abelian groups
#
# This file is part of Hecke.
#
# Copyright (c) 2015, 2016, 2017: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016, 2017 Tommy Hofmann, Claus Fieker
#
################################################################################

export AbelianGroup, DiagonalGroup, issnf, ngens, nrels, rels, snf, isfinite,
       isinfinite, rank, order, exponent, istrivial, isisomorphic,
       direct_product, istorsion, torsion_subgroup, sub, quo, iscyclic,
       psylow_subgroup

import Base.+, Nemo.snf, Nemo.parent, Base.rand, Nemo.issnf

################################################################################
#
#  Functions for some abstract interfaces
#
################################################################################

elem_type(G::GrpAbFinGen) = GrpAbFinGenElem

elem_type(::Type{GrpAbFinGen}) = GrpAbFinGenElem

parent_type(::Type{GrpAbFinGenElem}) = GrpAbFinGen

##############################################################################
#
#  Constructors
#
##############################################################################

# We do we have AbelianGroup and DiagonalGroup?
doc"""
***
    AbelianGroup(M::fmpz_mat) -> GrpAbFinGen

> Creates the abelian group with relation matrix `M`. That is, the group will
> have `cols(M)` generators and each row of `M` describes one relation.
"""
function AbelianGroup(M::fmpz_mat)
  return GrpAbFinGen(M)
end

doc"""
***
    AbelianGroup(M::Array{fmpz, 2}) -> GrpAbFinGen

> Creates the abelian group with relation matrix `M`. That is, the group will
> have `cols(M)` generators and each row of `M` describes one relation.
"""
function AbelianGroup(M::Array{fmpz, 2})
  return AbelianGroup(matrix(FlintZZ, M))
end

doc"""
***
    AbelianGroup(M::Array{Integer, 2}) -> GrpAbFinGen

> Creates the abelian group with relation matrix `M`. That is, the group will
> have `cols(M)` generators and each row of `M` describes one relation.
"""
function AbelianGroup(M::Array{T, 2}) where T <: Integer
  return AbelianGroup(matrix(FlintZZ, M))
end

doc"""
***
    AbelianGroup(M::Array{fmpz, 1}) -> GrpAbFinGen

> Creates the abelian group with relation matrix `M`. That is, the group will
> have `length(M)` generators and one relation.
"""
function AbelianGroup(M::Array{fmpz, 1})
  return AbelianGroup(matrix(FlintZZ, 1, length(M), M))
end

doc"""
***
    AbelianGroup(M::Array{Integer, 1}) -> GrpAbFinGen

> Creates the abelian group with relation matrix `M`. That is, the group will
> have `length(M)` generators and one relation.
"""
function AbelianGroup(M::Array{T, 1}) where T <: Integer
  return AbelianGroup(matrix(FlintZZ, 1, length(M), M))
end

doc"""
***
    DiagonalGroup(M::fmpz_mat) -> GrpAbFinGen

Assuming that $M$ has only one row, this function creates the direct product of
the cyclic groups $\mathbf{Z}/m_i$, where $m_i$ is the $i$th entry of `M`.
"""
function DiagonalGroup(M::fmpz_mat)
  if rows(M) != 1
    error("The argument must have only one row")
  end

  N = zero_matrix(FlintZZ, cols(M), cols(M))
  for i = 1:cols(M)
    N[i,i] = M[1, i]
  end
  if issnf(N)
    return GrpAbFinGen(fmpz[M[1, i] for i = 1:cols(M)])
  else
    return GrpAbFinGen(N)
  end
end

doc"""
***
    DiagonalGroup(M::Array{Union{fmpz, Integer}, 1}) -> GrpAbFinGen

Creates the direct product of the cyclic groups $\mathbf{Z}/m_i$,
where $m_i$ is the $i$th entry of `M`.
"""
function DiagonalGroup(M::Array{T, 1}) where T <: Union{Integer, fmpz}
  N = zero_matrix(FlintZZ, length(M), length(M))
  for i = 1:length(M)
    N[i,i] = M[i]
  end
  if issnf(N)
    return GrpAbFinGen(M)
  else
    return GrpAbFinGen(N)
  end
end

################################################################################
#
#  String I/O
#
################################################################################

function show(io::IO, A::GrpAbFinGen)
  if issnf(A)
    show_snf(io, A)
  else
    show_gen(io, A)
  end
end

function show_gen(io::IO, A::GrpAbFinGen)
  print(io, "(General) abelian group with relation matrix\n$(A.rels)")
  if isdefined(A, :snf_map)
    println(io, "\nwith structure of ", domain(A.snf_map))
  end
end

function show_snf(io::IO, A::GrpAbFinGen)
  len = length(A.snf)

  println(io, "Abelian group")

  if len == 0
    print(io, "Z/1")
    return
  end

  if A.snf[1] == 0
    if len == 1
      print(io, "Z")
    else
      print(io, "Z^$(len)")
    end
    return
  end

  i = 1
  while i <= len
    inv = A.snf[i]
    j = 1
    while i + j <= len && inv == A.snf[i + j]
      j += 1
    end
    if iszero(inv)
      print(io, "Z")
    else
      print(io, "Z/$(inv)")
    end
    if j > 1
      print(io, "^$(j)")
    end
    if i + j - 1 < len
      print(io, " x ")
    end
    i += j
  end
end

################################################################################
#
#  Field access
#
################################################################################

doc"""
***
    issnf(G::GrpAbFinGen) -> Bool

> Returns whether the current relation matrix of the group $G$ is in Smith
> normal form.
"""
issnf(A::GrpAbFinGen) = A.issnf

doc"""
***
    ngens(G::GrpAbFinGen) -> Int

> Returns the number of generators of $G$ in the current representation.
"""
function ngens(A::GrpAbFinGen)
  if issnf(A)
    return length(A.snf)
  else
    return cols(A.rels)
  end
end

doc"""
***
    nrels(G::GrpAbFinGen) -> Int

> Returns the number of relations of $G$ in the current representation.
"""
function nrels(A::GrpAbFinGen)
  if issnf(A)
    return length(A.snf)
  else
    return rows(A.rels)
  end
end

doc"""
***
    rels(A::GrpAbFinGen) -> fmpz_mat

> Returns the currently used relations of $G$ as a single matrix.
"""
function rels(A::GrpAbFinGen)
  if issnf(A)
    return rels_snf(A)
  else
    return rels_gen(A)
  end
end

function rels_gen(A::GrpAbFinGen)
  return A.rels
end

function rels_snf(A::GrpAbFinGen)
  M = zero_matrix(FlintZZ, ngens(A), ngens(A))
  for i = 1:ngens(A)
    M[i,i] = A.snf[i]
  end
  return M
end

################################################################################
#
#  Hermite normal form
#
################################################################################

function assure_has_hnf(A::GrpAbFinGen)
  isdefined(A, :hnf) && return
  A.hnf = hnf(A.rels)
end

################################################################################
#
#  Smith normal form
#
################################################################################

doc"""
***
    snf(A::GrpAbFinGen) -> GrpAbFinGen, Map

> Returns a pair $(G, f)$, where $G$ is an abelian group in canonical Smith
> normal form isomorphic to $G$ and an isomorphism $f : G \to A$.
"""
function snf(G::GrpAbFinGen)
  if isdefined(G, :snf_map)
    return domain(G.snf_map)::GrpAbFinGen, G.snf_map::GrpAbFinGenMap
  end

  if issnf(G)
    G.snf_map = GrpAbFinGenMap(G) # identity
    return G, G.snf_map::GrpAbFinGenMap
  end

  S, _, T = snf_with_transform(G.rels, false, true)
  d = fmpz[S[i,i] for i = 1:min(rows(S), cols(S))]

  while length(d) < ngens(G)
    push!(d, 0)
  end

  #s = Array{fmpz, 1}()
  s = fmpz[ d[i] for i in 1:length(d) if d[i] !=  1]
  #for i = 1:length(d)
  #  if d[i] != 1
  #    push!(s, d[i])
  #  end
  #end
  TT = zero_matrix(FlintZZ, rows(T), length(s))
  j = 1
  for i = 1:length(d)
    if d[i] != 1
      for k=1:rows(T)
        TT[k, j] = T[k, i]
      end
      j += 1
    end
  end

  TTi = zero_matrix(FlintZZ, length(s), rows(T))
  Ti = inv(T)

  j = 1
  for i = 1:length(d)
    if d[i] != 1
      for k=1:rows(T)
        TTi[j, k] = Ti[i, k]
      end
      j += 1
    end
  end

  H = GrpAbFinGen(s)
  mp = hom(H, G, TTi, TT, false)
  G.snf_map = mp
  return H, mp::GrpAbFinGenMap
end

################################################################################
#
#  Finiteness
#
################################################################################

doc"""
***
    isfinite(A::GrpAbFinGen) -> Bool

> Returns whether $A$ is finite.
"""
isfinite(A::GrpAbFinGen) = issnf(A) ? isfinite_snf(A) : isfinite_gen(A)

isfinite_snf(A::GrpAbFinGen) = length(A.snf) == 0 || !iszero(A.snf[end])

isfinite_gen(A::GrpAbFinGen) = isfinite(snf(A)[1])

doc"""
***
    isinfinite(A::GrpAbFinGen) -> Bool

> Returns whether $A$ is infinite.
"""
isinfinite(A::GrpAbFinGen) = !isfinite(A)

################################################################################
#
#  Rank
#
################################################################################

doc"""
***
    rank(A::GrpAbFinGenGen) -> Int

> Returns the rank of $A$, that is, the dimension of the
> $\mathbf{Q}$-vectorspace $A \otimes_{\mathbf Z} \mathbf Q$.
"""
rank(A::GrpAbFinGen) = issnf(A) ? rank_snf(A) : rank_gen(A)

rank_snf(A::GrpAbFinGen) = length(find(x -> x == 0, A.snf))

rank_gen(A::GrpAbFinGen) = rank(snf(A)[1])

################################################################################
#
#  Order
#
################################################################################

doc"""
***
    order(A::GrpAbFinGen) -> fmpz

> Returns the order of $A$. It is assumed that $A$ is finite.
"""
order(A::GrpAbFinGen) = issnf(A) ? order_snf(A) : order_gen(A)

function order_snf(A::GrpAbFinGen)
  isinfinite(A) && error("Group must be finite")
  return prod(A.snf)
end

order_gen(A::GrpAbFinGen) = order(snf(A)[1])

################################################################################
#
#  Exponent
#
################################################################################

doc"""
***
    exponent(A::GrpAbFinGen) -> fmpz

> Returns the exponent of $A$. It is assumed that $A$ is finite.
"""
exponent(A::GrpAbFinGen) = issnf(A) ? exponent_snf(A) : exponent_gen(A)

function exponent_snf(A::GrpAbFinGen)
  isinfinite(A) && error("Group must be finite")
  ngens(A)==0 && return 1
  return A.snf[end]
end

exponent_gen(A::GrpAbFinGen) = exponent(snf(A)[1])

################################################################################
#
#  Trivial
#
################################################################################

doc"""
***
    istrivial(A::GrpAbFinGen) -> Bool

> Checks if $A$ is the trivial group.
"""
istrivial(A::GrpAbFinGen) = isfinite(A) && order(A) == 1

################################################################################
#
#  Isomorphism
#
################################################################################

doc"""
***
    isisomorphic(G::GrpAbFinGen, H::GrpAbFinGen) -> Bool

> Checks if $G$ and $H$ are isomorphic.
"""
function isisomorphic(G::GrpAbFinGen, H::GrpAbFinGen)
  b = filter(x -> x != 1, snf(G)[1].snf) == filter(x -> x != 1, snf(H)[1].snf)
  return b
end

################################################################################
#
#  Direct product
#
################################################################################

doc"""
***
    direct_product(G::GrpAbFinGen, H::GrpAbFinGen) -> GrpAbFinGen

> Returns the abelian group $G\times H$.
"""
function direct_product(G::GrpAbFinGen, H::GrpAbFinGen)
  A = vcat(rels(G), zero_matrix(FlintZZ, rows(rels(H)), cols(rels(G))))
  B = vcat(zero_matrix(FlintZZ, rows(rels(G)), cols(rels(H))), rels(H))

  return AbelianGroup(hcat(A,B))
end

################################################################################
#
#  Torsion
#
################################################################################

istorsion(G::GrpAbFinGen) = isfinite(G)

function torsion_subgroup(G::GrpAbFinGen)
  S, mS = snf(G)
  subs = GrpAbFinGenElem[]
  for i in 1:ngens(S)
    if !iszero(S.snf[i])
      push!(subs, mS(S[i]))
    end
  end
  return sub(G, subs)
end

##############################################################################
#
#  Subgroups
#
##############################################################################

doc"""
***
    sub(G::GrpAbFinGen, s::Array{GrpAbFinGenElem, 1}) -> GrpAbFinGen, Map

> Create the subgroup $H$ of $G$ generated by the elements in `s` together
> with the injection $\iota : H \to G$.
"""
function sub(G::GrpAbFinGen, s::Array{GrpAbFinGenElem, 1},
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)

  if length(s) == 0
    S = GrpAbFinGen(fmpz[1])
    I = zero_matrix(FlintZZ, ngens(S), ngens(G))
    m = hom(S, G, I, false)
    if add_to_lattice
      append!(L, m)
    end
    return S, m
  end

  p = s[1].parent
  @assert G == p
  m = zero_matrix(FlintZZ, length(s) + nrels(p), ngens(p) + length(s))
  for i = 1:length(s)
    for j = 1:ngens(p)
      m[i + nrels(p), j] = s[i][j]
    end
    m[i + nrels(p), i + ngens(p)] = 1
  end
  if issnf(p)
    for i = 1:ngens(p)
      m[i, i] = p.snf[i]
    end
  else
    for i = 1:nrels(p)
      for j = 1:ngens(p)
        m[i, j] = p.rels[i, j]
      end
    end
  end
  h = hnf(m)
  fstWithoutOldGens = 1
  for i in rows(h):-1:1, j in ngens(p):-1:1
    if !iszero(h[i,j])
      fstWithoutOldGens = i + 1
      break
    end
  end
  r = sub(h, fstWithoutOldGens:rows(h), ngens(p) + 1:cols(h))
  S = AbelianGroup(r)

  mS = hom(S, p, sub(m, (nrels(p) + 1):rows(h), 1:ngens(p)), false)

  if add_to_lattice
    append!(L, mS)
  end
  return S, mS
end

doc"""
***
    sub(s::Array{GrpAbFinGenElem, 1}) -> GrpAbFinGen, Map

> Assuming that the non-empty array `s` contains elements of an abelian group
$G$, this functions returns the subgroup $H$ of $G$ generated by the elements
> in `s` together with the injection $\iota : H \to G$.
"""
function sub(s::Array{GrpAbFinGenElem, 1},
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  length(s) == 0 && error("Array must be non-empty")
  return sub(parent(s[1]), s, add_to_lattice, L)
end

doc"""
***
    sub(G::GrpAbFinGen, M::fmpz_mat) -> GrpAbFinGen, Map

> Create the subgroup $H$ of $G$ generated by the elements corresponding to the
> rows of $M$ together with the injection $\iota : H \to G$.
"""
function sub(G::GrpAbFinGen, M::fmpz_mat,
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  m = zero_matrix(FlintZZ, rows(M) + nrels(G), ngens(G) + rows(M))
  for i = 1:rows(M)
    for j = 1:ngens(G)
      m[i + nrels(G), j] = M[i,j]
    end
    m[i + nrels(G), i + ngens(G)] = 1
  end
  if issnf(G)
    for i = 1:ngens(G)
      m[i, i] = G.snf[i]
    end
  else
    for i = 1:nrels(G)
      for j = 1:ngens(G)
        m[i, j] = G.rels[i, j]
      end
    end
  end

  h = hnf(m)
  fstWithoutOldGens = 1

  for i in rows(h):-1:1, j in ngens(G):-1:1
    if !iszero(h[i,j])
      fstWithoutOldGens = i + 1
      break
    end
  end
  r = sub(h, fstWithoutOldGens:rows(h), ngens(G) + 1:cols(h))
  S = AbelianGroup(r)
  mS = hom(S, G, sub(m, (nrels(G) + 1):rows(h), 1:ngens(G)), false)

  if add_to_lattice
    append!(L, mS)
  end
  return S, mS
end

doc"""
***
    sub(G::GrpAbFinGen, n::fmpz) -> GrpAbFinGen, Map

> Create the subgroup $n \cdot G$ of $G$ together
> with the injection $\iota : n\cdot G \to G$.
"""
function sub(G::GrpAbFinGen, n::fmpz,
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  sg = [ n*g for g in gens(G) if !iszero(n*g)]
  return sub(G, sg, add_to_lattice, L)
end

doc"""
***
    sub(G::GrpAbFinGen, n::Integer) -> GrpAbFinGen, Map

> Create the subgroup $n \cdot G$ of $G$ together
> with the injection $\iota : n \cdot G \to G$.
"""
function sub(G::GrpAbFinGen, n::Integer,
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  return sub(G, fmpz(n), add_to_lattice, L)
end

################################################################################
#
#  Quotients
#
################################################################################

doc"""
***
  quo(G::GrpAbFinGen, s::Array{GrpAbFinGenElem, 1}) -> GrpAbFinGen, Map

> Create the quotient $H$ of $G$ by the subgroup generated by the elements in
> $s$, together with the projection $p : G \to H$.
"""
function quo(G::GrpAbFinGen, s::Array{GrpAbFinGenElem, 1},
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  if length(s) == 0
    I = identity_matrix(FlintZZ, ngens(G))
    m = hom(G, G, I, I, false)
    if add_to_lattice
      append!(L, m)
    end
    return G, m
  end

  p = s[1].parent
  @assert G == p
  m = zero_matrix(FlintZZ, length(s)+nrels(p), ngens(p))
  for i = 1:length(s)
    for j = 1:ngens(p)
      m[i + nrels(p),j] = s[i][j]
    end
  end
  if issnf(p)
    for i = 1:ngens(p)
      m[i, i] = p.snf[i]
    end
  else
    for i = 1:nrels(p)
      for j = 1:ngens(p)
        m[i, j] = p.rels[i, j]
      end
    end
  end

  Q = AbelianGroup(m)
  I = identity_matrix(FlintZZ, ngens(p))
  m = hom(p, Q, I, I, false)
  if add_to_lattice
    append!(L, m)
  end
  return Q, m
end

doc"""
***
  quo(G::GrpAbFinGen, M::fmpz_mat) -> GrpAbFinGen, Map

> Create the quotient $H$ of $G$ by the subgroup generated by the elements
> corresponding to the rows of $M$. together with the projection $p : G \to H$.
"""
function quo(G::GrpAbFinGen, M::fmpz_mat,
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  m = vcat(rels(G), M)
  Q = AbelianGroup(m)
  I = identity_matrix(FlintZZ, ngens(G))
  m = hom(G, Q, I, I, false)
  if add_to_lattice
    append!(L, m)
  end
  return Q, m
end

doc"""
***
    quo(G::GrpAbFinGen, n::Integer}) -> GrpAbFinGen, Map
    quo(G::GrpAbFinGen, n::fmpz}) -> GrpAbFinGen, Map

> Returns the quotient $H = G/nG$ together with the projection $p : G \to H$.
"""
function quo(G::GrpAbFinGen, n::Union{fmpz, Integer},
             add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  if issnf(G)
    return quo_snf(G, n, add_to_lattice, L)
  else
    quo_gen(G, n, add_to_lattice, L)
  end
end

function quo_snf(G::GrpAbFinGen, n::Union{fmpz, Integer},
                 add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  r = [gcd(x, n) for x = G.snf]
  I = identity_matrix(FlintZZ, ngens(G))
  Q = DiagonalGroup(r)
  m = hom(G, Q, I, I, false)
  if add_to_lattice
    append!(L, m)
  end
  return Q, m
end

function quo_gen(G::GrpAbFinGen, n::Union{fmpz, Integer},
                 add_to_lattice::Bool = true, L::GrpAbLattice = GroupLattice)
  m = vcat(G.rels, n*identity_matrix(FlintZZ, ngens(G)))
  Q = AbelianGroup(m)
  I = identity_matrix(FlintZZ, ngens(G))
  m = hom(G, Q, I, I, false)
  if add_to_lattice
    append!(L, m)
  end
  return Q, m
end

################################################################################
#
#  Make Smith normal form
#
################################################################################

# TODO: This is type unstable
function make_snf(m::Map{GrpAbFinGen, T}) where T
  G = domain(m)
  if issnf(G)
    return m
  end
  S, mS = snf(G)
  return _compose(mS, m)
end

################################################################################
#
#  Cyclic
#
################################################################################

doc"""
    iscyclic(G::GrpAbFinGen) -> Bool

> Returns whether $G$ is cyclic.
"""
function iscyclic(G::GrpAbFinGen)
  if !issnf(G)
    S = snf(G)[1]
    return ngens(S) == 1
  else
    return ngens(G) == 1
  end
end

################################################################################
#
#  p-Sylow subgroup
#
################################################################################

function _psylow_subgroup_gens(G::GrpAbFinGen, p::Union{fmpz, Integer})
  @assert issnf(G)
  z = GrpAbFinGenElem[]
  for i in 1:ngens(G)
    k, m = remove(G.snf[i], p)
    if k != 0
      push!(z, m*G[i])
    end
  end
  return z
end

function psylow_subgroup(G::GrpAbFinGen, p::Union{fmpz, Integer},
                         to_lattice::Bool = true)
  S, mS = snf(G)
  z = _psylow_subgroup_gens(S, p)
  zz = [ image(mS, x) for x in z ]
  return sub(G, zz, to_lattice)
end

################################################################################
#
#  Some special abelian groups
#
################################################################################

# TH: Isn't this the same as UnitsModM.jl?
# TODO: remove this from here. It does not belong here
doc"""
***
    multgrp_of_cyclic_grp(n::fmpz) -> GrpAbFinGen

> Returns the multiplicative group of the cyclic group with $n$ elements.
"""
function multgrp_of_cyclic_grp(n::fmpz)
  composition = Array{fmpz,1}()
  for (p,mp) in factor(n)
    if (p == 2) && (mp >= 2)
      push!(composition,2)
      push!(composition,fmpz(2)^(mp-2))
    else
      push!(composition,(p-1)*p^(mp-1))
    end
  end
  return DiagonalGroup(composition)
end

multgrp_of_cyclic_grp(n::Integer) = multgrp_of_cyclic_grp(fmpz(n))


text:
using Dates
using Flux
using Flux.Data: DataLoader

# Represent a single sequence with observation, target and mask
struct RepeatCopy
    obs
    target
    mask
end

_normalise(r, normval) = r/normval
_unnormalise(r, normval) = r*normval

function RepeatCopy(;
        nbits=4,
        minlength=1,
        maxlength=2,
        minrepeats=1,
        maxrepeats=2)
    # Reserve one channel for start marker and one for num-repeats marker
    fullobssize = nbits+2
    # Reserve one target channel for the end-marker
    fulltargsize = nbits+1
    timingmarker_idx = nbits+1
    numrepeatsmarker_idx = fullobssize

    # Sample variables for this sequence
    seqlength = rand(minlength:maxlength)
    numrepeats = rand(minrepeats:maxrepeats)
    # Get total size of sequence with padding
    totallength = (maxrepeats+1)*maxlength + 3

    obspattern = rand([0.f0,1.f0], nbits, seqlength) #Use Float32
    targpattern = repeat(obspattern, 1, numrepeats)

    # Attach start- and numrepeats markers to observation pattern
    startflag = Flux.onehot(timingmarker_idx, 1:fullobssize)
    repeatsflag = zeros(eltype(obspattern),  fullobssize)
    repeatsflag[numrepeatsmarker_idx] = _normalise(numrepeats, 10)
    obspattern = pad(obspattern, fullobssize, seqlength)
    obspattern = hcat(startflag, obspattern, repeatsflag)
    # Pad observation with space for the target
    obs = pad(obspattern, fullobssize, totallength)

    # Same for target pattern, but with only one marker
    endflag = Flux.onehot(timingmarker_idx, 1:fulltargsize)
    targpattern = pad(targpattern, fulltargsize, size(targpattern)[2])
    targpattern = hcat(targpattern, endflag)
    targ = pad(targpattern, fulltargsize, size(obspattern)[2]+size(targpattern)[2]; left=true)
    targ = pad(targ, fulltargsize, totallength)

    # Create mask
    mask = [zeros(Float32, size(obspattern)[2]); ones(Float32, size(targpattern)[2])]
    mask = [mask; zeros(Float32, totallength-length(mask))]
    RepeatCopy(obs, targ, mask)
end

import Base.show
function Base.show(io::IO, rp::RepeatCopy)
    print(io, "RepeatCopy(obssize: $(size(rp.obs)), targetsize $(size(rp.target))")
end

struct RepeatCopyBatchLoader
    dataloader::DataLoader
end

"""
    RepeatCopyBatchLoader(data...; batchsize=1, shuffle=false, partial=true)

An object that iterates over batches of RepeatCopy sequences.
The output of iteration are three list of matrices: observations, targets, and mask.
The lists represents the time steps, so the ``j``th column of observation matrix ``i`` is the ``i``th column of sequence ``j``.

See also: [`Flux.Data.DataLoader`](@ref), [`RepeatCopy`](@ref)

"""
RepeatCopyBatchLoader(a...; ka...) = RepeatCopyBatchLoader(DataLoader(a...; ka...))

function Base.iterate(d::RepeatCopyBatchLoader, i=0)
    i >= d.dataloader.imax && return nothing
    batch, nexti = Base.iterate(d.dataloader, i)
    batchsize = d.dataloader.batchsize
    obsbits, seqlength = size(batch[1].obs)
    targetbits = size(batch[1].target, 1)
    batchobsmatrix = [similar(batch[1].obs, eltype(batch[1].obs), (obsbits, batchsize))
                      for _ in 1:seqlength]
    @views for i in 1:seqlength
        for rc in 1:batchsize
            batchobsmatrix[i][:,rc] = batch[rc].obs[:, i]
        end
    end
    batchtargetmatrix = [similar(batch[1].target, eltype(batch[1].target), (targetbits, batchsize))
                      for _ in 1:seqlength]
    @views for i in 1:seqlength
        for rc in 1:batchsize
            batchtargetmatrix[i][:,rc] = batch[rc].target[:, i]
        end
    end
    batchmaskmatrix = [similar(batch[1].mask, eltype(batch[1].mask), (batchsize))
                      for _ in 1:seqlength]
    @views for i in 1:seqlength
        for rc in 1:batchsize
            batchmaskmatrix[i][rc] = batch[rc].mask[i]
        end
    end
    return ((batchobsmatrix, batchtargetmatrix, batchmaskmatrix), nexti)
end

Base.length(d::RepeatCopyBatchLoader) = length(d.dataloader)

function pad(arr::AbstractArray{T, 2}, dims::Integer...; left=false, padval=zero(T)) where T
    arrx, arry = size(arr)
    padlenx = dims[1] - arrx
    padleny = dims[2] - arry
    (padlenx < 0 || padleny < 0) && return arr
    hpadmat = fill(padval, (arrx, padleny))
    vpadmat = fill(padval, (padlenx, arry+padleny))
    if left
        return vcat(hcat(hpadmat, arr), vpadmat)
    else
        return vcat(hcat(arr, hpadmat), vpadmat)
    end
end

function printrow(row)
    row = map(row) do x
        x==0 && return "-"
        return "$x"
    end
    "+ $(join(row, " ")) +"
end

function prettyprint(data; isnormalized=false)
    r, c = size(data)
    for row in 1:r
        if isnormalized && row == r
            println(printrow(Int32.(_unnormalise.(data[row,:], 10))))
        else
            println(printrow(Int32.(data[row, :])))
        end
    end
    print("\n")
end

function prettyprint(rp::RepeatCopy)
    print("Observation:\n")
    prettyprint(rp.obs; isnormalized=true)
    print("\nTarget:\n")
    prettyprint(rp.target)
end
"""
    sigmoidlogitscrossentropy(logits, targets)

Formula taken from: https://www.tensorflow.org/api_docs/python/tf/nn/sigmoid_cross_entropy_with_logits
"""
function sigmoidlogitscrossentropy(logits, targets)
    x, z = logits, targets
    return max.(x, 0) .- x .* z + log.(ones(size(x)) + exp.(-abs.(x)))
end

function maskedsigmoidcrossentropy(logits, target, mask)
    l = sum(Flux.logitbinarycrossentropy.(logits, target); dims=1) * mask
    l[1]
end

function printmodeloutput(logits, mask)
    r, c  = size(logits)
    convertedmatrix = Int32.(repeat(mask', r) .* round.(Flux.σ.(logits)))
    print("Model output:\n")
    prettyprint(convertedmatrix)
end

function runmodel(model, rp::RepeatCopy)
    r, c = size(rp.obs)
    inputs = [rp.obs[:,col] for col in 1:c]
    logits = model.(inputs)
    logits = reshape(vcat(logits...), size(rp.target))
end

function loss(model, rp::RepeatCopy; printoutput=false)
    logits = runmodel(model, rp)
    l = maskedsigmoidcrossentropy(logits, rp.target, rp.mask)
    if printoutput
        prettyprint(rp)
        printmodeloutput(logits, rp.mask)
        println("Loss = $l")
    end
    l
end

function loss(model, batcheddata::Array{RepeatCopy, 1}; printoutput=false)
    batchloss = 0
    for rp in batcheddata
        l = loss(model, rp; printoutput=printoutput)
        printoutput=false # only print from first batch
        batchloss += l
    end
    batchloss
end

"""
This loss function is intented to handle the output of a RepeatCopyBatchLoader iteration.
"""
function loss(model, obs::AbstractArray{T, 1},
              target::AbstractArray{T, 1},
              mask::AbstractArray
              ; printoutput = false) where {T<:AbstractArray{E, 2}} where E
    length(obs) == length(target) == length(mask) || throw(ArgumentError("Observations, targets and mask needs to be of same length"))
    logits = model.(obs)
    l = maskedsigmoidcrossentropy.(logits, target, mask)
    if printoutput
        first_rc = reconstruct_rc(obs, target, mask, 1)
        prettyprint(first_rc)
        first_logits = getcolumn(logits, 1)
        printmodeloutput(first_logits, first_rc.mask)
        println("Average loss = ", sum(l)/size(obs[1], 2))
        println("Total batch loss = ", sum(l))
    end
    sum(l)
end

getcolumn(arr::AbstractArray{T, 1}, colno::Integer) where {T<:AbstractArray} =
    reshape(reduce(vcat, [arr[i][:, colno] for i in 1:length(arr)]), (size(arr[1], colno), length(arr)))

function reconstruct_rc(obs, target, mask, rcno::Integer)
    rcobs = getcolumn(obs, rcno)
    rctarget = getcolumn(target, rcno)
    rcmask = [mask[i][rcno] for i in 1:length(mask)]
    return RepeatCopy(rcobs, rctarget, rcmask)
end

using BSON: @save
using Flux: @progress, throttle
using Flux.Optimise: update!, runall, StopException
using Zygote: Params, gradient
using Dates

function mytrain!(loss, ps, data, opt; cb=()->())
    ps = Params(ps)
    cb = runall(cb)
    for d in data
        try
            gs = gradient(ps) do
                loss(d)
            end
            update!(opt, ps, gs)
            cb()
        catch ex
            if ex isa StopException
                break
            else
                rethrow(ex)
            end
        end
    end
    @save "$(today())-repeatcopy-$niter-$batchsize.bson" model
end

# Allow evalcb to report based on iteration number
mutable struct ThrottleIterations
    fn
    interval
    cnt
end

ThrottleIterations(fn, interval) = ThrottleIterations(fn, interval, 0)

function (fncint::ThrottleIterations)(a...;ka...)
    fncint.cnt += 1
    if fncint.cnt % fncint.interval == 0
        println(now(), ", iteration number ", fncint.cnt)
        return fncint.fn(a...,ka...)
    else
        return nothing
    end
end


text:
#Set Gunrock path
let
    global gunrock_lib
    succeeded = false
    if !isdefined(:gunrock_lib)
        @linux_only lib = "libgunrock.so"
        Libdl.dlopen(lib)
        succeeded = true
        succeeded || error("Gunrock library not found")
        @eval const gunrock_lib = $lib
    end
end


text:
# The uncertainty handling framework in this file will be added
# as part of a 1.X release. Can be ignored for now.

using CausalityTools, UncertainData

n = 20

# Mean values for two time series x and y, and standard deviations for those values
vals_x, stds_x = rand(n), rand(n) * 0.1
vals_y, stds_y = rand(n), rand(n) * 0.1
vals_z, stds_z = rand(n), rand(n) * 0.1

# Represent values as normal distributions
uvals_x = [UncertainValue(Normal, vals_x[i], stds_x[i]) for i = 1:n]
uvals_y = [UncertainValue(Normal, vals_y[i], stds_y[i]) for i = 1:n]
uvals_z = [UncertainValue(Normal, vals_z[i], stds_z[i]) for i = 1:n]

X = UncertainValueDataset(uvals_x)
Y = UncertainValueDataset(uvals_y)
Z = UncertainValueDataset(uvals_z)


text:
function construct_plankton(arch::Architecture, sp::Int64, params::Dict, maxN)
    rawdata = StructArray(x   = zeros(maxN), y   = zeros(maxN), z   = zeros(maxN),
                          xi  = zeros(Int,maxN), yi  = zeros(Int,maxN), zi  = zeros(Int,maxN),
                          iS  = zeros(maxN), Sz  = zeros(maxN), Bm  = zeros(maxN), Chl = zeros(maxN),
                          gen = zeros(maxN), age = zeros(maxN), ac  = zeros(maxN), idx = zeros(maxN),
                          PS  = zeros(maxN), resp= zeros(maxN), Th  = zeros(maxN),
                          graz= zeros(maxN), mort= zeros(maxN), dvid= zeros(maxN)
                          )
    data = replace_storage(array_type(arch), rawdata)

    param_names=(:Nsuper, :Cquota, :mean, :var, :α, :Φ, :T⁺, :Ea,
                 :PCmax, :PC_b, :Chl2C, :respir_a, :respir_b,
                 :grz_P, :dvid_type, :dvid_P, :dvid_stp, :dvid_reg, :dvid_stp2, :dvid_reg2,
                 :mort_P, :mort_reg, :grazFracC, :mortFracC, :ther_mort)

    pkeys = collect(keys(params))
    tmp = zeros(length(param_names))
    for i in 1:length(param_names)
        if length(findall(x->x==string(param_names[i]),pkeys))==0
            throw(ArgumentError("PARAM: parameter not found $(param_names[i])"))
        else
            tmp[i] = params[string(param_names[i])][sp]
        end
    end
    p = NamedTuple{param_names}(tmp)
    return plankton(data, p)
end

function generate_plankton!(plank, N::Int64, g::AbstractGrid, arch::Architecture)
    mean = plank.p.mean
    var = plank.p.var
    Cquota = plank.p.Cquota
    Nsuper = plank.p.Nsuper
    Chl2C = plank.p.Chl2C

    plank.data.ac[1:N]  .= 1.0                                                                             # activity
    plank.data.gen[1:N] .= 1.0                                                                             # generation
    plank.data.age[1:N] .= 0.0                                                                             # age

    randn!(rng_type(arch), plank.data.iS)
    rand!(rng_type(arch), plank.data.x)
    rand!(rng_type(arch), plank.data.y)
    rand!(rng_type(arch), plank.data.z)

    plank.data.x   .=(plank.data.x .* g.Nx) .* plank.data.ac                                         # x, unit: grid spacing, starting from 0
    plank.data.y   .=(plank.data.y .* g.Ny) .* plank.data.ac                                         # y, unit: grid spacing, starting from 0
    plank.data.z   .=(plank.data.z .* g.Nz) .* plank.data.ac                                         # z, unit: grid spacing, starting from 0
    plank.data.iS  .= max.(1.0, plank.data.iS .* var .+ mean) .* plank.data.ac                       # init_size
    plank.data.Sz  .= copy(plank.data.iS)                                                            # size
    plank.data.Bm  .= Cquota .* plank.data.Sz .* Nsuper                                              # Bm
    plank.data.Chl .= plank.data.Bm .* Chl2C                                                         # Chl

    mask_individuals!(plank.data, g, N, arch)
end


text:
using Test, DataFrames, CSV, BDisposal
println("Testing BDisposal on js-data...")

# oecd-data test with bad inputs separated...
data = CSV.read(joinpath(dirname(pathof(BDisposal)),"..","test","data","js-data","oecd.txt"),DataFrame; delim=' ',ignorerepeated=true,copycols=true,header=false)
rename!(data,[:ccid,:year,:gdp,:co2,:capital,:labour,:energy])
goodInputsLabels  = ["capital","labour"]
badInputsLabels   = ["energy"]
goodOutputsLabels = ["gdp"]
badOutputsLabels  = ["co2"]
sort!(data, [:year, :ccid]) # sort data by period and dmu
periods = unique(data.year)
dmus    = unique(data.ccid)

nGI, nBI, nGO, nBO, nPer, nDMUs,  = length(goodInputsLabels), length(badInputsLabels), length(goodOutputsLabels), length(badOutputsLabels), length(periods),length(dmus)

gI = Array{Float64}(undef, (nDMUs,nGI,nPer))
bI = Array{Float64}(undef, (nDMUs,nBI,nPer))
gO = Array{Float64}(undef, (nDMUs,nGO,nPer))
bO = Array{Float64}(undef, (nDMUs,nBO,nPer))

myeps = 0.000000001
for (p,period) in enumerate(periods)
    periodData = data[data.year .== period,:]
    gI[:,:,p] = Matrix{Float64}(periodData[:,goodInputsLabels])
    if nBI > 0
        bI[:,:,p] =
         Matrix{Float64}(periodData[:,badInputsLabels])
    end
    gO[:,:,p] = Matrix{Float64}(periodData[:,goodOutputsLabels])
    bO[:,:,p] = Matrix{Float64}(periodData[:,badOutputsLabels])
end

#=
using Distributions
ud = Uniform(1,100000)

gI = rand(ud,nDMUs,nGI,nPer)
bI = rand(ud,nDMUs,nBI,nPer)
gO = rand(ud,nDMUs,nGO,nPer)
bO = rand(ud,nDMUs,nBO,nPer)
=#

t = length(periods)-1
gIₜ = gI[:,:,t]
gIᵤ = gI[:,:,t+1]
bIₜ = bI[:,:,t]
bIᵤ = bI[:,:,t+1]
gOₜ = gO[:,:,t]
gOᵤ = gO[:,:,t+1]
bOₜ = bO[:,:,t]
bOᵤ = bO[:,:,t+1]

z = 2
gI₀ₜ = gI[z,:,t]
gI₀ᵤ = gI[z,:,t+1]
bI₀ₜ = bI[z,:,t]
bI₀ᵤ = bI[z,:,t+1]
gO₀ₜ = gO[z,:,t]
gO₀ᵤ = gO[z,:,t+1]
bO₀ₜ = bO[z,:,t]
bO₀ᵤ = bO[z,:,t+1]

convex = false

(dirGIm,dirGIa) = (-1,0,0,0) , (1,0,0,0)
(dirBIm,dirBIa) = (0,-1,0,0) , (0,1,0,0)
(dirGOm,dirGOa) = (0,0,1,0)  , (0,0,1,0)
(dirBOm,dirBOa) = (0,0,0,-1) , (0,0,0,1)

#=
gI_ut_mult = []
for z in 1:nDMUs

gI₀ₜ = gI[z,:,t]
gI₀ᵤ = gI[z,:,t+1]
bI₀ₜ = bI[z,:,t]
bI₀ᵤ = bI[z,:,t+1]
gO₀ₜ = gO[z,:,t]
gO₀ᵤ = gO[z,:,t+1]
bO₀ₜ = bO[z,:,t]
bO₀ᵤ = bO[z,:,t+1]

gI_ut_mult_z = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
  directions=dirGIm,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)
push!(gI_ut_mult,gI_ut_mult_z)
end

gI_ut_mult
=#

for z in 1:nDMUs

println(z)

gI₀ₜ = gI[z,:,t]
gI₀ᵤ = gI[z,:,t+1]
bI₀ₜ = bI[z,:,t]
bI₀ᵤ = bI[z,:,t+1]
gO₀ₜ = gO[z,:,t]
gO₀ᵤ = gO[z,:,t+1]
bO₀ₜ = bO[z,:,t]
bO₀ᵤ = bO[z,:,t+1]

 # Everything at time t....
gI_mult = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
  directions=dirGIm,startValues=(),forceLinearModel=true,convexAssumption=convex)
gI_add = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="additive",
  directions=dirGIa,startValues=(),forceLinearModel=true,convexAssumption=convex)
@test gI_mult >= 1.0 - myeps
@test begin (gI_mult ≈ 1.0) ?  (0.0 - myeps) < gI_add < (0.0 + myeps) : gI_add > 0.0 - myeps end

bI_mult = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
  directions=dirBIm,startValues=(),forceLinearModel=true,convexAssumption=convex)
bI_add = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="additive",
  directions=dirBIa,startValues=(),forceLinearModel=true,convexAssumption=convex)
@test bI_mult >= 1.0 - myeps
@test begin (bI_mult ≈ 1.0) ?  (0.0 - myeps) < bI_add < (0.0 + myeps) : bI_add > 0.0 - myeps end

gO_mult = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
 directions=dirGOm,startValues=(),forceLinearModel=true,convexAssumption=convex)
gO_add = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="additive",
 directions=dirGOa,startValues=(),forceLinearModel=true,convexAssumption=convex)
@test gO_mult >= 1.0 - myeps
@test begin (gO_mult ≈ 1.0) ?  (0.0 - myeps) < gO_add < (0.0 + myeps) : gO_add > 0.0 - myeps end

bO_mult = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
 directions=dirBOm,startValues=(),forceLinearModel=true,convexAssumption=convex)
bO_add = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="additive",
 directions=dirBOa,startValues=(),forceLinearModel=true,convexAssumption=convex)
@test bO_mult >= 1.0 - myeps
@test begin (bO_mult ≈ 1.0) ?  (0.0 - myeps) < bO_add < (0.0 + myeps) : bO_add > 0.0 - myeps end

# Observation at time t+1, all other dmu at time t....
gI_t̃_mult = problem(gI₀ᵤ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
 directions=dirGIm,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
gI_t̃_add = problem(gI₀ᵤ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="additive",
 directions=dirGIa,startValues=(),forceLinearModel=true,crossTime=true,convexAssumption=convex)
if convex
  @test  isapprox(gI_t̃_add,1 - 1/gI_t̃_mult, atol=0.000001)
else
  @test  isapprox(gI_t̃_add,1 - gI_t̃_mult, atol=0.000001)
end

bI_t̃_mult = problem(gI₀ₜ,bI₀ᵤ,gO₀ₜ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
 directions=dirBIm,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
bI_t̃_add = problem(gI₀ₜ,bI₀ᵤ,gO₀ₜ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="additive",
 directions=dirBIa,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
if convex
  @test bI_t̃_add ≈ 1 - 1/bI_t̃_mult
else
  @test bI_t̃_add ≈ 1 - bI_t̃_mult
end

gO_t̃_mult = problem(gI₀ₜ,bI₀ₜ,gO₀ᵤ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
 directions=dirGOm,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
gO_t̃_add = problem(gI₀ₜ,bI₀ₜ,gO₀ᵤ,bO₀ₜ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="additive",
 directions=dirGOa,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
@test isapprox(gO_t̃_add,(gO_t̃_mult -1);atol=0.0001)

bO_t̃_mult = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ᵤ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
 directions=dirBOm,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
bO_t̃_add = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ᵤ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="additive",
 directions=dirBOa,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
@test bO_t̃_add ≈ 1 - 1/bO_t̃_mult

# Observation at time t, all other dmu at time t+1....
gI_ũ_mult = problem(gI₀ₜ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="multiplicative",
 directions=dirGIm,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
gI_ũ_add = problem(gI₀ₜ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="additive",
 directions=dirGIa,startValues=(),forceLinearModel=true,crossTime=true,convexAssumption=convex)
if convex
  @test  isapprox(gI_ũ_add,1 - 1/gI_ũ_mult, atol=0.000001)
else
  @test  isapprox(gI_ũ_add,1 - gI_ũ_mult, atol=0.000001)
end

bI_ũ_mult = problem(gI₀ᵤ,bI₀ₜ,gO₀ᵤ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="multiplicative",
 directions=dirBIm,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
bI_ũ_add = problem(gI₀ᵤ,bI₀ₜ,gO₀ᵤ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="additive",
 directions=dirBIa,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
if convex
  @test isapprox(bI_ũ_add, 1 - 1/bI_ũ_mult, atol=0.0000001)
else
  @test isapprox(bI_ũ_add, 1 - bI_ũ_mult, atol=0.0000001)
end

gO_ũ_mult = problem(gI₀ᵤ,bI₀ᵤ,gO₀ₜ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="multiplicative",
 directions=dirGOm,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
gO_ũ_add = problem(gI₀ᵤ,bI₀ᵤ,gO₀ₜ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="additive",
 directions=dirGOa,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
@test isapprox(gO_ũ_add, (gO_ũ_mult -1), atol = 0.0001)

bO_ũ_mult = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ₜ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="multiplicative",
 directions=dirBOm,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
bO_ũ_add = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ₜ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="additive",
 directions=dirBOa,startValues=(),forceLinearModel=true, crossTime=true,convexAssumption=convex)
@test bO_ũ_add ≈ 1 - 1/bO_ũ_mult

# Everything at time t+1 (u)....

gI_u_mult = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="multiplicative",
  directions=dirGIm,startValues=(),forceLinearModel=true,convexAssumption=convex)
gI_u_add = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="additive",
  directions=dirGIa,startValues=(),forceLinearModel=true,convexAssumption=convex)

bI_u_mult = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="multiplicative",
  directions=dirBIm,startValues=(),forceLinearModel=true,convexAssumption=convex)
bI_u_add = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="additive",
  directions=dirBIa,startValues=(),forceLinearModel=true,convexAssumption=convex)

gO_u_mult = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="multiplicative",
 directions=dirGOm,startValues=(),forceLinearModel=true,convexAssumption=convex)
gO_u_add = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="additive",
 directions=dirGOa,startValues=(),forceLinearModel=true,convexAssumption=convex)

bO_u_mult = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="multiplicative",
 directions=dirBOm,startValues=(),forceLinearModel=true,convexAssumption=convex)
bO_u_add = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="additive",
 directions=dirBOa,startValues=(),forceLinearModel=true,convexAssumption=convex)

# DMU measures at time t+1 and frontier at time t...
gI_ut_mult = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
  directions=dirGIm,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)
gI_ut_add = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="additive",
  directions=dirGIa,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)

if convex
  @test  isapprox(gI_ut_add,1 - 1/gI_ut_mult, atol=0.000001)
else
  @test  isapprox(gI_ut_add,1 - gI_ut_mult, atol=0.000001)
end

bI_ut_mult = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
  directions=dirBIm,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)
bI_ut_add = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="additive",
  directions=dirBIa,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)

gO_ut_mult = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
 directions=dirGOm,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)
gO_ut_add = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="additive",
 directions=dirGOa,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)

bO_ut_mult = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="multiplicative",
 directions=dirBOm,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)
bO_ut_add = problem(gI₀ᵤ,bI₀ᵤ,gO₀ᵤ,bO₀ᵤ,gIₜ,bIₜ,gOₜ,bOₜ,retToScale="variable",prodStructure="additive",
 directions=dirBOa,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)

# DMU measures at time t and frontier at time t+1...
gI_tu_mult = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="multiplicative",
  directions=dirGIm,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)
gI_tu_add = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="additive",
  directions=dirGIa,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)

if convex
  @test  isapprox(gI_tu_add,1 - 1/gI_tu_mult, atol=0.000001)
else
  @test  isapprox(gI_tu_add,1 - gI_tu_mult, atol=0.000001)
end

bI_tu_mult = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="multiplicative",
  directions=dirBIm,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)
bI_tu_add = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="addittive",
  directions=dirBIa,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)

gO_tu_mult = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="multiplicative",
 directions=dirGOm,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)
gO_tu_add = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="additive",
 directions=dirGOa,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)

bO_tu_mult = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="multiplicative",
 directions=dirBOm,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)
bO_tu_add = problem(gI₀ₜ,bI₀ₜ,gO₀ₜ,bO₀ₜ,gIᵤ,bIᵤ,gOᵤ,bOᵤ,retToScale="variable",prodStructure="additive",
 directions=dirBOa,startValues=(),forceLinearModel=true,convexAssumption=convex,crossTime=true)

################################
############# M - MULTIPLICATIVE CASE

# M.1 Data copy
idx_gi_t = gI_mult
idx_bi_t = bI_mult
idx_go_t = gO_mult
idx_bo_t = bO_mult
idx_gi_t̃ = gI_t̃_mult
idx_bi_t̃ = bI_t̃_mult
idx_go_t̃ = gO_t̃_mult
idx_bo_t̃ = bO_t̃_mult
idx_gi_u = gI_u_mult
idx_bi_u = bI_u_mult
idx_go_u = gO_u_mult
idx_bo_u = bO_u_mult

idx_gi_ut,idx_bi_ut,idx_go_ut,idx_bo_ut = gI_ut_mult,bI_ut_mult,gO_ut_mult,bO_ut_mult
idx_gi_tu,idx_bi_tu,idx_go_tu,idx_bo_tu = gI_tu_mult,bI_tu_mult,gO_tu_mult,bO_tu_mult

idx_gi_ũ = gI_ũ_mult
idx_bi_ũ = bI_ũ_mult
idx_go_ũ = gO_ũ_mult
idx_bo_ũ = bO_ũ_mult

############################
# M.2 Convex aggregation/dissaggregations
# Convex case aggregation

idx_i_t = (idx_gi_t̃/idx_gi_t) * (idx_bi_t̃/idx_bi_t)
idx_o_t = (idx_go_t/idx_go_t̃) * (idx_bo_t/idx_bo_t̃)
idx_t   = idx_o_t/idx_i_t
idx_i_u = (idx_gi_u/idx_gi_ũ) * (idx_bi_u/idx_bi_ũ)
idx_o_u = (idx_go_ũ/idx_go_u) * (idx_bo_ũ/idx_bo_u)
idx_u   = idx_o_u/idx_i_u
idx     = (idx_t * idx_u)^(1/2)

idx_go_t/idx_go_t̃
idx_bo_t̃/idx_bo_t

1/idx_o_t

# Disaggregation  Good/ bads

idx_Gi_t = (idx_gi_t̃/idx_gi_t)
idx_Go_t = (idx_go_t/idx_go_t̃)
idx_Gt   = idx_Go_t/idx_Gi_t
idx_Gi_u = (idx_gi_u/idx_gi_ũ)
idx_Go_u = (idx_go_ũ/idx_go_u)
idx_Gu   = idx_Go_u/idx_Gi_u
idx_G    = (idx_Gt * idx_Gu)^(1/2)

idx_Bi_t = (idx_bi_t̃/idx_bi_t)
idx_Bo_t = (idx_bo_t̃/idx_bo_t)
idx_Bt   = 1/(idx_Bo_t * idx_Bi_t)
idx_Bi_u = (idx_bi_u/idx_bi_ũ)
idx_Bo_u = (idx_bo_u/idx_bo_ũ)
idx_Bu   = 1/(idx_Bo_u*idx_Bi_u)
idx_B    = (idx_Bt * idx_Bu)^(1/2)

# Disaggregation T/E/S...
idx_T_G_O = ((idx_go_t/idx_go_tu) * (idx_go_ut/idx_go_u) )^-(1/2)
idx_T_B_O = ((idx_bo_t/idx_bo_tu) * (idx_bo_ut/idx_bo_u) )^-(1/2)
idx_T_O   = idx_T_G_O * idx_T_B_O
idx_T_G_I = ((idx_gi_t/idx_gi_tu) * (idx_gi_ut/idx_gi_u) )^-(1/2)
idx_T_B_I = ((idx_bi_t/idx_bi_tu) * (idx_bi_ut/idx_bi_u) )^-(1/2)
idx_T_I   = idx_T_G_I * idx_T_B_I
idx_T     = (idx_T_O * idx_T_I)

idx_E_G_O = (idx_go_u/idx_go_t)^-1
idx_E_B_O = (idx_bo_u/idx_bo_t)^-1
idx_E_O   = idx_E_G_O * idx_E_B_O
idx_E_G_I = (idx_gi_u/idx_gi_t)^-1
idx_E_B_I = (idx_bi_u/idx_bi_t)^-1
idx_E_I   = idx_E_G_I * idx_E_B_I
idx_E     = (idx_E_O * idx_E_I)

idx_S_I   = idx / (idx_T_I * idx_E_I)
idx_S_O   = idx / (idx_T_O * idx_E_O)
idx_S_G_O = ((idx_go_t̃ / idx_go_ut) * (idx_gi_t̃ / idx_gi_t) * (idx_go_tu / idx_go_ũ) * (idx_gi_u/idx_gi_ũ)  )^-(1/2)
idx_S_B_O = ((idx_bo_t̃ / idx_bo_ut) * (idx_bi_t̃ / idx_bi_t) * (idx_bo_tu / idx_bo_ũ) * (idx_bi_u/idx_bi_ũ)  )^-(1/2)
@test idx_S_O ≈ idx_S_G_O * idx_S_B_O
@test idx ≈ idx_T_O * idx_E_O * idx_S_O

idx_S_G_I = ((idx_gi_t̃ / idx_gi_ut) * (idx_go_t̃ / idx_go_t) * (idx_gi_tu / idx_gi_ũ) * (idx_go_u/idx_go_ũ)  )^-(1/2)
idx_S_B_I = ((idx_bi_t̃ / idx_bi_ut) * (idx_bo_t̃ / idx_bo_t) * (idx_bi_tu / idx_bi_ũ) * (idx_bo_u/idx_bo_ũ)  )^-(1/2)
@test idx_S_I ≈ idx_S_G_I * idx_S_B_I
@test idx ≈ idx_T_I * idx_E_I * idx_S_I

idx_S_G_O = idx_G /  (idx_T_G_O * idx_E_G_O)
idx_S_B_O = idx_B /  (idx_T_B_O * idx_E_B_O)
idx_S_O   = idx_S_G_O * idx_S_B_O
idx_S_G_I = idx_G /  (idx_T_G_I * idx_E_G_I)
idx_S_B_I = idx_B /  (idx_T_B_I * idx_E_B_I)
idx_S_I   = idx_S_G_I * idx_S_B_I
idx_S     = (idx_S_O * idx_S_I)^(1/2)

############################
# M.3 Non-convex case aggregation
idx_i_t = (idx_gi_t/idx_gi_t̃) * (idx_bi_t/idx_bi_t̃)
idx_o_t = (idx_go_t/idx_go_t̃) / (idx_bo_t̃/idx_bo_t)
idx_t   = idx_o_t/idx_i_t
idx_i_u = (idx_gi_ũ/idx_gi_u) * (idx_bi_ũ/idx_bi_u)
idx_o_u = (idx_go_ũ/idx_go_u) / (idx_bo_u/idx_bo_ũ)
idx_u   = idx_o_u/idx_i_u
idx     = (idx_t * idx_u)^(1/2)

# Disaggregation  Good/ bads

idx_Gi_t = 1/(idx_gi_t̃*idx_gi_t)
idx_Go_t = (idx_go_t/idx_go_t̃)
idx_Gt   = idx_Go_t/idx_Gi_t
idx_Gi_u = 1/(idx_gi_u/idx_gi_ũ)
idx_Go_u = (idx_go_ũ/idx_go_u)
idx_Gu   = idx_Go_u/idx_Gi_u
idx_G    = (idx_Gt * idx_Gu)^(1/2)

idx_Bi_t = 1/(idx_bi_t̃/idx_bi_t)
idx_Bo_t = 1/(idx_bo_t̃/idx_bo_t)
idx_Bt   = (idx_Bo_t / idx_Bi_t)
idx_Bi_u = 1/(idx_bi_u/idx_bi_ũ)
idx_Bo_u = 1/(idx_bo_u/idx_bo_ũ)
idx_Bu   = (idx_Bo_u/idx_Bi_u)
idx_B    = (idx_Bt * idx_Bu)^(1/2)

# A - ADDITTIVE CASE ################################

# A.1 Data copy
(idx_gi_t, idx_bi_t, idx_go_t, idx_bo_t) = (gI_add, bI_add, gO_add, bO_add)
(idx_gi_t̃, idx_bi_t̃, idx_go_t̃, idx_bo_t̃) = (gI_t̃_add, bI_t̃_add, gO_t̃_add, bO_t̃_add)
(idx_gi_u, idx_bi_u, idx_go_u, idx_bo_u) = (gI_u_add, bI_u_add, gO_u_add, bO_u_add)
(idx_gi_ũ, idx_bi_ũ, idx_go_ũ, idx_bo_ũ) = (gI_ũ_add, bI_ũ_add, gO_ũ_add, bO_ũ_add)
idx_gi_ut,idx_bi_ut,idx_go_ut,idx_bo_ut = gI_ut_add,bI_ut_add,gO_ut_add,bO_ut_add
idx_gi_tu,idx_bi_tu,idx_go_tu,idx_bo_tu = gI_tu_add,bI_tu_add,gO_tu_add,bO_tu_add

############
# A.2 Convex case

idx_i_t = (idx_gi_t̃ - idx_gi_t) + (idx_bi_t̃ - idx_bi_t)
idx_o_t = (idx_go_t - idx_go_t̃) - (idx_bo_t̃ - idx_bo_t)
idx_t   = idx_o_t - idx_i_t
idx_i_u = (idx_gi_u - idx_gi_ũ) + (idx_bi_u - idx_bi_ũ)
idx_o_u = (idx_go_ũ - idx_go_u) - (idx_bo_u - idx_bo_ũ)
idx_u   = idx_o_u - idx_i_u
idx     = (idx_t + idx_u) / 2

idx_T_G_O = ((-idx_go_t+idx_go_tu) + (-idx_go_ut+idx_go_u) )/2
idx_T_B_O = ((-idx_bo_t-idx_bo_tu) + (+idx_bo_ut+idx_bo_u) )/2
idx_T_O   = idx_T_G_O + idx_T_B_O
idx_T_G_I = ((-idx_gi_t+idx_gi_tu) + (-idx_gi_ut+idx_gi_u) )/2
idx_T_B_I = ((-idx_bi_t+idx_bi_tu) + (-idx_bi_ut+idx_bi_u) )/2
idx_T_I   = idx_T_G_I + idx_T_B_I
idx_T     = (idx_T_O + idx_T_I)

idx_E_G_O = -(idx_go_u-idx_go_t)
idx_E_B_O = -(idx_bo_u-idx_bo_t)
idx_E_O   = idx_E_G_O + idx_E_B_O
idx_E_G_I = -(idx_gi_u-idx_gi_t)
idx_E_B_I = -(idx_bi_u-idx_bi_t)
idx_E_I   = idx_E_G_I + idx_E_B_I
idx_E     = (idx_E_O + idx_E_I)

idx_S_I   = idx - (idx_T_I + idx_E_I)
idx_S_O   = idx - (idx_T_O + idx_E_O)
idx_S_G_O = -((idx_go_t̃ - idx_go_ut) + (idx_gi_t̃ - idx_gi_t) + (idx_go_tu - idx_go_ũ) + (idx_gi_u-idx_gi_ũ)  )/2
idx_S_B_O = -((idx_bo_t̃ + idx_bo_ut) + (idx_bi_t̃ - idx_bi_t) + (-idx_bo_tu - idx_bo_ũ) + (idx_bi_u-idx_bi_ũ)  )/2
idx_S_G_I = -((idx_gi_t̃ - idx_gi_ut) + (idx_go_t̃ - idx_go_t) + (idx_gi_tu - idx_gi_ũ) + (idx_go_u-idx_go_ũ)  )/2
idx_S_B_I = -((idx_bi_t̃ - idx_bi_ut) + (idx_bo_t̃ - idx_bo_t) + (idx_bi_tu - idx_bi_ũ) + (idx_bo_u-idx_bo_ũ)  )/2

##################
# A.3 Non-convex case
idx_i_t = -(idx_gi_t-idx_gi_t̃) - (idx_bi_t-idx_bi_t̃)
idx_o_t = (idx_go_t-idx_go_t̃) - (idx_bo_t̃-idx_bo_t)
idx_t   = idx_o_t-idx_i_t
idx_i_u = -(idx_gi_ũ-idx_gi_u) - (idx_bi_ũ-idx_bi_u)
idx_o_u = (idx_go_ũ-idx_go_u) - (idx_bo_u-idx_bo_ũ)
idx_u   = idx_o_u-idx_i_u
idx     = (idx_t + idx_u)/2

# Disaggregation  Good/ bads

idx_Gi_t = (idx_gi_t̃-idx_gi_t)
idx_Go_t = (idx_go_t-idx_go_t̃)
idx_Gt   = idx_Go_t-idx_Gi_t
idx_Gi_u = (idx_gi_u-idx_gi_ũ)
idx_Go_u = (idx_go_ũ-idx_go_u)
idx_Gu   = idx_Go_u-idx_Gi_u
idx_G    = (idx_Gt + idx_Gu)/2

idx_Bi_t = (idx_bi_t̃-idx_bi_t)
idx_Bo_t = (idx_bo_t̃-idx_bo_t)
idx_Bt   = (-idx_Bo_t - idx_Bi_t)
idx_Bi_u = (idx_bi_u-idx_bi_ũ)
idx_Bo_u = (idx_bo_u-idx_bo_ũ)
idx_Bu   = (-idx_Bo_u-idx_Bi_u)
idx_B    = (idx_Bt + idx_Bu)/2

################################################################################
###### OVERALL ANALYSIS ########################################################

oecdAnalysis  = prodIndex(gI,gO,bO,bI;
                   retToScale="variable",prodStructure="multiplicative",convexAssumption=true)

isapprox(oecdAnalysis.prodIndexes_G .* oecdAnalysis.prodIndexes_B, oecdAnalysis.prodIndexes, atol=0.000001)
#isapprox(oecdAnalysis.prodIndexes_T .* oecdAnalysis.prodIndexes_E .* oecdAnalysis.prodIndexes_S, oecdAnalysis.prodIndexes, atol=0.000001)

oecdAnalysis.prodIndexes

oecdAnalysis.prodIndexes_T

oecdAnalysis.prodIndexes_E

oecdAnalysis.prodIndexes_S

oecdAnalysisA  = prodIndex(gI,gO,bO,bI;
                   retToScale="variable",prodStructure="addittive",convexAssumption=true)

oecdAnalysisA.prodIndexes_T
oecdAnalysisA.prodIndexes_E
oecdAnalysisA.prodIndexes_S

isapprox(oecdAnalysisA.prodIndexes_G .+ oecdAnalysisA.prodIndexes_B, oecdAnalysisA.prodIndexes, atol=0.000001)
#isapprox(oecdAnalysisA.prodIndexes_T .+ oecdAnalysisA.prodIndexes_E .+ oecdAnalysisA.prodIndexes_S, oecdAnalysisA.prodIndexes, atol=0.000001)

################################################################################
# Non convex test ##############################################################

oecdAnalysis_nc  = prodIndex(gI,gO,bO,bI;
                   retToScale="variable",prodStructure="multiplicative",convexAssumption=false)

isapprox(oecdAnalysis_nc.prodIndexes_G .* oecdAnalysis_nc.prodIndexes_B, oecdAnalysis_nc.prodIndexes, atol=0.000001)

oecdAnalysis_ncA  = prodIndex(gI,gO,bO,bI;
                   retToScale="variable",prodStructure="additive",convexAssumption=false)

isapprox(oecdAnalysis_ncA.prodIndexes_G .+ oecdAnalysis_ncA.prodIndexes_B, oecdAnalysis_ncA.prodIndexes, atol=0.000001)

log.(oecdAnalysis_nc.prodIndexes)

oecdAnalysis_ncA.prodIndexes

log.(oecdAnalysis.prodIndexes)

#=
mBool = mult .> (1 - myeps)
aBool = add .> (0 - myeps)
mBool == aBool
mBool .== aBool
=#


text:
@testset "reservoir simple" begin
    res = Wflow.SimpleReservoir{Float64}(
        demand = [52.523],
        maxrelease = [420.184],
        maxvolume = [25_000_000.0],
        volume = [1.925e7],
        area = [1885665.353626924],
        targetfullfrac = [0.8],
        targetminfrac = [0.2425554726620697],
        precipitation = [4.2],
        evaporation = [1.5],
    )

    Wflow.update(res, 1, 100.0, 86400.0)
    @test res.outflow[1] ≈ 91.3783714867453
    @test res.volume[1] ≈ 2.0e7
    @test res.percfull[1] ≈ 0.80
    @test res.demandrelease[1] ≈ 52.5229994727611
    @test res.precipitation[1] ≈ 4.2
    @test res.evaporation[1] ≈ 1.5
end

@testset "natural lake" begin
    lake = Wflow.NaturalLake{Float64}(
        loc_id = [1],
        lowerlake_ind = [0],
        area = [180510409.0],
        threshold = [0.0],
        storfunc = [1],
        outflowfunc = [3],
        b = [0.22],
        e = [2.0],
        sh = [Wflow.DataFrame()],
        hq = [Wflow.DataFrame()],
        avg_waterlevel = [18.5],
        precipitation = [20.0],
        evaporation = [3.2],
    )

    Wflow.update(lake, 1, 2500.0, 181, 86400.0)
    @test lake.outflow[1] ≈ 85.31903276150577
    @test lake.storage[1] ≈ 3.551103576940606e9
    @test lake.waterlevel[1] ≈ 19.672569557695734
    @test lake.precipitation[1] ≈ 20.0
    @test lake.evaporation[1] ≈ 3.2
end


text:
#################### Model Graph ####################

#################### Display ####################

function draw(m::Model; filename::AbstractString="")
  dot = graph2dot(m)
  if length(filename) == 0
    print(dot)
  else
    if search(filename, '.') == 0
      filename = string(filename, ".dot")
    end
    f = open(filename, "w")
    write(f, dot)
    close(f)
  end
end

function graph(m::Model)
  g = graph(KeyVertex{Symbol}[], Edge{KeyVertex{Symbol}}[])
  lookup = Dict{Symbol, Int}()
  for key in keys(m, :all)
    lookup[key] = length(lookup) + 1
    add_vertex!(g, KeyVertex(lookup[key], key))
  end
  V = vertices(g)
  for key in keys(m)
    node = m[key]
    if isa(node, AbstractDependent)
      for src in node.sources
        add_edge!(g, V[lookup[src]], V[lookup[key]])
      end
    end
  end
  g
end

function graph2dot(m::Model)
  g = graph(m)
  io = IOBuffer()
  write(io, "digraph MambaModel {\n")
  deps = keys(m, :dependent)
  for v in vertices(g)
    attr = Dict{AbstractString, AbstractString}()
    if v.key in deps
      node = m[v.key]
      if isa(node, AbstractLogical)
        attr["shape"] = "diamond"
      elseif isa(node, AbstractStochastic)
        attr["shape"] = "ellipse"
      end
      if isempty(node.monitor)
        attr["style"] = "filled"
        attr["fillcolor"] = "gray85"
      end
    else
      attr["shape"] = "box"
      attr["style"] = "filled"
      attr["fillcolor"] = "gray85"
    end
    write(io, "\t\"")
    write(io, v.key)
    write(io, "\" [")
    write(io, join(map(x -> "$(x[1])=\"$(x[2])\"", attr), ", "))
    write(io, "];\n")
    for e in out_edges(v, g)
      t = target(e, g)
      write(io, "\t\t\"")
      write(io, v.key)
      write(io, "\" -> \"")
      write(io, t.key)
      write(io, "\";\n")
     end
  end
  write(io, "}\n")
  bytestring(io)
end

#################### Auxiliary Functions ####################

function any_stochastic(v::KeyVertex{Symbol}, g::AbstractGraph, m::Model)
  found = false
  for v in out_neighbors(v, g)
    if isa(m[v.key], AbstractStochastic) || any_stochastic(v, g, m)
      found = true
      break
    end
  end
  found
end

function gettargets(v::KeyVertex{Symbol}, g::AbstractGraph, m::Model,
                    terminal::Vector{Symbol}=keys(m, :stochastic))
  values = Symbol[]
  for v in out_neighbors(v, g)
    push!(values, v.key)
    if !(v.key in terminal)
      values = union(values, gettargets(v, g, m, terminal))
    end
  end
  values
end

function tsort{T}(g::AbstractGraph{KeyVertex{T}, Edge{KeyVertex{T}}})
  V = topological_sort_by_dfs(g)
  map(v -> v.key, V)
end

function tsort(m::Model)
  tsort(graph(m))
end


text:
## Convert_s functions for the mdp
# Converts from vector to a Scene
function POMDPs.convert_s(::Type{Scene}, s::AbstractArray{Float64}, mdp::MDP{Scene, A}) where A
    entities = []

    # Loop through all the agents of the mdp
    index = 1
    for agent in agents(mdp)
        ent = agent.vec_to_entity(s[index : index + agent.entity_dim - 1], id(agent), mdp.roadway, agent.model)
        index += agent.entity_dim
        !end_of_road(ent, mdp.roadway, mdp.end_of_road) && push!(entities, ent)
    end
    isempty(entities) ? Scene(typeof(sut(mdp).get_initial_entity())) : Scene([entities...])
end

# Convert from Scene to a vector
function POMDPs.convert_s(::Type{AbstractArray}, state::Scene, mdp::MDP{Scene, A}) where A
    isempty(mdp.last_observation) && (mdp.last_observation = zeros(sum([a.entity_dim for a in agents(mdp)])))
    index = 1
    for a in agents(mdp)
        !(id(a) in state) && continue
        entity = get_by_id(state, id(a))
        mdp.last_observation[index : index + a.entity_dim - 1] .= a.entity_to_vec(entity)
        index += a.entity_dim
    end
    copy(mdp.last_observation)
end

POMDPs.convert_s(::Type{Array{Float64, 1}}, state::Scene, mdp::MDP{Scene, A}) where A = convert_s(AbstractArray, state, mdp)
POMDPs.convert_s(::Type{Array{Float32, 1}}, state::Scene, mdp::MDP{Scene, A}) where A = convert(Array{Float32,1}, convert_s(AbstractArray, state, mdp))
POMDPs.convert_s(::Type{Array{Float32}}, state::Scene, mdp::MDP{Scene, A}) where A = convert_s(Array{Float32,1}, state, mdp)

## Adversarial Pedestrians vehicles
const PEDESTRIAN_ENTITY_DIM = 4

# Converts state of a pedestrian to a vector
function NoisyPedestrian_to_vec(ped::Entity{NoisyPedState, VehicleDef, Int64})
    Float64[posg(ped)..., vel(ped)]
end

function vec_to_NoisyPedestrian_fn(crosswalk_id::Int)
    function vec_to_NoisyPedestrian(arr::AbstractArray, id, roadway::Roadway, model)
        @assert length(arr) == PEDESTRIAN_ENTITY_DIM

        pos = VecSE2(arr[1], arr[2], arr[3]) # Distance along the lane
        v = arr[4]

        vs = VehicleState(pos, roadway_lane(roadway, crosswalk_id), roadway, v)
        bv = Entity(NoisyPedState(vs, Noise()), PEDESTRIAN_DEF, id)
    end
end

## Blinker vehicles
const BLINKERVEHICLE_ENTITY_DIM = 4
const BLINKERVEHICLE_EXPANDED_ENTITY_DIM = 30

# Converts the state of a blinker vehicle to a vector
function BlinkerVehicle_to_vec(veh::Entity{BlinkerState, VehicleDef, Int64})
    Float64[posf(veh).s, vel(veh), laneid(veh), veh.state.blinker]
end

# Converts the BlinkerVehicle vector back to an agent
function vec_to_BlinkerVehicle(arr::AbstractArray, id, roadway::Roadway, model)
    @assert length(arr) == BLINKERVEHICLE_ENTITY_DIM
    s = arr[1] # Distance along the lane
    v = arr[2] # velocity
    g = Int(arr[3]) # Goal (lane id)
    b = Bool(arr[4]) # blinker

    vs = VehicleState(Frenet(roadway, g, s), roadway, v)
    bs = BlinkerState(vs, b, model.goals[g], Noise())
    bv = Entity(bs, VehicleDef(), id)
end

# Converts the state of a blinker vehicle to an expanded state space representation
function BlinkerVehicle_to_expanded_vec(veh::Entity{BlinkerState, VehicleDef, Int64})
    one_hot = zeros(6)
    one_hot[laneid(veh)] = 1
    s = posf(veh.state).s .* one_hot
    v = vel(veh.state) .* one_hot
    v2 = v.^2
    b = veh.state.blinker .* one_hot
    Float64[one_hot..., s..., v..., v2..., b...]
end


text:
using DeclarativePackages
using Base.Test

if VERSION < v"0.5.0"
    readstring = readall
end

test(name, r, value) = test(name, r, x->x==value)
function test(name, r, f::Function)
    a = filter(x->length(x)>7 && x[1:7]!="Cloning", split(r[1], '\n'))
    a = filter(x->length(x)>4 && x[1:4]!="http", a)
    a = join(a)
    r2 = (a,r[2])
    if !f(r2)
        println("Test '$name' failed")
        println("the output was: ---------")
        @show r
        println("------ end of output")
        exit(1)
    end
end

decdir = "$(homedir())/.julia/declarative/"
existinginstallation = exists(decdir)

pathof(a...) = joinpath(dirname(@__FILE__), "..",a...)
jdp = pathof("bin/jdp")
function runjdp(file)
    file = pathof("test",file)
    listpackages = pathof("test/listpackages.jl")

    println("Testing $file")
    tmp = tempname()
    cp(file, tmp)
    ENV["DECLARE"] = tmp
    ENV["DECLARE_VERBOSITY"] = 0
    r = (readstring(`$jdp $listpackages`), readstring(tmp))
    rm(tmp)
    r
end

ENV["DECLARE_INCLUDETEST"] = ""
test("empty", runjdp("DECLARE.empty"), x->isempty(x[1]))
test("METADATA1", runjdp("DECLARE.METADATA1"), x->isempty(x[1]))
test("METADATA2", runjdp("DECLARE.METADATA2"), x->isempty(x[1]))
test("JSON1", runjdp("DECLARE.JSON1"), x->contains(x[2],"JSON 0.3.9"))
test("JSON2", runjdp("DECLARE.JSON2"), x->contains(x[2],"JSON 0.3.7"))
test("HDF51", runjdp("DECLARE.HDF5_1"), x->contains(x[1],"HDF5 0.7.0") && !contains(x[2],"DataFrames"))
test("HDF52", runjdp("DECLARE.HDF5_2"), x->contains(x[1],"HDF5 0.7.0") && contains(x[2],(VERSION < v"0.5.0" ? "SHA 0.2.2" : "staticfloat/SHA.jl.git 0.2.2")))
test("HDF53", runjdp("DECLARE.HDF5_3"), x->contains(x[2],"HDF5.jl.git 0.7.0") && contains(x[2],"rened/HDF5"))
ENV["DECLARE_INCLUDETEST"] = "true"
test("HDF55_withtest", runjdp("DECLARE.HDF5_1"), x->contains(x[1],"HDF5 0.7.0"))

if !existinginstallation
    run(`chmod -R a+w $decdir`)
    run(`rm -rf $decdir`)
end


text:
import Dates
using AutoHashEquals

"""
    MonthIncrementPDOM()

When incrementing by months (or years) then preserve the day of month from
originally requested. Uses the "PDOM" shorthand
"""
struct MonthIncrementPDOM <: MonthIncrementConvention end
adjust(::MonthIncrementPDOM, date::Dates.Date, new_month, new_year, cal_mgr::CalendarManager) = (new_year, new_month, Dates.day(date))
Base.show(io::IO, ::MonthIncrementPDOM) = print(io, "PDOM")

"""
    MonthIncrementPDOMEOM()
    MonthIncrementPDOMEOM(calendars)

When incrementing by months (or years) then preserve the day of month from
originally requested, unless it's the last day of the month then maintain that.
Uses the "PDOMEOM" short hand.

To preserve the last business day of the month, then you can pass calendars as
well.
"""
@auto_hash_equals struct MonthIncrementPDOMEOM <: MonthIncrementConvention
    calendars::Union{Vector{String}, Nothing}

    MonthIncrementPDOMEOM() = new(nothing)
    MonthIncrementPDOMEOM(calendars) = new(calendars)
    MonthIncrementPDOMEOM(calendars::String) = new(split(calendars, "|"))
end

function adjust(mic::MonthIncrementPDOMEOM, date::Dates.Date, new_month, new_year, cal_mgr::CalendarManager)
    ldom_rdate = LDOM(mic.calendars)
    ldom = apply(ldom_rdate, date, cal_mgr)

    if ldom == date
        new_ldom = apply(ldom_rdate, Dates.Date(new_year, new_month), cal_mgr)
        return Dates.yearmonthday(new_ldom)
    else
        return (new_year, new_month, Dates.day(date))
    end
end

Base.show(io::IO, mic::MonthIncrementPDOMEOM) = mic.calendars === nothing ? print(io, "PDOMEOM") : print(io, "PDOMEOM@$(join(mic.calendars, "|"))")

const MONTH_INCREMENT_MAPPINGS = Dict(
    "PDOM" => MonthIncrementPDOM(),
    "PDOMEOM" => MonthIncrementPDOMEOM()
)


text:
const BITBUCKET_TEAM = "UnofficialJuliaMirror"
const BITBUCKET_BOT_USERNAME = "UnofficialJuliaMirrorBot"
const BITBUCKET_BOT_APP_PASSWORD = ENV["BITBUCKET_BOT_APP_PASSWORD"]


text:
# QuEST_jl/src/base/decoherence.jl
#

function mixDamping(qureg       :: QuEST_Types.Qureg,
                    targetQubit :: Integer,
                    prob        :: Qreal)   ::Nothing

    ccall(:mixDamping, Cvoid,
          (QuEST_Types.Qureg, Cint,        Qreal),
          qureg,              targetQubit, prob)

    return nothing
end

function mixDensityMatrix(combineQureg      :: QuEST_Types.Qureg,
                          prob              :: Qreal,
                          otherQureg        :: QuEST_Types.Qureg)   :: Nothing

    ccall(:mixDensityMatrix, Cvoid,
          (QuEST_Types.Qureg, Qreal, QuEST_Types.Qureg),
          combineQureg,       prob,  otherQureg)

    return nothing
end

function mixDephasing(qureg         :: QuEST_Types.Qureg,
                      targetQubit   :: Integer,
                      prob          :: Qreal)   ::Nothing

    ccall(:mixDephasing, Cvoid,
          (QuEST_Types.Qureg, Cint,        Qreal),
          qureg,              targetQubit, prob)

    return nothing
end

function mixDepolarising(qureg          :: QuEST_Types.Qureg,
                         targetQubit    :: Integer,
                         prob           :: Qreal)   ::Nothing

    ccall(:mixDepolarising, Cvoid,
          (QuEST_Types.Qureg, Cint,        Qreal),
          qureg,              targetQubit, prob)

    return nothing
end

function mixKrausMap(qureg          :: QuEST_Types.Qureg,
                     targetQubit    :: Integer,
                     ops            :: Vector{QuEST_Types.ComplexMatrix2}) ::Nothing

    ccall(:mixKrausMap, Cvoid,
          (QuEST_Types.Qureg, Cint,         Ptr{QuEST_Types.ComplexMatrix2}, Cint),
          qureg,              targetQubit,  ops,                             length(ops))

    return nothing
end

function mixMultiQubitKrausMap(qureg        :: QuEST_Types.Qureg,
                               targetQubits :: Vector{QubitIdx},
                               ops          :: Vector{QuEST_Types.ComplexMatrixN}) ::Nothing

    local numQubits ::Cint = length(targetQubits)

    @assert all( opidx -> ops[opidx].numQubits==numQubits ,
                 1:length(ops) )

    ccall(:mixMultiQubitKrausMap, Cvoid,
          (QuEST_Types.Qureg, Ptr{Cint},    Cint,      Ptr{QuEST_Types.ComplexMatrixN}, Cint),
          qureg,              targetQubits, numQubits, ops,                             length(ops))

    return nothing
end

function mixPauli(qureg         :: QuEST_Types.Qureg,
                  targetQubit   :: Integer,
                  probX         :: Qreal,
                  probY         :: Qreal,
                  probZ         :: Qreal)   ::Nothing

    ccall(:mixPauli, Cvoid,
          (QuEST_Types.Qureg, Cint,          Qreal, Qreal, Qreal),
          qureg,              targetQubit,   probX, probY, probZ)

    return nothing
end

function mixTwoQubitDephasing(qureg         :: QuEST_Types.Qureg,
                              qubit1        :: Integer,
                              qubit2        :: Integer,
                              prob          :: Qreal)     ::Nothing

    ccall(:mixTwoQubitDephasing, Cvoid,
          (QuEST_Types.Qureg, Cint,   Cint,   Qreal),
          qureg,              qubit1, qubit2, prob)

    return nothing
end

function mixTwoQubitDepolarising(qureg         :: QuEST_Types.Qureg,
                                 qubit1        :: Integer,
                                 qubit2        :: Integer,
                                 prob          :: Qreal)     ::Nothing

    ccall(:mixTwoQubitDepolarising, Cvoid,
    (QuEST_Types.Qureg, Cint,   Cint,    Qreal),
    qureg,              qubit1, qubit2,  prob)

    return nothing
end

function mixTwoQubitKrausMap(qureg         :: QuEST_Types.Qureg,
                             qubit1        :: Integer,
                             qubit2        :: Integer,
                             ops           :: Vector{QuEST_Types.ComplexMatrix4}) ::Nothing

    ccall(:mixTwoQubitKrausMap, Cvoid,
          (QuEST_Types.Qureg, Cint,   Cint,   Ptr{QuEST_Types.ComplexMatrix4}, Cint),
          qureg,              qubit1, qubit2, ops,                             length(ops))

    return nothing
end

#EOF


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/ModelReduction.jl/blob/master/LICENSE

using ModelReduction
using Base.Test

@testset "ModelReduction.jl" begin
    include("test_guyan_reduction.jl")
    include("test_craig_bampton.jl")
end

@testset "test model reduction using Craig-Bampton method" begin
    include("test_model_reduction_craig_bampton.jl")
end


text:
module Dataclass
export @dataclass

macro dataclass(expr::Expr)
    structName = expr.args[2]
    ismutable = expr.args[1]
    fields = []
    initvars = []
    raw_fields = expr.args[3]
    if raw_fields isa Expr && raw_fields.head == :block
        raw_fields = raw_fields.args
    end

    for field in raw_fields
        field isa LineNumberNode && continue

        if isa(field, Expr) &&
           length(field.args) == 2 &&
           isa(field.args[1], Symbol)
            name = field.args[1]
            type =
                (field.head == :(::)) ? Core.eval(__module__, field.args[2]) :
                (type = field.args[2])
            # save _initvars
            (name != :(_initvars)) ? push!(fields, :($(name)::$(type))) :
            initvars = type
        end

    end

    map(eval, initvars.args)
    fields = Tuple(fields)

    quote
        # If _frozen is true, should we not allow __init__(self) to be created and
        # remove mutable attribute from the struct?

        # Struct definition
        if $(esc(ismutable))
            mutable struct $(esc(structName))
                $(map(esc, fields)...)
            end
        else
            struct $(esc(structName))
                $(map(esc, fields)...)
            end
        end

        # Dataclass Functions
        # __init__
        if $(esc(_init))
            function $(esc(:(__init__)))(
                $(esc(:(self)))::$(esc(structName)),
                $(map(esc, fields)...),
            )
                fNames = fieldnames($(esc(structName)))
                for (field, fName) in zip([$(map(esc, fields)...)], fNames)
                    setfield!(self, fName, field)
                end
            end
        end

        # __repr__
        if $(esc(_repr))
            function $(esc(:(__repr__)))($(esc(:(self)))::$(esc(structName)))
                fNames = fieldnames($(esc(structName)))
                str_repr = string($(esc(structName)))
                field_arr = []
                for fName in fNames
                    push!(field_arr, getfield(self, fName))
                end
                field_str = join(field_arr, ", ")
                return str_repr * "(" * field_str * ")"
            end
        end

        # __eq__
        if $(esc(_eq))
            function $(esc(:(__eq__)))(
                $(esc(:(self)))::$(esc(structName)),
                $(esc(:(other)))::$(esc(structName)),
            )
                return $(esc(:__key))(self) == $(esc(:__key))(other)
            end
        end

        if $(esc(_order))
            # __lt__
            function $(esc(:(__lt__)))(
                $(esc(:(self)))::$(esc(structName)),
                $(esc(:(other)))::$(esc(structName)),
            )::Bool
                return $(esc(:__key))(self) < $(esc(:__key))(other)
            end

            # __le__
            function $(esc(:(__le__)))(
                $(esc(:(self)))::$(esc(structName)),
                $(esc(:(other)))::$(esc(structName)),
            )::Bool
                return $(esc(:__key))(self) <= $(esc(:__key))(other)
            end

            # __gt__
            function $(esc(:(__gt__)))(
                $(esc(:(self)))::$(esc(structName)),
                $(esc(:(other)))::$(esc(structName)),
            )::Bool
                return $(esc(:__key))(self) > $(esc(:__key))(other)
            end

            # __ge__
            function $(esc(:(__ge__)))(
                $(esc(:(self)))::$(esc(structName)),
                $(esc(:(other)))::$(esc(structName)),
            )::Bool
                return $(esc(:__key))(self) >= $(esc(:__key))(other)
            end
        end

        # __unsafe_hash__
        if !$(esc(_unsafe_hash))
            if $(esc(_eq)) && $(esc(ismutable))
                # __hash__
                # Rules:
                # - If a class does not define an __eq__() method it should not define a __hash__() operation either --> OK
                # - If it defines __eq__() but not __hash__(), its instances will not be usable as items in hashable collections --> Not covered here
                # - If a class defines mutable objects and implements an __eq__() method, it should not implement __hash__() --> OK
                function $(esc(:(__hash__)))(
                    $(esc(:(self)))::$(esc(structName)),
                )
                    return hash($(esc(:__key))(self))
                end
            end
        end

        # Additional Functions
        # __key
        function $(esc(:(__key)))($(esc(:(self)))::$(esc(structName)))
            fNames = fieldnames($(esc(structName)))
            fieldValues = []
            for fName in fNames
                push!(fieldValues, getfield(self, fName))
            end
            return Tuple(fieldValues)
        end
    end
end

end #module


text:
@testset  "Header" begin
    table = init_log_table(
        (id=:open_nodes, name="#Open"),
        (id=:closed_nodes, name="#Closed"),
    )
    header = TableLogger.get_header(table)
    lines = split(header, "\n")
    @test length(lines) == 2
    @test occursin("#Open", lines[1])
    @test occursin("#Closed", lines[1])
    @test length(lines[2]) == sum(table.setup.widths)
end


text:
    @testset "Test slice" begin

        x = randn((1,2,3,4))
        y = NDTools.slice(x, 2, 2)
        @test x[:, 2:2, :, :] == y

        x = randn((5,2,3,4))
        y = NDTools.slice(x, 1, 4)
        @test x[4:4, :, :, :] == y

        x = randn((5))
        y = NDTools.slice(x, 1, 5)
        @test x[5:5] == y

    end

    @testset "Test slice indices" begin
        x = randn((1,2,3))
        y = NDTools.slice_indices(axes(x), 1, 1)
        @test y == (1:1, 1:2, 1:3)

        x = randn((20,4,20, 1, 2))
        y = NDTools.slice_indices(axes(x), 2, 3)
        @test y == (1:20, 3:3, 1:20, 1:1, 1:2)
    end

    @testset "Flatten trailing dims" begin
        @test flatten_trailing_dims(ones((1, 1, 2, 2)), Val(1)) == [1.0, 1.0, 1.0, 1.0]
        @test flatten_trailing_dims(ones((1, 1, 2, 2)), Val(2)) == [1.0 1.0 1.0 1.0]
        @test flatten_trailing_dims(ones((1, 1, 2, 2)), Val(3)) == [1.0;;; 1.0;;; 1.0;;; 1.0]
        @test flatten_trailing_dims(ones((1, 1, 2, 2)), Val(4)) == [1.0;;; 1.0;;;; 1.0;;; 1.0]
    end

    @testset "Test expand_dims second" begin
        x = zeros(Int, 2,2,2,2,2,2)
        @test size(expand_dims(x,1,2,3)) == (1,1,1,2,2,2,2,2,2)
        @test size(expand_dims(x,3)) == (2,2,1,2,2,2,2)
        @test size(expand_dims(x,6)) == (2,2,2,2,2,1, 2)
        @test size(expand_dims(x,7)) == (2,2,2,2,2,2,1)
        @test size(expand_dims(x,1)) == (1, 2,2,2,2,2,2)
        x = zeros(Int, 2)
        @test size(expand_dims(x,1)) == (1,2)
        @test size(expand_dims(x,2)) == (2,1)

        @test expand_dims(zeros((2,2)), 1) |> size == (1, 2, 2)

        @test expand_dims(zeros((2,2)), 2) |> size == (2, 1, 2)

        @test expand_dims(zeros((2,2)), 3) |> size == (2, 2, 1)

        @test expand_dims(zeros((2,2)), 1,3,4) |> size == (1, 2, 1, 1, 2)
    end

    @testset "Test expand_dims" begin
        function f(s, N)
            @test expand_dims(randn(s), Val(N + length(s)))|> size == (s..., ones(Int,N)...)
        end
        f((1,2,3), 2)
        f((1,2,3,4,5), 8)
        f((1), 5)
    end

    @testset "Test select_region" begin
        @test select_region_view(ones(3,),new_size=(7,),center=(1,), pad_value=-1) == [-1.0;-1.0;-1.0;1.0;1.0;1.0;-1.0]
        for d=1:5
            for n=1:10
                sz = Tuple(rand(1:5) for q in 1:d)
                a = rand(sz...)
                nsz = Tuple(rand(1:5) for q in 1:d)
                nc = Tuple(rand(1:5) for q in 1:d)
                pad = rand()
                @test all(select_region(a,new_size=nsz, center=nc, pad_value=pad) .== select_region(a,new_size=nsz, center=nc, pad_value=pad))
            end
        end
        a = ones(10,10)
        @test all(select_region(a) .== 1) # simplest version
        nz = (20,20)
        @test all(select_region(a, new_size=nz, pad_value=1) .== 1) # with padding
        @test all(select_region(a, new_size=nz, center=(-100,100), pad_value=10) .== 10) # only pad values
        function f(a,b) a.+=1 end # user-defined function
        @test all(select_region!(2 .*a, a, operator! = f) .== 2) # let the operator add one to destination
    end

    @testset "Test assignment functions" begin
        dst = ones(10,10)
        src = ones(10,10)
        @test all(NDTools.assign_to!(dst,src) .== 1)
        dst = ones(10,10)
        @test all(NDTools.add_to!(dst,src) .== 2)
        dst = ones(10,10)
        @test all(NDTools.sub_to!(dst,src) .== 0)
        src = ones(10,10) .* 2
        dst = ones(10,10) .* 2
        @test all(NDTools.mul_to!(dst,src) .== 4)
        src = ones(10,10) .* 2
        dst = ones(10,10)
        @test all(NDTools.div_to!(dst,src) .== 0.5)
    end