text:
"""
    b, e10 = reduce_shortest(f[, maxsignif])

Reduce to shortest decimal representation where `abs(f) == b * 10^e10` and `b` is an
integer. If a `maxsignif` argument is provided, then `b < maxsignif`.
"""
@inline function reduce_shortest(f::T, maxsignif=nothing) where {T}
    U = uinttype(T)
    uf = reinterpret(U, f)
    m = uf & significand_mask(T)
    e = ((uf & exponent_mask(T)) >> significand_bits(T)) % Int

    ## Step 1
    #  mf * 2^ef == f
    mf = (one(U) << significand_bits(T)) | m
    ef = e - exponent_bias(T) - significand_bits(T)
    f_isinteger = mf & ((one(U) << -ef) - one(U)) == 0

    if ef > 0 || ef < -Base.significand_bits(T) || !f_isinteger
        # fixup subnormals
        if e == 0
            ef = 1 - exponent_bias(T) - significand_bits(T)
            mf = m
        end

        ## Step 2
        #  u * 2^e2 == (f + prevfloat(f))/2
        #  v * 2^e2 == f
        #  w * 2^e2 == (f + nextfloat(f))/2
        e2 = ef - 2
        mf_iseven = iseven(mf) # trailing bit of significand is zero

        v = U(4) * mf
        w = v + U(2)
        u_shift_half = m == 0 && e > 1 # if first element of binade, other than first normal one
        u = v - U(2) + u_shift_half

        ## Step 3
        #  a == floor(u * 2^e2 / 10^e10), exact if a_allzero
        #  b == floor(v * 2^e2 / 10^e10), exact if b_allzero
        #  c == floor(w * 2^e2 / 10^e10)
        a_allzero = false
        b_allzero = false
        b_lastdigit = 0x00
        if e2 >= 0
            q = log10pow2(e2) - (T == Float64 ? (e2 > 3) : 0)
            e10 = q
            k = pow5_inv_bitcount(T) + pow5bits(q) - 1
            i = -e2 + q + k
            a, b, c = mulshiftinvsplit(T, u, v, w, q, i)
            if T == Float32 || T == Float16
                if q != 0 && div(c - 1, 10) <= div(a, 10)
                    l = pow5_inv_bitcount(T) + pow5bits(q - 1) - 1
                    mul = pow5invsplit_lookup(T, q-1)
                    b_lastdigit = (mulshift(v, mul, -e2 + q - 1 + l) % 10) % UInt8
                end
            end
            if q <= qinvbound(T)
                if ((v % UInt32) - 5 * div(v, 5)) == 0
                    b_allzero = pow5(v, q)
                elseif mf_iseven
                    a_allzero = pow5(u, q)
                else
                    c -= pow5(w, q)
                end
            end
        else
            q = log10pow5(-e2) - (T == Float64 ? (-e2 > 1) : 0)
            e10 = q + e2
            i = -e2 - q
            k = pow5bits(i) - pow5_bitcount(T)
            j = q - k
            a, b, c = mulshiftsplit(T, u, v, w, i, j)
            if T == Float32 || T == Float16
                if q != 0 && div(c - 1, 10) <= div(a, 10)
                    j = q - 1 - (pow5bits(i + 1) - pow5_bitcount(T))
                    mul = pow5split_lookup(T, i+1)
                    b_lastdigit = (mulshift(v, mul, j) % 10) % UInt8
                end
            end
            if q <= 1
                b_allzero = true
                if mf_iseven
                    a_allzero = !u_shift_half
                else
                    c -= 1
                end
            elseif q < qbound(T)
                b_allzero = pow2(v, q - (T != Float64))
            end
        end

        ## Step 4: reduction
        if a_allzero || b_allzero
            # a) slow loop
            while true
                c_div10 = div(c, 10)
                a_div10 = div(a, 10)
                if c_div10 <= a_div10
                    break
                end
                a_mod10 = (a % UInt32) - UInt32(10) * (a_div10 % UInt32)
                b_div10 = div(b, 10)
                b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
                a_allzero &= a_mod10 == 0
                b_allzero &= b_lastdigit == 0
                b_lastdigit = b_mod10 % UInt8
                b = b_div10
                c = c_div10
                a = a_div10
                e10 += 1
            end
            if a_allzero
                while true
                    a_div10 = div(a, 10)
                    a_mod10 = (a % UInt32) - UInt32(10) * (a_div10 % UInt32)
                    if a_mod10 != 0 && (maxsignif === nothing || b < maxsignif)
                        break
                    end
                    c_div10 = div(c, 10)
                    b_div10 = div(b, 10)
                    b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
                    b_allzero &= b_lastdigit == 0
                    b_lastdigit = b_mod10 % UInt8
                    b = b_div10
                    c = c_div10
                    a = a_div10
                    e10 += 1
                end
            end
            if b_allzero && b_lastdigit == 5 && iseven(b)
                b_lastdigit = UInt8(4)
            end
            roundup = (b == a && (!mf_iseven || !a_allzero)) || b_lastdigit >= 5
        else
            # b) specialized for common case (99% Float64, 96% Float32)
            roundup = b_lastdigit >= 5
            c_div100 = div(c, 100)
            a_div100 = div(a, 100)
            if c_div100 > a_div100
                b_div100 = div(b, 100)
                b_mod100 = (b % UInt32) - UInt32(100) * (b_div100 % UInt32)
                roundup = b_mod100 >= 50
                b = b_div100
                c = c_div100
                a = a_div100
                e10 += 2
            end
            while true
                c_div10 = div(c, 10)
                a_div10 = div(a, 10)
                if c_div10 <= a_div10
                    break
                end
                b_div10 = div(b, 10)
                b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
                roundup = b_mod10 >= 5
                b = b_div10
                c = c_div10
                a = a_div10
                e10 += 1
            end
            roundup = (b == a || roundup)
        end
        if maxsignif !== nothing && b > maxsignif
            # reduce to max significant digits
            while true
                b_div10 = div(b, 10)
                b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
                if b <= maxsignif
                    break
                end
                b = b_div10
                roundup = (b_allzero && iseven(b)) ? b_mod10 > 5 : b_mod10 >= 5
                b_allzero &= b_mod10 == 0
                e10 += 1
            end
            b = b + roundup

            # remove trailing zeros
            while true
                b_div10 = div(b, 10)
                b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
                if b_mod10 != 0
                    break
                end
                b = b_div10
                e10 += 1
            end
        else
            b = b + roundup
        end
    else
        # c) specialized f an integer < 2^53
        b = mf >> -ef
        e10 = 0

        if maxsignif !== nothing && b > maxsignif
            b_allzero = true
            # reduce to max significant digits
            while true
                b_div10 = div(b, 10)
                b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
                if b <= maxsignif
                    break
                end
                b = b_div10
                roundup = (b_allzero && iseven(b)) ? b_mod10 > 5 : b_mod10 >= 5
                b_allzero &= b_mod10 == 0
                e10 += 1
            end
            b = b + roundup
        end
        while true
            b_div10 = div(b, 10)
            b_mod10 = (b % UInt32) - UInt32(10) * (b_div10 % UInt32)
            if b_mod10 != 0
                break
            end
            b = b_div10
            e10 += 1
        end
    end
    return b, e10
end

@inline function writeshortest(buf::Vector{UInt8}, pos, x::T,
                               plus=false, space=false, hash=true,
                               precision=-1, expchar=UInt8('e'), padexp=false, decchar=UInt8('.'),
                               typed=false, compact=false) where {T}
    @assert 0 < pos <= length(buf)
    neg = signbit(x)
    # special cases
    if x == 0
        if neg
            buf[pos] = UInt8('-')
            pos += 1
        elseif plus
            buf[pos] = UInt8('+')
            pos += 1
        elseif space
            buf[pos] = UInt8(' ')
            pos += 1
        end
        buf[pos] = UInt8('0')
        pos += 1
        if hash
            buf[pos] = decchar
            pos += 1
        end
        if precision == -1
            buf[pos] = UInt8('0')
            pos += 1
            if typed && x isa Float32
                buf[pos] = UInt8('f')
                buf[pos + 1] = UInt8('0')
                pos += 2
            end
            return pos
        end
        while precision > 1
            buf[pos] = UInt8('0')
            pos += 1
            precision -= 1
        end
        if typed && x isa Float32
            buf[pos] = UInt8('f')
            buf[pos + 1] = UInt8('0')
            pos += 2
        end
        return pos
    elseif isnan(x)
        buf[pos] = UInt8('N')
        buf[pos + 1] = UInt8('a')
        buf[pos + 2] = UInt8('N')
        if typed
            if x isa Float32
                buf[pos + 3] = UInt8('3')
                buf[pos + 4] = UInt8('2')
            elseif x isa Float16
                buf[pos + 3] = UInt8('1')
                buf[pos + 4] = UInt8('6')
            end
        end
        return pos + 3 + (typed && x isa Union{Float32, Float16} ? 2 : 0)
    elseif !isfinite(x)
        if neg
            buf[pos] = UInt8('-')
        end
        buf[pos + neg] = UInt8('I')
        buf[pos + neg + 1] = UInt8('n')
        buf[pos + neg + 2] = UInt8('f')
        if typed
            if x isa Float32
                buf[pos + neg + 3] = UInt8('3')
                buf[pos + neg + 4] = UInt8('2')
            elseif x isa Float16
                buf[pos + neg + 3] = UInt8('1')
                buf[pos + neg + 4] = UInt8('6')
            end
        end
        return pos + neg + 3 + (typed && x isa Union{Float32, Float16} ? 2 : 0)
    end

    output, nexp = reduce_shortest(x, compact ? 999_999 : nothing)

    if typed && x isa Float16
        buf[pos] = UInt8('F')
        buf[pos + 1] = UInt8('l')
        buf[pos + 2] = UInt8('o')
        buf[pos + 3] = UInt8('a')
        buf[pos + 4] = UInt8('t')
        buf[pos + 5] = UInt8('1')
        buf[pos + 6] = UInt8('6')
        buf[pos + 7] = UInt8('(')
        pos += 8
    end
    if neg
        buf[pos] = UInt8('-')
        pos += 1
    elseif plus
        buf[pos] = UInt8('+')
        pos += 1
    elseif space
        buf[pos] = UInt8(' ')
        pos += 1
    end

    olength = decimallength(output)
    exp_form = true
    pt = nexp + olength
    if -4 < pt <= (precision == -1 ? (T == Float16 ? 3 : 6) : precision)
        exp_form = false
        if pt <= 0
            buf[pos] = UInt8('0')
            pos += 1
            buf[pos] = decchar
            pos += 1
            for _ = 1:abs(pt)
                buf[pos] = UInt8('0')
                pos += 1
            end
            # elseif pt >= olength
            # nothing to do at this point
            # else
            # nothing to do at this point
        end
    else
        pos += 1
    end
    i = 0
    ptr = pointer(buf)
    ptr2 = pointer(DIGIT_TABLE)
    if (output >> 32) != 0
        q = output ÷ 100000000
        output2 = (output % UInt32) - UInt32(100000000) * (q % UInt32)
        output = q

        c = output2 % UInt32(10000)
        output2 = div(output2, UInt32(10000))
        d = output2 % UInt32(10000)
        c0 = (c % 100) << 1
        c1 = (c ÷ 100) << 1
        d0 = (d % 100) << 1
        d1 = (d ÷ 100) << 1
        memcpy(ptr, pos + olength - 2, ptr2, c0 + 1, 2)
        memcpy(ptr, pos + olength - 4, ptr2, c1 + 1, 2)
        memcpy(ptr, pos + olength - 6, ptr2, d0 + 1, 2)
        memcpy(ptr, pos + olength - 8, ptr2, d1 + 1, 2)
        i += 8
    end
    output2 = output % UInt32
    while output2 >= 10000
        c = output2 % UInt32(10000)
        output2 = div(output2, UInt32(10000))
        c0 = (c % 100) << 1
        c1 = (c ÷ 100) << 1
        memcpy(ptr, pos + olength - i - 2, ptr2, c0 + 1, 2)
        memcpy(ptr, pos + olength - i - 4, ptr2, c1 + 1, 2)
        i += 4
    end
    if output2 >= 100
        c = (output2 % UInt32(100)) << 1
        output2 = div(output2, UInt32(100))
        memcpy(ptr, pos + olength - i - 2, ptr2, c + 1, 2)
        i += 2
    end
    if output2 >= 10
        c = output2 << 1
        buf[pos + 1] = DIGIT_TABLE[c + 2]
        buf[pos - exp_form] = DIGIT_TABLE[c + 1]
    else
        buf[pos - exp_form] = UInt8('0') + (output2 % UInt8)
    end

    if !exp_form
        if pt <= 0
            pos += olength
            precision -= olength
            while hash && precision > 0
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
        elseif pt >= olength
            pos += olength
            precision -= olength
            for _ = 1:nexp
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
            if hash
                buf[pos] = decchar
                pos += 1
                if precision < 0
                    buf[pos] = UInt8('0')
                    pos += 1
                end
                while precision > 0
                    buf[pos] = UInt8('0')
                    pos += 1
                    precision -= 1
                end
            end
        else
            pointoff = olength - abs(nexp)
            memmove(ptr, pos + pointoff + 1, ptr, pos + pointoff, olength - pointoff + 1)
            buf[pos + pointoff] = decchar
            pos += olength + 1
            precision -= olength
            while hash && precision > 0
                buf[pos] = UInt8('0')
                pos += 1
                precision -= 1
            end
        end
        if typed && x isa Float32
            buf[pos] = UInt8('f')
            buf[pos + 1] = UInt8('0')
            pos += 2
        end
    else
        if olength > 1 || hash
            buf[pos] = decchar
            pos += olength
            precision -= olength
        end
        if hash && olength == 1
            buf[pos] = UInt8('0')
            pos += 1
        end
        while hash && precision > 0
            buf[pos] = UInt8('0')
            pos += 1
            precision -= 1
        end

        buf[pos] = expchar
        pos += 1
        exp2 = nexp + olength - 1
        if exp2 < 0
            buf[pos] = UInt8('-')
            pos += 1
            exp2 = -exp2
        elseif padexp
            buf[pos] = UInt8('+')
            pos += 1
        end

        if exp2 >= 100
            c = exp2 % 10
            memcpy(ptr, pos, ptr2, 2 * div(exp2, 10) + 1, 2)
            buf[pos + 2] = UInt8('0') + (c % UInt8)
            pos += 3
        elseif exp2 >= 10
            memcpy(ptr, pos, ptr2, 2 * exp2 + 1, 2)
            pos += 2
        else
            if padexp
                buf[pos] = UInt8('0')
                pos += 1
            end
            buf[pos] = UInt8('0') + (exp2 % UInt8)
            pos += 1
        end
    end
    if typed && x isa Float16
        buf[pos] = UInt8(')')
        pos += 1
    end

    return pos
end

max_stars_repo_name: TimoLarson/julia

max_issues_repo_name: TimoLarson/julia

max_forks_repo_name: TimoLarson/julia


text:
using DisplayAs
using Test

@testset "DisplayAs.jl" begin
    @test showable("text/plain", DisplayAs.Text(nothing))
    @test showable("text/html", DisplayAs.HTML(nothing))
    @test showable("text/markdown", DisplayAs.MD(nothing))
    @test showable("image/png", DisplayAs.PNG(nothing))
end

using Aqua
Aqua.test_all(DisplayAs)

ext: jl


text:
"""
    AbstractSampler

Abstract supertype for sampler objects.
"""
abstract type AbstractSampler{S <: AbstractScenario} end
"""
    sample(sampler::AbstractSampler{S})

Sample a scenario of type `S` using `sampler`.
"""
function sample(sampler::AbstractSampler)
    return sampler()
end
function Base.show(io::IO, sampler::AbstractSampler{S}) where S <: AbstractScenario
    print(io, "$(S.name.name) sampler")
    return io
end
function Base.show(io::IO, sampler::AbstractSampler{S}) where S <: Scenario
    print(io, "Scenario sampler")
    return io
end
"""
    sample(sampler::AbstractSampler{S}, π::AbstractSampler)

Sample a scenario of type `S` using `sampler` and set the probability of the sampled scenario to `π`.
"""
function sample(sampler::AbstractSampler, π::AbstractFloat)
    scenario = sampler()
    set_probability!(scenario, π)
    return scenario
end
"""
    Sampler

General purpose sampler object that samples Scenario.

See also: [`Scenario`](@ref), [`@sampler`](@ref)
"""
struct Sampler <: AbstractSampler{Scenario}
    sampler::Function
end

function (sampler::Sampler)()
    return sampler.sampler()
end


text:
module Sysimage
import Pkg, Libdl, PackageCompiler

function compile(sysimage_path = "JuliaSysimage.$(Libdl.dlext)")
    env_to_precompile = joinpath(@__DIR__, "..")
    precompile_execution_file = joinpath(env_to_precompile, "Precompile.jl")
    project_filename = joinpath(env_to_precompile, "Project.toml")
    project = Pkg.API.read_project(project_filename)
    used_packages = Symbol.(collect(keys(project.deps)))
    # don't need these ever after building this so no need to have them in the sysimage
    used_packages = filter(x -> x ∉ [:Libdl, :PackageCompiler, :Pkg], used_packages)
    if Libdl.dlext == "dll"
        @warn "Ignoring packages which use gl dlls on Windows as these cause build errors"
        # see https://github.com/JuliaLang/PackageCompiler.jl/issues/365
        # recently the build leaves a corrupt dll if we include these rather than the error in the issue above
        used_packages = filter(x -> x ∉ [:Makie, :ImageView], used_packages)
    end
    @info "Building a custom sysimage for OpticSim.jl."
    PackageCompiler.create_sysimage(used_packages, sysimage_path = sysimage_path, project = env_to_precompile, precompile_execution_file = precompile_execution_file)
end

end # module Sysimage


text:
#=
The main module for Finch.
=#
module Finch

export init_finch, set_language, set_custom_gen_target, dendro, set_solver, set_stepper, set_specified_steps, set_matrix_free,
        reformat_for_stepper, reformat_for_stepper_fv,
        add_mesh, output_mesh, add_boundary_ID, add_test_function,
        add_initial_condition, add_boundary_condition, add_reference_point, set_rhs, set_lhs, set_lhs_surface, set_rhs_surface,
        set_symexpressions, set_assembly_loops
export build_cache_level, build_cache, build_cache_auto
export sp_parse
export generate_code_layer, generate_code_layer_surface, generate_code_layer_fv
export Variable, add_variable, VariableTransform
export Coefficient, add_coefficient
export Parameter, add_parameter
export Indexer
export FVInfo

### Module's global variables ###
# configuration, output files, etc.
config = nothing;
prob = nothing;
project_name = "unnamedProject";
output_dir = pwd();
language = 0;
gen_framework = 0;
generate_external = false;
solver = nothing;
codegen_params = nothing;
#log
use_log = false;
log_file = "";
log_line_index = 1;
#mesh
mesh_data = nothing;    # The basic element information as read from a MSH file or generated here.
elemental_order = [];   # Determines the order of the elemental loops
needed_grid_types = [false, false, false]; # [CG, DG, FV] true if that type is needed
# FEM specific
grid_data = nothing;    # The full collection of nodes(including internal nodes) and other mesh info in the actual DOF ordering.
dg_grid_data = nothing; # This DG version is only made if using mixed CG/DG. Otherwise it is in grid_data.
geo_factors = nothing;  # Geometric factors
refel = nothing;        # Reference element(s?)
# FVM specific
fv_info = nothing;      # Finite volume info
fv_refel = nothing;     # Reference element for FV
parent_maps = nothing;  # For high order FV
fv_grid = nothing;      # Similar to grid_data, but only first order CG elements, and includes ghost info when partitioned.
fv_geo_factors = nothing;# Geometric factors for FV

#problem variables
var_count = 0;
variables = [];
coefficients = [];
parameters = [];
test_functions = [];
indexers = [];
variable_transforms = [];
#generated functions
genfunc_count = 0;
genfunctions = [];
callback_functions = [];
#rhs
linears = [];
face_linears = [];
#lhs
bilinears = [];
face_bilinears = [];
#assembly loop functions
assembly_loops = [];
#symbolic layer
symexpressions = [[],[],[],[]];
# string versions of generated code for printing
code_strings = [[],[],[],[],[]];
#time stepper
time_stepper = nothing;
specified_dt = 0;
specified_Nsteps = 0;
use_specified_steps = false;

use_cachesim = false;

#handles for custom code gen functions
custom_gen_funcs = [];

###############################################################
include("finch_includes.jl");
include("macros.jl");
include("finch_interface.jl"); # included here after globals are defined
###############################################################

config = Finch_config(); # These need to be initialized here
prob = Finch_prob();

# This has already been initialized on loading, so this is more of a reset function
function init_finch(name="unnamedProject")
    global config = Finch_config();
    global prob = Finch_prob();
    global project_name = name;
    global language = JULIA;
    global framework = 0;
    global generate_external = false;
    global codegen_params = nothing;
    global log_file = "";
    global use_log = false;
    global log_line_index = 1;
    global mesh_data = nothing;
    global needed_grid_types = [false, false, false];
    global grid_data = nothing;
    global dg_grid_data = nothing;
    global geo_factors = nothing;
    global fv_info = nothing;
    global refel = nothing;
    global fv_refel = nothing;
    global elemental_order = [];
    global parent_maps = nothing;
    global fv_grid = nothing;
    global fv_geo_factors = nothing;
    global var_count = 0;
    global variables = [];
    global coefficients = [];
    global parameters = [];
    global test_functions = [];
    global indexers = [];
    global variable_transforms = [];
    global genfunc_count = 0;
    global genfunctions = [];
    global callback_functions = [];
    global linears = [];
    global bilinears = [];
    global face_linears = [];
    global face_bilinears = [];
    global assembly_loops = [];
    global symexpressions = [[],[],[],[]];
    global code_strings = [[],[],[],[],[]];
    global time_stepper = nothing;
    global specified_dt = 0;
    global specified_Nsteps = 0;
    global use_specified_steps = false;
    global use_cachesim = false;

    # check for MPI and initialize
    if @isdefined(MPI)
        MPI.Init();
        config.num_procs = MPI.Comm_size(MPI.COMM_WORLD);
        if config.num_procs < 2
            # If only one process, ignore MPI.
            # MPI.Finalize();
        else
            config.use_mpi = true;
            config.proc_rank = MPI.Comm_rank(MPI.COMM_WORLD);
        end
    end

    # check for thread availability
    config.num_threads = Threads.nthreads();
    if (config.num_threads > 1 || config.num_procs > 1) && config.proc_rank == 0
        println("Initialized with "*string(config.num_procs)*" processes and "*string(config.num_threads)*" threads per proc.");
    end
end

# Sets the code generation target
function set_included_gen_target(lang, framework, dirpath, name; head="")
    #TODO
    println("Premade generation targets currently under construction. Sorry. Reverting to Julia target.");
    return;

    global language = lang;
    global gen_framework = framework;
    global output_dir = dirpath;
    global project_name = name;
    global generate_external = true;
    init_codegenerator(dirpath, name, head);

    # Need to set these three functions

    set_generation_target(get_external_language_elements, generate_external_code_layer, generate_external_files);
end

# Setting a custom target requires three functions
# 1. get_external_language_elements() - file extensions, comment chars etc.
# 2. generate_external_code_layer(var, entities, terms, lorr, vors) - Turns symbolic expressions into code
# 3. generate_external_files(var, lhs_vol, lhs_surf, rhs_vol, rhs_surf) - Writes all files based on generated code
function set_custom_gen_target(lang_elements, code_layer, file_maker, dirpath, name; head="", params=nothing)
    global language = CUSTOM_GEN_TARGET;
    global gen_framework = CUSTOM_GEN_TARGET;
    global output_dir = dirpath;
    global project_name = name;
    global generate_external = true;
    init_code_generator(dirpath, name, head);
    set_generation_target(lang_elements, code_layer, file_maker);
end

# Set parameters used by code generation. This is specific to the target.
function set_codegen_parameters(params)
    global codegen_params = params;
end

# Set the needed discretizations and backend
function set_solver(stype, backend)
    config.linalg_backend = backend;
    if typeof(stype) <: Array
        config.solver_type = MIXED;
        global solver = MixedSolver;
        for s in stype
            if s == CG
                global needed_grid_types[1] = true;
            elseif s == DG
                global needed_grid_types[2] = true;
            elseif s == FV
                global needed_grid_types[3] = true;
            end
        end
    elseif stype == DG
        config.solver_type = stype;
        global solver = DGSolver;
        global needed_grid_types[2] = true;
    elseif stype == CG
        config.solver_type = stype;
        global solver = CGSolver;
        global needed_grid_types[1] = true;
    elseif stype == FV
        config.solver_type = stype;
        global solver = FVSolver;
        global needed_grid_types[3] = true;
    end
    solver.init_solver();

    if backend == PETSC_SOLVER
        # Initialize using the first available library.
        # This should have been set up beforehand.
        if length(PETSc.petsclibs) == 0
            printerr(
"No PETSc library found. Please build PETSc first.
If you have a preferred PETSc library or the one supplied with PETSc.jl
is causing trouble, do this:
\$ export JULIA_PETSC_LIBRARY=/path/to/your/petsc_lib.so
julia> ]build PETSc", fatal=true);
        end
        petsclib = PETSc.petsclibs[1];
        PETSc.initialize(petsclib)
    end
end

# Sets the time stepper
function set_stepper(type, cfl)
    if typeof(type) <: Array
        first_stepper = Stepper(type[1], cfl);
        second_stepper = Stepper(type[2], cfl);

        global time_stepper = [first_stepper, second_stepper];
        global config.stepper = MIXED_STEPPER;
        log_entry("Set mixed time stepper to ["*type[1]*", "*type[2]*"]", 1);
    else
        global time_stepper = Stepper(type, cfl);
        global config.stepper = type;
        log_entry("Set time stepper to "*type, 1);
    end
end

# Time steps can be chosen automatically or specified. This is for manual specifying.
function set_specified_steps(dt, steps)
    global specified_dt = dt;
    global specified_Nsteps = steps;
    global use_specified_steps = true;
    prob.time_dependent = true;
    prob.end_time = dt*steps;
    log_entry("Set time stepper values to dt="*string(dt)*", Nsteps="*string(steps), 1);
end

# Uses matrix free iterative solver with the given max iterations and error tolerance.
function set_matrix_free(max, tol)
    config.linalg_matrixfree = true;
    config.linalg_matfree_max = max;
    config.linalg_matfree_tol = tol;
end

# Adds a mesh and builds the full grid and reference elements.
function add_mesh(mesh; partitions=0)
    global mesh_data = mesh;
    global refel;
    global fv_refel;
    global grid_data;
    global dg_grid_data;
    global fv_grid;
    global geo_factors;
    global fv_geo_factors;
    global fv_info;
    global elemental_order;

    # If no method has been specified, assume CG
    if !(needed_grid_types[1] || needed_grid_types[2] || needed_grid_types[3])
        set_solver(CG, DEFAULT_SOLVER);
    end

    if partitions==0
        np = config.num_procs; # By default each process gets a partition
        config.num_partitions = np;
        config.partition_index = config.proc_rank;
    else
        # If other partitioning strategies are desired.
        # More than one proc may be assigned to each mesh partition.
        # They are grouped by partition number: 0,0,0,1,1,1,...,p,p,p
        np = partitions;
        config.num_partitions = np;
        config.partition_index = Int(floor(((config.proc_rank+0.5) * np) / config.num_procs));
    end
    if config.use_mpi && np > 1
        if config.proc_rank == 0
            # Partition the mesh and send the partition info to each proc to build their subgrids.
            # For now, partition into the number of procs. This will be modified in the future.
            epart = get_element_partitions(mesh, np);

        else # other ranks recieve the partition info
            epart = fill(Cint(-1), mesh.nel);

        end
        # Broadcast partition info to all
        MPI.Bcast!(epart, 0, MPI.COMM_WORLD);

        # Each proc will build a subgrid containing only their elements and ghost neighbors.
        # Make a Grid struct for each needed type
        if needed_grid_types[1] # CG
            (refel, grid_data) = partitioned_grid_from_mesh(mesh_data, epart, grid_type=CG, order=config.basis_order_min);
        end
        if needed_grid_types[2] # DG
            if needed_grid_types[1] # CG also needed
                (refel, dg_grid_data) = partitioned_grid_from_mesh(mesh_data, epart, grid_type=DG, order=config.basis_order_min);
            else
                (refel, grid_data) = partitioned_grid_from_mesh(mesh_data, epart, grid_type=DG, order=config.basis_order_min);
            end
        end
        if needed_grid_types[3] # FV
            (fv_refel, fv_grid) = partitioned_grid_from_mesh(mesh_data, epart, grid_type=FV, order=1);
            # If only FV is used, also set grid_data and refel to this. Just in case.
            if !(needed_grid_types[1] || needed_grid_types[2])
                grid_data = fv_grid;
                refel = fv_refel;
            end
        end

    else
        # Make a Grid struct for each needed type
        if needed_grid_types[1] # CG
            (refel, grid_data) = grid_from_mesh(mesh_data, grid_type=CG, order=config.basis_order_min);
        end
        if needed_grid_types[2] # DG
            if needed_grid_types[1] # CG also needed
                (refel, dg_grid_data) = grid_from_mesh(mesh_data, grid_type=DG, order=config.basis_order_min);
            else
                (refel, grid_data) = grid_from_mesh(mesh_data, grid_type=DG, order=config.basis_order_min);
            end
        end
        if needed_grid_types[3] # FV
            (fv_refel, fv_grid) = grid_from_mesh(mesh_data, grid_type=FV, order=1);
            # If only FV is used, also set grid_data and refel to this. Just in case.
            if !(needed_grid_types[1] || needed_grid_types[2])
                grid_data = fv_grid;
                refel = fv_refel;
            end
        end
    end

    # regular parallel sided elements or simplexes have constant Jacobians, so only store one value per element.
    if ((config.geometry == SQUARE && config.mesh_type == UNIFORM_GRID) ||
        config.dimension == 1 ||
        (config.dimension == 2 && refel.Nfaces == 3) ||
        (config.dimension == 3 && refel.Nfaces == 4) )
        constantJ = true;
    else
        constantJ = false;
    end
    do_faces = false;
    do_vol = false;
    if config.solver_type == DG || config.solver_type == FV || config.solver_type == MIXED
        do_faces = true;
    end
    if config.solver_type == FV || config.solver_type == MIXED
        do_vol = true;
    end

    # FE version is always made?
    geo_factors = build_geometric_factors(refel, grid_data, do_face_detj=do_faces, do_vol_area=do_vol, constant_jacobian=constantJ);
    if needed_grid_types[3] # FV
        fv_geo_factors = build_geometric_factors(fv_refel, fv_grid, do_face_detj=do_faces, do_vol_area=do_vol, constant_jacobian=constantJ);
        fv_info = build_FV_info(fv_grid);
    end

    log_entry("Added mesh with "*string(mesh_data.nx)*" vertices and "*string(mesh_data.nel)*" elements.", 1);
    if np > 1
        e_count = zeros(Int, np);
        for ei=1:length(epart)
            e_count[epart[ei]+1] += 1;
        end
        log_entry("Number of elements in each partition: "*string(e_count), 2);
        log_entry("Full grid has "*string(length(grid_data.global_bdry_index))*" nodes.", 2);

    else
        log_entry("Full grid has "*string(size(grid_data.allnodes,2))*" nodes.", 2);
    end

    # Set elemental loop ordering to match the order from the grid for now.
    elemental_order = 1:grid_data.nel_owned;
end

# Write the mesh to a MSH file
function output_mesh(file, format)
    if config.proc_rank == 0
        write_mesh(file, format, mesh_data);
        log_entry("Wrote mesh data to file.", 1);
    end
end

function set_parent_and_child(p_maps, c_grid, order)
    global parent_maps = p_maps;
    global fv_grid = c_grid;
    dim = config.dimension;
    nfaces = size(fv_grid.element2face,1);
    global fv_refel = build_refel(dim, 1, nfaces, config.elemental_nodes);
    global fv_geo_factors = build_geometric_factors(fv_refel, fv_grid, do_face_detj=true, do_vol_area=true, constant_jacobian=true);
    global fv_info = build_FV_info(fv_grid, order);
    log_entry("Created parent-child grid with "*string(size(c_grid.allnodes,2))*" nodes and "*string(size(c_grid.loc2glb,2))*" elements.", 2);

    # If CELL variables exist, resize their values
    N = size(fv_grid.loc2glb, 2);
    for i=1:length(variables)
        if variables[i].location == CELL
            if variables[i].type == SCALAR
                val_size = (1, N);
            elseif variables[i].type == VECTOR
                val_size = (config.dimension, N);
            elseif variables[i].type == TENSOR
                val_size = (config.dimension*config.dimension, N);
            elseif variables[i].type == SYM_TENSOR
                val_size = (Int((config.dimension*(config.dimension+1))/2), N);
            elseif variables[i].type == VAR_ARRAY
                if typeof(variables[i].indexer) <: Array
                    comps = 1;
                    for j=1:length(variables[i].indexer)
                        comps *= length(variables[i].indexer[j].range);
                    end
                elseif typeof(variables[i].indexer) == Indexer
                    comps = length(variables[i].indexer.range);
                else
                    comps = 1;
                end
                val_size = (comps, N);
            end

            variables[i].values = zeros(val_size);
        end
    end
end

# Maybe remove. This is in grid.jl
function add_boundary_ID(bid, on_bdry)
    add_boundary_ID_to_grid(bid, on_bdry, grid_data);
end

# Defines a test function symbol by creating a special coefficient object.
function add_test_function(v, type)
    varind = length(test_functions) + 1;
    # make SymType
    dim = config.dimension;
    components = 1;
    if type == SCALAR
        components = 1;
    elseif type == VECTOR
        components = dim;
    elseif type == TENSOR
        components = dim*dim;
    elseif type == SYM_TENSOR
        if dim == 1
            components = 1;
        elseif dim == 2
            components = 3;
        elseif dim == 3
            components = 6;;
        elseif dim == 4
            components = 10;
        end
    elseif type == VAR_ARRAY
        components = 1;
    end
    symvar = sym_var(string(v), type, components);

    push!(test_functions, Finch.Coefficient(v, symvar, varind, type, NODAL, [], false););
    log_entry("Set test function symbol: "*string(v)*" of type: "*type, 2);
end

# Adds a variable and allocates everything associated with it.
function add_variable(var)
    global var_count += 1;

    # adjust values arrays
    if grid_data === nothing
        N = 1;
    else
        if var.location == CELL
            if fv_grid === nothing
                N = size(grid_data.loc2glb, 2);
            else
                N = size(fv_grid.loc2glb, 2);
            end
        else
            N = size(grid_data.allnodes,2);
        end
    end

    if var.type == SCALAR
        val_size = (1, N);
        comps = 1;
    elseif var.type == VECTOR
        val_size = (config.dimension, N);
        comps = config.dimension;
    elseif var.type == TENSOR
        val_size = (config.dimension*config.dimension, N);
        comps = config.dimension*config.dimension;
    elseif var.type == SYM_TENSOR
        val_size = (Int((config.dimension*(config.dimension+1))/2), N);
        comps = val_size[1];
    elseif var.type == VAR_ARRAY
        if typeof(var.indexer) <: Array
            comps = 1;
            for i=1:length(var.indexer)
                comps *= length(var.indexer[i].range);
            end
        elseif typeof(var.indexer) == Indexer
            comps = length(var.indexer.range);
        else
            comps = 1;
        end
        val_size = (comps, N);
    end
    var.total_components = comps;

    if language == JULIA || language == 0
        var.values = zeros(val_size);
    end

    # make symbolic layer variable symbols
    var.symvar = sym_var(string(var.symbol), var.type, var.total_components);

    global variables = [variables; var];

    global linears = [linears; nothing];
    global bilinears = [bilinears; nothing];
    global face_linears = [face_linears; nothing];
    global face_bilinears = [face_bilinears; nothing];
    global assembly_loops = [assembly_loops; nothing];
    global symexpressions[1] = [symexpressions[1]; nothing];
    global symexpressions[2] = [symexpressions[2]; nothing];
    global symexpressions[3] = [symexpressions[3]; nothing];
    global symexpressions[4] = [symexpressions[4]; nothing];
    global code_strings[1] = [code_strings[1]; ""];
    global code_strings[2] = [code_strings[2]; ""];
    global code_strings[3] = [code_strings[3]; ""];
    global code_strings[4] = [code_strings[4]; ""];
    global code_strings[5] = [code_strings[5]; ""];

    log_entry("Added variable: "*string(var.symbol)*" of type: "*var.type*", location: "*var.location, 2);
end

# Adds a coefficient with either constant value or some generated function of (x,y,z,t)
function add_coefficient(c, type, location, val, nfuns, element_array=false)
    global coefficients;
    # The values of c will have the same array structure as val
    if typeof(val) <: Array
        vals = Array{Any}(undef,size(val));
    else
        vals = [];
    end
    if nfuns == 0 # constant values
        vals = val;
        if length(vals) == 1 && !(typeof(vals) <: Array)
            vals = [val];
        end

    else # genfunction values
        if typeof(val) <: Array
            ind = length(genfunctions) - nfuns + 1;
            for i=1:length(val)
                if typeof(val[i]) == String
                    vals[i] = genfunctions[ind];
                    ind += 1;
                else
                    vals[i] = val[i];
                end
            end
        else
            push!(vals, genfunctions[end]);
        end
    end
    # If element_array, the last index in size(vals) should be per element.
    # c[n1, nel] for n1 components by nel elements.
    if element_array
        sz = size(vals);
        components = 0;
        for ind=1:length(sz)-1
            components += sz[ind];
        end
        if components == 0
            components = 1;
        end
    else
        components = length(vals);
    end

    symvar = sym_var(string(c), type, components);

    index = length(coefficients) + 1;
    push!(coefficients, Coefficient(c, symvar, index, type, location, vals, element_array));

    if element_array
        log_entry("Added coefficient "*string(c)*" : (array of elemental values)", 2);
    else
        log_entry("Added coefficient "*string(c)*" : "*string(val), 2);
    end

    return coefficients[end];
end

# Adds a parameter entity.
function add_parameter(p, type, val)
    index = length(parameters);
    push!(parameters, Parameter(p, index, type, val));

    log_entry("Added parameter "*string(p)*" : "*string(val), 2);

    return parameters[end];
end

# Adds an Indexer entity
function add_indexer(indexer)
    push!(indexers, indexer);
    if length(indexer.range) < 3
        log_entry("Added indexer "*string(indexer.symbol)*" : "*string(indexer.range));
    else
        log_entry("Added indexer "*string(indexer.symbol)*" : ["*string(indexer.range[1])*", ... "*string(indexer.range[end])*"]");
    end

end

# Defines a variable transform
function add_variable_transform(var1, var2, func)
    push!(variable_transforms, VariableTransform(var1, var2, func));
    log_entry("Added transform "*string(func)*": "*string(var1)*" -> "*string(var2));
    return variable_transforms[end];
end

# Needed to insert parameter expressions into the weak form expressions.
function swap_parameter_xyzt(ex)
    if typeof(ex) == Symbol
        if ex === :x
            return :parameterCoefficientForx ;
        elseif ex === :y
            return :parameterCoefficientFory ;
        elseif ex === :z
            return :parameterCoefficientForz ;
        elseif ex === :t
            return :parameterCoefficientFort ;
        end
    elseif typeof(ex) == Expr
        for i=1:length(ex.args)
            ex.args[i] = swap_parameter_xyzt(ex.args[i]); # Recursively swap
        end
    end
    return ex;
end

# Sets initial condition for the given variable, but does not initialize.
function add_initial_condition(varindex, ex, nfuns)
    global prob;
    while length(prob.initial) < varindex
        prob.initial = [prob.initial; nothing];
    end
    if typeof(ex) <: Array
        vals = [];
        ind = length(genfunctions) - nfuns + 1;
        for i=1:length(ex)
            if typeof(ex[i]) == String
                push!(vals, genfunctions[ind]);
                ind += 1;
            else
                push!(vals, ex[i]);
            end
        end
        prob.initial[varindex] = vals;
    else
        var_components = size(variables[varindex].values,1);
        if typeof(ex) == String
            prob.initial[varindex] = fill(genfunctions[end], var_components);
        else
            prob.initial[varindex] = fill(ex, var_components);
        end
    end

    if length(prob.initial[varindex]) < 10
        log_entry("Initial condition for "*string(variables[varindex].symbol)*" : "*string(prob.initial[varindex]), 2);
    else
        log_entry("Initial condition for "*string(variables[varindex].symbol)*" : "*string(prob.initial[varindex][1:6])*" (truncated for printing)", 2);
    end
    # hold off on initializing till solve or generate is determined.
end

# Evaluates all available intitial conditions setting the values for variables.
function eval_initial_conditions()
    dim = config.dimension;

    # build initial conditions
    for vind=1:length(variables)
        if !(variables[vind].ready) && vind <= length(prob.initial)
            if !(prob.initial[vind] === nothing)
                if variables[vind].location == CELL && !(fv_grid === nothing)
                    # Need to use the fv_grid instead of grid_data
                    this_grid_data = fv_grid;
                    this_geo_factors = fv_geo_factors;
                    this_refel = fv_refel;
                else
                    this_grid_data = grid_data;
                    this_geo_factors = geo_factors;
                    this_refel = refel;
                end

                # Evaluate at nodes
                nodal_values = zeros(length(prob.initial[vind]), size(this_grid_data.allnodes,2));
                for ci=1:length(prob.initial[vind])
                    for ni=1:size(this_grid_data.allnodes,2)
                        if typeof(prob.initial[vind][ci]) <: Number
                            nodal_values[ci,ni] = prob.initial[vind][ci];
                        elseif dim == 1
                            nodal_values[ci,ni] = prob.initial[vind][ci].func(this_grid_data.allnodes[ni],0,0,0);
                        elseif dim == 2
                            nodal_values[ci,ni] = prob.initial[vind][ci].func(this_grid_data.allnodes[1,ni],this_grid_data.allnodes[2,ni],0,0);
                        elseif dim == 3
                            nodal_values[ci,ni] = prob.initial[vind][ci].func(this_grid_data.allnodes[1,ni],this_grid_data.allnodes[2,ni],this_grid_data.allnodes[3,ni],0);
                        end
                    end
                end

                # compute cell averages using nodal values if needed
                if variables[vind].location == CELL
                    nel = size(this_grid_data.loc2glb, 2);
                    for ei=1:nel
                        # e = elemental_order[ei];
                        e = ei;
                        glb = this_grid_data.loc2glb[:,e];
                        vol = this_geo_factors.volume[e];
                        detj = this_geo_factors.detJ[e];

                        for ci=1:length(prob.initial[vind])
                            variables[vind].values[ci,e] = detj / vol * (this_refel.wg' * this_refel.Q * (nodal_values[ci,glb][:]))[1];
                        end
                    end
                else
                    variables[vind].values = nodal_values;
                end

                variables[vind].ready = true;
                log_entry("Built initial conditions for: "*string(variables[vind].symbol));
            end
        end
    end
end

# Sets boundary condition for a variable and BID.
function add_boundary_condition(var, bid, type, ex, nfuns)
    global prob;
    # make sure the arrays are big enough
    if size(prob.bc_func, 1) < var_count || size(prob.bc_func, 2) < bid
        if !(grid_data===nothing)
            nbid = length(grid_data.bids);
        else
            # For some targets no grid is created
            nbid = 1;
        end

        tmp1 = fill("", var_count, nbid);
        tmp2 = Array{Any,2}(undef, (var_count, nbid)); # need to keep this array open to any type
        fill!(tmp2, 0);
        tmp3 = zeros(Int, (var_count, nbid));

        for i=1:size(prob.bc_func,1)
            for j=1:size(prob.bc_func,2)
                tmp1[i,j] = prob.bc_type[i,j];
                tmp2[i,j] = prob.bc_func[i,j];
                tmp3[i,j] = prob.bid[i,j];
            end
        end
        prob.bc_type = tmp1;
        prob.bc_func = tmp2;
        prob.bid = tmp3;
    end

    # Add this boundary condition to the struct
    valstr = "";
    if typeof(ex) <: Array
        vals = [];
        valstr = "[";
        ind = length(genfunctions) - nfuns + 1;
        for i=1:length(ex)
            if typeof(ex[i]) <: Number
                push!(vals, ex[i]);
                valstr *= string(ex[i]);
            elseif typeof(ex[i]) == String
                push!(vals, genfunctions[ind]);
                valstr *= genfunctions[ind].name;
                ind += 1;
            else # callback
                push!(vals, ex[i]);
                valstr *= ex[i].name;
            end
            if i < length(ex)
                valstr *= ", ";
            end
        end
        valstr *= "]";
        prob.bc_func[var.index, bid] = vals;
    else
        var_components = size(var.values,1);
        if typeof(ex) <: Number
            prob.bc_func[var.index, bid] = fill(ex, var_components);
            valstr = string(ex);
        elseif typeof(ex) == String
            prob.bc_func[var.index, bid] = fill(genfunctions[end], var_components);
            valstr = genfunctions[end].name;
        else # callback
            prob.bc_func[var.index, bid] = fill(ex, var_components);
            valstr = ex.name;
        end
    end
    prob.bc_type[var.index, bid] = type;
    prob.bid[var.index, bid] = bid;

    log_entry("Boundary condition: var="*string(var.symbol)*" bid="*string(bid)*" type="*type*" val="*valstr, 2);
end

# A reference point is a single point with a defined value.
# It's treated like a Dirichlet boundary with one boundary node.
function add_reference_point(var, pos, val)
    global prob;
    # make sure the array is big enough
    if size(prob.ref_point,1) < var_count
        tmp = Array{Any,2}(undef, (var_count, 3));
        for i=1:size(prob.ref_point,1)
            tmp[i,:] = prob.ref_point[i,:];
        end
        for i=(size(prob.ref_point,1)+1):var_count
            tmp[i,1] = false;
            tmp[i,2] = [0,0];
            tmp[i,3] = [0];
        end
        prob.ref_point = tmp;
    end
    if typeof(pos) <: Number
        pos = pos*ones(config.dimension);
    end
    if typeof(val) <: Number
        val = val*ones(length(var.symvar));
    end

    # Find the closest vertex to pos
    # The stored pos is actually the index into glbvertex pointing to the closest vertex
    ind = [1,1];
    mindist = 12345;
    nel = size(grid_data.loc2glb,2);
    for ei=1:nel
        for i=1:size(grid_data.glbvertex,1)
            d = 0;
            for comp=1:length(pos)
                d = d + abs(grid_data.allnodes[comp, grid_data.glbvertex[i, ei]] - pos[comp]);
            end
            if d<mindist
                ind = [i,ei];
                mindist = d;
            end
        end
    end

    prob.ref_point[var.index, 1] = true;
    prob.ref_point[var.index, 2] = ind;
    prob.ref_point[var.index, 3] = val;

    log_entry("Reference point: var="*string(var.symbol)*" position="*string(pos)*" value="*string(val), 2);
end

# adds a callback function to the
function add_callback_function(f)
    push!(callback_functions, f);
end

# Sets the RHS(linear) function for the given variable(s).
function set_rhs(var, code="")
    global linears;
    if language == 0 || language == JULIA
        if typeof(var) <:Array
            for i=1:length(var)
                linears[var[i].index] = genfunctions[end];
                code_strings[2][var[i].index] = code;
            end
        else
            linears[var.index] = genfunctions[end];
            code_strings[2][var.index] = code;
        end

    else # external generation
        if typeof(var) <:Array
            for i=1:length(var)
                linears[var[i].index] = code;
            end
        else
            linears[var.index] = code;
        end
    end
end

# Sets the LHS(bilinear) function for the given variable(s).
function set_lhs(var, code="")
    global bilinears;
    if language == 0 || language == JULIA
        if typeof(var) <:Array
            for i=1:length(var)
                bilinears[var[i].index] = genfunctions[end];
                code_strings[1][var[i].index] = code;
            end
        else
            bilinears[var.index] = genfunctions[end];
            code_strings[1][var.index] = code;
        end

    else # external generation
        if typeof(var) <:Array
            for i=1:length(var)
                bilinears[var[i].index] = code;
            end
        else
            bilinears[var.index] = code;
        end
    end
end

# Sets the surface RHS(linear) function for the given variable(s).
function set_rhs_surface(var, code="")
    global face_linears;
    if language == 0 || language == JULIA
        if typeof(var) <:Array
            for i=1:length(var)
                face_linears[var[i].index] = genfunctions[end];
                code_strings[4][var[i].index] = code;
            end
        else
            face_linears[var.index] = genfunctions[end];
            code_strings[4][var.index] = code;
        end

    else # external generation
        if typeof(var) <:Array
            for i=1:length(var)
                face_linears[var[i].index] = code;
            end
        else
            face_linears[var.index] = code;
        end
    end
end

# Sets the surface LHS(bilinear) function for the given variable(s).
function set_lhs_surface(var, code="")
    global face_bilinears;
    if language == 0 || language == JULIA
        if typeof(var) <:Array
            for i=1:length(var)
                face_bilinears[var[i].index] = genfunctions[end];
                code_strings[3][var[i].index] = code;
            end
        else
            face_bilinears[var.index] = genfunctions[end];
            code_strings[3][var.index] = code;
        end

    else # external generation
        if typeof(var) <:Array
            for i=1:length(var)
                face_bilinears[var[i].index] = code;
            end
        else
            face_bilinears[var.index] = code;
        end
    end
end

function set_symexpressions(var, ex, lorr, vors)
    # If the ex is empty, don't set anything
    if ex == [] || ex == [[]]
        return;
    end

    if typeof(var) <:Array
        for i=1:length(var)
            set_symexpressions(var[i], ex, lorr, vors);
        end

    else
        global symexpressions;
        if lorr == LHS
            if vors == "volume"
                ind = 1;
            else
                ind = 2;
            end
        else # rhs
            if vors == "volume"
                ind = 3;
            else
                ind = 4;
            end
        end
        symexpressions[ind][var.index] = ex;
    end
end

function set_assembly_loops(var, code="")
    global assembly_loops;
    if language == 0 || language == JULIA
        if typeof(var) <:Array
            for i=1:length(var)
                assembly_loops[var[i].index] = genfunctions[end];
                code_strings[5][var[i].index] = code;
            end
        else
            assembly_loops[var.index] = genfunctions[end];
            code_strings[5][var.index] = code;
        end

    else # external generation
        if typeof(var) <:Array
            for i=1:length(var)
                assembly_loops[var[i].index] = code;
            end
        else
            assembly_loops[var.index] = code;
        end
    end
end

function finalize()
    finalize_finch();
end

end # module

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
@testset "acceptance tracker" begin
    at = BSI.AccptTracker(0)

    @testset "initialisation" begin
        @test at.accpt == 0
        @test at.prop == 0
    end

    accept_reject = [true, false, false, true, true, false]
    for ar in accept_reject BSI.register_accpt!(at, ar) end

    @testset "testing update!" begin
        @test BSI.acceptance_rate(at) == sum(accept_reject)/length(accept_reject)
    end
    BSI.reset!(at)
    @testset "after reset" begin
        @test at.accpt == 0
        @test at.prop == 0
    end

    at_vec = [BSI.AccptTracker(0) for i in 1:6]
    BSI.register_accpt!(at_vec, accept_reject)
    @testset "testing registration for vector of accpt tracker!" begin
        @test all([a.prop == 1 for a in at_vec])
        @test all([at_vec[i].accpt == 1*accept_reject[i] for i in 1:length(at_vec)])
    end
    BSI.reset!(at_vec)
    @testset "vector after reset" begin
        @test all([a.prop == 0 for a in at_vec])
        @test all([a.accpt == 0 for a in at_vec])
    end
end

@testset "definition of parameter update" begin
    θ = fill(0.0,5) # just for the dimension
    pu1 = BSI.ParamUpdate(BSI.MetropolisHastingsUpdt(), (1,), θ,
                          BSI.UniformRandomWalk(1.0,true), BSI.ImproperPosPrior(),
                          BSI.UpdtAuxiliary(BSI.Vern7(), true))
    mvn_prior = MvNormal([0.0,0.0], diagm(0=>fill(1000.0, 2)))
    pu2 = BSI.ParamUpdate(BSI.ConjugateUpdt(), (2,3), θ, nothing, mvn_prior,
                          BSI.UpdtAuxiliary(BSI.Vern7(), true))
    @testset "initialisation" begin
        @test typeof(pu1) <: BSI.ParamUpdate{BSI.MetropolisHastingsUpdt}
        @test typeof(pu2) <: BSI.ParamUpdate{BSI.ConjugateUpdt}
        @test pu1.updt_coord == Val((true,false,false,false,false))
        @test pu2.updt_coord == Val((false,true,true,false,false))
        @test pu1.t_kernel.ϵ == 1.0
        @test pu1.t_kernel.pos
        @test pu2.t_kernel === nothing
        @test typeof(pu1.priors) <: BSI.Priors
        @test typeof(pu2.priors) <: BSI.Priors
        @test typeof(pu1.priors.priors) <: Tuple{BSI.ImproperPosPrior}
        @test pu2.priors.priors == (mvn_prior,)
    end
end


text:
import ForwardDiff.Dual

"""
    ivol(core::CoreModel, c::Contract, value, TianModel, nsteps=100; ...)

Compute the implied volatility of `c` using the Tian binomial model with `nsteps` steps.

This uses Newton's method, the settings of which are controlled by the keyword arguments:

 - `maxiter`: maximum number of iterations (default = 20).
 - `rtol`: relative tolerance for the implied volatility (default = 1e-6).
 - `σ0`: initial guess for volatility (default = 0.1).
"""

function ivol(core::CoreModel{T}, c::Contract, val, ::Type{S}, nsteps::Int=100;
                maxiter::Int=20, rtol::Float64=1e-6, σ0::Float64=0.1) where {T,S<:BinomialModel}

    startdate = core.carrycurve.reference_date
    enddate = maturitydate(c)
    startvalue = core.startprice
    interestrate = core.yieldcurve.rate
    carryrate = core.carrycurve.rate

    f(σ) = value(S(startdate, enddate, nsteps, startvalue,
                    interestrate, carryrate, σ), c) - val

    σ = σ0
    for i = 1:maxiter
        dσ = Dual(σ, 1.0)
        dv = f(dσ)

        v = ForwardDiff.value(dv)
        Δ = extract_derivative(dv)
        ϵ = v / Δ
        σ -= ϵ

        if abs(ϵ) < rtol*σ
            return σ
        end
    end
    error("Failed to converge")
end


text:
# function Base.:+(sys1::NLStateSpace, sys2::NLStateSpace)
#     @assert(size(sys1) == size(sys2), "Systems have different shapes")
#     @assert(sys1.Ts == sys2.Ts, "Sampling time mismatch")
#
#     # f(x) = sys1.f(x[1:nstates(sys1)])
# end

function combine(f2::StateFunction, f1::StateFunction, h1::OutputFunction{true})
    function f!(dx, x, u, t)
        dx1, dx2, x1, x2 = @views dx[1], dx[2], x[1], x[2]
        y = h1(x1, u, t)
        f1(dx1, x1, u, t)
        f2(dx2, x2, y, t)
        return nothing
    end
end
combine(h2::OutputFunction{true}, h1::OutputFunction{true}) =
    OutputFunction((x, u, t) -> h2(x[2], h1(x[1], u, t), t))
combine(h2::OutputFunction{true}, h1::OutputFunction{false}) =
    OutputFunction((x, t) -> h2(x[2], h1(x[1], t), t))
combine(h2::OutputFunction{false}, h1) = OutputFunction((x, t) -> h2(x[2], t))

# Note reverse order of sys1 and sys2 because of right-to-left operation
function Base.:*(sys2::NLStateSpace, sys1::NLStateSpace)
    # Check that input/output dimensions and sampling times are compatable
    @assert(CS.ninputs(sys1) == CS.noutputs(sys2),
        "sys1 must have same number of inputs as sys2 has outputs"
    )
    @assert(sys1.Ts == sys2.Ts, "Sampling time mismatch")

    # Input/output dimensions
    nu = sys1.nu
    ny = sys2.ny

    # Initial conditions
    init = initial_conditions.([sys1, sys2])

    # New state function
    function _f!(dx, x, u, t)
        dx1, dx2, x1, x2 = @views dx[1], dx[2], x[1], x[2]
        y = sys1.h(x1, u, t)
        sys1.f(dx1, x1, u, t)
        sys2.f(dx2, x2, y, t)
        return nothing
    end
    f! = StateFunction(_f!, nu, init)

    # New output function
    h = combine(sys2.h, sys1.h)

    return NLStateSpace(f!, h, sys1.Ts)
end

max_stars_repo_name: jonniedie/ComposableDiffEqs

max_issues_repo_name: jonniedie/ComposableDiffEqs

max_forks_repo_name: jonniedie/ComposableDiffEqs


text:
"""
RCall.jl's type `Sxp` mirrors the R symbolic expression record `SEXPREC` in R API.
These are represented by a pointer `Ptr{S<:Sxp}` (which is called `SEXP` in R API).
"""
abstract type Sxp end # SEXPREC
const SxpPtrInfo = UInt32 # sxpinfo_struct

"R Sxp header: a pointer to this is used for unknown types."
struct SxpHead <: Sxp # SEXPREC_HEADER
    info::SxpPtrInfo
    attrib::Ptr{Cvoid}
    gc_next::Ptr{Cvoid}
    gc_prev::Ptr{Cvoid}
end
const UnknownSxp = SxpHead

abstract type VectorSxp <: Sxp end
abstract type VectorAtomicSxp <: VectorSxp end
abstract type VectorNumericSxp <: VectorAtomicSxp end
abstract type VectorListSxp <: VectorSxp end
abstract type PairListSxp <: Sxp end
abstract type FunctionSxp <: Sxp end

"R NULL value"
struct NilSxp <: PairListSxp   # type tag 0
    head::SxpHead
end

"R pairs (cons) list cell"
struct ListSxp <: PairListSxp  # type tag 2
    head::SxpHead
    car::Ptr{UnknownSxp}
    cdr::Ptr{UnknownSxp}
    tag::Ptr{UnknownSxp}
end

"R function closure"
struct ClosSxp <: FunctionSxp  # type tag 3
    head::SxpHead
    formals::Ptr{ListSxp}
    body::Ptr{UnknownSxp}
    env::Ptr{UnknownSxp}
end

"R environment"
struct EnvSxp <: Sxp  # type tag 4
    head::SxpHead
    frame::Ptr{UnknownSxp}
    enclos::Ptr{UnknownSxp}
    hashtab::Ptr{UnknownSxp}
end

"R promise"
struct PromSxp <: Sxp  # type tag 5
    head::SxpHead
    value::Ptr{UnknownSxp}
    expr::Ptr{UnknownSxp}
    env::Ptr{UnknownSxp}
end

"R function call"
struct LangSxp <: PairListSxp  # type tag 6
    head::SxpHead
    car::Ptr{UnknownSxp}
    cdr::Ptr{UnknownSxp}
    tag::Ptr{UnknownSxp}
end

"R special function"
struct SpecialSxp <: FunctionSxp  # type tag 7
    head::SxpHead
end

"R built-in function"
struct BuiltinSxp <: FunctionSxp  # type tag 8
    head::SxpHead
end

"R character string"
struct CharSxp <: VectorAtomicSxp     # type tag 9
    head::SxpHead
    length::Cint
    truelength::Cint
end

"R symbol"
struct SymSxp <: Sxp   # type tag 1
    head::SxpHead
    name::Ptr{CharSxp}
    value::Ptr{UnknownSxp}
    internal::Ptr{UnknownSxp}
end

"R logical vector"
struct LglSxp <: VectorNumericSxp     # type tag 10
    head::SxpHead
    length::Cint
    truelength::Cint
end

"R integer vector"
struct IntSxp <: VectorNumericSxp     # type tag 13
    head::SxpHead
    length::Cint
    truelength::Cint
end

"R real vector"
struct RealSxp <: VectorNumericSxp    # type tag 14
    head::SxpHead
    length::Cint
    truelength::Cint
end

"R complex vector"
struct CplxSxp <: VectorNumericSxp    # type tag 15
    head::SxpHead
    length::Cint
    truelength::Cint
end

"R vector of character strings"
struct StrSxp <: VectorListSxp     # type tag 16
    head::SxpHead
    length::Cint
    truelength::Cint
end

"R dot-dot-dot object"
struct DotSxp <: Sxp     # type tag 17
    head::SxpHead
end

"R \"any\" object"
struct AnySxp <: Sxp     # type tag 18
    head::SxpHead
end

"R list (i.e. Array{Any,1})"
struct VecSxp <: VectorListSxp     # type tag 19
    head::SxpHead
    length::Cint
    truelength::Cint
end

"R expression vector"
struct ExprSxp <: VectorListSxp    # type tag 20
    head::SxpHead
    length::Cint
    truelength::Cint
end

"R byte code"
struct BcodeSxp <: Sxp   # type tag 21
    head::SxpHead
end

"R external pointer"
struct ExtPtrSxp <: Sxp  # type tag 22
    head::SxpHead
    ptr::Ptr{Cvoid}
    prot::Ptr{Cvoid}
    tag::Ptr{UnknownSxp}
end

"R weak reference"
struct WeakRefSxp <: Sxp  # type tag 23
    head::SxpHead
end

"R byte vector"
struct RawSxp <: VectorAtomicSxp      # type tag 24
    head::SxpHead
    length::Cint
    truelength::Cint
end

"R S4 object"
struct S4Sxp <: Sxp      # type tag 25
    head::SxpHead
end

"""
An `RObject` is a Julia wrapper for an R object (known as an "S-expression" or "SEXP"). It is stored as a pointer which is protected from the R garbage collector, until the `RObject` itself is finalized by Julia. The parameter is the type of the S-expression.

When called with a Julia object as an argument, a corresponding R object is constructed.

```julia_skip
julia> RObject(1)
RObject{IntSxp}
[1] 1

julia> RObject(1:3)
RObject{IntSxp}
[1] 1 2 3

julia> RObject(1.0:3.0)
RObject{RealSxp}
[1] 1 2 3
```
"""
mutable struct RObject{S<:Sxp}
    p::Ptr{S}
    # used for pre-defined constants
    function RObject{S}() where S
        new{S}(C_NULL)
    end
    function RObject{S}(p::Ptr{S}) where S
        preserve(p)
        r = new{S}(p)
        finalizer(release, r)
        r
    end
    # SymSxps are not garbage collected, so preserve not necessary.
    RObject{S}(p::Ptr{SymSxp}) where S = new{S}(p)
end

RObject(p::Ptr{S}) where S<:Sxp = RObject{S}(p)
RObject(x::RObject) = x

"""
R Class
"""
struct RClass{Symbol} end

"""
Element types of R vectors.
"""
eltype(::Type{LglSxp}) = Cint
eltype(::Type{IntSxp}) = Cint
eltype(::Type{RealSxp}) = Float64
eltype(::Type{CplxSxp}) = ComplexF64
eltype(::Type{CharSxp}) = UInt8
eltype(::Type{RawSxp}) = UInt8

eltype(::Type{StrSxp}) = Ptr{CharSxp}
eltype(::Type{VecSxp}) = Ptr{UnknownSxp}
eltype(::Type{ExprSxp}) = Ptr{UnknownSxp}

eltype(s::Ptr{S}) where S<:Sxp = eltype(S)
eltype(s::RObject{S}) where S<:Sxp = eltype(S)

"""
Prevent garbage collection of an R object. Object can be released via `release`.

This is slower than `protect`, as it requires searching an internal list, but
more flexible.
"""
preserve(p::Ptr{S}) where S<:Sxp = ccall((:R_PreserveObject,libR), Nothing, (Ptr{S},), p)

"""
Release object that has been gc protected by `preserve`.
"""
release(p::Ptr{S}) where S<:Sxp = ccall((:R_ReleaseObject,libR),Nothing,(Ptr{S},),p)
release(r::RObject{S}) where S<:Sxp = release(r.p)

"""
Stack-based protection of garbage collection of R objects. Objects are
released via `unprotect`. Returns the same pointer, allowing inline use.

This is faster than `preserve`, but more restrictive. Really only useful
inside functions.
"""
protect(p::Ptr{S}) where S<:Sxp = ccall((:Rf_protect,libR), Ptr{S}, (Ptr{S},), p)

"""
Release last `n` objects gc-protected by `protect`.
"""
unprotect(n::Integer) = ccall((:Rf_unprotect,libR), Nothing, (Cint,), n)

"""
The SEXPTYPE number of a `Sxp`

Determined from the trailing 5 bits of the first 32-bit word. Is
a 0-based index into the `info` field of a `SxpHead`.
"""
sexpnum(h::SxpHead) = h.info & 0x1f
sexpnum(p::Ptr{S}) where S<:Sxp = sexpnum(unsafe_load(p))

"vector of R Sxp types"
const typs = [NilSxp,SymSxp,ListSxp,ClosSxp,EnvSxp,
              PromSxp,LangSxp,SpecialSxp,BuiltinSxp,CharSxp,
              LglSxp,Nothing,Nothing,IntSxp,RealSxp,
              CplxSxp,StrSxp,DotSxp,AnySxp,VecSxp,
              ExprSxp,BcodeSxp,ExtPtrSxp,WeakRefSxp,RawSxp,
              S4Sxp]

for (i,T) in enumerate(typs)
    if T != Nothing
        @eval sexpnum(::Type{$T}) = $(i-1)
    end
end

"""
Convert a `Ptr{UnknownSxp}` to an approptiate `Ptr{S<:Sxp}`.
"""
function sexp(p::Ptr{UnknownSxp})
    typ = sexpnum(p)
    0 ≤ typ ≤ 10 || 13 ≤ typ ≤ 25 || error("Unknown SEXPTYPE $typ")
    styp = typs[typ+1]
    Ptr{styp}(p)
end
sexp(s::Ptr{S}) where S<:Sxp = s
sexp(r::RObject) = r.p

sexp(::Type{S},s::Ptr{S}) where S<:Sxp = s
sexp(::Type{S},r::RObject{S}) where S<:Sxp = r.p


text:
using FEASTSolver
using LinearAlgebra
using Printf
using Random
using MatrixMarket
using NonlinearEigenproblems: nep_gallery, compute_Mder, polyeig, PEP
using DelimitedFiles

function bf_arrays()
    N = diagm(-1 => [1,1,1,1,1,1,1])
    Mh0 = 1/6*(4*I + N + N')
    Mh1 = N - N'
    Mh2 = -1*(2*I - N - N')
    Mh3 = Mh1
    Mh4 = -Mh2
    c = [0.6 1.3; 1.3 0.1; 0.1 1.2; 1.0 1.0; 1.2 1.0]
    I8 = Matrix(I, 8, 8)
    M0 = c[1,1] * kron(I8, Mh0) + c[1,2] * kron(Mh0, I8)
    M1 = c[2,1] * kron(I8, Mh1) + c[2,2] * kron(Mh1, I8)
    M2 = c[3,1] * kron(I8, Mh2) + c[3,2] * kron(Mh2, I8)
    M3 = c[4,1] * kron(I8, Mh3) + c[4,2] * kron(Mh3, I8)
    M4 = c[5,1] * kron(I8, Mh4) + c[5,2] * kron(Mh4, I8)
    return [M0, M1, M2, M3, M4]
end

function deficient_arrays()
	# Random.seed!(1234)
	# A0 = complex.(rand(15,15))
	# A1 = complex.(rand(15,15))
	# A0[:,1] .= 0
	A0 = mmread("data/quadraticM0.mtx")
	A1 = mmread("data/quadraticM1.mtx")

	return [(A0 .- 0.02.*A1), (0.1.*A1), A1]
end

function quadratic()
	A0 = mmread("data/system5A0.mtx")
	A1 = mmread("data/system5A1.mtx")
	A2 = mmread("data/system5A2.mtx")

	function T(z::ComplexF64)
		return (A2 .* (z*z)) .+ (A1 .* z) .+ A0
	end

	return T
end

function gun()
	nep_gun = nep_gallery("nlevp_native_gun")
	T(z) = compute_Mder(nep_gun, z)
	return T
end

function hadeler()
	nep_hadeler = nep_gallery("nlevp_native_hadeler", 100, 200)
	T(z) = compute_Mder(nep_hadeler, z)
	return T
end

# for i=2:8
# 	println("\nBeyn N=$(2^i)")
# 	e, v, res = nlfeast!(bf(), rand(ComplexF64,64,30), 2^i, 0, c=complex(1.0,1.0), r=0.5)
# 	info(e, v, res, complex(1.0,1.0), 0.5)
# end

# Random.seed!(1729)

println("\nButterfly Problem")
nep = PEP(bf_arrays())
e,v = polyeig(nep)
e_in = e[abs.(e .- (1.0+1.0im)) .<= 0.5]
bf_eig = hcat(real.(e), imag.(e))
bf_eig_in = hcat(real.(e_in), imag.(e_in))
writedlm("test/data/bf_eig.dat", bf_eig)
writedlm("test/data/bf_eig_in.dat", bf_eig_in)

println("\nDeficient Quadratic Problem")
nep = PEP(deficient_arrays())
e,v = polyeig(nep)
e_in = e[abs.(e) .<= 0.25]
def_eig = hcat(real.(e), imag.(e))
def_eig_in = hcat(real.(e_in), imag.(e_in))
writedlm("test/data/deficient_eig.dat", def_eig)
writedlm("test/data/deficient_eig_in.dat", def_eig_in)

# println("\nLarger Quadratic Problem")
# quad_feast_call(X, N, iter, err_tol) = nlfeast_moments!(quadratic(), X, N, iter, c=complex(-1.55,0.0), r=0.05, ϵ=err_tol, moments=2, store=true)
# quad_res = res_table(quad_feast_call, (1000,70), max_iter_all, 6, 1e-6, complex(-1.55,0.0), 0.05, err_tol_all)
# writedlm("test/data/quad.dat", zip(1:max_iter_all, eachcol(quad_res)...))

println("\nGun Cavity Problem")
# gun_feast_call(X, N, iter, err_tol) = nlfeast!(gun(), X, N, iter, c=complex(140000.0, 0.0), r=30000, ϵ=err_tol, store=false)
# gun_res = res_table(gun_feast_call, (9956,32), max_iter_all, 6, 1e-6, complex(140000.0, 0.0), 30000, err_tol_all)
e, v, res = nlfeast!(gun(), rand(ComplexF64,9956,84), 2^4, 10, c=complex(140000.0, 0.0), r=2*30000, debug=true, ϵ=10e-16, store=false)
e_conv = e[res .<= 1e-10]
e_in = e_conv[abs.(e_conv .- 140000) .<= 30000]
gun_eig = hcat(real.(e_conv), imag.(e_conv))
gun_eig_in = hcat(real.(e_in), imag.(e_in))
writedlm("test/data/gun_eig.dat", gun_eig)
writedlm("test/data/gun_eig_in.dat", gun_eig_in)

println("\nHadeler Problem")
# gun_feast_call(X, N, iter, err_tol) = nlfeast!(gun(), X, N, iter, c=complex(140000.0, 0.0), r=30000, ϵ=err_tol, store=false)
# gun_res = res_table(gun_feast_call, (9956,32), max_iter_all, 6, 1e-6, complex(140000.0, 0.0), 30000, err_tol_all)
e, v, res = nlfeast!(hadeler(), rand(ComplexF64,200,100), 2^5, 10, c=complex(-30.0, 0.0), r=20, debug=true, ϵ=10e-16, store=true)
e_conv = e[res .<= 1e-10]
e_in = e_conv[abs.(e_conv .+ 30) .<= 10]
hadeler_eig = hcat(real.(e_conv), imag.(e_conv))
hadeler_eig_in = hcat(real.(e_in), imag.(e_in))
writedlm("test/data/hadeler_eig.dat", hadeler_eig)
writedlm("test/data/hadeler_eig_in.dat", hadeler_eig_in)


text:
@require DataStreams begin
using DataStreams
using WeakRefStrings

immutable EnumerableDataStream{T, S<:DataStreams.Data.Source, TC} <: Enumerable{T}
    source::S
    schema::DataStreams.Data.Schema
end

function query{S<:DataStreams.Data.Source}(source::S)
	schema = Data.schema(source)

    col_expressions = Array{Expr,1}()
    columns_tuple_type = Expr(:curly, :Tuple)

    for i in 1:schema.cols
        if schema.types[i] <: WeakRefString
            col_type = String
        elseif schema.types[i] <: Nullable && schema.types[i].parameters[1] <: WeakRefString
            col_type = Nullable{String}
        else
            col_type = schema.types[i]
        end

        push!(col_expressions, Expr(:(::), schema.header[i], col_type))
        push!(columns_tuple_type.args, col_type)
    end
    t_expr = NamedTuples.make_tuple(col_expressions)

    t2 = :(Query.EnumerableDataStream{Float64,Float64,Float64})
    t2.args[2] = t_expr
    t2.args[3] = typeof(source)
    t2.args[4] = columns_tuple_type

    eval(NamedTuples, :(import Query))
    t = eval(NamedTuples, t2)

    e_df = t(source, schema)

    return e_df
end

function length{T, S<:DataStreams.Data.Source, TC}(iter::EnumerableDataStream{T,S,TC})
    return iter.schema.rows
end

function eltype{T, S<:DataStreams.Data.Source, TC}(iter::EnumerableDataStream{T,S,TC})
    return T
end

function start{T, S<:DataStreams.Data.Source, TC}(iter::EnumerableDataStream{T,S,TC})
    return 1
end

@generated function next{T, S<:DataStreams.Data.Source, TC}(iter::EnumerableDataStream{T,S,TC}, state)
    constructor_call = Expr(:call, :($T))
    for i in 1:length(TC.types)
    	col_type = TC.types[i] <: WeakRefString ? String : TC.types[i]
        push!(constructor_call.args, :(Data.getfield(source, $col_type, row, $i)))
    end

    quote
    	source = iter.source
        row = state
        a = $constructor_call
        return a, state+1
    end
end

function done{T, S<:DataStreams.Data.Source, TC}(iter::EnumerableDataStream{T,S,TC}, state)
    return Data.isdone(iter.source,state,1)
end

end


text:
using Plots
using IndexedContour

x = collect(-4:0.01:3) # hide
y = collect(-5:0.01:5) # hide
z = [Float64((xi^2 + yi^2)) for xi in x, yi in y] # hide
w = [Float64((xi^2 - yi^2)) for xi in x, yi in y] # hide
p = plot([0], [0], lw = 0.001, linecolor=:white, label = "", title="Contours of x^2+y^2 colored by x^2-y^2")
pw = heatmap(x,y,w', title="x^2 - y^2")
display(pw)

for cl in levels(contours(x,y,z))
    lvl = level(cl) # the z-value of this contour level
    for line in lines(cl)
        xs, ys = coordinates(line) # coordinates of this line segment
        xi, yi, zi = indices(line)
        xss = x[xi]
        yss = y[yi]
        zss = w[zi]
        #plot!(p, xs, ys, label = "", line_z = w[zi], lw = 3)
        scatter!(p, xss, yss, label = "", markersize = 1.5, marker_z = zss, msw=0.0001)
        #plot!(p, xss, yss, label = "", linecolor = :red, lw=3)
    end
end
display(p)


text:
@testset "Bold x*" begin # issue 223
    h = raw"**x\***" |> seval
    @test h == "
x*
\n"

    h = raw"_x\__" |> seval
    @test h == "
x_
\n"
end

@testset "Bold code" begin # issue 222
    h = raw"""A **`master`** B.""" |> seval
    @test h == "
A master B.
\n"
end

@testset "Tickssss" begin # issue 219
    st = raw"""A `B` C"""
    tokens = F.find_tokens(st, F.MD_TOKENS, F.MD_1C_TOKENS)
    @test tokens[1].name == :CODE_SINGLE
    @test tokens[2].name == :CODE_SINGLE

    st = raw"""A ``B`` C"""
    tokens = F.find_tokens(st, F.MD_TOKENS, F.MD_1C_TOKENS)
    @test tokens[1].name == :CODE_DOUBLE
    @test tokens[2].name == :CODE_DOUBLE

    st = raw"""A ``` B ``` C"""
    tokens = F.find_tokens(st, F.MD_TOKENS, F.MD_1C_TOKENS)
    @test tokens[1].name == :CODE_TRIPLE
    @test tokens[2].name == :CODE_TRIPLE

    st = raw"""A ````` B ````` C"""
    tokens = F.find_tokens(st, F.MD_TOKENS, F.MD_1C_TOKENS)
    @test tokens[1].name == :CODE_PENTA
    @test tokens[2].name == :CODE_PENTA

    st = raw"""A ```b B ``` C"""
    tokens = F.find_tokens(st, F.MD_TOKENS, F.MD_1C_TOKENS)
    @test tokens[1].name == :CODE_LANG
    @test tokens[2].name == :CODE_TRIPLE

    st = raw"""A `````b B ````` C"""
    tokens = F.find_tokens(st, F.MD_TOKENS, F.MD_1C_TOKENS)
    @test tokens[1].name == :CODE_LANG2
    @test tokens[2].name == :CODE_PENTA

    h = raw"""
        A
        `````markdown
        B
        `````
        C
        """ |> fd2html_td

    @test isapproxstr(h, raw"""

A

B

C

            """)

    h = raw"""
        A
        `````markdown
        ```julia
        B
        ```
        `````
        C
        """ |> fd2html_td

    @test isapproxstr(h, raw"""

A

```julia
            B
            ```

C

            """)
end

@testset "Nested ind" begin # issue 285
    h = raw"""
    \newcommand{\hello}{
        yaya
        bar bar
    }
    \hello
    """ |> fd2html_td
    @test isapproxstr(h, raw"""yaya  bar bar""")

    h = raw"""
    @@da
        @@db
            @@dc
                blah
            @@
        @@
    @@
    """ |> fd2html_td
    @test isapproxstr(h, raw"""

                        blah

            """)
    h = raw"""
    \newcommand{\hello}[1]{#1}
    \hello{
        good lord
    }
    """ |> fd2html_td
    @test isapproxstr(h, "good lord")
end

@testset "Double brace" begin
    s = """
        @def title = "hello"
        {{title}}{{title}}
        """ |> fd2html_td
    @test isapproxstr(s, "
hellohello
")
    s = """
        @def a_b = "hello"
        @def c_d = "goodbye"
        {{a_b}}{{c_d}}
        """ |> fd2html_td
    @test isapproxstr(s, "
hellogoodbye
")
end

# issue 424 with double braces
@testset "Double brace2" begin
    s = raw"""
        @def title = "hello"
        {{title}}
        $\rho=\frac{e^{-\beta \mathcal{E}_{s}}} {\mathcal{Z}} $
        """ |> fd2html_td
    @test s // raw"""

hello \(\rho=\frac{e^{-\beta \mathcal{E}_{s}}} {\mathcal{Z}} \)
"""
end

# issue 432 and consequences
@testset "Hz rule" begin
    # issue 432
    s = raw"""
        hello[^a]

        [^a]: world

        ---
        """ |> fd2html_td
    @test isapproxstr(s, """

hello[1]

                [1]
                world

        """)
    s = raw"""
        A
        ---
        """ |> fd2html_td
    @test isapproxstr(s, """

A
""")
    s = raw"""
        A
        ----
        ****
        """ |> fd2html_td
    @test isapproxstr(s, "
A
\n
")

    # cases where nothing should happen
    s = raw"""
        A
        ---*
        ***_
        """ |> fd2html_td
    @test isapproxstr(s, "
A –-* ***_
") # note double -- is transformed -
end

# issue 439 and consequences
@testset "mathunicode" begin
    s = raw"""
        \newcommand{\u}{1}
        $$ φφφ\u abcdef $$
        """ |> fd2html_td
    @test isapproxstr(s, raw"\[ φφφ1 abcdef \]")
end

@testset "emoji-basics" begin
    s = "abc :+1: def" |> fd2html
    @test s // "
abc 👍 def
"
    s = "abc :+10: def" |> fd2html
    @test s // "
abc :+10: def
"
    s = "abc :joy: def" |> fd2html
    @test s // "
abc 😂 def
"
    s = "abc :joy def:" |> fd2html
    @test s // "
abc :joy def:
"
    s = "abc : joy: def" |> fd2html
    @test s // "
abc : joy: def
"
    s = "abc :joi: def" |> fd2html
    @test s // "
abc :joi: def
"
end


text:
@graphs begin
	attack_threat_cond = """
	A *mobile* army can attack *threat*.
	"""

	move_adjacent_cond = """
	Army/armies can move so that a *mobile* army is created adjacent to *threat*.
	"""

	move_stronghold_cond = """
	Army/armies can move such that the *value* at a stronghold under *threat* increases.
	"""

	move_toward_cond = """
	*Mobile* army/armies can move towards their *target*/*targets* and reduce their distance to the *threat*.
	"""

	move_exposed_cond = """
	Army/armies can move towards an *exposed* region.
	"""

	muster_cond = """
	Possible to muster at a stronghold under *threat*.
	"""

	character_cond = """
	A *threat* is besieging a Shadow stronghold whose army leadership is less than 5 and less than the number of army units.
	"""

	move_remain_cond = """
	The Army Die have one move remaining.
	"""

	attack_text = """
	Attack according to the latest statement. Select army at random if several can perform such an attack.
	"""

	move_text = """
	Move according to the latest statement. Select army at random if several can perform such a move.
	"""

	muster_text = """
	*Focus* priority:
	1. Stronghold under *threat*
	2. Random

	Muster:
	*Primary*: Elite
	*Secondary*: Regular

	If unit is unavailable, rotate as:
	Elite -> Regular -> Nazgûl -> Elite
	"""

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

	@node threat_exposed = Start() -> tx_t
	@node tx_t = BinaryCondition("A *threat* exist.") -> [n_true = tx_1, n_false = tx_5]

	########################################
	@node tx_1 = SetActiveDie('C') -> [next = tx_1_char_cond, no_die = tx_1_army]

	@node tx_1_char_cond = BinaryCondition(attack_threat_cond) -> [n_true = tx_1_use_die, n_false = tx_1_army]
	@node tx_1_army = SetActiveDie('A', may_use_ring = true) -> [next = tx_1_army_cond, no_die = tx_2]

	@node tx_1_army_cond = BinaryCondition(attack_threat_cond) -> [n_true = tx_1_use_die, n_false = tx_2]
	@node tx_1_use_die = UseActiveDie() -> tx_1_action
	@node tx_1_action = PerformAction(attack_text) -> tx_1_end
	@node tx_1_end = End() -> []

	########################################
	@node tx_2 = SetActiveDie('C') -> [next = tx_2_char_cond, no_die = tx_2_army]

	@node tx_2_char_cond = BinaryCondition(move_adjacent_cond) -> [n_true = tx_move_use_die, n_false = tx_2_army]
	@node tx_2_army = SetActiveDie('A', may_use_ring = true) -> [next = tx_2_army_cond, no_die = tx_3]

	@node tx_2_army_cond = BinaryCondition(move_adjacent_cond) -> [n_true = tx_move_use_die, n_false = tx_3]

	########################################
	@node tx_3 = SetActiveDie('C') -> [next = tx_3_char_cond, no_die = tx_3_army]

	@node tx_3_char_cond = BinaryCondition(move_stronghold_cond) -> [n_true = tx_move_use_die, n_false = tx_3_army]
	@node tx_3_army = SetActiveDie('A', may_use_ring = true) -> [next = tx_3_army_cond, no_die = tx_4]

	@node tx_3_army_cond = BinaryCondition(move_stronghold_cond) -> [n_true = tx_move_use_die, n_false = tx_4]

	########################################
	@node tx_4 = SetActiveDie('C') -> [next = tx_4_char_cond, no_die = tx_4_army]

	@node tx_4_char_cond = BinaryCondition(move_toward_cond) -> [n_true = tx_move_use_die, n_false = tx_4_army]
	@node tx_4_army = SetActiveDie('A', may_use_ring = true) -> [next = tx_4_army_cond, no_die = tx_5]

	@node tx_4_army_cond = BinaryCondition(move_toward_cond) -> [n_true = tx_move_use_die, n_false = tx_5]

	########################################
	@node tx_5 = SetActiveDie('C') -> [next = tx_5_char_cond, no_die = tx_5_army]

	@node tx_5_char_cond = BinaryCondition(move_exposed_cond) -> [n_true = tx_move_use_die, n_false = tx_5_army]
	@node tx_5_army = SetActiveDie('A', may_use_ring = true) -> [next = tx_5_army_cond, no_die = tx_m]

	@node tx_5_army_cond = BinaryCondition(move_exposed_cond) -> [n_true = tx_move_use_die, n_false = tx_m]

	########################################
	@node tx_move_use_die = UseActiveDie() -> tx_move_action
	@node tx_move_action = PerformAction(move_text) -> tx_move_army_die_to_move
	@node tx_move_army_die_to_move = CheckActiveDie('A') -> [n_true=tx_move_movement_remains, n_false=tx_move_end]
	@node tx_move_movement_remains = BinaryCondition(move_remain_cond) -> [n_true = tx_use_remaining_movement, n_false = tx_move_end]
	@node tx_use_remaining_movement = JumpToGraph("movement_attack_corr") -> tx_m
	@node tx_move_end = End() -> []

	########################################
	@node tx_m = SetActiveDie('M', may_use_ring = true) -> [next = tx_m_cond, no_die = tx_c]

	@node tx_m_cond = BinaryCondition(muster_cond) -> [n_true = tx_m_die, n_false = tx_c]
	@node tx_m_die = UseActiveDie() -> tx_m_action
	@node tx_m_action = PerformAction(muster_text) -> tx_m_end
	@node tx_m_end = End() -> []

	########################################
	@node tx_c = SetActiveDie('C', may_use_ring = true) -> [next = tx_c_cond, no_die = tx_return]

	@node tx_c_cond = BinaryCondition(character_cond) -> [n_true = tx_c_1, n_false = tx_return]
	@node tx_c_1 = JumpToGraph("character_move") -> tx_return

	########################################
	@node tx_return = ReturnFromGraph() -> []
end

max_stars_repo_name: mvmorin/queller-bot

max_issues_repo_name: mvmorin/queller-bot

max_forks_repo_name: mvmorin/queller-bot


text:
module Lowering
export lower_args_default

using ExprParsers
using Setfield
using WhereTraits.Utils

lower_args_default(expr::Expr) = lower_args_default(parse_expr(EP.Function(), expr))
function lower_args_default(func_parsed::EP.Function_Parsed)
  args_parsed = [parse_expr(EP.Arg(), a) for a in func_parsed.args]
  defaults = [a.default for a in args_parsed]
  @assert _nodefault_only_at_front(defaults) "no normal positional argument allowed after positional argument with default value"
  count_defaults = sum(d isa EP.NoDefault ? 0 : 1 for d in defaults)

  # collect several expressions

  # the whole function is returned without defaults
  args = [to_expr(@set a.default = EP.nodefault) for a in args_parsed]
  func_without_defaults = @set func_parsed.args = args
  # collect all lowerings
  lowerings = []
  for i in 1:count_defaults  # 1:0 is empty
    symbols = [a.name for a in args_parsed[1:end-i]]
    dropped_symbols = [a.name for a in args_parsed[end-i+1:end]]
    defaults = [a.default for a in args_parsed[end-i+1:end]]
    drop_unapplicable_traits(wheres) = filter(wheres) do expr
      !depends_on(expr, dropped_symbols)
    end

    lowering = EP.Function_Parsed(
      name = func_parsed.name,
      curlies = copy(func_parsed.curlies),
      args = [to_expr(@set a.default = EP.nodefault) for a in args_parsed[1:end-i]],
      kwargs = [],
      # the traits normalization kicks out all unused wheres where TypeVariables matter
      # but dependencies on args still need to be handled
      wheres = drop_unapplicable_traits(func_parsed.wheres),
      body = Expr(:call, func_parsed.name, symbols..., defaults...),
    )
    push!(lowerings, lowering)
  end
  func_without_defaults, lowerings
end

function _nodefault_only_at_front(defaults::Vector)
  if isempty(defaults)
    true
  else
    first = defaults[1]
    if first isa EP.NoDefault
      _nodefault_only_at_front(defaults[2:end])
    else
      _no_nodefault_any_more(defaults[2:end])
    end
  end
end

function _no_nodefault_any_more(defaults::Vector)
  if isempty(defaults)
    true
  else
    first = defaults[1]
    if first isa EP.NoDefault
      false
    else
      _no_nodefault_any_more(defaults[2:end])
    end
  end
end

end # module


text:
r=0.1

#theta(x::FloatOrVector,r::Float64)
#theta's vanishe in 0, are bounded by 1, increasing, concave
# positive for positive values, and negative for negative values.
@test MPCCSolver.theta(0.0, r) .== 0.0
@test MPCCSolver.theta(1.0+1e6, r) .<= 1
@test MPCCSolver.theta(-1e-6, r) .< 0.0
@test MPCCSolver.theta(1.0, r) .>= MPCCSolver.theta(0.5, r)

#dtheta(x::FloatOrVector,r::Float64)
#la dérivée des fonctions theta doit être positive
@test MPCCSolver.dtheta(1.0+1e6, r) >= 0.0
@test MPCCSolver.dtheta(-0.5, r)    >= 0.0

#ddtheta(x::FloatOrVector,r::Float64)
#la dérivée seconde des fonctions theta doit être négative
@test MPCCSolver.ddtheta(1.0+1e6, r) <= 0.0
@test MPCCSolver.ddtheta(-0.5, r)    <= 0.0

#invtheta(x::Vector,r::Float64)
@test abs(MPCCSolver.invtheta(MPCCSolver.theta(1.,r),r)-1.)      <= eps(Float64)*10.
@test abs(MPCCSolver.invtheta(MPCCSolver.theta(-1e-6,r),r)+1e-6) <= eps(Float64)*10.

#_poly_degree_2_positif
#Compute the smallest positive root of a second order polynomial
#or a negative/complex value if no positive real root.
#If no solution or infinite nb of solutions, return Inf.
@test abs(MPCCSolver._poly_degree_2_positif(0.5,-1.,0.25) - 0.2928932188134524) <= eps(Float64)
@test abs(MPCCSolver._poly_degree_2_positif(0.,1.,-0.5) - 0.5) <= eps(Float64)
@test !isfinite(MPCCSolver._poly_degree_2_positif(0.,0.,-0.))
@test !isfinite(MPCCSolver._poly_degree_2_positif(1.,1.,1.))

#alpha_theta_max(x::Float64,dx::Float64,y::Float64,dy::Float64,r::Float64,s::Float64,t::Float64)
#alpha_theta_max : résoud l'équation en alpha où psi(x,r,s,t)=s+t*theta(x,r)
r=1.0;s=1.0;t=0.0 #résoud l'éq.

x=1.0;y=2.0;dx=1.0;dy=-1.0;
#doit rendre (0,1)
@test isequal(MPCCSolver.alpha_theta_max(x,dx,y,dy,r,s,t), [0.0,1.0])

x=0.0;y=2.0;dx=1.0;dy=-2.0;
#solution : (1,0.5)

@test isequal(MPCCSolver.alpha_theta_max(x,dx,y,dy,r,s,t), [1.0,0.5])

r=1.0;s=1.0;t=1.0
#on tape la contrainte :
x=1.0;y=2.0;dx=1.0;dy=0.0;

@test minimum(MPCCSolver.alpha_theta_max(x,dx,y,dy,r,s,t))>=0.0
@test maximum(MPCCSolver.alpha_theta_max(x,dx,y,dy,r,s,t))==Inf

x=2.0;y=1.0;dx=-1.0;dy=2.0;

@test minimum(MPCCSolver.alpha_theta_max(x,dx,y,dy,r,s,t))>=0.0
@test maximum(MPCCSolver.alpha_theta_max(x,dx,y,dy,r,s,t))<Inf

#on est sous (s,s) et donc on ne tape pas la contrainte virtuelle
x=0.0;y=0.0;dx=0.0;dy=-2.0;

@test isequal(maximum(MPCCSolver.alpha_theta_max(x,dx,y,dy,r,s,t)),Inf)


text:
# Quadratic Discriminant Analysis

mutable struct QuadraticDiscriminantModel{T} <: DiscriminantModel{T}
    "Standard discriminant model parameters"
    Θ::DiscriminantParameters{T}
    "Shrinkage parameter towards the common covariance matrix"
    λ::Union{Nothing,T}
    "Covariance matrix for each class"
    Σₘ::Union{Nothing,Vector{AbstractMatrix{T}}}
    "Whitening transformation by class"
    Wₘ::Vector{AbstractMatrix{T}}
    "Discriminant function intercept for each class"
    δ₀::Vector{T}
    function QuadraticDiscriminantModel{T}() where T
        new{T}(DiscriminantParameters{T}(), nothing, Array{T}(undef,0,0,0),
               Vector{T}(undef,0))
    end
end

function fit!(QDA::QuadraticDiscriminantModel{T},
              y::Vector{<:Integer},
              X::Matrix{T};
              dims::Integer=1,
              canonical::Bool=false,
              store_cov::Bool=false,
              store_class_covs::Bool=false,
              centroids::Union{Nothing,AbstractMatrix}=nothing,
              priors::Union{Nothing,AbstractVector}=nothing,
              gamma::Union{Nothing,Real}=nothing,
              lambda::Union{Nothing,Real}=nothing) where T
    Θ = QDA.Θ

    if lambda !== nothing
        0 ≤ lambda ≤ 1 || throw(DomainError(lambda, "λ must be in the interval [0,1]"))
    end
    QDA.λ = lambda

    set_dimensions!(Θ, y, X, dims)
    set_gamma!(Θ, gamma)
    set_statistics!(Θ, y, X, centroids)
    set_priors!(Θ, priors)

    # Compute overall degrees of freedom
    df = size(X, dims) - Θ.m

    if dims == 1
        # Overall centroid is prior-weighted average of class centroids
        Θ.μ = vec(transpose(Θ.π)*Θ.M)
        # Center the data matrix with respect to classes to compute whitening matrix
        X .-= view(Θ.M, y, :)
        # Compute scatter matrix if requested
        S = (store_cov || QDA.λ !== nothing) ? transpose(X)*X : nothing
    else
        Θ.μ = Θ.M*Θ.π
        X .-= view(Θ.M, :, y)
        S = (store_cov || QDA.λ !== nothing) ? X*transpose(X) : nothing
    end

    # Set up some references
    Σₘ = QDA.Σₘ
    Wₘ = QDA.Wₘ
    δ₀ = QDA.δ₀

    # Allocate space for each covariance matrix if required
    QDA.Σₘ = store_class_covs ? [zeros(T, 0, 0) for k=1:QDA.m] : nothing

    # QDA - find whitening matrix for each class
    for (k, nₖ) in enumerate(Θ.nₘ)
        Xₖ = dims == 1 ? X[y .== k, :] : X[:, y .== k]
        dfₖ = nₖ - 1

        if store_class_covs || QDA.λ !== nothing
            Sₖ = dims == 1 ? transpose(Xₖ)*Xₖ : Xₖ*transpose(Xₖ)

            if QDA.λ === nothing
                # Compute covariance matrix using within-class degrees of freedom
                Σₖ = lmul!(one(T)/dfₖ, Sₖ)
                Σₘ[k] = deepcopy(Σₖ)

                # Compute whitening
                Wₘ[k], δ₀[k] = whiten_cov_chol!(Σₖ, Θ.γ, dims=dims)
            else
                if store_class_covs
                    QDA.Σₘ[k] = lmul!(one(T)/dfₖ, deepcopy(Sₖ))
                end

                # Regularize covariance matrix
                scale = one(T)/regularize(dfₖ, df, QDA.λ)  # denominator is separate
                Σₖ = lmul!(scale, broadcast!(regularize, Sₖ, Sₖ, S, QDA.λ))

                # Compute whitening
                Wₘ[k], δ₀[k] = whiten_cov_chol!(Σₖ, Θ.γ, dims=dims)
            end
        else
            if Θ.γ === nothing
                Wₘ[k], δ₀[k] = whiten_data_chol!(Xₖ, dims=dims, df=dfₖ)
            else
                Wₘ[k], δ₀[k] = whiten_data_svd!(Xₖ, Θ.γ, dims=dims, df=dfₖ)
            end
        end
    end

    if store_cov
        Θ.Σ = lmul!(one(T)/df, S)
    end

    Θ.fit = true

    return LDA
end

# discriminants

#==========================================================================
  Regularized Quadratic Discriminant Analysis Solver
==========================================================================#

#immutable ModelQDA{T<:BlasReal}
#    order::Union{Type{Val{:row}},Type{Val{:col}}}
#    W_k::Vector{AbstractMatrix{T}}  # Vector of class whitening matrices
#    M::Matrix{T}                    # Matrix of class means (one per row)
#    priors::Vector{T}               # Vector of class priors
#    gamma::Nullable{T}
#    lambda::Nullable{T}
#end
#
#function show(io::IO, model::ModelQDA)
#    println(io, "Quadratic Discriminant Model")
#    println(io, "\nRegularization Parameters:")
#    print(io, "    γ = ")
#    isnull(model.gamma) ? print(io, "N/A") : showcompact(get(model.gamma))
#    print(io, "\n    λ = ")
#    isnull(model.lambda) ? print(io, "N/A") : showcompact(get(model.lambda))
#    println(io, "\n\nClass Priors:")
#    for i in eachindex(model.priors)
#        print(io, "    Class ", i, " Probability: ")
#        showcompact(io, model.priors[i]*100)
#        println(io, "%")
#    end
#    println(io, "\nClass Means:")
#    for j in eachindex(model.priors)
#        print(io, "    Class ", j, " Mean: [")
#        print(io, join([sprint(showcompact, v) for v in subvector(model.order, model.M, j)], ", "))
#        println(io, "]")
#    end
#    print(io, "\n")
#end
#
#function covmatrix{T<:BlasReal}(::Type{Val{:row}}, H::Matrix{T})
#    n = size(H,1)
#    Σ = H'H
#    broadcast!(/, Σ, Σ, n-1)
#end
#
#function covmatrix{T<:BlasReal}(::Type{Val{:col}}, H::Matrix{T})
#    n = size(H,2)
#    Σ = A_mul_Bt(H, H)
#    broadcast!(/, Σ, Σ, n-1)
#end
#
#for (scheme, dim_obs) in ((:(:row), 1), (:(:col), 2))
#    isrowmajor = dim_obs == 1
#    dim_param = isrowmajor ? 2 : 1
#
#    H_i      = isrowmajor ? :(H[y .== i, :]) : :(H[:, y .== i])
#    _ij, _ji = isrowmajor ? (:i, :j) : (:j, :i)  # Swapped variables for row and column ordering
#    _nk, _kn = isrowmajor ? (:n, :k) : (:k, :n)
#    _np, _pn = isrowmajor ? (:n, :p) : (:p, :n)
#    W_j, H   = isrowmajor ? (:(W_k[j]), :H) : (:H, :(W_k[j]))
#
#    @eval begin
#        # λ-regularized QDA - require full covariance matrices
#        function classwhiteners!{T<:BlasReal,U}(
#                 ::Type{Val{$scheme}},
#                H::Matrix{T},
#                y::RefVector{U},
#                γ::Nullable{T},
#                λ::T
#            )
#            k = convert(Int64, y.k)
#            Σ_k = AbstractMatrix{T}[covmatrix(Val{$scheme}, $H_i) for i = 1:k]
#            Σ   = covmatrix(Val{$scheme}, H)
#            for i in eachindex(Σ_k)
#                regularize!(Σ_k[i], λ, Σ)
#                Σ_k[i] = whitencov_chol!(Val{$scheme}, Σ_k[i], γ)
#            end
#            Σ_k
#        end
#
#        # No λ-regularization - only need data matrices
#        function classwhiteners!{T<:BlasReal,U}(
#                 ::Type{Val{$scheme}},
#                H::Matrix{T},
#                y::RefVector{U},
#                γ::Nullable{T}
#            )
#            k = convert(Int64, y.k)
#            if isnull(γ)
#                AbstractMatrix{T}[whitendata_qr!(Val{$scheme}, $H_i) for i = 1:k]
#            else
#                AbstractMatrix{T}[whitendata_svd!(Val{$scheme}, $H_i, get(γ)) for i = 1:k]
#            end
#        end
#
#        function qda!{T<:BlasReal,U}(
#                 ::Type{Val{$scheme}},
#                X::Matrix{T},
#                M::Matrix{T},
#                y::RefVector{U},
#                γ::Nullable{T},
#                λ::Nullable{T}
#            )
#            H = centerclasses!(Val{$scheme}, X, M, y)
#            if isnull(λ)
#                classwhiteners!(Val{$scheme}, H, y, γ)
#            else
#                classwhiteners!(Val{$scheme}, H, y, γ, get(λ))
#            end
#        end
#
#        function discriminants_qda{T<:BlasReal}(
#                 ::Type{Val{$scheme}},
#                W_k::Vector{AbstractMatrix{T}},
#                M::Matrix{T},
#                priors::Vector{T},
#                Z::Matrix{T}
#            )
#            n = size(Z, $dim_obs)
#            p = size(Z, $dim_param)
#            k = length(priors)
#            D   = Array(T, $_nk, $_kn)  # discriminant function values
#            H   = Array(T, $_np, $_pn)  # temporary array to prevent re-allocation k times
#            hᵀh = Array(T, n)           # diag(H'H)
#            for j = 1:k
#                broadcast!(-, H, Z, subvector(Val{$scheme}, M, j))
#                dotvectors!(Val{$scheme}, $H * $W_j, hᵀh)
#                for i = 1:n
#                    D[$_ij, $_ji] = -hᵀh[i]/2 + log(priors[j])
#                end
#            end
#            D
#        end
#    end
#end
#
#doc"`qda(X, y; order, M, lambda, gamma, priors)` Fits a regularized quadratic discriminant model to
#the data in `X` based on class identifier `y`."
#function qda{T<:BlasReal,U<:Integer}(
#        X::Matrix{T},
#        y::AbstractVector{U};
#        order::Union{Type{Val{:row}},Type{Val{:col}}} = Val{:row},
#        M::Matrix{T} = classmeans(order, X, RefVector(y)),
#        gamma::Union{T,Nullable{T}}  = Nullable{T}(),
#        lambda::Union{T,Nullable{T}} = Nullable{T}(),
#        priors::Vector{T} = ones(T,maximum(y))/maximum(y)
#    )
#    all(priors .> 0) || error("Argument priors must have positive values only")
#    isapprox(sum(priors), one(T)) || error("Argument priors must sum to 1")
#    yref = isa(y,RefVector) ? y : RefVector(y)
#    length(priors) == yref.k || throw(DimensionMismatch("Prior length must match class count"))
#    γ = isa(gamma,  Nullable) ? gamma  : Nullable(gamma)
#    λ = isa(lambda, Nullable) ? lambda : Nullable(lambda)
#    W_k = qda!(order, copy(X), M, yref, γ, λ)
#    ModelQDA{T}(order, W_k, M, priors, γ, λ)
#end
#
#doc"`discriminants(Model, Z)` Uses `Model` on input `Z` to product the class discriminants."
#function discriminants{T<:BlasReal}(mod::ModelQDA{T}, Z::Matrix{T})
#    discriminants_qda(mod.order, mod.W_k, mod.M, mod.priors, Z)
#end
#
#doc"`classify(Model, Z)` Uses `Model` on input `Z`."
#function classify{T<:BlasReal}(mod::ModelQDA{T}, Z::Matrix{T})
#    vec(mapslices(indmax, discriminants(mod, Z), mod.order == Val{:row} ? 2 : 1))
#end


text:
function re_smp_deci(x::Neuro1DRealSignal; rate::Real)
    r, c = size(x.signal)

    resampled_data::Vector{Vector{Float64}} = []

    for i = 1:c
        append!(resampled_data, [resample(x.signal[:, i], rate)])
    end

    new_r = minimum(length.(resampled_data))
    result = zeros(new_r, c)

    for i = 1:c
        result[:, c] = resampled_data[i][1:new_r]
    end

    Neuro1DRealSignal(
        signal = result,
        sample_rate = x.sample_rate * rate
    )
end

function re_smp(
    x::Neuro1DRealSignal,
    method = :decimate;
    kwargs...)

    if method == :decimate
        re_smp_deci(x; kwargs...)
    else
        throw(ArgumentError("invalid method"))
    end
end