text:
using MacroTools

 # splitarg produces a tuple for each arg of the form (arg_name, arg_type, slurp, default)
_arg_name(arg_tup) = arg_tup[1]
_arg_type(arg_tup) = arg_tup[2]
_arg_slurp(arg_tup) = arg_tup[3]
_arg_default(arg_tup) = arg_tup[4]

const NumericArray = Array{T, N} where {T <: Number, N}

function _collect_bindings(exprs)
    bindings = []
    # @info "_collect_bindings: $exprs"

    for expr in exprs
        if @capture(expr, name_ => val_) && name isa Symbol &&
            (val isa Symbol || val isa Number || val.head in (:vcat, :hcat, :vect))
            push!(bindings, name => val)
        else
            error("Elements of bindings list must Pair{Symbol, Symbol} or Pair{Symbol, Number or Array of Number} got $expr")
        end
    end

    # @info "returning $bindings"
    return bindings
end

function _subcomp(calling_module, args, kwargs)
    # splitarg produces a tuple for each arg of the form (arg_name, arg_type, slurp, default)
    arg_tups = map(splitarg, args)

    if kwargs === nothing
        # If a ";" was not used to separate kwargs, move any kwargs from args.
        kwarg_tups = filter(tup -> _arg_default(tup) !== nothing, arg_tups)
        arg_tups   = filter(tup -> _arg_default(tup) === nothing, arg_tups)
    else
        kwarg_tups = map(splitarg, kwargs)
    end

    if 1 > length(arg_tups) > 2
        @error "component() must have one or two non-keyword values"
    end

    arg1 = _arg_name(arg_tups[1])
    alias = length(arg_tups) == 2 ? _arg_name(args_tups[2]) : nothing

    cmodule = nothing
    if ! (@capture(arg1, cmodule_.cname_) || @capture(arg1, cname_Symbol))
        error("Component name must be a Module.name expression or a symbol, got $arg1")
    end

    valid_kws = (:bindings,)    # valid keyword args to the component() psuedo-function
    kw = Dict([key => [] for key in valid_kws])

    for (arg_name, arg_type, slurp, default) in kwarg_tups
        if arg_name in valid_kws
            if default isa Expr && hasmethod(Base.iterate, (typeof(default.args),))
                append!(kw[arg_name], default.args)
            else
                @error "Value of $arg_name argument must be iterable"
            end
        else
            @error "Unknown keyword $arg_name; valid keywords are $valid_kws"
        end
    end

    bindings = _collect_bindings(kw[:bindings])
    module_obj = (cmodule === nothing ? calling_module : getfield(calling_module, cmodule))
    return SubComponent(module_obj, cname, alias, bindings)
end

# Convert an expr like `a.b.c.d` to `[:a, :b, :c, :d]`
function parse_dotted_symbols(expr)
    global Args = expr
    syms = Symbol[]

    ex = expr
    while @capture(ex, left_.right_) && right isa Symbol
        push!(syms, right)
        ex = left
    end

    if ex isa Symbol
        push!(syms, ex)
    else
        # @warn "Expected Symbol or Symbol.Symbol..., got $expr"
        return nothing
    end

    syms = reverse(syms)
    var_or_par = pop!(syms)
    return ComponentPath(syms), var_or_par
end

"""
    defcomposite(cc_name::Symbol, ex::Expr)

Define a Mimi CompositeComponent `cc_name` with the expressions in `ex`.  Expressions
are all variations on `component(...)`, which adds a component to the composite. The
calling signature for `component()` processed herein is:

    component(comp_name, local_name;
              bindings=[list Pair{Symbol, Symbol or Number or Array of Numbers}])

Bindings are expressed as a vector of `Pair` objects, where the first element of the
pair is the name (again, without the `:` prefix) representing a parameter in the component
being added, and the second element is either a numeric constant, a matrix of the
appropriate shape, or the name of a variable in another component. The variable name
is expressed as the component id (which may be prefixed by a module, e.g., `Mimi.adder`)
followed by a `.` and the variable name in that component. So the form is either
`modname.compname.varname` or `compname.varname`, which must be known in the current module.

Unlike leaf components, composite components do not have user-defined `init` or `run_timestep`
functions; these are defined internally to iterate over constituent components and call the
associated method on each.
"""
macro defcomposite(cc_name, ex)
    # @info "defining composite $cc_name in module $(fullname(__module__))"

    @capture(ex, elements__)
    comps = SubComponent[]
    imports = []
    conns = []

    calling_module = __module__
    # @info "defcomposite calling module: $calling_module"

    for elt in elements
        # @info "parsing $elt"; dump(elt)

        if @capture(elt, (component(args__; kwargs__) | component(args__)))
            push!(comps, _subcomp(calling_module, args, kwargs))

        # distinguish imports, e.g., :(EXP_VAR = CHILD_COMP1.COMP2.VAR3),
        #    from connections, e.g., :(COMP1.PAR2 = COMP2.COMP5.VAR2)

        # elseif elt.head == :tuple && length(elt.args) > 0 && @capture(elt.args[1], left_ = right_) && left isa Symbol
        #     # Aliasing a local name to several parameters at once is possible using an expr like
        #     # :(EXP_PAR1 = CHILD_COMP1.PAR2, CHILD_COMP2.PAR2, CHILD_COMP3.PAR5, CHILD_COMP3.PAR6)
        #     # Note that this parses as a tuple expression with first element being `EXP_PAR1 = CHILD_COMP1`.
        #     # Here we parse everything on the right side, at once using broadcasting and add the initial
        #     # component (immediately after "=") to the list, and then store a Vector of param refs.
        #     args = [right, elt.args[2:end]...]
        #     vars_pars = parse_dotted_symbols.(args)
        #     @info "import as $left = $vars_pars"
        #     push!(imports, (left, vars_pars))

        elseif @capture(elt, left_ = right_)

            if left isa Symbol # simple import case
                # Save a singletons as a 1-element Vector for consistency with multiple linked params
                var_par = right.head == :tuple ? parse_dotted_symbols.(right.args) : [parse_dotted_symbols(right)]
                push!(imports, (left, var_par))
                # @info "import as $left = $var_par"

            # note that `comp_Symbol.name_Symbol` failed; bug in MacroTools?
            elseif @capture(left, comp_.name_) # simple connection case
                dst = parse_dotted_symbols(left)
                dst === nothing && error("Expected dot-delimited sequence of symbols, got $left")

                src = parse_dotted_symbols(right)
                src === nothing && error("Expected dot-delimited sequence of symbols, got $right")

                push!(conns, (dst, src))
                # @info "connection: $dst = $src"

            else
                error("Unrecognized expression on left hand side of '=' in @defcomposite: $elt")
            end
        else
            error("Unrecognized element in @defcomposite: $elt")
        end
    end

    # @info "imports: $imports"
    # @info "  $(length(imports)) elements"
    # global IMP = imports

    result = :(
        let conns = $conns,
            imports = $imports,

            cc_id = Mimi.ComponentId($calling_module, $(QuoteNode(cc_name)))

            global $cc_name = Mimi.CompositeComponentDef(cc_id, $(QuoteNode(cc_name)), $comps, $__module__)

            # @info "Defining composite $cc_id"

            function _store_in_ns(refs, local_name)
                isempty(refs) && return

                if length(refs) == 1
                    $cc_name[local_name] = refs[1]
                else
                    # We will eventually allow linking parameters, but not variables. For now, neither.
                    error("Variables and parameters may only be aliased individually: $refs")
                end
            end

            # This is more complicated than needed for now since we're leaving in place some of
            # the structure to accommodate linking multiple parameters to a single imported name.
            # We're postponing this feature to accelerate merging the component branch and will
            # return to this later.
            for (local_name, item) in imports
                var_refs = []
                par_refs = []

                for (src_path, src_name) in item
                    dr = Mimi.DatumReference(src_name, $cc_name, src_path)
                    if Mimi.is_parameter(dr)
                        push!(par_refs, Mimi.ParameterDefReference(dr))
                    else
                        push!(var_refs, Mimi.VariableDefReference(dr))
                    end
                end

                _store_in_ns(var_refs, local_name)
                _store_in_ns(par_refs, local_name)
            end

            for ((dst_path, dst_name), (src_path, src_name)) in conns
                # @info "connect_param!($(nameof($cc_name)), $dst_path, :$dst_name, $src_path, :$src_name)"
                Mimi.connect_param!($cc_name, dst_path, dst_name, src_path, src_name)
            end

            # propagates unconnected parameters
            Mimi.propagate_params!($cc_name)

            $cc_name
        end
    )

    # @info "defcomposite:\n$result"
    return esc(result)
end

nothing


text:
using Test
using Cleaner: CleanTable, get_all_repeated, level_distribution, compare_table_columns

@testset "get_all_repeated is working as expected" begin
    testNT = (; A=["y", "x", "y"], B=["x", "x", "x"])

    result = get_all_repeated(testNT, [:A])
    @test result.row_index == [1, 3]
    @test result.A == ["y", "y"]

    result = get_all_repeated(testNT, [:A, :B])
    @test result.row_index == [1, 3]
    @test result.A == ["y", "y"]
    @test result.B == ["x", "x"]

    let err = nothing
        try
            get_all_repeated(testNT, [:C])
        catch err
        end

        @test err isa Exception
        @test sprint(showerror, err) == "All column names specified must exist in the table"
    end

    if Threads.nthreads() > 1
        testCT = CleanTable(
            [:A, :B],
            [
                append!(repeat([1, 2], 2), collect(5:1_000_000)),
                append!(repeat([1, 2], 2), collect(5:1_000_000)),
            ],
        )
        result = get_all_repeated(testCT, [:A, :B])

        @test result.row_index == [2, 4, 1, 3]
        @test result.A == [2, 2, 1, 1]
        @test result.B == [2, 2, 1, 1]
    end
end

@testset "level_distribution is working as expected" begin
    testNT = (; A=["y", "x", "y"], B=["x", "x", "x"])

    result = level_distribution(testNT, [:A, :B])
    @test result.value == [Any["y", "x"], Any["x", "x"]]
    @test result.percent == [0.667, 0.333]

    result = level_distribution(testNT, [:A, :B]; round_digits=2)
    @test result.value == [Any["y", "x"], Any["x", "x"]]
    @test result.percent == [0.67, 0.33]

    let err = nothing
        try
            level_distribution(testNT, [:C])
        catch err
        end

        @test err isa Exception
        @test sprint(showerror, err) == "All column names specified must exist in the table"
    end

    if Threads.nthreads() > 1
        testCT = CleanTable(
            [:A, :B],
            [
                append!(repeat([1, 2], 200_000), repeat([3, 4], 300_000)),
                append!(repeat([1, 2], 200_000), repeat([3, 4], 300_000)),
            ],
        )
        result = level_distribution(testCT, [:A, :B])
        @test length(result.value) == 4
        @test result.percent == [0.3, 0.3, 0.2, 0.2]
    end
end

@testset "compare_table_columns is working as expected" begin
    testCT = CleanTable([:A, :A, :B], [[1, 2, 3], [4, 5, 6], String["7", "8", "9"]])

    testCT2 = CleanTable(
        [:A, :A, :C], [[missing, missing, missing], [1, missing, 3], ["x", "", "z"]]
    )

    result = compare_table_columns(testCT, testCT2)
    @test result.column_name == [:A, :A_1, :B, :C]
    @test result.tbl1 == [Int, Int, String, Nothing]
    @test result.tbl2 == [Missing, Union{Missing, Int}, Nothing, String]

    let err = nothing
        try
            compare_table_columns(testCT, testCT2; dupe_sanitize=false)
        catch err
        end

        @test err isa Exception
        @test sprint(showerror, err) == "Duplicated column name 'A' found on table number 1 passed. \nPlease use `dupe_sanitize=true`, `polish_names` or sanitize duplicated column names manually"
    end
end

ext: jl

max_stars_repo_name: TheRoniOne/Cleaner

max_issues_repo_name: TheRoniOne/Cleaner

max_forks_repo_name: TheRoniOne/Cleaner


text:
function nlp_cvx_202_010(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    # Test Goals:
    # - linear objective
    # - intersection convex quadratic constraints
    # Variants
    #   010 - one binding constraint (inflection point)
    #   011 - one binding constraint (inflection point)
    #   012 - one binding constraint (non-inflection point)
    #   013 - one binding constraint (non-inflection point)
    #   014 - intersection set

    model = Model(optimizer)

    @variable(model, x)
    @variable(model, y)
    @variable(model, z)

    @objective(model, Min, -z)
    @NLconstraint(model, x^2 + y^2 <= z)
    @NLconstraint(model, x^2 + y^2 <= -z + 1)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, -1, tol = objective_tol)
    return check_solution([x, y, z], [0, 0, 1], tol = primal_tol)
end


text:
using Revise
using VolterraGP
using Flux
using Flux:@adjoint

pars = (threeEQs, collect(1.:1.:5.), 0.1 * collect(1.:1.:5.), [0.01])

function fill_test(f, ts, Gs, u)
	reduce(hcat,
				[[f(ts[j], ts[i], Gs[j], Gs[i], u)
			for i in 1:size(ts)[1]]
		for	j in 1:size(ts)[1]])
end

function fill_phi_zip(f, ts, Gs, u)
	reduce(hcat,
				[[f(ti, tpi, pi, ppi, u)
			for (tpi, ppi) in zip(ts, Gs)]
		for	(ti, pi) in zip(ts, Gs)])
end

@adjoint function fill_phi_zip(f, ts, Gs, u)

    Gadj = zip_adj_G(f, ts, Gs, u)
    uadj = zip_adj_u(f, ts, Gs, u)

    fill_phi_zip(f, ts, Gs, u), Δ -> (nothing, nothing,   Δ .* Gadj, Δ .* uadj)
end

function zip_adj_G(f, ts, Gs, u)
    N = size(ts)[1]
    Gadj = zeros((N, N, N))
    for n = 1:N
        for i = 1:N
            # only do upper triangle cos of symmetry
            for j = i:N
                if i == n && j == n
                    g = gradient(x -> f(ts[i], ts[j], x, x, u), Gs[i])[1]
                    Gadj[i, j, n], Gadj[j, i, n] = g, g
                elseif i == n
                    g = gradient(x -> f(ts[i], ts[j], x, Gs[j], u), Gs[i])[1]
                    Gadj[i, j, n], Gadj[j, i, n] = g, g
                elseif j == n
                    g = gradient(x -> f(ts[i], ts[j], Gs[i], x, u), Gs[j])[1]
                    Gadj[i, j, n], Gadj[j, i, n] = g, g
                end
            end
        end
    end
    Gadj
end

function zip_adj_u(f, ts, Gs, u)
    N = size(ts)[1]
    uadj = zeros((N, N, size(u)[1]))
    for i = 1:N
        for j = 1:N
            uadj[i, j, :] = gradient(x -> f(ts[i], ts[j], Gs[i], Gs[j], x), u)[1]
        end
    end
    uadj
end

# @btime Gadj = zip_adj_G(pars...)
# @btime uadj = zip_adj_u(pars...)
ftest(t, tp, Ga, Gb, u) = (Ga + Gb) * (t + tp) + Gb * Ga * (t * tp)^2 + u[1]

# pars = (, [1., 2.], [2., 3.], [5.])
# a = fill_phi_zip(pars...)
@btime gradient(x -> sum(fill_phi_zip(ftest, collect(0.0:0.1:1), x, [0.1])), collect(0.0:0.1:1))

@btime gradient(x -> sum(fill_test(ftest, collect(0.0:0.1:1), x, [0.1])), collect(0.0:0.1:1))

y, back = Flux.pullback(fill_test, ftest, [1., 2.], [1., 1.], [0.1])


text:
module Verifier

using Random
using ..Utils
using ..Optimizer

using ProgressMeter: @showprogress
using nGraph

export verify

# We need to be able to verify that the output of modified graphs matches that of a standard
# graph.
#
# To do that, we use an auxiliary function that accepts the constructor function for a model
# as well as any passes to be performed on the model and checks that the base computation
# matches the modifed computation.

# Get all of the input and output parameters of an executable.
io_arrays(fex) = vflatten(nGraph.splat_inputs(fex), nGraph.splat_outputs(fex))

function verify(backend::nGraph.Backend{nGraph.CPU}, f, opt, cache;
        seed = 8086,
        inner_iterations = 5,
        outer_iterations = 2,
    )

    success = true
    for i in 1:outer_iterations
        println("Working on iteration $i")

        # Create a new seed for each outer iteration.
        outer_seed = rand(UInt64)

        # Wrapped the callable in a seed generator
        seeded_f = @closure begin
            Random.seed!(outer_seed)
            return f()
        end

        # Create a baseline copy and an optimized copy of the function.
        fex_baseline = actualize(backend, seeded_f)
        fex_optimized = first(Optimizer.factory(
            backend,
            seeded_f,
            opt;
            cache = cache,
            search_ratio = false
        ))

        # Sanity check - make sure all input tensors are equivalent.
        # Because some output tensors may be undefined - don't check those until the
        # executables are run for the first time.
        for (a, b) in zip(nGraph.splat_inputs(fex_baseline), nGraph.splat_inputs(fex_optimized))
            if parent(a) != parent(b)
                equal_error(i, 0, seed, outer_seed)
                return false
            end
        end

        for j in 1:inner_iterations
            println("    Inner Iteration $j")

            # Run both functions - make sure ALL input and output tensors match each itertion.
            @time fex_baseline()
            @time fex_optimized()
            for (a,b) in zip(io_arrays(fex_baseline), io_arrays(fex_optimized))
                # Check for NaN's
                if any(isnan, parent(a)) || any(isnan, parent(b))
                    nan_error(i, j, seed, outer_seed)
                    return false
                end

                # Check for Subnormal
                if any(issubnormal, parent(a)) || any(issubnormal, parent(b))
                    subnormal_error(i, j, seed, outer_seed)
                    return false
                end

                # Check equality
                #
                # Use `isapprox` because of non-determinism in GPU computations.
                if !isapprox(parent(a), parent(b))
                    equal_error(i, j, seed, outer_seed)
                    return (a, b)
                end
            end
        end
    end
    return true
end

# For the GPU, nondeterminism means we have to trace over a period of time.
function verify(backend::nGraph.Backend{nGraph.GPU}, f, opt, cache;
        seed = 8086,
        inner_iterations = 100,
        outer_iterations = 2,
    )

    success = true
    for i in 1:outer_iterations
        println("Working on iteration $i")

        # Create a new seed for each outer iteration.
        outer_seed = rand(UInt64)

        # Wrapped the callable in a seed generator
        seeded_f = @closure begin
            Random.seed!(outer_seed)
            return f()
        end

        # GC fence everything to ensure we don't run out of memory.
        GC.gc()
        baseline = run_baseline(backend, seeded_f, opt, cache, inner_iterations)
        GC.gc()
        optimized = run_optimized(backend, seeded_f, opt, cache, inner_iterations)
        GC.gc()
        return baseline, optimized
    end
end

# Inner functions for GC purposes.
function run_baseline(backend, f, opt, cache, inner_iterations)
    fex = actualize(backend, f)
    vals = []
    @showprogress 1 for j in 1:inner_iterations
        push!(vals, collect(parent(fex())))
    end
    return vals
end

function run_optimized(backend, f, opt, cache, inner_iterations)
    #fex = actualize(backend, f)
    fex = first(Optimizer.factory(
        backend,
        f,
        opt;
        cache = cache,
    ))

    vals = []
    @showprogress 1 for j in 1:inner_iterations
        push!(vals, collect(parent(fex())))
    end
    return vals
end

#####
##### Error Generation.
#####

nan_error(x...) = generic_error("NaN Values", x...)
subnormal_error(x...) = generic_error("Subnormal Values", x...)
equal_error(x...) = generic_error("Match Error", x...)

function generic_error(prefix, outer_iteration, inner_iteration, seed, outer_seed)
    message = """
    $prefix
    Outer Iteration: $outer_iteration
    Inner Iteration: $inner_iteration
    Master seed: $seed
    Seed for this round: $outer_seed
    """
    printstyled(message; color = :red)
    println()
    return nothing
end

end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
abstract type TPs1a0s0 end # Testimate that uses probabilities (Ps1a0s0)

abstract type TNs1a0s0 end # Testimate that uses counts (Ns1a0s0)


text:
module QNaN

import Base:convert

export qnan

convert(::Type{Signed}, x::Type{Float16}) = abs(x)<=typemax(Int16) ? trunc(Int16, x) : throw(ArgumentError("too large"))
convert(::Type{Signed}, x::Type{Float32}) = abs(x)<=typemax(Int32) ? trunc(Int32, x) : throw(ArgumentError("too large"))
convert(::Type{Signed}, x::Type{Float64}) = abs(x)<=typemax(Int64) ? trunc(Int64, x) : throw(ArgumentError("too large"))

for (FL, I, UI, UPos, UNeg) in [(:Float64, :Int64, :UInt64, :0x7ff8000000000000, :0xfff8000000000000),
                                (:Float32, :Int32, :UInt32, :0x7fc00000, :0xffc00000),
                                (:Float16, :Int16, :UInt16, :0x7e00, :0xfe00) ]
  @eval begin
      @inline isqnan(x::$(UI))    = (x & $(UPos)) == $(UPos)
      @inline isqnan(x::$(FL))    = isqnan(reinterpret($(UI),x))

      function qnan(si::$(I))
          u = reinterpret($(UI), abs(si))
          if (u > ~$(UNeg)) # 2^51-1
              throw(ArgumentError("The magnitude of n $(si) exceeds available QNaN."))
          end
          u |= ((si > 0) ? $(UPos) : $(UNeg))
          reinterpret($(FL),u)
      end

      function qnan(fp::$(FL))
          u = reinterpret($(UI), fp)
          if !isqnan(u)
              if fp == trunc(fp)
                  return qnan(convert(Signed,fp))
              else
                  throw(ArgumentError("The value $(fp) is not a valid QNaN payload."))
              end
          end
          a = u & ~$(UNeg)
          b =  reinterpret($(I),a)
          (((u & ($(UPos) $ $(UNeg))) == 0) ? b : -b)
      end
  end
end

@doc """
Module QNaN
  qnan(si::{Int16,Int32,Int64}) generates a quiet NaN with a payload of `si`
  qnan(fp::FloatingPoint) recovers the Signed payload from `fp`
""" -> QNaN

@doc """
qnan(`si`::{Int16|32|64}) generates a quiet NaN using `si` for the payload
qnan(`fp`::FloatingPoint) recovers the signed payload from a quiet NaN `fp`
""" -> qnan

end # module


text:
#################################################################
###### Functionality related to observation dimensionality ######
#################################################################
		# Note: The functions are designed to work on two
		# dimensional data containers (i.e. matrices)

# Opposide dimension function
"""
	oppdim(::LearnBase.ObsDimension)

Returns the other dimension for a matrix i.e. if provided `ObsDim.Constant{1}`
returns `ObsDim.Constant{2}`.
"""
oppdim(::ObsDim.First) = ObsDim.Constant{2}()
oppdim(::ObsDim.Last) = ObsDim.First()
oppdim(::ObsDim.Constant{2}) = ObsDim.Constant{1}()

# get integer dimension
"""
	intdim(::LearnBase.ObsDimension)

Returns the integer associated to a dimension object
i.e. `intdim(ObsDim.Constant{3})`  returns `3`.
The function is designed to work on matices so
`intdim(::ObsDim.Last)` will return `2`.
"""
intdim(::ObsDim.First) = 1
intdim(::ObsDim.Last) = 2
intdim(::ObsDim.Constant{N}) where N = N

# nvars function
"""
Returns the number of variables given a data object which
supports the `nobs` function. The data object must ideally
present two dimensions i.e. matrix.
"""
nvars(X, arg) = nobs(X,oppdim(arg))

# pre-allocation
"""
	matrix_prealloc(no, nv, obsdim, val)

Returns a `Matrix{T}` filled with values equal to `val::T`, having
the size `no` (number of observations) on dimension `obsdim` and
`nv` (number of variables) in the other dimension.
"""
function matrix_prealloc(no::Int, nv::Int, obsdim::O, val::T=zero(T)) where {
		T,O <:LearnBase.ObsDimension}
	# get dimensions (inner function)
	_getdims_(no, nv, ::ObsDim.First) = no, nv
	_getdims_(no, nv, ::ObsDim.Last) = nv, no
	_getdims_(no, nv, ::ObsDim.Constant{2}) = nv, no
	_getdims_(no, nv, ::ObsDim.Undefined) = error("Undefined observation dimension")

	m,n = _getdims_(no, nv, obsdim)
	M = Array{T}(undef, m,n)
	fill!(M,val)

	return M
end

#################################################################
################# Additional utility functions ##################
#################################################################

# Function that calculates the number of  relational variables / each adjacency structure
get_size_out(y::AbstractVector{T}) where T<:Float64 = 1			# regression case
get_size_out(y::AbstractVector{T}) where T = length(unique(y))::Int	# classification case
get_size_out(y::AbstractArray) = error("Only vectors supported as targets in relational learning.")

# Function that calculates the priors of the dataset
getpriors(y::AbstractVector{T}) where T<:Float64 = [1.0]
getpriors(y::AbstractVector{T}) where T = [sum(yi.==y)/length(y) for yi in sort(unique(y))]
getpriors(y::AbstractArray) = error("Only vectors supported as targets in relational learning.")

encode_targets(labels::T where T<:AbstractVector{S}) where S = begin
	ulabels::Vector{S} = sort(unique(labels))
	enc = LabelEnc.NativeLabels(ulabels)
	return (enc, label2ind.(labels,enc))
end

encode_targets(labels::T where T<:AbstractVector{S}) where S<:AbstractFloat = begin
	return (nothing, labels)
end

encode_targets(labels::T where T<:AbstractArray{S}) where S = begin
	error("Targets must be in vector form, other arrays not supported.")
end

read_citation_data(content_file::String, cites_file::String) = begin

	# Read files
	content = readdlm(content_file,'\t')
	cites = readdlm(cites_file,'\t')

	# Construct datasets
	labels = content[:,end]
	paper_ids = Int.(content[:,1])
	data = Float64.(content[:,2:end-1]')
	content_data = (data, labels)

	# Construct citing/cited paper indices
	cited_papers = indexin(Int.(cites[:,1]), paper_ids)
	citing_papers = indexin(Int.(cites[:,2]), paper_ids)

	return content_data, cited_papers, citing_papers

end

# Function that grabs the Cora dataset
grab_cora_data(tmpdir::String="/tmp") = begin

	DATA_FILE = download("https://linqs-data.soe.ucsc.edu/public/lbc/cora.tgz")
	run(`tar zxvf $(DATA_FILE) -C $tmpdir`)

	cora_data = read_citation_data("$tmpdir/cora/cora.content","$tmpdir/cora/cora.cites")

	run(`rm -rf $tmpdir/cora`)
	run(`rm $DATA_FILE`)

	return cora_data
end

# Function that generates an adjacency matrix based on the citing and cited paper indices as well as
# the indices in these vectors of the citations that are to be considered
function generate_partial_adjacency(cited::T, citing::T, useidx::S) where {T<:AbstractVector, S<:AbstractVector}

	# Search which citing/cited papers appear in the subset indices
	citing_idx = indexin(cited, useidx)
	cited_idx = indexin(citing, useidx)

	# Construct adjacency
	n = length(useidx)
	W = zeros(n,n)
	for i in 1:length(citing)
	  if citing_idx[i] != nothing && cited_idx[i] != nothing
	      W[citing_idx[i],cited_idx[i]] += 1
	      W[cited_idx[i],citing_idx[i]] += 1
	  end
	end
	return W
end

macro print_verbose(level, message)
	esc( :(
	 	NetworkLearning.VERBOSE >= 0 &&
	 	NetworkLearning.VERBOSE >= $level &&
	 	print_with_color(:cyan, "[Network Learning] $($message)\n")
		)
	)
end


text:
using Test
using OptionalArgChecks

function f(x)
    @argcheck x<2 && if x == 1
        return "foo"
    else
        error("Test")
    end
    return x + 3
end

for x in 0:2
    @test @skipargcheck f(x) == x + 3
end

using Documenter
DocMeta.setdocmeta!(OptionalArgChecks, :DocTestSetup, :(using OptionalArgChecks); recursive=true)
doctest(OptionalArgChecks; manual = false)


text:
println("... basic tests: processes ...")
simex = Simulate.SimException(Simulate.Stop())
@test simex.ev == Simulate.Stop()
@test isnothing(simex.value)

# ===== test process registration
ch1 = Channel(32)
ch2 = Channel(32)

reset!(𝐶)
incr(c1::Channel, c2::Channel, a) = (a+1, yield())
a = [1,1,3.0,3.0,"A","A","A","A"]
b = [1,2,3.0,nextfloat(3.0),"A","A#1","A#2","A#3"]
for i in 1:8
    @test process!(SP(a[i], incr, ch1, ch2, 1)) == b[i]
end
for i in 1:8
    @test 𝐶.processes[b[i]].id == b[i]
end
@test process!(SP((1,2), incr, ch1, ch2, 1)) == (1,2)
@test_throws ArgumentError process!(SP((1,2), incr, ch1, ch2, 1))

println("... test channel 4 example ...")
A = []

function simple(input::Channel, output::Channel, name, id, op)
    token = take!(input)         # take something from the input
    push!(A, (τ(), name, id, token))
    d = delay!(rand())           # after a delay
    put!(output, op(token, id))  # put it out with some op applied
end

reset!(𝐶)
Random.seed!(123)

for i in 1:2:8    # create, register and start 8 SimProcesses SP
    process!(SP(i, simple, ch1, ch2, "foo", i, +))
    process!(SP(i+1, simple, ch2, ch1, "bar", i+1, *))
end

@test length(𝐶.processes) == 8
for p in values(𝐶.processes)
    @test p.state == Simulate.Idle()
    @test istaskstarted(p.task)
end

put!(ch1, 1)
sleep(0.01)
run!(𝐶, 10)

@test length(A) > 20
p = [i[3] for i in A]
for i in 1:8
    @test i ∈ p  # all processes did something
end

for p in values(𝐶.processes)
    @test istaskstarted(p.task)
    stop!(p, Simulate.Stop())
    @test istaskdone(p.task)
end

println("... test wait! ...")

ch1 = Channel(1)
ch2 = Channel(1)

a = 1
b = 1
res = []
incb() = global b +=1
checktime(x) = τ() ≥ x
checka(x) = a == x
checkb(x) = b ≥ x

function testwait(c1::Channel, c2::Channel)
    wait!((SF(checktime, 2), SF(checka, 1)))
    push!(res, (τ(), 1, a, b))
    wait!(SF(isa, a, Int)) # must return immediately
    push!(res, (τ(), 2, a, b))
    sample!(SF(incb))
    wait!(SF(checkb, 201))
    push!(res, (τ(), 3, a, b))
    take!(c1)
end

reset!(𝐶)
process!(SP(1, testwait, ch1, ch2))

run!(𝐶, 10)
r = [i[1] for i in res]
@test r[1] ≈ 2
@test r[2] ≈ 2
@test r[3] ≈ 4
@test res[3][4] == 201
@test b == 801


text:
# 3D FWI of the Overthrust model with stochastic SP
# Author: Philipp Witte
# Date: September 2017
#

Pkg.checkout("opesciSLIM","development")
using PyCall, PyPlot, JLD, HDF5, opesciSLIM.TimeModeling, opesciSLIM.SLIM_optim, SeisIO

# Load velocity model
n,d,o,m0 = h5open("../data/overthrust_3D_initial_model.h5","r") do file
	read(file, "n", "d", "o", "m0")
end

# Set up model structure
model0 = Model((n[1],n[2],n[3]), (d[1],d[2],d[3]), (o[1],o[2],o[3]), m0)

# Bound constraints
vmin = ones(Float32,model0.n) + 0.4f0
vmax = ones(Float32,model0.n) + 5.5f0
mmin = vec((1f0./vmax).^2); vmax=[]
mmax = vec((1f0./vmin).^2); vmin=[]

# Load data container containing filenames
container = segy_scan("../data","overthrust_3D_shot",["GroupX","GroupY","dt","RecGroupElevation","SourceSurfaceElevation"])
dobs = joData(container)

# Set up source
src_geometry = Geometry(container; key="source")
wavelet = source(src_geometry.t[1], src_geometry.dt[1], 0.008)	# 8 Hz peak frequency
q = joData(src_geometry, wavelet)

############################### FWI ###########################################

# Optimization parameters
srand(1)	# set seed of random number generator
fevals = 15
batchsize = 1080
count = 6
fvals = zeros(50)
opt = Options(limit_m=true,
			  save_rate=14.,		# sampling rate of gradient
			  buffer_size=500.
			  )

function objective_function(x)

	# update model and save snapshot
	model0.m = reshape(x,model0.n);
	save(join(["../data/nshot_1080_snapshot_minConf",string(count),".jld"]),"v",sqrt.(1./model0.m))

	# select batch
	idx = randperm(dobs.nsrc)[1:batchsize]
	dsub = subsample(dobs,idx)
	qsub = subsample(q,idx)

	# fwi function value and gradient
	fval, grad = fwi_objective(model0,qsub,dsub;options=opt); gc()
	grad = reshape(grad,model0.n); grad[:,:,1:21] = 0.f0	# reset gradient in water column to 0.
	grad = .125f0*grad/maximum(abs(grad))	# scale gradient to help line search
	global count; count += 1; fvals[count] = fval

    return fval, vec(grad)
end

# Bound projection
ProjBound(x) = median([mmin x mmax],2)

# FWI with SPG
spg_opt = spg_options(verbose=3, maxIter=fevals, memory=1)
x, fsave, funEvals= minConf_SPG(objective_function, vec(model0.m), ProjBound, spg_opt)

save("../data/results_stochastic_nshot_1080.jld","x",x,"fsave",fsave,"funEvals",funEvals)

max_stars_repo_name: SINBADconsortium/SLIM-release-jlapps

max_issues_repo_name: SINBADconsortium/SLIM-release-jlapps

max_forks_repo_name: SINBADconsortium/SLIM-release-jlapps


text:
# This file holds the definition of the functions pertaining to the
#   PersonnelAttribute type.

# The functions of the PersonnelAttribute type require no additional types.
requiredTypes = [ "personnelAttribute" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include( joinpath( typePath, reqType * ".jl" ) )
    end  # if !isdefined( Symbol( ...
end  # for reqType in requiredTypes

export setOrdinal!,
       setPossibleValues!,
       setAttrValues!,
       setFixed!,
       addValueToAttr!,
       removeValueFromAttr!,
       generateAttrValue

"""
```
setOrdinal!( attr::PersonnelAttribute,
             isOrdinal::Bool )
```
This function sets the isOrdinal flag of the personnel attribute `attr` to
`isOrdinal`.

This function returns `nothing`.
"""
function setOrdinal!( attr::PersonnelAttribute, isOrdinal::Bool )::Void

    attr.isOrdinal = isOrdinal
    return

end  # setOrdinal!( attr, isOrdinal )

"""
```
setPossibleValues!( attr::PersonnelAttribute,
                    vals::Vector{String} )
```
This function sets the list of possible values of the personnel attribute `attr`
to the list given in `vals`.

This function returns `nothing`.
"""
function setPossibleValues!( attr::PersonnelAttribute,
    vals::Vector{String} )::Void

    attr.possibleValues = unique( vals )
    return

end  # setPossibleValues!( attr, vals )

"""
```
setAttrValues!( attr::PersonnelAttribute,
                vals::Dict{String, Float64} )
```
This function sets the list of possible values of the personnel attribute `attr`
to the list given in `vals`. The function ignores values which have an
accompanying probability weight ⩽ 0, and will reweight the valid entries in the
list such that their sum of probabilities equals 1. If no valid entries remain,
the list is wiped, and a warning is issued.

This function returns `nothing`.
"""
function setAttrValues!( attr::PersonnelAttribute,
    vals::Dict{String, Float64} )::Void

    # Remove the entries with a probability weight <= 0.
    tmpVals = Dict{String, Float64}()
    foreach( val -> if vals[ val ] > 0 tmpVals[ val ] = vals[ val ] end,
        keys( vals ) )
        # XXX Actually removing them can give trouble.

    # Normalise the probability weights of the entries to 1.
    if isempty( tmpVals )
        warn( "No valid entries in list of possible values for the attribute. Wiping the list." )
    else
        probMass = sum( map( val -> tmpVals[ val ], keys( tmpVals ) ) )
        foreach( val -> tmpVals[ val ] /= probMass, keys( tmpVals ) )
    end

    attr.values = tmpVals
    return

end  # setAttrValues!( attr, vals )

"""
```
setFixed!( attr::PersonnelAttribute,
           isFixed::Bool )
```
This function sets the isFixed flag of the personnel attribute `attr` to
`isFixed`.

This function returns `nothing`.
"""
function setFixed!( attr::PersonnelAttribute, isFixed::Bool )::Void

    attr.isFixed = isFixed
    return

end  # setFixed!( attr, isFixed )

"""
```
addValueToAttr!( attr::PersonnelAttribute,
                 val::String,
                 prob::T )
    where T <: Real
```
This function adds a possible initialisation value to the personnel attribute
`attr`. The new value is `val` and will be generated with probability `prob`
which must lie between 0.0 and 1.0 (endpoints exclusive). The generation
probabilities of the other possible values are adjusted accordingly by a factor
`1 - prob`.

If an impossible probability is entered, the function does nothing and returns a
warning. If the value to be added is already a possibility, its probability will
be adjusted, as well as the probability of the other possible values. If the
list was empty, the new value will have probability 100%.

The function returns `nothing`.
"""
function addValueToAttr!( attr::PersonnelAttribute, val::String, prob::T )::Void where T <: Real

    # Check if the entered probabilty makes sense.
    if ( prob <= 0 ) || ( prob >= 1 )
        warn( "Probability of value occurring must be between 0 and 1. Not adding value." )
        return
    end  # if ( prob <= 0 ) || ( prob >= 1 )

    # Check if there are already entries in the list.
    if isempty( attr.values )
        attr.values[ val ] = 1
        return
    end  # if isempty( attr.values )

    # Check if the new value is already in the list..
    if haskey( attr.values, val )
        warn( "Value is already in the list, changing probabilities." )
    end  # if haskey( attr.values, val )

    # Set cumulative probabilty of existing entries to 1 - prob
    adjFactor = ( 1 - prob ) / ( 1 - get( attr.values, val, 0 ) )
    foreach( tmpVal -> attr.values[ tmpVal ] *= adjFactor, keys( attr.values ) )
    attr.values[ val ] = prob
    return

end  # addValueToAttr!( attr, val, prob )

"""
```
removeValueFromAttr!( attr::PersonnelAttribute,
                      val::String )
```
This function removes the value `val` from the possible values of the personnel
attribute `attr`, and adjusts the probabilities accordingly.

This function returns `nothing`.
"""
function removeValueFromAttr!( attr::PersonnelAttribute, val::String )::Void

    # Do nothing if the value isn't in the list of possible values.
    if !haskey( attr.values, val )
        return
    end  # if !haskey( attr.values, val )

    prob = attr.values[ val ]
    foreach( tmpVal -> attr.values[ tmpVal ] /= 1 - prob, keys( attr.values ) )
    delete!( attr.values, val )
    return

end  # removeValueFromAttr!( attr, val )

function Base.show( io::IO, attr::PersonnelAttribute )

    print( io, "  Attribute: $(attr.name) (" )
    print( io, attr.isFixed ? "fixed" : "variable" )
    print( io, ")" )

    if isempty( attr.possibleValues )
        print( io, "\n    Attribute can't take any values." )
        return
    end  # if isempty( attr.possibleValues )

    print( io, "\n    Possible values: " )
    print( io, join( attr.possibleValues, ", " ) )

    if isempty( attr.values )
        print( io, "\n    No possible initial values" )
        return
    end  # if isempty( attr.values )

    print( io, "\n    Initial values:" )
    foreach( val -> print( io, "\n      $val ($(signif(attr.values[ val ] * 100, 3))%)" ),
        keys( attr.values ) )

end  # show( io, attr )

"""
```
generateAttrValue( attr::PersonnelAttribute )
```
This function generates a value for the personnel attribute `attr` and returns
it as a `String`. If the attribute has no possible initialisation values, the
function returns `nothing`.
"""
function generateAttrValue( attr::PersonnelAttribute )

    if isempty( attr.values )
        return
    end  # if isempty( attr.values )

    vals = collect( keys( attr.values ) )
    return vals[ rand( Categorical( map( val -> attr.values[ val ], vals ) ) ) ]

end  # generateAttrValue( attr )

# ==============================================================================
# Non-exported methods.
# ==============================================================================

function readAttribute( sheet::XLSX.Worksheet, attrCat::XLSX.Worksheet,
    sLine::Integer )::PersonnelAttribute

    newAttr = PersonnelAttribute( string( sheet[ "A$sLine" ] ) )
    catLine = sheet[ "B$sLine" ]

    if isa( catLine, Missings.Missing )  # XXX Check what happens with N/A
        error( "Attribute '$(newAttr.name)' not defined in catalogue." )
    end  # if isa( catLine, Void )

    catLine = catLine + 1
    setFixed!( newAttr, attrCat[ "B$catLine" ] == "YES" )
    setOrdinal!( newAttr, attrCat[ "C$catLine" ] == "YES" )
    nVals = attrCat[ "E$catLine" ]
    vals = Vector{String}( nVals )

    for ii in 1:nVals
        vals[ ii ] = strip( attrCat[ XLSX.CellRef( catLine, 5 + ii ) ] )
    end  # for ii in 1:nVals

    setPossibleValues!( newAttr, vals )
    ii = 1
    nInitVals = sheet[ "B$(sLine + 2)" ]
    vals = Dict{String, Float64}()

    for ii in (1:nInitVals) + 2
        val = sheet[ XLSX.CellRef( sLine, ii ) ]
        weight = sheet[ XLSX.CellRef( sLine + 1, ii ) ]

        if isa( weight, Real )
            vals[ val ] = weight + get( vals, val, 0.0 )
        end  # if isa( weight, Real )
    end  # for ii in (1:nInitVals) + 2

    setAttrValues!( newAttr, vals )
    return newAttr

end  # readAttribute( sheet, attrCat, sLine )


text:
# Times out with the default limit of 30 seconds (package load time, precompilation, etc.).
const JUPYTER_BIN = readstring(open(joinpath(Pkg.dir("IJulia"), "deps", "JUPYTER")))
const CMD = `$JUPYTER_BIN nbconvert --ExecutePreprocessor.timeout=600 --to html --output-dir=. --stdout --execute ../src/DistributionsPlot.jl.ipynb`
run(CMD)


text:
let part1 = 0, part2 = 1
    # This reads the input in a transposed format — rows here are columns in
    # the file.
    chart = hcat(map(collect, readlines("input.txt"))...)

    nrow, ncol = size(chart)

    for (run, rise) in [(1, 1), (3, 1), (5, 1), (7, 1), (1, 2)]
        trees, row = 0, 1

        for col in 1:rise:ncol
            if chart[row, col] == '#'
                # only count the tree for part1 if the slope is (3, 1)
                part1 += (run == 3 && rise == 1)

                # For part2
                trees += 1
            end

            # # We could use the following if statement which doesn't require
            # # as much mental overhead.
            # row += run
            #
            # # This is the only place where we have to think about the indexing.
            # # If we were using 0-indexing, we'd need `row ≥ nrow`.
            # if row > nrow
            #     row -= nrow
            # end

            # By default, Julia arrays are 1-indexed, so we have to subtract 1
            # before we `%` and add one after. I like this better because its
            # short, fairly clear if you know how `%` works, and doesn't
            # require the overhead of an if-statement.
            row = (row + run - 1) % nrow + 1
        end

        part2 *= trees
    end

    @info "Solution" part1 part2
end

max_stars_repo_name: dglmoore/advent-of-code

max_issues_repo_name: dglmoore/advent-of-code

max_forks_repo_name: dglmoore/advent-of-code


text:
############################################################
# joData # #################################################
############################################################

# Authors: Philipp Witte (pwitte@eos.ubc.ca), Henryk Modzelewski (hmodzelewski@eos.ubc.ca)
# Date: January 2017

export joData, joDataException, subsample, joData_to_SeisBlock

############################################################

# structure for seismic data as an abstract vector
mutable struct joData{vDT<:Number} <: joAbstractLinearOperator{vDT,vDT}
	name::String
	m::Integer
	n::Integer
	nsrc::Integer
	geometry::Geometry
	data
end

mutable struct joDataException <: Exception
	msg :: String
end

############################################################

## outer constructors

"""
    joData
	    name::String
        m::Integer
        n::Integer
        nsrc::Integer
        geometry::Geometry
        data

Abstract vector for seismic data. This vector-like structure contains the geometry and data for either\\
receiver data (shot records) or source data (wavelets).

Constructors
============

Construct vector from `Geometry` structure and cell array of shot records or wavelets. The `data` keyword\\
can also be a single (non-cell) array, in which case the data is the same for all source positions:

    joData(geometry, data)

Construct vector for observed data from `SeisIO.SeisBlock`. `segy_depth_key` is the `SeisIO` keyword \\
that contains the receiver depth coordinate:

    joData(SeisIO.SeisBlock; segy_depth_key="RecGroupElevation")

Construct vector for observed data from out-of-core data container of type `SeisIO.SeisCon`:

    joData(SeisIO.SeisCon; segy_depth_key="RecGroupElevation")

Examples
========

(1) Construct data vector from `Geometry` structure and a cell array of shot records:

    dobs = joData(rec_geometry, shot_records)

(2) Construct data vector for a seismic wavelet (can be either cell arrays of individual\\
wavelets or a single wavelet as an array):

    q = joData(src_geometry, wavelet)

(3) Construct data vector from `SeisIO.SeisBlock` object:

    using SeisIO
    seis_block = segy_read("test_file.segy")
	dobs = joData(seis_block; segy_depth_key="RecGroupElevation")

(4) Construct out-of-core data vector from `SeisIO.SeisCon` object (for large SEG-Y files):

    using SeisIO
    seis_container = segy_scan("/path/to/data/directory","filenames",["GroupX","GroupY","RecGroupElevation","SourceDepth","dt"])
    dobs = joData(seis_container; segy_depth_key="RecGroupElevation")

"""
function joData(geometry::Geometry,data::Array; vDT::DataType=Float32)
	vDT == Float32 || throw(joDataException("Domain type not supported"))

	# length of vector
	n = 1
	if typeof(geometry) == GeometryOOC
		nsrc = length(geometry.container)
		m = sum(geometry.nsamples)
	else
		nsrc = length(geometry.xloc)
		m = 0
		for j=1:nsrc
			m += length(geometry.xloc[j])*geometry.nt[j]
		end
	end
	dataCell = Array{Array}(nsrc)
	for j=1:nsrc
		dataCell[j] = data
	end
	return joData{Float32}("Seismic data vector",m,n,nsrc,geometry,dataCell)
end

# constructor if data is passed as a cell array
function joData(geometry::Geometry,data::Union{Array{Any,1},Array{Array,1}}; vDT::DataType=Float32)
	vDT == Float32 || throw(joDataException("Domain and range types not supported"))

	# length of vector
	if typeof(geometry) == GeometryOOC
		nsrc = length(geometry.container)
		m = sum(geometry.nsamples)
	else
		nsrc = length(geometry.xloc)
		m = 0
		for j=1:nsrc
			m += length(geometry.xloc[j])*geometry.nt[j]
		end
	end
	n = 1
	return joData{Float32}("Seismic data vector",m,n,nsrc,geometry,data)
end

############################################################
# constructors from SEGY files or out-of-core containers

# contructor for in-core data container
function joData(data::SeisIO.SeisBlock; segy_depth_key="RecGroupElevation", vDT::DataType=Float32)
	vDT == Float32 || throw(joDataException("Domain and range types not supported"))

	# length of data vector
	src = get_header(data,"FieldRecord")
 	nsrc = length(unique(src))

	numTraces = get_header(data,"TraceNumber")[end] - get_header(data,"TraceNumber")[1] + 1
	numSamples = get_header(data,"ns")[end]
	m = numTraces*numSamples
	n = 1

	# extract geometry from data container
	geometry = Geometry(data; key="receiver", segy_depth_key=segy_depth_key)

	# fill data vector with pointers to data location
	dataCell = Array{Array}(nsrc)
	full_data = convert(Array{Float32,2},data.data)
	for j=1:nsrc
		traces = find(src .== unique(src)[j])
		dataCell[j] = full_data[:,traces]
	end

	return joData{Float32}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for in-core data container and given geometry
function joData(geometry::Geometry, data::SeisIO.SeisBlock; vDT::DataType=Float32)
	vDT == Float32 || throw(joDataException("Domain and range types not supported"))

	# length of data vector
	src = get_header(data,"FieldRecord")
 	nsrc = length(unique(src))

	numTraces = get_header(data,"TraceNumber")[end] - get_header(data,"TraceNumber")[1] + 1
	numSamples = get_header(data,"ns")[end]
	m = numTraces*numSamples
	n = 1

	# fill data vector with pointers to data location
	dataCell = Array{Array}(nsrc)
	for j=1:nsrc
		traces = find(src .== unique(src)[j])
		dataCell[j] = convert(Array{Float32,2},data.data[:,traces])
	end

	return joData{Float32}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for out-of-core data container from single container
function joData(data::SeisIO.SeisCon; segy_depth_key="RecGroupElevation", vDT::DataType=Float32)
	vDT == Float32 || throw(joDataException("Domain and range types not supported"))

	# length of data vector
	nsrc = length(data)
	numTraces = 0
	for j=1:nsrc
		numTraces += Int((data.blocks[j].endbyte - data.blocks[j].startbyte)/(240 + data.ns*4))
	end
	m = numTraces*data.ns
	n = 1

	# extract geometry from data container
	geometry = Geometry(data; key="receiver", segy_depth_key=segy_depth_key)

	# fill data vector with pointers to data location
	dataCell = Array{SeisIO.SeisCon}(nsrc)
	for j=1:nsrc
		dataCell[j] = split(data,j)
	end

	return joData{Float32}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for single out-of-core data container and given geometry
function joData(geometry::Geometry, data::SeisIO.SeisCon; vDT::DataType=Float32)
	vDT == Float32 || throw(joDataException("Domain and range types not supported"))

	# length of data vector
	nsrc = length(data)
	numTraces = 0
	for j=1:nsrc
		numTraces += Int((data.blocks[j].endbyte - data.blocks[j].startbyte)/(240 + data.ns*4))
	end
	m = numTraces*data.ns
	n = 1

	# fill data vector with pointers to data location
	dataCell = Array{SeisIO.SeisCon}(nsrc)
	for j=1:nsrc
		dataCell[j] = split(data,j)
	end

	return joData{Float32}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for out-of-core data container from cell array of containers
function joData(data::Array{SeisIO.SeisCon,1}; segy_depth_key="RecGroupElevation", vDT::DataType=Float32)
	vDT == Float32 || throw(joDataException("Domain and range types not supported"))

	# length of data vector
	nsrc = length(data)
	numTraces = 0
	for j=1:nsrc
		numTraces += Int((data[j].blocks[1].endbyte - data[j].blocks[1].startbyte)/(240 + data.ns*4))
	end
	m = numTraces*data.ns
	n = 1

	# extract geometry from data container
	geometry = Geometry(data; key="receiver", segy_depth_key=segy_depth_key)

	# fill data vector with pointers to data location
	dataCell = Array{SeisIO.SeisCon}(nsrc)
	for j=1:nsrc
		dataCell[j] = data[j]
	end

	return joData{Float32}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

# contructor for out-of-core data container from cell array of containers and given geometry
function joData(geometry::Geometry, data::Array{SeisIO.SeisCon}; vDT::DataType=Float32)
	vDT == Float32 || throw(joDataException("Domain and range types not supported"))

	# length of data vector
	nsrc = length(data)
	numTraces = 0
	for j=1:nsrc
		numTraces += Int((data[j].blocks[1].endbyte - data[j].blocks[1].startbyte)/(240 + data[j].ns*4))
	end
	m = numTraces*data[1].ns
	n = 1

	# fill data vector with pointers to data location
	dataCell = Array{SeisIO.SeisCon}(nsrc)
	for j=1:nsrc
		dataCell[j] = data[j]
	end

	return joData{Float32}("Julia seismic data container",m,n,nsrc,geometry,dataCell)
end

############################################################
## overloaded Base functions

# conj(jo)
conj{vDT}(a::joData{vDT}) =
	joData{vDT}("conj("*a.name*")",a.m,a.n,a.nsrc,a.geometry,a.data)

# transpose(jo)
transpose{vDT}(a::joData{vDT}) =
	joData{vDT}(""*a.name*".'",a.n,a.m,a.nsrc,a.geometry,a.data)

# ctranspose(jo)
ctranspose{vDT}(a::joData{vDT}) =
	joData{vDT}(""*a.name*"'",a.n,a.m,a.nsrc,a.geometry,a.data)

##########################################################

# +(joData, joData)
function +{avDT,bvDT}(a::joData{avDT}, b::joData{bvDT})
	size(a) == size(b) || throw(joDataException("dimension mismatch"))
	compareGeometry(a.geometry, b.geometry) == 1 || throw(joDataException("geometry mismatch"))
	c = deepcopy(a)
	c.data = a.data + b.data
	return c
end

# -(joData, joData)
function -{avDT,bvDT}(a::joData{avDT}, b::joData{bvDT})
	size(a) == size(b) || throw(joDataException("dimension mismatch"))
	compareGeometry(a.geometry, b.geometry) == 1 || throw(joDataException("geometry mismatch"))
	c = deepcopy(a)
	c.data = a.data - b.data
	return c
end

# *(joData, number)
function *{avDT}(a::joData{avDT},b::Number)
	c = deepcopy(a)
	c.data = c.data*b
	return c
end

# *(number, joData)
function *{bvDT}(a::Number,b::joData{bvDT})
	c = deepcopy(b)
	c.data = a*c.data
	return c
end

# *(joLinearFunction, joData)
function *{ADDT,ARDT,avDT}(A::joLinearFunction{ADDT,ARDT},v::joData{avDT})
	A.n == size(v,1) || throw(joDataException("shape mismatch"))
	jo_check_type_match(ADDT,avDT,join(["DDT for *(joLinearFunction,joData):",A.name,typeof(A),avDT]," / "))
	V = A.fop(v)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joLinearFunction,joData):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

# *(joLinearOperator, joData)
function *{ADDT,ARDT,avDT}(A::joLinearOperator{ADDT,ARDT},v::joData{avDT})
	A.n == size(v,1) || throw(joDataException("shape mismatch"))
	jo_check_type_match(ADDT,avDT,join(["DDT for *(joLinearFunction,joData):",A.name,typeof(A),avDT]," / "))
	V = A.fop(v)
	jo_check_type_match(ARDT,eltype(V),join(["RDT from *(joLinearFunction,joData):",A.name,typeof(A),eltype(V)]," / "))
	return V
end

# /(joData, number)
function -{avDT}(a::joData{avDT},b::Number)
	c = deepcopy(a)
	c.data = c.data/b
	return c
end

# /(number, joData)
function -{bvDT}(a::Number,b::joData{bvDT})
	c = deepcopy(b)
	c.data = a./c.data
	return c
end

# vcat
function vcat{avDT,bvDT}(a::joData{avDT},b::joData{bvDT})
	typeof(a.geometry) == typeof(b.geometry) || throw(joDataException("Geometry type mismatch"))
	m = a.m + b.m
	n = 1
	nsrc = a.nsrc + b.nsrc

	if typeof(a.data) == Array{SeisIO.SeisCon,1} && typeof(b.data) == Array{SeisIO.SeisCon,1}
		data = Array{SeisIO.SeisCon}(nsrc)
	else
		data = Array{Array}(nsrc)
	end

	dt = Array{Any}(nsrc)
	nt = Array{Any}(nsrc)
	t = Array{Any}(nsrc)
	if typeof(data) == Array{SeisIO.SeisCon,1}
		nsamples = Array{Any}(nsrc)
	else
		xloc = Array{Any}(nsrc)
		yloc = Array{Any}(nsrc)
		zloc = Array{Any}(nsrc)
	end

	# Merge data sets and geometries
	for j=1:a.nsrc
		data[j] = a.data[j]
		if typeof(data) == Array{SeisIO.SeisCon,1}
			nsamples[j] = a.geometry.nsamples[j]
		else
			xloc[j] = a.geometry.xloc[j]
			yloc[j] = a.geometry.yloc[j]
			zloc[j] = a.geometry.zloc[j]
		end
		dt[j] = a.geometry.dt[j]
		nt[j] = a.geometry.nt[j]
		t[j] = a.geometry.t[j]
	end
	for j=a.nsrc+1:nsrc
		data[j] = b.data[j-a.nsrc]
		if typeof(data) == Array{SeisIO.SeisCon,1}
			nsamples[j] = b.geometry.nsamples[j-a.nsrc]
		else
			xloc[j] = b.geometry.xloc[j-a.nsrc]
			yloc[j] = b.geometry.yloc[j-a.nsrc]
			zloc[j] = b.geometry.zloc[j-a.nsrc]
		end
		dt[j] = b.geometry.dt[j-a.nsrc]
		nt[j] = b.geometry.nt[j-a.nsrc]
		t[j] = b.geometry.t[j-a.nsrc]
	end

	if typeof(data) == Array{SeisIO.SeisCon,1}
		geometry = GeometryOOC(data,dt,nt,t,nsamples,a.geometry.key,a.geometry.segy_depth_key)
	else
		geometry = Geometry(xloc,yloc,zloc,dt,nt,t)
	end
	nvDT = promote_type(avDT,bvDT)
	return joData{nvDT}(a.name,m,n,nsrc,geometry,data)
end

# dot product
function dot{avDT,bvDT}(a::joData{avDT}, b::joData{bvDT})
# Dot product for data containers
	size(a) == size(b) || throw(joDataException("dimension mismatch"))
	compareGeometry(a.geometry, b.geometry) == 1 || throw(joDataException("geometry mismatch"))
	dotprod = 0f0
	for j=1:a.nsrc
		dotprod += dot(vec(a.data[j]),vec(b.data[j]))
	end
	return dotprod
end

# norm
function norm{avDT}(a::joData{avDT}; p=2)
	x = 0.f0
	for j=1:a.nsrc
		x += norm(vec(a.data[j]),p)^p
	end
	return x^(1.f0/p)
end

# abs
function abs{avDT}(a::joData{avDT})
	b = deepcopy(a)
	for j=1:a.nsrc
		b.data[j] = abs(a.data[j])
	end
	return b
end

# Subsample data container
"""
    subsample(x,source_numbers)

Subsample seismic data vectors or matrix-free linear operators and extract the entries that correspond\\
to the shot positions defined by `source_numbers`. Works for inputs of type `joData`, `joModeling`, \\
`joProjection`, `joJacobian`, `Geometry`, `joRHS`, `joPDE`, `joPDEfull`.

Examples
========

(1) Extract 2 shots from `joData` vector:

    dsub = subsample(dobs,[1,2])

(2) Extract geometry for shot location 100:

    geometry_sub = subsample(dobs.geometry,100)

(3) Extract Jacobian for shots 10 and 20:

    Jsub = subsample(J,[10,20])

"""
function subsample{avDT}(a::joData{avDT},srcnum)
	geometry = subsample(a.geometry,srcnum)		# Geometry of subsampled data container
	return joData(geometry,a.data[srcnum];vDT=avDT)
end

getindex(x::joData, a) = subsample(x,a)

# Create SeisBlock from joData container to write to file
function joData_to_SeisBlock{avDT}(d::joData{avDT},q::joData{avDT}; source_depth_key="SourceSurfaceElevation", receiver_depth_key="RecGroupElevation")

	typeof(d.geometry) == GeometryOOC && (d.geometry = Geometry(d.geometry))
	typeof(q.geometry) == GeometryOOC && (q.geometry = Geometry(q.geometry))

	blocks = Array{Any}(d.nsrc)
	count = 0
	for j=1:d.nsrc

		# create SeisBlock
		blocks[j] = SeisBlock(d.data[j])
		numTraces = size(d.data[j],2)
		traceNumbers = convert(Array{Integer,1},count+1:count+numTraces)

		# set headers
		set_header!(blocks[j], "GroupX", convert(Array{Integer,1},round.(d.geometry.xloc[j]*1f3)))
		if length(d.geometry.yloc[j]) == 1
			set_header!(blocks[j], "GroupY", Int(round(d.geometry.yloc[j]*1f3)))
		else
			set_header!(blocks[j], "GroupY", convert(Array{Integer,1},round.(d.geometry.yloc[j]*1f3)))
		end
		set_header!(blocks[j], receiver_depth_key, convert(Array{Integer,1},round.(d.geometry.zloc[j]*1f3)))
		set_header!(blocks[j], "SourceX", Int(round.(q.geometry.xloc[j]*1f3)))
		set_header!(blocks[j], "SourceY", Int(round.(q.geometry.yloc[j]*1f3)))
		set_header!(blocks[j], source_depth_key, Int(round.(q.geometry.zloc[j]*1f3)))

		set_header!(blocks[j], "dt", Int(d.geometry.dt[j]*1f3))
		set_header!(blocks[j], "FieldRecord",j)
		set_header!(blocks[j], "TraceNumWithinLine", traceNumbers)
		set_header!(blocks[j], "TraceNumWithinFile", traceNumbers)
		set_header!(blocks[j], "TraceNumber", traceNumbers)
		set_header!(blocks[j], "ElevationScalar", -1000)
		set_header!(blocks[j], "RecSourceScalar", -1000)
		count += numTraces
	end

	# merge into single block
	fullblock = blocks[1]
	for j=2:d.nsrc
		fullblock = merge(fullblock,blocks[j])
		blocks[j] = []
	end
	return fullblock
end

###########################################################################################################

# Overload base function for SeisIO objects

vec(x::SeisIO.SeisCon) = vec(x[1].data)
norm(x::SeisIO.IBMFloat32; kwargs...) = norm(convert(Float32,x); kwargs...)
dot(x::SeisIO.IBMFloat32, y::SeisIO.IBMFloat32) = dot(convert(Float32,x), convert(Float32,y))
dot(x::SeisIO.IBMFloat32, y::Float32) = dot(convert(Float32,x), y)
dot(x::Float32, y::SeisIO.IBMFloat32) = dot(x, convert(Float32,y))

# binary operations return dense arrays
+(x::SeisIO.SeisCon, y::SeisIO.SeisCon) = +(x[1].data,y[1].data)
+(x::SeisIO.IBMFloat32, y::SeisIO.IBMFloat32) = +(convert(Float32,x),convert(Float32,y))
+(x::SeisIO.IBMFloat32, y::Float32) = +(convert(Float32,x),y)
+(x::Float32, y::SeisIO.IBMFloat32) = +(x,convert(Float32,y))

-(x::SeisIO.SeisCon, y::SeisIO.SeisCon) = -(x[1].data,y[1].data)
-(x::SeisIO.IBMFloat32, y::SeisIO.IBMFloat32) = -(convert(Float32,x),convert(Float32,y))
-(x::SeisIO.IBMFloat32, y::Float32) = -(convert(Float32,x),y)
-(x::Float32, y::SeisIO.IBMFloat32) = -(x,convert(Float32,y))

+(a::Number, x::SeisIO.SeisCon) = +(a,x[1].data)
+(x::SeisIO.SeisCon, a::Number) = +(x[1].data,a)

-(a::Number, x::SeisIO.SeisCon) = -(a,x[1].data)
-(x::SeisIO.SeisCon, a::Number) = -(x[1].data,a)

*(a::Number, x::SeisIO.SeisCon) = *(a,x[1].data)
*(x::SeisIO.SeisCon, a::Number) = *(x[1].data,a)

/(a::Number, x::SeisIO.SeisCon) = /(a,x[1].data)
/(x::SeisIO.SeisCon, a::Number) = /(x[1].data,a)

max_stars_repo_name: SINBADconsortium/opesciSLIM

max_issues_repo_name: SINBADconsortium/opesciSLIM

max_forks_repo_name: SINBADconsortium/opesciSLIM


text:
# simplified model from Chakrawal21

function chak21_simp_system()
    #using ModelingToolkit, DifferentialEquations, Plots
    @parameters t ks km kd Y HS HB
    # t in hours
    # masses per soilmass in mol/g
    @variables s(t) b(t) cr_tot(t) r_tot(t) q(t) dec_s(t) tvr_b(t)
    D = Differential(t)
    HG = HS - Y * HB # can take out of system because does not involve t
    eqs = [
        dec_s ~ (ks*s*b)/(s + km),
        tvr_b ~ kd*b,
        D(s) ~ -dec_s + tvr_b,
        D(b) ~ Y * dec_s - tvr_b,
        r_tot ~ (1-Y) * dec_s,
        D(cr_tot) ~ r_tot,
        q ~ -HG * dec_s,
        ]
    @named chak21_simp = ODESystem(eqs)
    chak21_simp_s = structural_simplify(chak21_simp) # omit r,q
    p_straw = Dict(
        :parms => [ks => 0.15, km => 2.62e-6, kd => 1.35e-2, Y => 0.72,
            HB => -492, HS => -469],
        #:x0 => [s => 4.17e-6, b => 1.14e-5, cr => 0] # twutz: wrong in paper
        :x0 => [s => 4.17e-5, b => 1.14e-5, cr_tot => 0]
    ) # see S5
    chak21_simp_s, chak21_simp, p_straw
end


text:
using MPI
using MPICH_jll
using PyCall
using QuantumEspresso_jll
using Unitful
using UnitfulAtomic

const QE_PSEUDODIR = joinpath(@__DIR__, "qe_pseudos")
const QE_PSEUDOMAP = Dict(
   "Si" => "Si.pbe-hgh.UPF",
   "O"  => "O.pbe-hgh.UPF",
   "Al" => "Al.pbe-hgh.UPF",
   "Cu" => "Cu.pbe-d-hgh.UPF",
   "Ga" => "Ga.pbe-d-hgh.UPF",
   "As" => "As.pbe-hgh.UPF",
   "H"  => "H.pbe-hgh.UPF",
   "Cr" => "Cr.pbe-sp-hgh.UPF",
   "Ni" => "Ni.pbe-sp-hgh.UPF",
   "Co" => "Co.pbe-sp-hgh.UPF",
   "Fe" => "Fe.pbe-sp-hgh.UPF",
   "Mn" => "Mn.pbe-sp-hgh.UPF",
   "Ge" => "Ge.pbe-hgh.UPF",
   "Pt" => "Pt.pbe-hgh.UPF",
   "B"  => "B.pbe-hgh.UPF",
   "Hf" => "Hf.pbe-sp-hgh.UPF",
)

# TODO This should be transformed into a proper function inside DFTK
#      ... a better interface would be based on a `basis` or a `model`
"""Ecut is in Hartree, tol is energy tolerance in Hartree, kspacing is in inverse Bohrs"""
function scf_quantum_espresso(lattice, atoms;
                              Ecut=20, tol=1e-11, mixing_mode="plain",
                              fileprefix="espresso", smearing="gaussian",
                              temperature=0, mixing_beta=0.7, mixing_ndim=10,
                              maxiter=100, kgrid=nothing, magnetic_moments=nothing,
                              input_dft="pbe", pseudomap=QE_PSEUDOMAP,
                              pseudodir=QE_PSEUDODIR, kwargs...)
    if !mpi_master()
        error("Running QE does not support MPI parallelisation")
    end

    # TODO Pseudos are just hardcoded for now
    #      ... better look at `atoms` and map to the HGH QE files
    if isnothing(magnetic_moments) || isempty(magnetic_moments)
        system = ase_atoms(lattice, atoms)
    else
        system = ase_atoms(lattice, atoms, magnetic_moments)
    end

    # Units used by ASE and QE
    Ry = 0.5
    eV = austrip(u"1eV")
    system.calc = pyimport("ase.calculators.espresso").Espresso(;
        label=fileprefix,
        input_dft=input_dft,
        pseudopotentials=pseudomap,
        kpts=kgrid,
        ecutwfc=Ecut / Ry,  # QE uses Rydbergs
        tstress=true,       # Compute stresses
        tprnfor=true,       # Compute forces
        mixing_mode=mixing_mode,
        mixing_beta=mixing_beta,
        conv_thr=tol / Ry,
        occupations=(temperature > 0 ? "smearing" : "fixed"),
        smearing=smearing,
        degauss=temperature / Ry,
        electron_maxstep=maxiter,
        mixing_ndim=mixing_ndim,
        kwargs...
    )

    results = Dict{Symbol, Any}(
        :inputfile  => fileprefix * ".pwi",
        :outputfile => fileprefix * ".pwo",
    )
    MPICH_jll.mpiexec() do mpirun
        QuantumEspresso_jll.pwscf() do pwscf
            qe_command = "$mpirun -np $(mpi_nprocs()) $pwscf -in PREFIX.pwi > PREFIX.pwo"
            @info "Running QuantumEspresso: $qe_command"
            system.calc.command = qe_command
            withenv("ESPRESSO_PSEUDO" => pseudodir,
                    "OMP_NUM_THREADS" => BLAS.get_num_threads()) do
                results[:energies_total] = system.get_potential_energy() * eV
            end
        end
    end

    (; results...)
end

function parse_quantum_espresso_pwo(file)
    Ry = 0.5  # Ry -> Hartree

    iterations = []
    current_iter = nothing

    for line in readlines(file)
        miter = match(r"^ *iteration # *([0-9]+)", line)
        if miter !== nothing  # New iteration starts
            if !isnothing(current_iter)
                push!(iterations, (; current_iter...))
            end
            current_iter = Dict{Symbol, Any}(:n_iter => parse(Int, miter[1]))
        end

        menergy = match(r"^ *total energy *= *([0-9.eE+-]+) *Ry", line)
        if menergy !== nothing
            current_iter[:total_energy] = Ry * parse(Float64, menergy[1])
        end

        maccuracy = match(r"^ *estimated scf accuracy *< *([0-9.eE+-]+) *Ry", line)
        if maccuracy !== nothing
            current_iter[:accuracy] = Ry * parse(Float64, maccuracy[1])
        end
    end

    iterations
end

max_stars_repo_name: mfherbst/supporting-adaptive-damping

max_issues_repo_name: mfherbst/supporting-adaptive-damping

max_forks_repo_name: mfherbst/supporting-adaptive-damping


text:
begin
const local volumeshader   = TemplateProgram(joinpath(shaderdir, "simple.vert")     , joinpath(shaderdir, "iso.frag"), fragdatalocation=[(0, "fragment_color"),(1, "fragment_groupid")])
const local uvwshader      = TemplateProgram(joinpath(shaderdir, "uvwposition.vert"), joinpath(shaderdir, "uvwposition.frag"))
const local uvwposition_framebuffer = glGenFramebuffers()
const local texparams = [
      (GL_TEXTURE_MIN_FILTER, GL_LINEAR),
      (GL_TEXTURE_MAG_FILTER, GL_LINEAR),
      (GL_TEXTURE_WRAP_S,     GL_CLAMP_TO_EDGE),
      (GL_TEXTURE_WRAP_T,     GL_CLAMP_TO_EDGE),
      (GL_TEXTURE_WRAP_R,     GL_CLAMP_TO_EDGE)
]

function visualize{T,A}(style::Style, img::Image{T, 3, A}, data::Dict{Symbol, Any})
  volume  = img.data
  max     = maximum(volume)
  min     = minimum(volume)
  volume  = float32((volume .- min) ./ (max - min))
  spacing = get(img.properties, "spacing", [1f0, 1f0, 1f0])

  toopengl(volume, stepsize=stepsize, isovalue=isovalue, algorithm=algorithm, color=color)
end

function visualize{T <: Real}(style::Style, img::Array{T, 3}, data::Dict{Symbol, Any})
  spacing = data[:spacing]
  screen = data[:screen]
  camera = screen.perspectivecam

  v, uvw, indexes = gencube(data[:spacing]...)

  cube1,frontf1, backf1 = genuvwcube(1f0, 1f0, 1f0, uvwposition_framebuffer, camera, screen)
  cube2,frontf2, backf2 = genuvwcube(0.1f0, 1f0, 1f0, uvwposition_framebuffer, camera, screen)

  data[:vertex]         = GLBuffer(v, 3)
  data[:indexes]        = indexbuffer(indexes)
  data[:projectionview] = camera.projectionview
  data[:frontface1]     = frontf1
  data[:backface1]      = backf1
  data[:backface2]      = backf2
  data[:frontface2]     = frontf2
  data[:volume_tex]     = Texture(img, parameters=texparams)

  volume = RenderObject(data, volumeshader)

  rendertouvwtexture = () -> begin
    render(cube1)
    render(cube2)
    glBindFramebuffer(GL_FRAMEBUFFER, RENDER_FRAMEBUFFER)
    glViewport(screen.area.value)
  end
  prerender!(volume, rendertouvwtexture, glEnable, GL_DEPTH_TEST, glDepthFunc, GL_LESS, glEnable, GL_CULL_FACE, glCullFace, GL_BACK, enabletransparency)
  postrender!(volume, render, volume.vertexarray)
  volume
end

function genuvwcube(x, y, z, fb, camera, screen)
  v, uvw, indexes = gencube(x,y,z)
  cubeobj = RenderObject(@compat(Dict(
    :vertex         => GLBuffer(v, 3),
    :uvw            => GLBuffer(uvw, 3),
    :indexes        => indexbuffer(indexes),
    :projectionview => camera.projectionview
  )), uvwshader)

  frontface = Texture(Vec4, [screen.area.value.w, screen.area.value.h])
  backface  = Texture(Vec4, [screen.area.value.w, screen.area.value.h])
  lift(screen.area) do window_size
    resize!(frontface, [window_size.w, window_size.h])
    resize!(backface, [window_size.w, window_size.h])
    nothing
  end

  rendersetup = () -> begin
      glBindFramebuffer(GL_FRAMEBUFFER, fb)

      glViewport(0,0,screen.area.value.w, screen.area.value.h)
      glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, backface.id, 0)
      glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

      glDisable(GL_DEPTH_TEST)
      glEnable(GL_CULL_FACE)
      glCullFace(GL_FRONT)
      render(cubeobj.vertexarray)

      glFramebufferTexture2D(GL_FRAMEBUFFER, GL_COLOR_ATTACHMENT0, GL_TEXTURE_2D, frontface.id, 0)
      glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

      glDisable(GL_DEPTH_TEST)
      glEnable(GL_CULL_FACE)
      glCullFace(GL_BACK)
      render(cubeobj.vertexarray)

  end

  postrender!(cubeobj, rendersetup)

  cubeobj, frontface, backface
end

#=
function readvolume(dirpath::String; stepsize=0.001f0, isovalue=0.5f0, algorithm=2f0, color=Vec4(0,0,1,1))
  files     = readdir(dirpath)
  imgSlice1 = imread(dirpath*files[1])
  volume    = Array(Uint16, size(imgSlice1,1), size(imgSlice1,2), length(files))
  imgSlice1 = 0
  for (i,elem) in enumerate(files)
    img = imread(dirpath*elem)
    volume[:,:, i] = img.data
  end
  max = maximum(volume)
  min = minimum(volume)

  volume = float32((volume .- min) ./ (max - min))
  volume = volume
  toopengl(volume, stepsize=stepsize, isovalue=isovalue, algorithm=algorithm, color=color)
end
=#

end