text:
# Functions to automatically remove
# nodes from graphs

export compress_graph_dangling!
export compress_graph_zero_coeff!
export compress_graph_output_cleaning!
export has_identity_lincomb
export has_trivial_nodes
export compress_graph_trivial!
export compress_graph_redundant!
export compress_graph_passthrough!
export compress_graph!

function conditional_println(out_string, verbose)
    return verbose && println(out_string)
end

function delete_crefs!(cref, n)
    II = [] # II will be list of nodes not to remove
    for (i, c) in enumerate(cref)
        if (c[1] != n)
            push!(II, i)
        end
    end
    cref_org = copy(cref)
    empty!(cref)
    return append!(cref, cref_org[II])
end

"""
    compress_graph_dangling!(graph,cref=[];verbose=false)

Removes dangling nodes in the graph. If a node is not used anywhere and it is
not in the output, it can safely be removed.
"""
function compress_graph_dangling!(graph, cref = []; verbose = false)
    # Remove dangling nodes
    modified = true
    while modified
        modified = false
        for (key, parents) in graph.parents
            nof_children = 0
            for (pkey, pparents) in graph.parents
                if (pparents[1] == key || pparents[2] == key)
                    nof_children += 1
                end
            end

            if (nof_children == 0 && !(key in graph.outputs))
                conditional_println("Delete dangling: $key", verbose)
                delete!(graph.parents, key)
                delete!(graph.operations, key)
                if (haskey(graph.coeffs, key))
                    delete!(graph.coeffs, key)
                end
                delete_crefs!(cref, key)
                modified = true
                break
            end
        end
    end
end

"""
    compress_graph_zero_coeff!(graph,cref=[];droptol=0;verbose=false)

Searches for linear combinations with zero coeff value and tries to compress by
redirecting node references.
"""
function compress_graph_zero_coeff!(
    graph,
    cref = [];
    droptol = 0,
    verbose = false,
)
    modified = true
    while (modified)
        modified = false
        for (key, parents) in graph.parents
            big_break = false
            if (graph.operations[key] == :lincomb && !(key in graph.outputs)) # don't remove outputs
                v = graph.coeffs[key]
                for s = 1:2
                    other_s = 3 - s
                    if (abs(v[s]) <= droptol)
                        conditional_println(
                            "Remove node: $key because coeff $s is " *
                            string(v[s]) *
                            " ≈ 0",
                            verbose,
                        )

                        # Redirect all lincombs we can find
                        for (pkey, op) in graph.operations
                            pp = [graph.parents[pkey]...]
                            if (
                                op == :lincomb &&
                                (pp[1] == key || pp[2] == key)
                            )
                                # It can be merged

                                conditional_println(
                                    "Delete connection: $key -> $pkey",
                                    verbose,
                                )
                                # z: parent connected to the node to be removed
                                z = 1
                                if (pp[2] == key)
                                    z = 2
                                end

                                cc = [graph.coeffs[pkey]...]
                                cc[z] = cc[z] * v[other_s]

                                pp[z] = graph.parents[key][other_s]
                                # Overwrite
                                add_lincomb!(
                                    graph,
                                    pkey,
                                    cc[1],
                                    pp[1],
                                    cc[2],
                                    pp[2],
                                )

                                modified = true
                                big_break = true
                            end
                        end
                        if (big_break)
                            break
                        end
                    end
                end
            end
            if (big_break)
                break
            end
        end
    end
end
"""
    compress_graph_output_cleaning!(graph,cref=[];verbose=false)

Checks if the output is computed from the linear combination that can be
compressed.
"""
function compress_graph_output_cleaning!(graph, cref = []; verbose = false)
    ismodified = false
    while ismodified
        ismodified = false
        for (i, n) in enumerate(graph.outputs)
            if (graph.operations[n] == :lincomb)
                for s = 1:2
                    other_s = 3 - s
                    if (
                        graph.coeffs[n][s] == 1.0 &&
                        graph.coeffs[n][other_s] == 0.0
                    )
                        # Redirect since it is just a copy of
                        # another node
                        graph.outputs[i] = graph.parents[n][s]
                        ismodified = true
                    end
                end
            end
        end
    end
end

### Remove trivial nodes:
# 1) :mult with identity parent;
# 2) :ldiv with identity as first (left) parent.

# Check if node is trivial because of an identity on the left (:mult or :ldiv).
function is_trivial_left(graph, node)
    return (
        graph.parents[node][1] == :I &&
        (graph.operations[node] == :mult || graph.operations[node] == :ldiv)
    )
end

# Check if node is trivial because of an identity on the right (:mult).
function is_trivial_right(graph, node)
    return (graph.parents[node][2] == :I && graph.operations[node] == :mult)
end

# Return a parent that can replace a node.
# The empty array Any[] is returned if the node is not trivial.
function find_replacement_node(graph, node)
    if is_trivial_left(graph, node)
        return graph.parents[node][2]
    elseif is_trivial_right(graph, node)
        return graph.parents[node][1]
    else
        return []
    end
end

# Replace node with replacement node.
function replace_node!(graph, node, replacement, cref)
    # Make replacement an output node if current node was.
    if (node in graph.outputs)
        deleteat!(graph.outputs, graph.outputs .== node)
        add_output!(graph, replacement)
    end
    # Update children of node to point to replacement.
    for child in get_children(graph, node)
        graph.parents[child] =
            map(x -> (x == node ? replacement : x), graph.parents[child])
    end
    # Delete node from graph.
    del_node!(graph, node)
    return delete_crefs!(cref, node)
end

"""
    has_identity_lincomb(graph) -> Bool

Checks whether the graph has a node which is a linear combination of identity
matrices.
"""
function has_identity_lincomb(graph)
    has_identity_lincomb = false
    for (key, parents) in graph.parents
        if graph.operations[key] == :lincomb &&
           parents[1] == :I &&
           parents[2] == :I
            has_identity_lincomb = true
            break
        end
    end
    return has_identity_lincomb
end

"""
    has_trivial_node(graph) -> Bool

Checks whether the graph has trivial nodes, that is, multiplications by the
identity or linear systems whose coefficient is the identity matrix.
"""
function has_trivial_nodes(graph)
    has_trivial_nodes = false
    for (key, parents) in graph.parents
        if is_trivial_left(graph, key) || is_trivial_right(graph, key)
            has_trivial_nodes = true
            break
        end
    end
    return has_trivial_nodes
end

"""
    compress_graph_trivial!(graph,cref=[];verbose=false)

Removes from the graph trivial nodes, that is, multiplications by the identity
or linear systems whose coefficient is the identity matrix.
"""
function compress_graph_trivial!(graph, cref = []; verbose = false)
    # cref is not used as this function does not remove lincomb nodes.
    ismodified = true
    while ismodified
        ismodified = false
        for (key, parents) in graph.parents
            replacement = find_replacement_node(graph, key)
            if replacement != []
                conditional_println(
                    "Replace node: $key by $replacement",
                    verbose,
                )
                replace_node!(graph, key, replacement, cref)
                ismodified = true
            end
        end
    end
end

### Remove redundant nodes:
# 1) :mult with same coefficients;
# 2) :ldiv with same terms;
# 3) :lincomb with same terms and coefficients (not necessarily)

# Return list of nodes that can be merged with node.
function find_redundant_nodes(graph, node, compress_lincomb)
    # Find nodes with same operation and same parents.
    operation = graph.operations[node]
    if !compress_lincomb && operation == :lincomb
        return Array{Symbol}([])
    end
    parents = graph.parents[node]
    same_operation = keys(
        filter(
            kv -> kv.first != node && operation == kv.second,
            graph.operations,
        ),
    )
    if operation == :ldiv # Order of parents matter.
        redundant_candidates = keys(
            filter(
                kv -> kv.first in same_operation && kv.second == parents,
                graph.parents,
            ),
        )
    elseif operation == :mult # Order of parents doesn't matter.
        redundant_candidates = keys(
            filter(
                kv ->
                    kv.first in same_operation && (
                        kv.second == parents || kv.second == reverse(parents)
                    ),
                graph.parents,
            ),
        )
    elseif operation == :lincomb
        # For :lincomb nodes, check that coefficients are the same.
        coefficients = graph.coeffs[node]
        # Parents in the same order.
        same_parents_order = keys(
            filter(
                kv -> kv.first in same_operation && parents == kv.second,
                graph.parents,
            ),
        )
        redundant_candidates_order = keys(
            filter(
                kv ->
                    kv.first in same_parents_order &&
                        kv.second == coefficients,
                graph.coeffs,
            ),
        )
        # Parents in reverse order.
        same_parents_reverse = keys(
            filter(
                kv ->
                    kv.first in same_operation &&
                        parents == reverse(kv.second),
                graph.parents,
            ),
        )
        redundant_candidates_reverse = keys(
            filter(
                kv ->
                    kv.first in same_parents_reverse &&
                        kv.second == reverse(coefficients),
                graph.coeffs,
            ),
        )
        redundant_candidates =
            union(redundant_candidates_order, redundant_candidates_reverse)
    end

    # Rendundants nodes must both or neither be output nodes.
    is_output_node = node in graph.outputs
    redundant_nodes = collect(
        filter(
            x ->
                x in redundant_candidates &&
                    !xor(x in graph.outputs, is_output_node),
            redundant_candidates,
        ),
    )

    return redundant_nodes
end

# Merge redundant nodes into node.
function merge_redundant_nodes!(graph, node, redundant_nodes, cref, verbose)
    for key in redundant_nodes
        conditional_println("Merge node: $key with $node", verbose)
        replace_node!(graph, key, node, cref)
    end
end

"""
    compress_graph_redundant!(
        graph,
        cref = [];
        compress_lincomb = true,
        verbose = false,
    )

Removes from the graph redundant nodes, that is, nodes that repeat a computation
that is already present in the graph. Nodes corresponding to a linear
combination are removed only if the coefficients are the same and
`compress_lincomb` is set to `true`.
"""
function compress_graph_redundant!(
    graph,
    cref = [];
    compress_lincomb = true,
    verbose = false,
)
    # cref is not used as this function does not remove lincomb nodes.
    ismodified = true
    while ismodified
        ismodified = false
        for key in get_sorted_keys(graph) # No need to check input nodes.
            if key in get_sorted_keys(graph) # Check that node is not removed
                redundant_nodes =
                    find_redundant_nodes(graph, key, compress_lincomb)
                if !isempty(redundant_nodes)
                    merge_redundant_nodes!(
                        graph,
                        key,
                        redundant_nodes,
                        cref,
                        verbose,
                    )
                    ismodified = true
                end
            end
        end
    end
end

"""
    compress_graph_passthrough!(graph,cref=[];verbose=false);

Identifies lincombs that have coefficients (0 1) or (1 0) which correspond to
identity operations. It redirect appropriately.
"""
function compress_graph_passthrough!(graph, cref = []; verbose = false)
    ismodified = true
    while ismodified
        ismodified = false
        for (key, coeffs) in graph.coeffs
            if coeffs == (0, 1)
                k_one = 2
                k_zero = 1
            elseif coeffs == (1, 0)
                k_one = 1
                k_zero = 2
            else
                continue
            end
            # k_one is to keep
            # k_zero connection not important
            conditional_println("Redirect node: $key to $k_one.", verbose)

            p_passthrough = graph.parents[key][k_one]

            # Find all nodes that use key and redirect them
            for (key2, parents) in graph.parents
                if (parents[1] == key)
                    parents = (p_passthrough, parents[2])
                    ismodified = true
                end
                if (parents[2] == key)
                    parents = (parents[1], p_passthrough)
                    ismodified = true
                end
                graph.parents[key2] = parents
            end
        end
    end
end

"""
    compress_graph!(graph,cref=[];verbose=false)

Searches for nodes which can be removed or reorganized without changing the
function it represents. Corresponding references in the `cref`-vector are
removed.
"""
function compress_graph!(graph, cref = []; verbose = false)
    compress_graph_output_cleaning!(graph, cref, verbose = verbose)
    compress_graph_zero_coeff!(graph, cref, verbose = verbose)
    compress_graph_trivial!(graph, cref, verbose = verbose)
    compress_graph_passthrough!(graph, cref, verbose = verbose)
    compress_graph_dangling!(graph, cref, verbose = verbose)
    return graph
end


text:
function MOI.optimize!(
   model::Optimizer
)
   start_time = time()

   # create GAMS Workspace if not done so far
   if isnothing(model.gamswork)
      if isnothing(model.sysdir) && isnothing(model.workdir)
         model.gamswork = GAMSWorkspace()
      elseif isnothing(model.sysdir)
         model.gamswork = GAMSWorkspace()
         set_working_dir(model.gamswork, model.workdir)
      elseif isnothing(model.workdir)
         model.gamswork = GAMSWorkspace(model.sysdir)
      else
         model.gamswork = GAMSWorkspace(model.sysdir, model.workdir)
      end
      model.sysdir = nothing
      model.workdir = nothing
   end

   # check solver options
   if length(model.solver_options) > 0 && ! haskey(model.gams_options, "solver")
      error("No GAMS solver selected (attribute 'solver') but solver options specified: ",
         model.solver_options)
   end

   # dimensions
   model.n_binary = 0
   model.n_integer = 0
   model.n_semiint = 0
   model.n_semicont = 0
   for var in model.variable_info
      if var.type == VARTYPE_BINARY
         model.n_binary += 1
      end
      if var.type == VARTYPE_INTEGER
         model.n_integer += 1
      end
      if var.type == VARTYPE_SEMICONT
         model.n_semicont += 1
      end
      if var.type == VARTYPE_SEMIINT
         model.n_semiint += 1
      end
   end
   model.m_lin = length(model.linear_le_constraints) +
      length(model.linear_ge_constraints) + length(model.linear_eq_constraints)
   model.m_quad = length(model.quadratic_le_constraints) +
      length(model.quadratic_ge_constraints) + length(model.quadratic_eq_constraints)
   model.m_nonlin = !isnothing(model.nlp_data) ? length(model.nlp_data.constraint_bounds) : 0
   model.m = offset_nonlin(model) + model.m_nonlin

   # choose model type
   is_discrete = model.n_binary + model.n_integer + model.n_semicont + model.n_semiint > 0
   is_discrete |= length(model.sos1_constraints) + length(model.sos2_constraints) > 0
   is_nonlinear = model.m_nonlin > 0 || (!isnothing(model.nlp_data) && model.nlp_data.has_objective)
   is_quadratic = model.m_quad > 0 || isa(model.objective, MOI.ScalarQuadraticFunction{Float64})
   is_complementarity = length(model.complementarity_constraints) > 0
   model.model_type = auto_model_type(model.user_model_type, is_quadratic, is_nonlinear, is_discrete, is_complementarity)
   if model.model_type != model.user_model_type && ! MOI.get(model, MOI.Silent())
      @info "Updated GAMS model type: " * label(model.user_model_type) * " -> " * label(model.model_type)
   end

   # use additional objective variable?
   model.objvar = true
   if typeof(model.objective) == MOI.SingleVariable && model.m > 0
      model.objvar = false
   end
   if model.model_type == GAMS.MODEL_TYPE_MCP || model.model_type == GAMS.MODEL_TYPE_CNS
      model.objvar = false
   end

   # write GMS file
   filename = joinpath(model.gamswork.working_dir, "moi.gms")
   open(filename, "w") do fio
      io = GAMSTranslateStream(fio)
      translate_header(io)
      translate_defsets(io, model)
      translate_defvars(io, model)
      translate_defequs(io, model)
      translate_objective(io, model)
      translate_equations(io, model)
      translate_vardata(io, model)
      translate_solve(io, model, "moi")
   end

   # run GAMS
   job = GAMSJob(model.gamswork, filename, "moi")
   model.sol, stats = run(job, options=model.gams_options, solver_options=model.solver_options)

   # process optimal objective
   if model.model_type != GAMS.MODEL_TYPE_MCP && model.model_type != GAMS.MODEL_TYPE_CNS
      if model.objvar
         objvar_name = "objvar"
      else
         objvar_name = "x$(model.objective.variable.value)"
      end
      if haskey(model.sol.var, objvar_name)
         model.obj = model.sol.var[objvar_name].level[1]
      end
   end

   # process solution statistics
   model.solve_status = stats["solveStat"]
   model.model_status = stats["modelStat"]
   if haskey(stats, "objEst")
      model.obj_est = stats["objEst"]
   end

   if haskey(stats, "resUsd")
      model.solve_time = stats["resUsd"]
   else
      model.solve_time = 0.0
   end
   return
end

function MOI.get(
   model::Optimizer,
   ::MOI.SolveTime
)
   return model.solve_time
end

function MOI.get(
   model::Optimizer,
   ::MOI.TerminationStatus
)
   # not called
   if model.solve_status == SOLVE_STATUS_UNDEFINED && model.model_status == MODEL_STATUS_UNDEFINED
      return MOI.OPTIMIZE_NOT_CALLED
   end

   # (more or less) good outcomes
   if model.solve_status == SOLVE_STATUS_NORMAL
      if model.model_status == MODEL_STATUS_OPTIMAL_GLOBAL
         return MOI.OPTIMAL
      elseif model.model_status == MODEL_STATUS_OPTIMAL_LOCAL
         return MOI.LOCALLY_SOLVED
      elseif model.model_status == MODEL_STATUS_UNBOUNDED
         return MOI.DUAL_INFEASIBLE
      elseif model.model_status == MODEL_STATUS_INFEASIBLE_GLOBAL
         return MOI.INFEASIBLE
      elseif model.model_status == MODEL_STATUS_INFEASIBLE_LOCAL
         return MOI.LOCALLY_INFEASIBLE
      elseif model.model_status == MODEL_STATUS_SOLVED_UNIQUE
         return MOI.OPTIMAL
      elseif model.model_status == MODEL_STATUS_SOLVED
         return MOI.LOCALLY_SOLVED
      elseif model.model_status == MODEL_STATUS_SOLVED_SINGULAR
         return MOI.LOCALLY_SOLVED
      elseif model.model_status == MODEL_STATUS_UNBOUNDED_NO_SOLUTION
         return MOI.DUAL_INFEASIBLE
      elseif model.model_status == MODEL_STATUS_INFEASIBLE_NO_SOLUTION
         return MOI.INFEASIBLE
      elseif model.model_status == MODEL_STATUS_INTEGER
         return MOI.OPTIMAL
      elseif model.model_status == MODEL_STATUS_INTEGER_INFEASIBLE
         return MOI.INFEASIBLE
      elseif model.model_status == MODEL_STATUS_FEASIBLE
         if ! MOI.get(model, MOI.Silent())
            @info "Solver returned feasible solution but failed to prove optimality"
         end
         return MOI.OTHER_ERROR
      else
         error("Unsupported termination: " * string(model.solve_status) * " / " *
            string(model.model_status))
      end

   # bad outcomes
   elseif model.solve_status == SOLVE_STATUS_ITERATION
      return MOI.ITERATION_LIMIT
   elseif model.solve_status == SOLVE_STATUS_RESOURCE
      return MOI.TIME_LIMIT
   elseif model.solve_status == SOLVE_STATUS_SOLVER
      return MOI.OTHER_LIMIT
   elseif model.solve_status == SOLVE_STATUS_EVAL_ERROR
      return MOI.OTHER_ERROR
   elseif model.solve_status == SOLVE_STATUS_CAPABILITY
      return MOI.INVALID_MODEL
   elseif model.solve_status == SOLVE_STATUS_LICENSE
      return MOI.OTHER_ERROR
   elseif model.solve_status == SOLVE_STATUS_SETUP_ERROR
      return MOI.INVALID_MODEL
   elseif model.solve_status == SOLVE_STATUS_SOLVE_ERROR
      return MOI.OTHER_ERROR
   elseif model.solve_status == SOLVE_STATUS_INTERNAL_ERROR
      return MOI.OTHER_ERROR
   elseif model.solve_status == SOLVE_STATUS_SYSTEM_ERROR
      return MOI.OTHER_ERROR
   else
      error("Unsupported termination: " * string(model.solve_status))
   end
end

function MOI.get(
   model::Optimizer,
   ::MOI.RawStatusString
)
   return string(model.solve_status) * " / " * string(model.model_status)
end

function MOI.get(
   model::Optimizer,
   ::MOI.ResultCount
)
   if typeof(model.sol) == GAMSSolution
      return 1
   end
   return 0
end

function MOI.get(
   model::Optimizer,
   attr::MOI.PrimalStatus
)
   if ! (1 <= attr.N <= MOI.get(model, MOI.ResultCount()))
      return MOI.NO_SOLUTION
   end

   if model.solve_status == SOLVE_STATUS_NORMAL
      if model.model_status == MODEL_STATUS_OPTIMAL_GLOBAL ||
         model.model_status == MODEL_STATUS_OPTIMAL_LOCAL ||
         model.model_status == MODEL_STATUS_SOLVED_UNIQUE ||
         model.model_status == MODEL_STATUS_SOLVED ||
         model.model_status == MODEL_STATUS_SOLVED_SINGULAR ||
         model.model_status == MODEL_STATUS_INTEGER
         return MOI.FEASIBLE_POINT
      elseif model.model_status == MODEL_STATUS_INFEASIBLE_GLOBAL ||
         model.model_status == MODEL_STATUS_INFEASIBLE_NO_SOLUTION ||
         model.model_status == MODEL_STATUS_UNBOUNDED ||
         model.model_status == MODEL_STATUS_UNBOUNDED_NO_SOLUTION ||
         model.model_status == MODEL_STATUS_INTEGER_INFEASIBLE
         return MOI.INFEASIBILITY_CERTIFICATE
      elseif model.model_status == MODEL_STATUS_INFEASIBLE_LOCAL
         return MOI.INFEASIBLE_POINT
      end
   end
   return MOI.UNKNOWN_RESULT_STATUS
end

function MOI.get(
   model::Optimizer,
   attr::MOI.DualStatus
)
   if ! (1 <= attr.N <= MOI.get(model, MOI.ResultCount()))
      return MOI.NO_SOLUTION
   end

   if model.solve_status == SOLVE_STATUS_NORMAL
      if model.model_status == MODEL_STATUS_OPTIMAL_GLOBAL ||
         model.model_status == MODEL_STATUS_OPTIMAL_LOCAL ||
         model.model_status == MODEL_STATUS_SOLVED_UNIQUE ||
         model.model_status == MODEL_STATUS_SOLVED ||
         model.model_status == MODEL_STATUS_SOLVED_SINGULAR ||
         model.model_status == MODEL_STATUS_INTEGER
         return MOI.FEASIBLE_POINT
      elseif model.model_status == MODEL_STATUS_INFEASIBLE_GLOBAL ||
         model.model_status == MODEL_STATUS_INFEASIBLE_NO_SOLUTION ||
         model.model_status == MODEL_STATUS_UNBOUNDED ||
         model.model_status == MODEL_STATUS_UNBOUNDED_NO_SOLUTION ||
         model.model_status == MODEL_STATUS_INTEGER_INFEASIBLE
         return MOI.INFEASIBILITY_CERTIFICATE
      elseif model.model_status == MODEL_STATUS_INFEASIBLE_LOCAL
         return MOI.INFEASIBLE_POINT
      end
   end
   return MOI.UNKNOWN_RESULT_STATUS
end

ext: jl


text:
function n_tupleh!(bs::BlueStyle, fst::FST, s::State)
    line_margin = s.line_offset + length(fst) + fst.extra_margin
    lidx = findlast(n -> n.typ === PLACEHOLDER, fst.nodes)
    fidx = findfirst(n -> n.typ === PLACEHOLDER, fst.nodes)
    opener = findfirst(is_opener, fst.nodes) !== nothing
    multiline_arg = findfirst(is_block, fst.nodes) !== nothing
    multiline_arg && (fst.nest_behavior = AlwaysNest)

    if lidx !== nothing && (line_margin > s.opts.margin || must_nest(fst))
        args_range = fidx+1:lidx-1
        args_margin = sum(length.(fst[args_range]))

        nest_to_oneline = if can_nest(fst)
            (fst.indent + s.opts.indent + args_margin <= s.opts.margin) &&
                !contains_comment(fst.nodes[args_range])
        else
            false
        end

        # @info "" nest_to_oneline fst.indent fst.indent + s.opts.indent + args_margin args_margin

        line_offset = s.line_offset
        if opener
            fst[end].indent = fst.indent
        end
        if fst.typ !== CSTParser.TupleH || opener
            fst.indent += s.opts.indent
        end

        for (i, n) in enumerate(fst.nodes)
            if n.typ === NEWLINE
                s.line_offset = fst.indent
            elseif nest_to_oneline && (i == fidx || i == lidx)
                fst[i] = Newline(length = n.len)
                s.line_offset = fst.indent
            elseif !nest_to_oneline && n.typ === PLACEHOLDER
                fst[i] = Newline(length = n.len)
                s.line_offset = fst.indent
            elseif n.typ === TRAILINGCOMMA
                if !nest_to_oneline
                    n.val = ","
                    n.len = 1
                end
                nest!(bs, n, s)
            elseif n.typ === TRAILINGSEMICOLON
                n.val = ";"
                n.len = 1
                nest!(bs, n, s)
            elseif n.typ === INVERSETRAILINGSEMICOLON
                n.val = ""
                n.len = 0
                nest!(bs, n, s)
            elseif opener && (i == 1 || i == length(fst.nodes))
                nest!(bs, n, s)
            else
                diff = fst.indent - fst[i].indent
                add_indent!(n, s, diff)
                n.extra_margin = 1
                nest!(bs, n, s)
            end
        end

        if opener
            s.line_offset = fst[end].indent + 1
        end

    else
        extra_margin = fst.extra_margin
        opener && (extra_margin += 1)
        nest!(bs, fst.nodes, s, fst.indent, extra_margin = extra_margin)
    end
end
@inline n_call!(bs::BlueStyle, fst::FST, s::State) = n_tupleh!(bs, fst, s)
@inline n_curly!(bs::BlueStyle, fst::FST, s::State) = n_tupleh!(bs, fst, s)
@inline n_macrocall!(bs::BlueStyle, fst::FST, s::State) = n_tupleh!(bs, fst, s)
@inline n_ref!(bs::BlueStyle, fst::FST, s::State) = n_tupleh!(bs, fst, s)
@inline n_braces!(bs::BlueStyle, fst::FST, s::State) = n_tupleh!(bs, fst, s)
@inline n_vect!(bs::BlueStyle, fst::FST, s::State) = n_tupleh!(bs, fst, s)
@inline n_parameters!(bs::BlueStyle, fst::FST, s::State) = n_tupleh!(bs, fst, s)
@inline n_invisbrackets!(bs::BlueStyle, fst::FST, s::State) = n_tupleh!(bs, fst, s)
@inline n_bracescat!(bs::BlueStyle, fst::FST, s::State) = n_tupleh!(bs, fst, s)
# @inline n_vcat!(bs::BlueStyle, fst::FST, s::State) = n_tupleh!(bs, fst, s)

function n_conditionalopcall!(bs::BlueStyle, fst::FST, s::State)
    if fst[end].typ === CSTParser.ConditionalOpCall
        conditional_to_if_block!(fst, s)
        nest!(bs, fst, s)
    else
        n_conditionalopcall!(DefaultStyle(bs), fst, s)
    end
    return
end


text:
##################################
# ALGORITMO GENÉTICO ESTACIONARIO
##################################
using Random

include("generar-cromosoma.jl")
include("torneo-binario.jl")
include("../utils/distancias.jl")
include("../utils/funcion-objetivo.jl")
"""
    AuxCrucePorIndices(i,j,g, FunciónCruce)
Función auxiliar para ahorrarse el cruce si i==j
"""
function AuxCrucePorIndices(i,j,g, FunciónCruce)
    if i == j
        return [g[i], g[i]]
    else
        return FunciónCruce(g[i], g[j])
    end
end
"""
    AGE(   evaluaciones_máximas_FE,
            numero_cromosomas_por_generación,
            probabilidad_mutación,
            tamaño_cromosoma,
            función_evaluación,
            función_cruce,
            función_mutación)
    Algoritmo genético generacional, tras las iteraciones convenientes devuelve el cromosoma mejor y su función de valuación
"""
function AGE(   evaluaciones_máximas_FE,
                numero_cromosomas_por_generación,
                # probabilidad_cruce, siempre es 1
                probabilidad_mutación,
                tamaño_cromosoma,
                función_evaluación,
                función_cruce,
                función_mutación)

    # Generamos la primera generación
    generación = [GeneraCromosoma(tamaño_cromosoma, 0.1) for i in 1:numero_cromosomas_por_generación]
    evaluaciones_función_evaluación = map(x -> función_evaluación(x), generación)

    # Variable que acumula el número de evaluaciones de la función de evaluación
    evaluaciones = numero_cromosomas_por_generación
    # Array que contendrá a los elementos seleccionados
    Seleccionados = generación
    while evaluaciones < evaluaciones_máximas_FE

        # Paso 1: Selección (torneo binario)
        índice_seleccionado1 = TorneoBinarioEntreDos(evaluaciones_función_evaluación, numero_cromosomas_por_generación)
        índice_seleccionado2 = TorneoBinarioEntreDos(evaluaciones_función_evaluación, numero_cromosomas_por_generación)

        # Paso 2: Cruce
        Seleccionados =
            AuxCrucePorIndices(
                índice_seleccionado1,
                índice_seleccionado2,
                generación,
                función_cruce
        )
        # Paso 3: Mutación
        for i in 1:2
            if(rand()< probabilidad_mutación)
                Seleccionados[i] = función_mutación(Seleccionados[i])
            end
        end
        # Paso 4: Reemplazo
        # Sustituimos los dos peores cromosomas por los dos nuevos
        índice_peor_1 = argmin(evaluaciones_función_evaluación)
        # para que no vuelva a cogerlo, lo ponemos al valor máximo
        evaluaciones_función_evaluación[índice_peor_1] = 100
        índice_peor_2 = argmin(evaluaciones_función_evaluación)

        # actualizamos los cromosomas
        generación[índice_peor_1] = Seleccionados[1]
        generación[índice_peor_2] = Seleccionados[2]
        # y su valor de evaluación
        evaluaciones_función_evaluación[índice_peor_1] = función_evaluación(generación[índice_peor_1])
        evaluaciones_función_evaluación[índice_peor_2] = función_evaluación(generación[índice_peor_2])

        evaluaciones += 2
    end

    # Devolvemos el mejor cromosoma encontrado
    índice_mejor = argmax(evaluaciones_función_evaluación)
    evaluación_mejor = evaluaciones_función_evaluación[índice_mejor]
    mejor_cromosa = generación[índice_mejor]
    return mejor_cromosa, evaluación_mejor
end

function AGE_LearnerOneNN(data::Matrix{<:Real},labels,
        evaluaciones_máximas_FE,
        numero_cromosomas_por_generación,
        probabilidad_cruce,
        tamaño_cromosoma,
        función_cruce,
        función_mutación)
     # Creamos función objetivo
    a = 0.5
    umbral_tasa_reduccion = 0.1
    función_evaluación = CrearFuncionObjetivo(data, labels, a, umbral_tasa_reduccion )

    # Buscamos pesos de acorde al algoritmo AGG
    w, f_w = AGE(evaluaciones_máximas_FE,
        numero_cromosomas_por_generación,
        probabilidad_cruce,
        tamaño_cromosoma,
        función_evaluación,
        función_cruce,
        función_mutación)
    return WeightedLearnerEuclideanOneNN(w, data, labels), f_w, w

end


text:
@info "Running Git tests"

@testset "Git repositories" begin
    @testset "Does not create Git repo" begin
        t = tpl(; plugins=[!Git])
        with_pkg(t) do pkg
            pkg_dir = joinpath(t.dir, pkg)
            @test !isdir(joinpath(pkg_dir, ".git"))
        end
    end

    @testset "Creates Git repo" begin
        t = tpl(; plugins=[Git()])
        with_pkg(t) do pkg
            pkg_dir = joinpath(t.dir, pkg)
            @test isdir(joinpath(pkg_dir, ".git"))
        end
    end

    @testset "With HTTPS" begin
        t = tpl(; plugins=[Git(; ssh=false)])
        with_pkg(t) do pkg
            LibGit2.with(GitRepo(joinpath(t.dir, pkg))) do repo
                remote = LibGit2.get(GitRemote, repo, "origin")
                @test startswith(LibGit2.url(remote), "https://")
            end
        end
    end

    @testset "With SSH" begin
        t = tpl(; plugins=[Git(; ssh=true)])
        with_pkg(t) do pkg
            LibGit2.with(GitRepo(joinpath(t.dir, pkg))) do repo
                remote = LibGit2.get(GitRemote, repo, "origin")
                @test startswith(LibGit2.url(remote), "git@")
            end
        end
    end

    @testset "Without .jl suffix" begin
        t = tpl(; plugins=[Git(; jl=false)])
        with_pkg(t) do pkg
            LibGit2.with(GitRepo(joinpath(t.dir, pkg))) do repo
                remote = LibGit2.get(GitRemote, repo, "origin")
                @test !occursin(".jl", LibGit2.url(remote))
            end
        end
    end

    @testset "With custom name/email" begin
        t = tpl(; plugins=[Git(; name="me", email="a@b.c")])
        with_pkg(t) do pkg
            LibGit2.with(GitRepo(joinpath(t.dir, pkg))) do repo
                @test LibGit2.getconfig(repo, "user.name", "") == "me"
                @test LibGit2.getconfig(repo, "user.email", "") == "a@b.c"
            end
        end
    end

    @testset "With custom default branch" begin
        t = tpl(; plugins=[Git(; branch="main")])
        with_pkg(t) do pkg
            LibGit2.with(GitRepo(joinpath(t.dir, pkg))) do repo
                @test LibGit2.branch(repo) == "main"
            end
        end
    end

    @testset "Adds version to commit message" begin
        # We're careful to avoid a Pkg.update as it triggers Cassette#130.
        t = tpl(; plugins=[Git(), !Tests])
        mock(PT.version_of => _p -> v"1.2.3") do _i
            with_pkg(t) do pkg
                pkg_dir = joinpath(t.dir, pkg)
                LibGit2.with(GitRepo(pkg_dir)) do repo
                    commit = GitCommit(repo, "HEAD")
                    @test occursin("PkgTemplates version: 1.2.3", LibGit2.message(commit))
                end
            end
        end
    end
end


text:
module MicroCollections

# Use README as the docstring of the module:
@doc let path = joinpath(dirname(@__DIR__), "README.md")
    include_dependency(path)
    replace(read(path, String), r"^```julia"m => "```jldoctest README")
end MicroCollections

export EmptyDict,
    EmptySet,
    EmptyVector,
    OneHotArray,
    OneHotVector,
    SingletonDict,
    SingletonSet,
    SingletonVector,
    UndefArray,
    UndefVector,
    emptyshim,
    singletonshim,
    vec0,
    vec1

using Setfield: @lens, Setfield

include("core.jl")
include("vectors.jl")
include("dicts.jl")
include("sets.jl")
include("arrays.jl")
include("onehot.jl")
include("bangbang.jl")

end # module


text:
module ClusterManagers

export launch, manage, kill, init_worker, connect
import Base: launch, manage, kill, init_worker, connect

if VERSION >= v"0.5.0-dev+4047"
    worker_arg = `--worker $(Base.cluster_cookie())`
else
    worker_arg = `--worker`
end

# PBS doesn't have the same semantics as SGE wrt to file accumulate,
# a different solution will have to be found
include("qsub.jl")
include("scyld.jl")
include("condor.jl")
include("slurm.jl")
include("affinity.jl")

end


text:
using ADCME

if Sys.iswindows()
    global save_tensor = load_op_and_grad("./build/Release/libSaver","save_tensor")
    global get_tensor = load_op_and_grad("./build/Release/libSaver","get_tensor")
    global delete_tensor = load_op_and_grad("./build/Release/libSaver","delete_tensor")
else
    global save_tensor = load_op_and_grad("./build/libSaveTensor","save_tensor")
    global get_tensor = load_op_and_grad("./build/libGetTensor","get_tensor")
    global delete_tensor = load_op_and_grad("./build/libDeleteTensor","delete_tensor")
end

val = constant(rand(10))
t1 = constant("tensor1")
t2 = constant("tensor2")
t3 = constant("tensor3")
u1 = save_tensor(t1,val)
u2 = save_tensor(t2,2*val)
u3 = save_tensor(t3,3*val)

z1 = get_tensor(t1);
z2 = get_tensor(t2);
z3 = get_tensor(t3);

d1 = delete_tensor(t1);
d2 = delete_tensor(t2);
d3 = delete_tensor(t3);
sess = Session();
run(sess, [u1,u2,u3]) # add all the keys

# get the keys one by one
run(sess, z1)
run(sess, z2)
run(sess, z3)

# delete 2nd key
run(sess, d2)


text:
include("test_core.jl")
include("perf.jl")


text:
using VertexStaticGraphs: VStaticGraph

using Test
import Random

using LightGraphs
const LG = LightGraphs
import LightGraphs.SimpleGraphs: SimpleEdge

@testset "VStaticGraph constructor" begin
    g = VStaticGraph{Int,6,false}()
    @test length(g.fadjlist) == 6
    @test g.ne == 0
    @test all(isempty.(g.fadjlist))
    g2 = VStaticGraph{Int}(Val{6}())
    @test length(g2.fadjlist) == 6
    @test g2.ne == 0
    @test all(isempty.(g2.fadjlist))
end

@testset "LG interface" begin
    g = VStaticGraph{Int,6,false}()
    dg = VStaticGraph{Int,6,true}()
    @test ne(g) == ne(dg) == 0
    @test nv(g) == nv(dg) == 6
    @test !is_directed(g)
    @test is_directed(dg)
end

@testset "badjlist and neighbors" begin
    g = VStaticGraph{Int,6,false}()
    dg = VStaticGraph{Int,6,true}()
    @test g.badjlist === nothing
    @test length(dg.badjlist) == 6
    @test LG.SimpleGraphs.badj(g) === g.fadjlist
    @test LG.SimpleGraphs.badj(dg) !== dg.fadjlist
end

@testset "Edge addition and presence" begin
    g = VStaticGraph{Int,6,false}()
    dg = VStaticGraph{Int,6,true}()
    @test add_edge!(g, 1, 2)
    @test add_edge!(dg, 1, 3)
    @test !add_edge!(g, 1, 2)
    @test !add_edge!(dg, 1, 3)
    @test ne(g) == ne(dg) == 1
    @test length(outneighbors(g, 1)) == length(outneighbors(dg, 1)) == 1
    @test inneighbors(dg, 3) == [1]
    @test isempty(outneighbors(dg, 3))
    @test_broken add_edge!(g, LightGraphs.Edge(1,4))
    @test_broken add_edge!(dg, LightGraphs.Edge(1,4))
    # checking vertices bounds
    @test !add_edge!(g, 1, 7)
    @test !add_edge!(dg, 1, 10)
end

@testset "Edge addition and presence" begin
    g = VStaticGraph{Int,6,false}()
    dg = VStaticGraph{Int,6,true}()
    @test add_edge!(g, 1, 2)
    @test add_edge!(dg, 1, 3)
    @test ne(g) == ne(dg) == 1
    @test rem_edge!(g, 1, 2)
    @test rem_edge!(dg, 1, 3)
    @test ne(g) == ne(dg) == 0
    @test !rem_edge!(g, 1, 2)
    @test !rem_edge!(dg, 1, 3)

    # checking vertices bounds
    @test !rem_edge!(g, 1, 7)
    @test !rem_edge!(dg, 1, 10)
end

@testset "Edge iterator" begin
    g = VStaticGraph{Int,6,false}()
    dg = VStaticGraph{Int,6,true}()

    !in(Edge(1, 2), edges(g))
    add_edge!(g, 1, 2)
    add_edge!(dg, 1, 3)
    @test Edge(1, 2) in edges(g)
    @test Edge(1, 3) in edges(dg)
    @test Edge(2, 1) in edges(g)
    @test !in(Edge(3, 1), edges(dg))
end

@testset "SimpleGraphs test suite" begin
    @testset "collection operations" begin
        Random.seed!(33)
        ga = VStaticGraph{Int,10,false}()
        gb = VStaticGraph{Int,10,false}()
        dga = VStaticGraph{Int,10,true}()
        dgb = VStaticGraph{Int,10,true}()
        for _ in 1:20
            s = rand(1:10)
            d = rand(1:10)
            add_edge!(ga, s, d)
            s = rand(1:10)
            d = rand(1:10)
            add_edge!(dga, s, d)
        end
        Random.seed!(33)
        for _ in 1:20
            s = rand(1:10)
            d = rand(1:10)
            add_edge!(gb, s, d)
            s = rand(1:10)
            d = rand(1:10)
            add_edge!(dgb, s, d)
        end
        @test @inferred(edges(ga)) == edges(gb)
        @test @inferred(edges(ga)) == collect(Edge, edges(gb))
        @test edges(ga) != collect(Edge, Base.Iterators.take(edges(gb), 5))
        @test collect(Edge, edges(gb)) == edges(ga)
        @test Set{Edge}(collect(Edge, edges(gb))) == edges(ga)
        @test @inferred(edges(ga)) == Set{Edge}(collect(Edge, edges(gb)))

        @test @inferred(edges(dga)) == edges(dgb)
        @test @inferred(edges(dga)) == collect(SimpleEdge, edges(dgb))
        @test edges(dga) != collect(Edge, Base.Iterators.take(edges(dgb), 5))
        @test collect(SimpleEdge, edges(dgb)) == edges(dga)
        @test Set{Edge}(collect(SimpleEdge, edges(dgb))) == edges(dga)
        @test @inferred(edges(dga)) == Set{SimpleEdge}(collect(SimpleEdge, edges(dgb)))
    end

    ga = VStaticGraph{Int,10,false}()
    add_edge!(ga, 3, 2)
    add_edge!(ga, 3, 10)
    add_edge!(ga, 5, 10)
    add_edge!(ga, 10, 3)

    dga = VStaticGraph{Int,10,true}()
    add_edge!(dga, 3, 2)
    add_edge!(dga, 3, 10)
    add_edge!(dga, 5, 10)
    add_edge!(dga, 10, 3)

    e1 = Edge(3, 10)
    e2 = (3, 10)
    @testset "membership operations" begin
        @test e1 ∈ edges(ga)
        @test e2 ∈ edges(ga)
        @test (3, 9) ∉ edges(ga)

        for u in 1:12, v in 1:12
            b = has_edge(ga, u, v)
            @test b == @inferred (u, v) ∈ edges(ga)
            @test b == @inferred (u => v) ∈ edges(ga)
            @test b == @inferred Edge(u, v) ∈ edges(ga)

            db = has_edge(dga, u, v)
            @test db == @inferred (u, v) ∈ edges(dga)
            @test db == @inferred (u => v) ∈ edges(dga)
            @test db == @inferred Edge(u, v) ∈ edges(dga)
        end
    end

    @testset "iterator protocol" begin
        eit = edges(ga)
        # @inferred not valid for new interface anymore (return type is a Union)
        @test collect(eit) == [Edge(2, 3), Edge(3, 10), Edge(5, 10)]

        eit = @inferred(edges(dga))
        @test collect(eit) == [
            SimpleEdge(3, 2), SimpleEdge(3, 10),
            SimpleEdge(5, 10), SimpleEdge(10, 3)
        ]
    end

end


text:
"""
    CountingTropical{T,CT} <: Number

Counting tropical number type is also a semiring algebra.
It is tropical algebra with one extra field for counting, it is introduced in [arXiv:2008.06888](https://arxiv.org/abs/2008.06888).

Example
-------------------------
```jldoctest; setup=:(using TropicalNumbers)
julia> CountingTropical(1.0, 5.0) + CountingTropical(3.0, 2.0)
(3.0, 2.0)ₜ

julia> CountingTropical(1.0, 5.0) * CountingTropical(3.0, 2.0)
(4.0, 10.0)ₜ

julia> one(CountingTropicalF64)
(0.0, 1.0)ₜ

julia> zero(CountingTropicalF64)
(-Inf, 0.0)ₜ
```
"""
struct CountingTropical{T,CT} <: Number
    n::T
    c::CT
end
CountingTropical(x::T) where T<:Real = CountingTropical(x, one(T))
CountingTropical{T1}(x) where {T1} = CountingTropical{T1,T1}(x)

CountingTropical{T1,CT1}(x) where {T1, CT1} = CountingTropical{T1,CT1}(T1(x), one(CT1))
CountingTropical{T1,CT1}(x::CountingTropical{T1,CT1}) where {T1,CT1} = x
CountingTropical{T1,CT1}(x::CountingTropical) where {T1,CT1} = CountingTropical(T1(x.n), CT1(x.c))
CountingTropical{T1,CT1}(x::Tropical) where {T1,CT1} = CountingTropical{T1,CT1}(T1(TropicalNumbers.content(x)), one(CT1))

Base.:*(a::CountingTropical, b::CountingTropical) = CountingTropical(a.n + b.n, a.c * b.c)
Base.:^(a::CountingTropical, b::Real) = CountingTropical(a.n * b, a.c ^ b)
Base.:^(a::CountingTropical, b::Integer) = CountingTropical(a.n * b, a.c ^ b)
function Base.:+(a::CountingTropical, b::CountingTropical)
    n = max(a.n, b.n)
    if a.n > b.n
        c = a.c
    elseif a.n == b.n
        c =  a.c + b.c
    else
        c = b.c
    end
    CountingTropical(n, c)
end
# inverse and division
Base.inv(x::CountingTropical) = CountingTropical(-x.n, x.c)

Base.zero(::Type{CountingTropical{T}}) where T = zero(CountingTropical{T,T})
Base.zero(::Type{CountingTropical{T,CT}}) where {T<:Integer,CT} = CountingTropical(T(-999999), zero(CT))
Base.zero(::Type{CountingTropical{T,CT}}) where {T<:AbstractFloat,CT} = CountingTropical(typemin(T), zero(CT))
Base.zero(::T) where T<:CountingTropical = zero(T)
Base.one(::Type{CountingTropical{T}}) where T = one(CountingTropical{T,T})
Base.one(::Type{CountingTropical{T,CT}}) where {T<:Integer,CT} = CountingTropical(zero(T), one(CT))
Base.one(::Type{CountingTropical{T,CT}}) where {T<:AbstractFloat,CT} = CountingTropical(zero(T), one(CT))
Base.one(::T) where T<:CountingTropical = one(T)
Base.isapprox(a::CountingTropical, b::CountingTropical; kwargs...) = isapprox(a.n, b.n; kwargs...) && isapprox(a.c, b.c; kwargs...)

Base.show(io::IO, t::CountingTropical) = Base.print(io, "$((t.n, t.c))ₜ")

Base.promote_type(::Type{CountingTropical{T1,CT1}}, b::Type{CountingTropical{T2,CT2}}) where {T1,T2,CT1,CT2} = CountingTropical{promote_type(T1,T2), promote_type(CT1,CT2)}


text:
#This example has been adopted from https://github.com/IntelLabs/ParallelAccelerator.jl/blob/master/examples/black-scholes/black-scholes.jl

using ArrayFire, Printf, Random

function blackscholes_serial(sptprice, strike, rate, volatility, time)
    logterm = log10( sptprice / strike)
    powterm = .5f0 * volatility * volatility
    den = volatility * sqrt(time)
    d1 = (((rate + powterm) * time) + logterm) / den
    d2 = d1 - den
    NofXd1 = cndf2(d1)
    NofXd2 = cndf2(d2)
    futureValue = strike * exp(- rate * time)
    c1 = futureValue * NofXd2
    call_ = sptprice * NofXd1 - c1
    put  = call_ - futureValue + sptprice
end

@inline function cndf2(in)
    out = 0.5f0 + 0.5f0 * erf(0.707106781f0 * in)
    return out
end

function runs(iterations)
    sptprice   = Float32[ 42.0 for i = 1:iterations ]
    initStrike = Float32[ 40.0 + (i / iterations) for i = 1:iterations ]
    rate       = Float32[ 0.5 for i = 1:iterations ]
    volatility = Float32[ 0.2 for i = 1:iterations ]
    time       = Float32[ 0.5 for i = 1:iterations ]

    sptprice_gpu = AFArray(sptprice)
    initStrike_gpu = AFArray(initStrike)
    rate_gpu = AFArray(rate)
    volatility_gpu = AFArray(volatility)
    time_gpu = AFArray(time)

    t1 = @elapsed put1 = blackscholes_serial.(sptprice, initStrike, rate, volatility, time)
    #    println("Serial checksum: ", sum(put1))
    t2 = @elapsed put2 = sync(blackscholes_serial.(sptprice_gpu, initStrike_gpu, rate_gpu, volatility_gpu, time_gpu))
    #    println("Parallel checksum: ", sum(put2))
    @assert sum(put1) ≈ sum(put2)
    return t1, t2
end

function driver(pwr)
    iterations = 10^pwr
    Random.seed!(0)
    blackscholes_serial.(Float32[], Float32[], Float32[], Float32[], Float32[])
    blackscholes_serial.(AFArray(Float32[1., 2.]), AFArray(Float32[1., 2.]),
                         AFArray(Float32[1., 2.]), AFArray(Float32[1., 2.]), AFArray(Float32[1., 2.]))
    tserial, tp1 = runs(iterations)
    tserial, tp2 = runs(iterations)
    tserial, tp3 = runs(iterations)
    tparallel = (tp1 + tp2 + tp3) / 3
    #    println("Time taken for CPU = $tserial")
    #    println("Time taken for GPU = $tparallel")
    println("10^$(pwr) options:")
    @printf("    Speedup = %.2f\n", tserial / tparallel)
    @printf("    CPU rate = 10^%.2f opts/sec\n", log10(iterations / tserial))
    @printf("    GPU rate = 10^%.2f opts/sec\n", log10(iterations / tparallel))
end

for k = 3:8
    driver(k)
end


text:
using BayesianDataFusion
using Base.Test

## testing column dot
X1 = rand(5, 10)
X2 = rand(5, 3)
@test_approx_eq BayesianDataFusion.column_dot(X1, X2, 8, 1) dot(X1[:,8], X2[:,1])

Y  = sprand(15,10, 0.2)
rd = RelationData(Y, class_cut = 0.5)
assignToTest!(rd.relations[1], 2)

# running the data
result = macau_hmc(rd, burnin = 10, psamples = 10, verbose = false)


text:
# Bivariate Std. Normal distribution with correlation ρ
# pdf and cdf, joint cdf for P(X ⩽ x, Y ⩽ y)

# Source: https://github.com/mschauer/GaussianDistributions.jl/blob/master/src/bivariate.jl

# 100 Gauss Legendre points via "FastGaussQuadrature.jl"
const lnodes = [-0.9997137267734413,-0.9984919506395958,-0.9962951347331251,-0.9931249370374434,-0.9889843952429918,-0.983877540706057,-0.9778093584869183,-0.9707857757637063,-0.9628136542558156,-0.9539007829254917,-0.944055870136256,-0.9332885350430795,-0.921609298145334,-0.9090295709825297,-0.895561644970727,-0.8812186793850184,-0.8660146884971647,-0.8499645278795913,-0.8330838798884008,-0.8153892383391763,-0.7968978923903145,-0.7776279096494956,-0.7575981185197073,-0.7368280898020207,-0.7153381175730565,-0.693149199355802,-0.6702830156031411,-0.6467619085141293,-0.6226088602037078,-0.5978474702471789,-0.5725019326213813,-0.5465970120650943,-0.5201580198817632,-0.493210789208191,-0.4657816497733582,-0.43789740217203155,-0.40958529167830166,-0.38087298162462996,-0.3517885263724217,-0.32236034390052926,-0.292617188038472,-0.26258812037150336,-0.23230248184497404,-0.20178986409573646,-0.1710800805386034,-0.1402031372361141,-0.10918920358006115,-0.07806858281343654,-0.046871682421591974,-0.015628984421543188,0.015628984421543188,0.046871682421591974,0.07806858281343654,0.10918920358006115,0.1402031372361141,0.1710800805386034,0.20178986409573646,0.23230248184497404,0.26258812037150336,0.292617188038472,0.32236034390052926,0.3517885263724217,0.38087298162462996,0.40958529167830166,0.43789740217203155,0.4657816497733582,0.493210789208191,0.5201580198817632,0.5465970120650943,0.5725019326213813,0.5978474702471789,0.6226088602037078,0.6467619085141293,0.6702830156031411,0.693149199355802,0.7153381175730565,0.7368280898020207,0.7575981185197073,0.7776279096494956,0.7968978923903145,0.8153892383391763,0.8330838798884008,0.8499645278795913,0.8660146884971647,0.8812186793850184,0.895561644970727,0.9090295709825297,0.921609298145334,0.9332885350430795,0.944055870136256,0.9539007829254917,0.9628136542558156,0.9707857757637063,0.9778093584869183,0.983877540706057,0.9889843952429918,0.9931249370374434,0.9962951347331251,0.9984919506395958,0.9997137267734413,]
const lweights = [0.0007346344905056717,0.001709392653518105,0.0026839253715534818,0.003655961201326376,0.004624450063422119,0.005588428003865517,0.006546948450845323,0.007499073255464713,0.008443871469668972,0.009380419653694457,0.010307802574868971,0.011225114023185977,0.012131457662979496,0.013025947892971542,0.013907710703718773,0.014775884527441305,0.015629621077546,0.016468086176145213,0.01729046056832358,0.018095940722128112,0.018883739613374903,0.019653087494435305,0.02040323264620943,0.021133442112527635,0.021843002416247394,0.02253122025633627,0.02319742318525412,0.02384096026596821,0.024461202707957052,0.025057544481579586,0.025629402910208116,0.026176219239545672,0.02669745918357096,0.02719261344657688,0.027661198220792382,0.02810275565910117,0.028516854322395098,0.028903089601125212,0.029261084110638276,0.029590488059912642,0.029890979593332836,0.030162265105169145,0.030404079526454818,0.030616186583980444,0.03079837903115259,0.03095047885049098,0.03107233742756652,0.031163835696209907,0.031224884254849355,0.03125542345386336,0.03125542345386336,0.031224884254849355,0.031163835696209907,0.03107233742756652,0.03095047885049098,0.03079837903115259,0.030616186583980444,0.030404079526454818,0.030162265105169145,0.029890979593332836,0.029590488059912642,0.029261084110638276,0.028903089601125212,0.028516854322395098,0.02810275565910117,0.027661198220792382,0.02719261344657688,0.02669745918357096,0.026176219239545672,0.025629402910208116,0.025057544481579586,0.024461202707957052,0.02384096026596821,0.02319742318525412,0.02253122025633627,0.021843002416247394,0.021133442112527635,0.02040323264620943,0.019653087494435305,0.018883739613374903,0.018095940722128112,0.01729046056832358,0.016468086176145213,0.015629621077546,0.014775884527441305,0.013907710703718773,0.013025947892971542,0.012131457662979496,0.011225114023185977,0.010307802574868971,0.009380419653694457,0.008443871469668972,0.007499073255464713,0.006546948450845323,0.005588428003865517,0.004624450063422119,0.003655961201326376,0.0026839253715534818,0.001709392653518105,0.0007346344905056717,]

struct BiNormal <: ContinuousMultivariateDistribution
	ρ::Float64
end
import Distributions: rand
rand(P::BiNormal) = let x = randn(); x, P.ρ*x + P.ρ * randn(); end

# Distribution function
Phi(x) = Distributions.cdf(Normal(), x)

# densities
phi(x) = Distributions.pdf(Normal(), x)

# bivariate density
phi(x, y, rho) =  1/(2*pi*sqrt(1-rho^2))*exp(-0.5*(x^2 + y^2 - 2x*y*rho)/(1-rho^2))
# transformed to the interval [-1,1]
phigauss(s, x, y, rho) = (rho)/2 * phi(x, y, 0.5*rho*(s + 1))
# substitute r = sqrt(1-rho^2) for backward integration
function phiback(x, y, r)
	r̄ = sqrt(1-r^2)
	(1/(2pi*r̄))*exp(-0.5*(x^2 + y^2 - 2x*y*r̄)/r^2)
end
# transformed to the interval [-1,1]
phibackgauss(s, x, y, rho) = 0.5*sqrt(1-rho^2)*phiback(x, y, (sqrt(1-rho^2))*0.5*(s + 1))

function Phi(x, y, ρ)
	if ρ == 1
		return Phi(min(x, y))
	end

	if x == Inf || y == Inf
		Phi(min(x, y))
	elseif x == -Inf || y == -Inf
		0.0
	elseif ρ < -0.95
		Phi(x) - Phi(x, -y, -ρ)
	elseif ρ > 0.95
		Phi(min(x,y)) - sum( lweights[i] * phibackgauss(lnodes[i], x, y, ρ) for i in 1:length(lnodes) )
	else
		Phi(x)*Phi(y) + sum( lweights[i] * phigauss(lnodes[i], x, y, ρ) for i in 1:length(lnodes) )
	end
end

import Distributions: pdf,cdf
pdf(P::BiNormal, x::Array{T} where T) = 1/(2*pi*P.ρ)*exp(-0.5*(x[1]^2 + x[2]^2 - 2x[1]*x[2]*P.ρ)/(1-P.ρ^2))
cdf(P::BiNormal, x::Array{T} where T) = Phi(x[1], x[2], P.ρ)


text:
# SPDX-License-Identifier: MIT

struct Gamma{T1, T2} <: ExponentialFamilyDistribution
    α::T1
	β::T2
end

η(p::Gamma) = [-p.β, p.α]
ξ(p::Gamma, η) = vcat(η[2], η[1])
unpack(p::Gamma, μ) = (x=μ[1], lnx=μ[2])

function A(p::Gamma, η)
    D = length(η)
    α, β = η[2], -η[1]
    loggamma(α) - α*log(β)
end

function sample(p::Gamma)
	α, β = ξ(p)
	rand(Distributions.Gamma(α, 1/β))
end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
# Function for running simplified manifold Metropolis adjusted Langevin
# algorithm (SMMALA)
function smmala(model::Model, opts::SmmalaOpts)
  mcmc = Array(Float64, opts.mcmc.nPostBurnin, model.nPars)
  z = Array(Float64, opts.mcmc.nPostBurnin, model.nPars)
  proposed, accepted = 0., 0.

  println("Running burn-in iterations...")

  currentPars = model.randPrior()
  currentLogPosterior = model.logPosterior(currentPars)
  currentGradLogPosterior = model.gradLogPosterior(currentPars)
  currentG = model.tensor(currentPars)
  currentInvG = inv(currentG)
  currentFirstTerm = currentInvG*currentGradLogPosterior
  driftStep = opts.setDriftStep(1, 0., opts.mcmc.n, opts.mcmc.nBurnin, 0.)

  for i = 1:opts.mcmc.n
    proposed += 1

    proposedPars = copy(currentPars)
    parameterMeans = proposedPars+(driftStep/2)*currentFirstTerm
    cholCurrentInvG = chol(driftStep*currentInvG)
    proposedPars = parameterMeans+cholCurrentInvG'*randn(model.nPars)
    proposedLogPosterior = model.logPosterior(proposedPars)
    probNewGivenOld = (-sum(log(diag(cholCurrentInvG)))
      -(0.5*(parameterMeans-proposedPars)'
      *(currentG/driftStep)*(parameterMeans-proposedPars))[1])
    proposedG = model.tensor(proposedPars)
    proposedInvG = inv(proposedG)
    proposedGradLogPosterior = model.gradLogPosterior(proposedPars)
    proposedFirstTerm = proposedG\proposedGradLogPosterior
    parameterMeans = proposedPars+(driftStep/2)*proposedFirstTerm
    probOldGivenNew =
      (-sum(log(diag(chol(driftStep*eye(model.nPars)*proposedInvG))))
      -(0.5*(parameterMeans-currentPars)'
      *(proposedG/driftStep)*(parameterMeans-currentPars))[1])

   ratio =
     proposedLogPosterior+probOldGivenNew-currentLogPosterior-probNewGivenOld

    if ratio > 0 || (ratio > log(rand()))
      accepted += 1

      currentPars = copy(proposedPars)
      currentLogPosterior = copy(proposedLogPosterior)
      currentGradLogPosterior = copy(proposedGradLogPosterior)
      currentG = copy(proposedG)
      currentInvG = copy(proposedInvG)
      currentFirstTerm  = copy(proposedFirstTerm)
    end

    if i > opts.mcmc.nBurnin
      mcmc[i-opts.mcmc.nBurnin, :] = currentPars;
      z[i-opts.mcmc.nBurnin, :] = -currentGradLogPosterior/2;
    end

    if mod(i, opts.mcmc.monitorRate) == 0
      acceptanceRatio = accepted/proposed

    driftStep = opts.setDriftStep(i, acceptanceRatio, opts.mcmc.n,
      opts.mcmc.nBurnin, driftStep)

      println("Iteration $i of $(opts.mcmc.n): ", round(100*acceptanceRatio, 2),
        " % acceptance ratio")

      proposed, accepted = 0., 0.
    end

    if i == opts.mcmc.nBurnin
      println("Burn-in completed...\n\nRunning post burn-in MCMC...");
    end
  end

  return mcmc, z
end


text:
"""
    init_particle!(model, de, p, n_iter)

Initializes values for a particle.

# Arguments

- `model`: model containing a likelihood function with data and priors
- `de`: differential evolution object
- `p`: a particle
- `n_iter`: the number of iterations
"""
function init_particle(de, model, id, n_iter)
    @unpack n_initial = de
    Θ = de.n_initial > 0 ? de.samples[1,:,id] : model.sample_prior()
    p = Particle(;Θ, id)
    N = n_iter + n_initial
    p.accept = fill(false, N)
    de.evaluate_fitness!(de, model, p)
    p.lp = fill(0.0, N)
    return p
end

function initialize_samples(sample_prior)
    s = sample_prior()
    return Array{eltype(s), 3}(undef,0,0,0)
end

function initialize_samples(de::DE, model, n_iter)
    n = n_iter + de.n_initial
    s = model.sample_prior()
    n_parms = length(s)
    n_p = de.Np * de.n_groups
    samples = Array{eltype(de.samples), 3}(undef,n,n_parms,n_p)
    for p in 1:n_p
        for i in 1:de.n_initial
            samples[i,:, p] = model.sample_prior()
        end
    end
    return samples
end

"""
    accept(proposal, current, log_adj=0.0)

Metropolis-Hastings proposal selection
Note: assumes weights are posterior log likelihoods

# Arguments

- `proposal`: weight of proposal e.g. posterior log likelihood
- `current`: weight of current value e.g. posterior log likelihood
- `adj`: an adjustment term for the snooker update
"""
function accept(proposal, current, log_adj=0.0)
    p = min(1.0, exp(proposal - current + log_adj))
    return rand() <= p ? true : false
end

"""
    in_bounds(b, θ::Real)

Checks whether parameter is within lower and upper bounds.

# Arguments

- `b`: boundary (lowerbound,upperbound)
- `θ`: a parameter value
"""
in_bounds(b, θ::Real) = θ >= b[1] && θ <= b[2]
in_bounds(b, θ::Array{<:Real,N}) where {N} = all(x -> in_bounds(b, x), θ)

function in_bounds(de::DE, proposal)
    for (b,θ) in zip(de.bounds, proposal.Θ)
        !in_bounds(b, θ) ? (return false) : nothing
    end
    return true
end

"""
    compute_posterior!(de, model, proposal)

Computes posterior log likelihood of proposal particle. The value -Inf
is returned if the proposal is out of bounds.

# Arguments

- `de`: differential evolution object
- `model`: model containing a likelihood function with data and priors
- `proposal`: the proposed particle
"""
function compute_posterior!(de, model, proposal)
    if in_bounds(de, proposal)
        proposal.weight = model.prior_loglike(proposal.Θ) + model.loglike(proposal.Θ)
    else
        proposal.weight = -Inf
    end
    return nothing
end

"""
    evaluate_fun!(de, model, proposal))

Evaluates the fitness of an arbitrary function called `loglike`. This is used for
point estimation as it does not use a prior distribution.

# Arguments

- `de`: differential evolution object
- `model`: model containing a likelihood function with data and priors
- `proposal`: the proposed particle
"""
function evaluate_fun!(de, model, proposal)
    if in_bounds(de, proposal)
        proposal.weight = model.loglike(proposal.Θ)
    else
        proposal.weight = de.update_particle! == maximize! ? -Inf : Inf
    end
    return nothing
end

"""
    get_names(model, p)

Returns parameters names.

# Arguments
- `p`: a particle
- `model`: model containing a likelihood function with data and priors
"""
function get_names(model, p)
    N = size.(p.Θ)
    n_parms = length.(p.Θ)
    parm_names = fill("", sum(n_parms))
    cnt = 0
    for (k,n) in zip(model.names, N)
        if isempty(n)
            cnt += 1
            parm_names[cnt] = string(k)
        else
            for i in CartesianIndices(n)
                cnt += 1
                parm_names[cnt] = string(k, "[", join([i.I...], ","), "]")
            end
        end
    end
    push!(parm_names, "acceptance", "lp")
    return parm_names
end

"""
    store_samples!(de, groups)

Store samples after burnin period.

# Arguments

- `de`: differential evolution object
- `groups`: groups of particles
"""
function store_samples!(de, groups)
    for group in groups
        for p in group
            add_sample!(de, p)
        end
    end
    return nothing
end

function add_sample!(de, p::Particle{T}) where {T <: Real}
    i = de.iter
    de.samples[i,:,p.id] = p.Θ'
    return nothing
end

function add_sample!(de, p)
    i = de.iter
    de.samples[i,:,p.id] = p.Θ
    return nothing
end

function as_union(p)
    T = find_type(p)
    return Array{T}(p)
end

find_type(p) = Union{unique(typeof.(p))...}

"""
    mh_update!(de, current, proposal, log_adj=0.0)

Update particle based on Metropolis-Hastings rule.

# Arguments

- `de`: differential evolution object
- `current`: current particle
- `proposal`: proposal particle
- `log_adj`: an adjustment term for snooker update
"""
function mh_update!(de, current, proposal, log_adj=0.0)
    accepted = accept(proposal.weight, current.weight, log_adj)
    if accepted
        current.Θ = proposal.Θ
        current.weight = proposal.weight
    end
    current.accept[de.iter] = accepted
    current.lp[de.iter] = current.weight
    return nothing
end

function maximize!(de, current, proposal)
    if proposal.weight > current.weight
        current.Θ = proposal.Θ
        current.weight = proposal.weight
     end
    return nothing
end

function minimize!(de, current, proposal)
    if proposal.weight < current.weight
        current.Θ = proposal.Θ
        current.weight = proposal.weight
     end
    return nothing
end

"""
    project(p1::Particle, p2::Particle)

Returns the projection of particle `p1` onto particle `p2`. This method
works for vectors and nested arrays.

# Arguments

- `p1::Particle`: base particle
- `p2::Particle`: particle projected onto `p1`
"""
function project(p1::Particle, p2::Particle)
    v1,v2 = (0.0,0.0)
    for (Θ1,Θ2) in zip(p1.Θ, p2.Θ)
        v1 += sum(Θ1 .* Θ2)
        v2 += sum(Θ2.^2)
    end
    return p2 * (v1 / v2)
end

norm(p::Particle) = norm(p.Θ)

function best_particle(particles, fun)
    mx = particles[1]
    for p in particles
        if fun(p.weight, mx.weight)
            mx = p
        end
    end
    return mx
end

function get_optimal(de, model, particles)
    fun = de.update_particle!  == maximize! ? (>) : (<)
    mxp = best_particle(particles, fun)
    Θ = NamedTuple{Symbol.(model.names)}(mxp.Θ)
    max_val = mxp.weight
    return Θ,max_val
end

# Type-stable arithmatic operations for Union{Array{T,1},T} types (which return Any otherwise)
import Base: +, - ,*

function +(x::Particle, y::Particle)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] + y.Θ[i]
    end
    return Particle(Θ=z)
end

+(x::Real, y::Particle) = +(y, x)

function +(x::Particle, y::Real)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] .+ y
    end
    return Particle(Θ=z)
end

function +(x::Particle, d::Distribution)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] .+′ draw(d, x.Θ[i])
    end
    return Particle(Θ=z)
end

function draw(d, v::Float64)
    return rand(d)
end

function draw(d, v)
    return rand(d, size(v))
end

function *(x::Particle, y::Particle)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] .*′ y.Θ[i]
    end
    return Particle(Θ=z)
end

*(x::Real, y::Particle) = *(y, x)

function *(x::Particle, y::Real)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] .*′ y
    end
    return Particle(Θ=z)
end

*(x::Array{<:Real,1}, y::Particle) = *(y, x)

function *(x::Particle, y::Array{<:Real,1})
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] .* y[i]
    end
    return Particle(Θ=z)
end

function -(x::Particle, y::Particle)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] - y.Θ[i]
    end
    return Particle(Θ=z)
end

-(x::Real, y::Particle) = -(y, x)

function -(x::Particle, y::Real)
    N = length(x.Θ)
    z = similar(x.Θ)
    for i in 1:N
        z[i] = x.Θ[i] .- y
    end
    return Particle(Θ=z)
end

# arithmetic methods for hanlding discrete parameters
*′(x, y) = x * y
*′(x::Int64, y::Float64) = Int(round(x * y))
*′(x::Float64, y::Int64) = Int(round(x * y))
*′(x::Array{Int64,N}, y::Float64) where {N} = @. Int(round(x * y))
*′(x::Float64, y::Array{Int64,N}) where {N} = @. Int(round(x * y))
+′(x, y) = x + y
+′(x::Int64, y::Float64) = Int(round(x + y))
+′(x::Float64, y::Int64) = Int(round(x + y))
+′(x::Array{Int64,N}, y::Float64) where {N} = @. Int(round(x + y))
+′(x::Float64, y::Array{Int64,N}) where {N} = @. Int(round(x + y))


text:
module PDiscRefFEsTests

using Test
using Gridap.TensorValues
using Gridap.ReferenceFEs

T = VectorValue{2,Float64}
reffe = LagrangianRefFE(T,QUAD,2,space=:P)
@test Conformity(reffe) == L2Conformity()
test_lagrangian_reference_fe(reffe)
@test is_n_cube(get_polytope(reffe))

T = VectorValue{2,Float64}
reffe = LagrangianRefFE(T,HEX,3,space=:P)
@test Conformity(reffe) == L2Conformity()
test_lagrangian_reference_fe(reffe)
@test is_n_cube(get_polytope(reffe))

#using Gridap.Visualization
#writevtk(reffe,"reffe")

end # module


text:
"""
    KeplerianElements

Struct representing the standard classical orbital elements.
"""
struct KeplerianElements <: AbstractElements
    "Semi-Major Axis"
    a::Float64
    "Eccentricity"
    e::Float64
    "Inclination"
    i::Float64
    "Argument of Periapsis"
    ω::Float64
    "Right Ascension of the Ascending Node"
    Ω::Float64
    "True Anomaly"
    ν::Float64
end

# Helper validation functions for keplerian elements
_valid_sma(sma) = sma != 0
_valid_ecc(ecc) = ecc >= 0
_valid_inc(inc) = 0 ≤ inc ≤ π
_valid_aop(aop) = true
_valid_raan(raan) = true
_valid_ta(ta) = true

function _validate(ke::KeplerianElements)
    a = semi_major_axis(ke)
    e = eccentricity(ke)
    i = inclination(ke)
    ω = argument_of_periapsis(ke)
    Ω = right_ascension(ke)
    ν = true_anomaly(ke)
    _valid_sma(a) || throw(DomainError(a, "Invalid semi-major axis"))
    _valid_ecc(e) || throw(DomainError(e, "Invalid eccentricity"))
    _valid_inc(i) || throw(DomainError(i, "Invalid inclination"))
    _valid_aop(ω) || throw(DomainError(ω, "Invalid argument of periapsis"))
    _valid_raan(Ω) || throw(DomainError(raan, "Invalid right ascension"))
    _valid_ta(ν) || throw(DomainError(ta, "Invalid true anomaly"))

    asign = a > 0
    ic = isclosed(ke)

    if (ishyperbolic(ke) && asign) || isparabolic(ke) && !isinf(a) || (ic && !asign)
        throw(DomainError((a=a, e=e), "Invalid eccentricity and semi-major axis pair"))
    end

    return ke
end

"""
    KeplerianElements(sma; ecc=0, inc=0, aop=0, raan=0, ta=0)

Contruct a new [`KeplerianElements`](@ref) struct specifying optional orbit elements.

See Also:
[`semi_major_axis`](@ref),
[`eccentricity`](@ref)
[`inclination`](@ref),
[`argument_of_periapsis`](@ref),
[`right_ascension`](@ref),
[`true_anomaly`](@ref)
"""
function KeplerianElements(sma; ecc=0.0, inc=0.0, aop=0.0, raan=0.0, ta=0.0)
    return KeplerianElements(
        sma,
        ecc,
        inc,
        aop,
        raan,
        ta
    ) |> _validate
end

semi_major_axis(k::KeplerianElements) = k.a
eccentricity(k::KeplerianElements) = k.e
inclination(k::KeplerianElements) = k.i
argument_of_periapsis(k::KeplerianElements) = k.ω
right_ascension(k::KeplerianElements) = k.Ω
true_anomaly(k::KeplerianElements) = k.ν


text:
# An implementation of the baconian cipher in Julia

function baconian()

    # The key for each letter
    keys = ["aaaaa", "aaaab", "aaaba", "aaabb", "aabaa", "aabab", "aabba", "aabbb", "abaaa","abaaa", "abaab", "ababa", "ababb", "abbaa", "abbab", "abbba", "abbbb", "baaaa", "baaab", "baaba", "baabb","baabb", "babaa", "babab", "babba", "babbb"]
#=  A = "aaaaa"    I = "abaaa"    R = "baaaa"
    B = "aaaab"    K = "abaab"    S = "baaab"
    C = "aaaba"    L = "ababa"    T = "baaba"
    D = "aaabb"    M = "ababb"    U = "baabb"
    E = "aabaa"    N = "abbaa"    W = "babaa"
    F = "aabab"    O = "abbab"    X = "babab"
    G = "aabba"    P = "abbba"    Y = "babba"
    H = "aabbb"    Q = "abbbb"    Z = "babbb" =#

    # Take input

    println("Please input the secret message")
    secret = readline()
    secret = uppercase(secret)
    secret = replace(secret, "V" => "U")
    secret = replace(secret, "J" => "I")
    secret = replace(secret, " " => "")

    println("Please input the normal text (5 times longer than the secret)")
    text = readline()
    text = replace(text, " " => "")
    while length(text) != 5 * length(secret)
        println("Not 5 times the length of the secret. Try again")
        text = readline()
    end
    text = lowercase(text)

    cipherarray = []

    # Encrypt
    for i in 0:(length(secret)-1)
        letter = secret[i+1]
        code = keys['Z' - letter + 1]
        for j in 1:length(code)
            if code[j] == 'b'
                push!(cipherarray, uppercase(text[5i+j]))
            else
                push!(cipherarray, lowercase(text[5i+j]))
            end
        end
    end
    # Covert array to String
    ciphertext = join(cipherarray)
    println(ciphertext)
end

baconian()


text:
using NCDatasets
using DataStructures

using Random
rng = MersenneTwister(1234)

# Create mock data in a temporary directory
mockdir = mktempdir(homedir())
mockdata_path = joinpath(mockdir, "GEOTRACES_mockdata.nc")
ENV["GEOTRACES_IDP21_PATH"] = mockdata_path

# The following was generated using something like below
# (Giving up on making this reproducible because of GEOTRACES own agreements/rules.)
#=
]dev GEOTRACES
]add NCDatasets
using GEOTRACES, NCDatasets
mockdir = mktempdir(homedir())
ncfile = GEOTRACES.GEOTRACES_IDP21_DiscreteSamples_path()
ncgen(ncfile, joinpath(mockdir, "mockdata.jl"))
=#

ds = NCDataset(mockdata_path, "c", attrib = OrderedDict(
    "Conventions"               => "CF-1.7",
    "comment"                   => "ODV NetCDF Export File V2.0",
    "Creator"                   => "rschlitz@BGEO04M097-1",
    "CreateTime"                => "2021-11-11T08:45:54",
    "Software"                  => "Ocean Data View 5.5.2 - 64 bit (Windows)",
    "source"                    => "C:/GEOTRACES/IDP2021/output/data/seawater/GEOTRACES_IDP2021_Seawater_Discrete_Sample_Data_v1.odv",
    "SourceLastModified"        => "2021-11-09T15:31:14",
    "DataField"                 => "Ocean",
    "DataType"                  => "Profiles",
    "Description"               => "GEOTRACES IDP2021 seawater discrete sample hydrographic and tracer data",
    "featureType"               => "profile",
    "AdditionnalComment"        => "This is a mock dataset for CI of GEOTRACES.jl"
))

# Dimensions

ds.dim["N_STATIONS"] = 3149
ds.dim["N_SAMPLES"] = 698
ds.dim["STRING26"] = 26
ds.dim["STRING6"] = 6
ds.dim["STRING20"] = 20
ds.dim["STRING23"] = 23
ds.dim["STRING13"] = 13
ds.dim["STRING25"] = 25
ds.dim["STRING75"] = 75
ds.dim["STRING31"] = 31
ds.dim["STRING77"] = 77
ds.dim["STRING7"] = 7
ds.dim["STRING4"] = 4
ds.dim["STRING12"] = 12
ds.dim["STRING22"] = 22
ds.dim["STRING56"] = 56
ds.dim["STRING36"] = 36
ds.dim["STRING10"] = 10

# Declare variables
# _FillValue modified by BP from "" to '\0'

ncmetavar1 = defVar(ds,"metavar1", Char, ("STRING6", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Cruise",
    "units"                     => "",
    "comment"                   => "",
))

ncmetavar2 = defVar(ds,"metavar2", Char, ("STRING26", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Station",
    "units"                     => "",
    "comment"                   => "",
))

nclongitude = defVar(ds,"longitude", Float32, ("N_STATIONS",), attrib = OrderedDict(
    "long_name"                 => "Longitude",
    "standard_name"             => "longitude",
    "units"                     => "degrees_east",
    "comment"                   => "",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

nclatitude = defVar(ds,"latitude", Float32, ("N_STATIONS",), attrib = OrderedDict(
    "long_name"                 => "Latitude",
    "standard_name"             => "latitude",
    "units"                     => "degrees_north",
    "comment"                   => "",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncdate_time = defVar(ds,"date_time", Float64, ("N_STATIONS",), attrib = OrderedDict(
    "long_name"                 => "Decimal Gregorian Days of the station",
    "standard_name"             => "time",
    "units"                     => "days since 2006-01-01 00:00:00 UTC",
    "comment"                   => "Relative Gregorian Days with decimal part",
    "C_format"                  => "%.5f",
    "FORTRAN_format"            => "F12.5",
    "_FillValue"                => -1.0e10,
))

ncvar1 = defVar(ds,"var1", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "CTDPRS_T_VALUE_SENSOR",
    "units"                     => "dbar",
    "comment"                   => "Pressure from CTD sensor",
    "ancillary_variables"       => "var1_qc",
    "C_format"                  => "%.0f",
    "FORTRAN_format"            => "F12.0",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar1_qc = defVar(ds,"var1_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of CTDPRS_T_VALUE_SENSOR",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar2 = defVar(ds,"var2", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "DEPTH",
    "units"                     => "m",
    "comment"                   => "Depth below sea surface calculated from pressure",
    "ancillary_variables"       => "var2_qc",
    "C_format"                  => "%.0f",
    "FORTRAN_format"            => "F12.0",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar2_qc = defVar(ds,"var2_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of DEPTH",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar15 = defVar(ds,"var15", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "CTDTMP_T_VALUE_SENSOR",
    "units"                     => "deg C",
    "comment"                   => "Temperature from CTD sensor in the ITS-90 convention. Metadata must include make and model numbers and recent calibration information.",
    "ancillary_variables"       => "var15_qc var15_err",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar15_qc = defVar(ds,"var15_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of CTDTMP_T_VALUE_SENSOR",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar16 = defVar(ds,"var16", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "CTDSAL_D_CONC_SENSOR",
    "units"                     => "pss-78",
    "comment"                   => "Practical salinity from CTD sensor on the PSS-1978 scale. Metadata must include make and model numbers and recent calibration information.",
    "ancillary_variables"       => "var16_qc var16_err",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar16_qc = defVar(ds,"var16_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of CTDSAL_D_CONC_SENSOR",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar22 = defVar(ds,"var22", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "He_D_CONC_BOTTLE",
    "units"                     => "nmol/kg",
    "comment"                   => "Concentration of dissolved Helium",
    "ancillary_variables"       => "var22_qc var22_err",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar22_qc = defVar(ds,"var22_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of He_D_CONC_BOTTLE",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar22_err = defVar(ds,"var22_err", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Error of He_D_CONC_BOTTLE",
    "units"                     => "nmol/kg",
    "comment"                   => "Concentration of dissolved Helium",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar35 = defVar(ds,"var35", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "OXYGEN_D_CONC_BOTTLE",
    "units"                     => "umol/kg",
    "comment"                   => "Concentration of dissolved oxygen from a bottle sample",
    "ancillary_variables"       => "var35_qc",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar35_qc = defVar(ds,"var35_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of OXYGEN_D_CONC_BOTTLE",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar37 = defVar(ds,"var37", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "PHOSPHATE_D_CONC_BOTTLE",
    "units"                     => "umol/kg",
    "comment"                   => "Concentration of dissolved phosphate, samples may or may not have been filtered",
    "ancillary_variables"       => "var37_qc var37_err",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar37_qc = defVar(ds,"var37_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of PHOSPHATE_D_CONC_BOTTLE",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar37_err = defVar(ds,"var37_err", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Error of PHOSPHATE_D_CONC_BOTTLE",
    "units"                     => "umol/kg",
    "comment"                   => "Concentration of dissolved phosphate, samples may or may not have been filtered",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar39 = defVar(ds,"var39", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "SILICATE_D_CONC_BOTTLE",
    "units"                     => "umol/kg",
    "comment"                   => "Concentration of dissolved silicate (silicic acid), samples may or may not have been filtered",
    "ancillary_variables"       => "var39_qc var39_err",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar39_qc = defVar(ds,"var39_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of SILICATE_D_CONC_BOTTLE",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar39_err = defVar(ds,"var39_err", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Error of SILICATE_D_CONC_BOTTLE",
    "units"                     => "umol/kg",
    "comment"                   => "Concentration of dissolved silicate (silicic acid), samples may or may not have been filtered",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar40 = defVar(ds,"var40", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "NITRATE_D_CONC_BOTTLE",
    "units"                     => "umol/kg",
    "comment"                   => "Concentration of dissolved nitrate, samples may or may not have been filtered",
    "ancillary_variables"       => "var40_qc var40_err",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar40_qc = defVar(ds,"var40_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of NITRATE_D_CONC_BOTTLE",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar40_err = defVar(ds,"var40_err", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Error of NITRATE_D_CONC_BOTTLE",
    "units"                     => "umol/kg",
    "comment"                   => "Concentration of dissolved nitrate, samples may or may not have been filtered",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar85 = defVar(ds,"var85", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Cd_D_CONC_BOTTLE",
    "units"                     => "nmol/kg",
    "comment"                   => "Concentration of dissolved Cd",
    "ancillary_variables"       => "var85_qc var85_err",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar85_qc = defVar(ds,"var85_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of Cd_D_CONC_BOTTLE",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar85_err = defVar(ds,"var85_err", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Error of Cd_D_CONC_BOTTLE",
    "units"                     => "nmol/kg",
    "comment"                   => "Concentration of dissolved Cd",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar88 = defVar(ds,"var88", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Fe_D_CONC_BOTTLE",
    "units"                     => "nmol/kg",
    "comment"                   => "Concentration of dissolved Fe",
    "ancillary_variables"       => "var88_qc var88_err",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar88_qc = defVar(ds,"var88_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of Fe_D_CONC_BOTTLE",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar88_err = defVar(ds,"var88_err", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Error of Fe_D_CONC_BOTTLE",
    "units"                     => "nmol/kg",
    "comment"                   => "Concentration of dissolved Fe",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar104 = defVar(ds,"var104", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Ni_D_CONC_BOTTLE",
    "units"                     => "nmol/kg",
    "comment"                   => "Concentration of dissolved Ni",
    "ancillary_variables"       => "var104_qc var104_err",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar104_qc = defVar(ds,"var104_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of Ni_D_CONC_BOTTLE",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar104_err = defVar(ds,"var104_err", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Error of Ni_D_CONC_BOTTLE",
    "units"                     => "nmol/kg",
    "comment"                   => "Concentration of dissolved Ni",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar158 = defVar(ds,"var158", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Cd_114_110_D_DELTA_BOTTLE",
    "units"                     => "per 10^3",
    "comment"                   => "Atom ratio of dissolved Cd isotopes expressed in conventional DELTA notation referenced to {NIST3108}",
    "ancillary_variables"       => "var158_qc var158_err",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar158_qc = defVar(ds,"var158_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of Cd_114_110_D_DELTA_BOTTLE",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar158_err = defVar(ds,"var158_err", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Error of Cd_114_110_D_DELTA_BOTTLE",
    "units"                     => "per 10^3",
    "comment"                   => "Atom ratio of dissolved Cd isotopes expressed in conventional DELTA notation referenced to {NIST3108}",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar160 = defVar(ds,"var160", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Fe_56_54_D_DELTA_BOTTLE",
    "units"                     => "per 10^3",
    "comment"                   => "Atom ratio of dissolved Fe isotopes expressed in conventional DELTA notation referenced to {IRMM-14}",
    "ancillary_variables"       => "var160_qc var160_err",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar160_qc = defVar(ds,"var160_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of Fe_56_54_D_DELTA_BOTTLE",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

ncvar160_err = defVar(ds,"var160_err", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Error of Fe_56_54_D_DELTA_BOTTLE",
    "units"                     => "per 10^3",
    "comment"                   => "Atom ratio of dissolved Fe isotopes expressed in conventional DELTA notation referenced to {IRMM-14}",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar109 = defVar(ds,"var109", Float32, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "V_D_CONC_BOTTLE",
    "units"                     => "nmol/kg",
    "comment"                   => "Concentration of dissolved V",
    "ancillary_variables"       => "var109_qc",
    "C_format"                  => "%.3f",
    "FORTRAN_format"            => "F12.3",
    "_FillValue"                => Float32(-1.0e10),
))

ncvar109_qc = defVar(ds,"var109_qc", Int8, ("N_SAMPLES", "N_STATIONS"), attrib = OrderedDict(
    "long_name"                 => "Quality flag of V_D_CONC_BOTTLE",
    "standard_name"             => "status_flag",
    "comment"                   => "SEADATANET - SeaDataNet quality codes",
    "flag_values"               => Int8[48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 65, 66, 81],
    "flag_meanings"             => "no_quality_control good_value probably_good_value probably_bad_value bad_value changed_value value_below_detection value_in_excess interpolated_value missing_value value_phenomenon_uncertain nominal_value value_below_limit_of_quantification",
    "_FillValue"                => Int8(57),
))

# Define variables

# metavars \
# Recreate cruises
ranges = [1:79, 80:298, 299:403, 404:471, 472:504, 505:740, 741:798, 799:991, 992:1014, 1015:1053, 1054:1080, 1081:1087, 1088:1089, 1090:1215, 1216:1228, 1229:1241, 1242:1251, 1252:1266, 1267:1281, 1282:1299, 1300:1320, 1321:1370, 1371:1394, 1395:1473, 1474:1636, 1637:1723, 1724:1848, 1849:1862, 1863:1908, 1909:1937, 1938:1984, 1985:2051, 2052:2082, 2083:2092, 2093:2187, 2188:2360, 2361:2403, 2404:2436, 2437:2520, 2521:2643, 2644:2702, 2703:2768, 2769:2786, 2787:2802, 2803:2809, 2810:2820, 2821:2835, 2836:2849, 2850:2863, 2864:2890, 2891:2906, 2907:2921, 2922:2927, 2928:2942, 2943:2985, 2986:2986, 2987:3013, 3014:3115, 3116:3137, 3138:3149]
cruises = ["GA01", "GA02", "GA03", "GA04N", "GA04S", "GA06", "GA08", "GA10", "GA11", "GA13", "GAc01", "GAc02", "GAc03", "GApr08", "GApr09", "GI01", "GI02", "GI03", "GI04", "GI05", "GI06", "GIPY01", "GIPY02", "GIPY04", "GIPY05", "GIPY06", "GIPY11", "GIPY13", "GIpr01", "GIpr05", "GIpr06", "GN01", "GN02", "GN03", "GN04", "GN05", "GP02", "GP06", "GP12", "GP13", "GP15", "GP16", "GP18", "GP19", "GPc01", "GPc02", "GPc03", "GPc05", "GPc06", "GPpr01", "GPpr02", "GPpr04", "GPpr05", "GPpr07", "GPpr08", "GPpr10", "GPpr11", "GS01", "GSc01", "GSc02"]
for (c, r) in zip(cruises, ranges), i in 1:6
    ncmetavar1[i,r] = string(c, repeat('\0', 6-length(c)))[i]
end
# Fake station names going from 1 to 3149
for i in 1:3149
    s = string(i)
    for j in 1:20
        ncmetavar2[j,i] = j≤length(s) ? s[j] : '\0'
    end
end
# number of (depth) indices to mock-fill
idx = rand(reduce(vcat, repeat(0:n, Int(ceil(round(100/n^2)))) for n in 1:698), 3149)
function myfill(v)
    out = fill(fillvalue(v.var), size(v)...)
    for i in eachindex(idx)
        out[1:idx[i], i] .= rand(rng, eltype(v.var), idx[i])
    end
    out
end
function myfill_qc(v)
    out = Array{Int8,2}(undef, size(v)...)
    out .= '9' |> Int8
    for i in eachindex(idx)
        out[1:idx[i], i] .= rand(rng, QC, idx[i])
    end
    out
end
# random lat lon date
nclongitude[:] = rand(eltype(nclongitude.var), size(nclongitude)...)
nclatitude[:] = rand(eltype(nclatitude.var), size(nclatitude)...)
ncdate_time[:] = rand(eltype(ncdate_time.var), size(ncdate_time)...)
# QC characters 1:5
QC = [s[1] |> Int8 for s in string.([1,2,3,4,5])]
# pressure
ncvar1[:] = rand(eltype(ncvar1.var), size(ncvar1)...)
ncvar1_qc[:] = rand(QC, size(ncvar1_qc)...)
# depth
ncvar2[:] = rand(eltype(ncvar2.var), size(ncvar2)...)
ncvar2_qc[:] = rand(QC, size(ncvar2_qc)...)
# temp
ncvar15[:] = rand(eltype(ncvar15.var), size(ncvar15)...)
ncvar15_qc[:] = rand(QC, size(ncvar15_qc)...)
# salinity
ncvar16[:] = rand(eltype(ncvar16.var), size(ncvar16)...)
ncvar16_qc[:] = rand(QC, size(ncvar16_qc)...)
# some tracers
for v in [ncvar22, ncvar35, ncvar37, ncvar39, ncvar40, ncvar85, ncvar88, ncvar104, ncvar109, ncvar158, ncvar160,
          ncvar22_err, ncvar37_err, ncvar39_err, ncvar40_err, ncvar85_err, ncvar88_err, ncvar104_err, ncvar158_err, ncvar160_err]
    v[:] = myfill(v)
end
for v in [ncvar22_qc, ncvar35_qc, ncvar37_qc, ncvar39_qc, ncvar40_qc, ncvar85_qc, ncvar88_qc, ncvar104_qc, ncvar109_qc, ncvar158_qc, ncvar160_qc]
    v[:] = myfill_qc(v)
end

close(ds)