text:
function ticks_and_labels end

module Formatters

    using Showoff

    function scientific(ticks::AbstractVector)
        Showoff.showoff(ticks, :scientific)
    end

    function plain(ticks::AbstractVector)
        try
            Showoff.showoff(ticks, :plain)
        catch e
            Base.showerror(stderr, e)
            println("with ticks: ", ticks)
            String["-Inf", "Inf"]
        end
    end

end
using .Formatters

"""
    $(SIGNATURES)

Plots a 2-dimensional axis.

## Theme
$(ATTRIBUTES)
"""
@recipe(Axis2D) do scene
    Theme(
        visible = true,

        showgrid = true,
        showticks = true,
		showtickmarks = true,
        padding = 0.1,

        ticks = Theme(

            ranges_labels = (automatic, automatic),
            formatter = Formatters.plain,

            gap = 3,
            title_gap = 3,

            linewidth = (1, 1),
            linecolor = ((:black, 0.4), (:black, 0.4)),
            linestyle = (nothing, nothing),

            textcolor = (:black, :black),
            textsize = (5, 5),
            rotation = (0.0, 0.0),
            align = ((:center, :top), (:right, :center)),
            font = lift(dim2, theme(scene, :font)),
        ),
		tickmarks = Theme(
			length = (3.0, 3.0),
			linewidth = (1,1),
			linecolor = ((:black, 0.4), (:black, 0.2)),
			linestyle = (nothing, nothing)
		),

        grid = Theme(
            linewidth = (0.5, 0.5),
            linecolor = ((:black, 0.3), (:black, 0.3)),
            linestyle = (nothing, nothing),
        ),

        frame = Theme(
            linewidth = 1.0,
            linecolor = :black,
            linestyle = nothing,
            axis_position = nothing,
            axis_arrow = false,
            arrow_size = 2.5,
            frames = ((false, false), (false, false)),
        ),

        names = Theme(
            axisnames = ("x", "y"),
            title = nothing,
            textcolor = (:black, :black),
            textsize = (6, 6),
            rotation = (0.0, -1.5pi),
            align = ((:center, :top), (:center, :bottom)),
            font = lift(dim2, theme(scene, :font)),
        )
    )
end

"""
    $(SIGNATURES)

Plots a 3-dimensional Axis.

## Theme
$(ATTRIBUTES)
"""
@recipe(Axis3D) do scene

    q1 = qrotation(Vec3f0(1, 0, 0), -0.5f0*pi)
    q2 = qrotation(Vec3f0(0, 0, 1), 1f0*pi)
    tickrotations3d = (
        qrotation(Vec3f0(0,0,1), -1.5pi),
        q2,
        qrotation(Vec3f0(1, 0, 0), -0.5pi) * q2
    )
    axisnames_rotation3d = tickrotations3d
    tickalign3d = (
        (:left, :center), # x axis
        (:right, :center), # y axis
        (:right, :center), # z axis
    )
    axisnames_align3d = tickalign3d
    tick_color = RGBAf0(0.5, 0.5, 0.5, 0.6)
    grid_color = RGBAf0(0.5, 0.5, 0.5, 0.4)
    grid_thickness = 1
    gridthickness = ntuple(x-> 1f0, Val(3))
    tsize = 5 # in percent
    Theme(
        visible = true,
        showticks = (true, true, true),
        showaxis = (true, true, true),
        showgrid = (true, true, true),
        scale = Vec3f0(1),
        padding = 0.1,
        names = Theme(
            axisnames = ("x", "y", "z"),
            textcolor = (:black, :black, :black),
            rotation = axisnames_rotation3d,
            textsize = (6.0, 6.0, 6.0),
            align = axisnames_align3d,
            font = lift(dim3, theme(scene, :font)),
            gap = 1
        ),

        ticks = Theme(
            ranges_labels = (automatic, automatic),
            formatter = Formatters.plain,

            textcolor = (tick_color, tick_color, tick_color),

            rotation = tickrotations3d,
            textsize =  (tsize, tsize, tsize),
            align = tickalign3d,
            gap = 1,
            font = lift(dim3, theme(scene, :font)),
        ),

        frame = Theme(
            linecolor = (grid_color, grid_color, grid_color),
            linewidth = (grid_thickness, grid_thickness, grid_thickness),
            axiscolor = (:black, :black, :black),
        )
    )
end

isaxis(x) = false
isaxis(x::Union{Axis2D, Axis3D}) = true

const Limits{N} = NTuple{N, Tuple{Number, Number}}

default_ticks(limits::Limits, ticks, scale_func::Function) = default_ticks.(limits, (ticks,), scale_func)
default_ticks(limits::Tuple{Number, Number}, ticks, scale_func::Function) = default_ticks(limits..., ticks, scale_func)

function default_ticks(lmin::Number, lmax::Number, ticks::AbstractVector{<: Number}, scale_func::Function)
    scale_func.((filter(t -> lmin <= t <= lmax, ticks)))
end
function default_ticks(lmin::Number, lmax::Number, ::Automatic, scale_func::Function)
    # scale the limits
    scaled_ticks, mini, maxi = optimize_ticks(
        scale_func(lmin),
        scale_func(lmax);
        k_min = 4, # minimum number of ticks
        k_max = 8, # maximum number of ticks
    )
    length(scaled_ticks) == 1 && isnan(scaled_ticks[1]) && return [-Inf, Inf]
    scaled_ticks
end

function default_ticks(lmin::Number, lmax::Number, n::Integer, scale_func = identity)
    scaled_ticks, mini, maxi = optimize_ticks(
        scale_func(lmin),
        scale_func(lmax);
        k_min = ticks, # minimum number of ticks
        k_max = ticks, # maximum number of ticks
        k_ideal = ticks,
        # `strict_span = false` rewards cases where the span of the
        # chosen  ticks is not too much bigger than amin - amax:
        strict_span = false,
    )
    scaled_ticks
end

default_ticks(x::Automatic, limits::Tuple, n) = default_ticks(limits, n, identity)
function default_ticks(ticks::Tuple, limits::Tuple, n::Tuple)
    default_ticks.(ticks, (limits,), n)
end

default_ticks(ticks::Tuple, limits::Limits, n) = default_ticks.(ticks, limits, (n,))

default_ticks(ticks::Tuple, limits::Limits, n::Tuple) = default_ticks.(ticks, limits, n)

default_ticks(ticks::AbstractVector{<: Number}, limits, n) = ticks

function default_labels(x::NTuple{N, Any}, formatter::Function) where N
    default_labels.(x, formatter)
end

function default_labels(x::AbstractVector, y::AbstractVector, formatter::Function = Formatters.plain)
    default_labels.((x, y), formatter)
end

function default_labels(ticks::AbstractVector, formatter::Function = Formatters.plain)
    if applicable(formatter, ticks)
        return formatter(ticks) # takes the whole array
    elseif applicable(formatter, first(ticks))
        return formatter.(ticks)
    else
        error("Formatting function $(formatter) is neither applicable to $(typeof(ticks)) nor $(eltype(ticks)).")
    end
end

default_labels(x::Automatic, ranges, formatter) = default_labels(ranges, formatter)
default_labels(x::Tuple, ranges::Tuple, formatter) = default_labels.(x, ranges, (formatter,))
default_labels(x::Tuple, ranges, formatter) = default_labels.(x, (ranges,), (formatter,))
default_labels(x::AbstractVector{<: AbstractString}, ranges, formatter::Function) = x
default_labels(x::AbstractVector{<: AbstractString}, ranges::AbstractVector, formatter::Function) = x

function convert_arguments(::Type{<: Axis2D}, limits::Rect)
    e = (minimum(limits), maximum(limits))
    (((e[1][1], e[2][1]), (e[1][2], e[2][2])),)
end
function convert_arguments(::Type{<: Axis3D}, limits::Rect)
    e = (minimum(limits), maximum(limits))
    (((e[1][1], e[2][1]), (e[1][2], e[2][2]), (e[1][3], e[2][3])),)
end
a_length(x::AbstractVector) = length(x)
a_length(x::Automatic) = x

function calculated_attributes!(::Type{<: Union{Axis2D, Axis3D}}, plot)
    ticks = plot.ticks
    args = (plot.padding, plot[1], ticks.ranges, ticks.labels, ticks.formatter)
    ticks[:ranges_labels] = lift(args...) do pad, lims, ranges, labels, formatter
        limit_widths = map(x-> x[2] - x[1], lims)
        pad = (limit_widths .* pad)
        # pad the drawn limits and use them as the ranges
        lim_pad = map((lim, p)-> (lim[1] - p, lim[2] + p), lims, pad)
        num_ticks = labels === automatic ? automatic : a_length.(labels)
        ranges = default_ticks(ranges, lim_pad, num_ticks)
        labels = default_labels(labels, ranges, formatter)
        (ranges, labels)
    end
    return
end

function draw_ticks(
        textbuffer, dim, origin, ticks,
        linewidth, linecolor, linestyle,
        textcolor, textsize, rotation, align, font
    )
    for (tick, str) in ticks
        pos = ntuple(i-> i != dim ? origin[i] : tick, Val(2))
        push!(
            textbuffer,
            str, pos,
            rotation = rotation[dim], textsize = textsize[dim],
            align = align[dim], color = textcolor[dim], font = font[dim]
        )
    end
end

function draw_tickmarks(
						linebuffer, dim, origin, ticks,dir::NTuple{N},
						linewidth, linecolor, linestyle,
						textcolor, textsize, rotation,align,font
						) where N
	for (tick, str) in ticks
		pos = ntuple(i-> i != dim ? origin[i] : tick, Val(2))
		posf0 = Point2f0(pos)
		dirf0  = Pointf0{N}(dir)
		append!(linebuffer,
		      [posf0, posf0 .+ dirf0],
		      color = linecolor[dim], linewidth = linewidth[dim]
		      )
	end
end

function draw_grid(
        linebuffer, dim, origin, ticks, dir::NTuple{N},
        linewidth, linecolor, linestyle
    ) where N
    dirf0 = Pointf0{N}(dir)
    for (tick, str) in ticks
        tup = ntuple(i-> i != dim ? origin[i] : tick, Val(N))
        posf0 = Pointf0{N}(tup)
        append!(
            linebuffer,
            [posf0, posf0 .+ dirf0],
            color = linecolor[dim], linewidth = linewidth[dim]#, linestyle = linestyle[dim]
        )
    end
end

function draw_frame(
        linebuffer, limits::NTuple{N, Any},
        linewidth, linecolor, linestyle,
        axis_position, axis_arrow, arrow_size, frames
    ) where N

    mini = minimum.(limits)
    maxi = maximum.(limits)
    rect = HyperRectangle(Vec(mini), Vec(maxi .- mini))
    origin = Vec{N}(0.0)
    if (origin in rect) && axis_position == :origin
        for i = 1:N
            start = unit(Point{N, Float32}, i) * Float32(mini[i])
            to = unit(Point{N, Float32}, i) * Float32(maxi[i])
            if false #axis_arrow
                arrows(
                    scene, [start, to],
                    linewidth = linewidth, linecolor = linecolor, linestyle = linestyle,
                    arrowsize = arrow_size
                )
            else
                append!(
                    linebuffer, [start, to],
                    linewidth = linewidth, color = linecolor #linestyle = linestyle,
                )
            end
        end
    end
    limit_widths = maxi .- mini
    frames = convert_attribute(frames, key"frames"())
    for side in 1:2
        from = Point{N, Float32}(getindex.(limits, side))
        # if axis is drawn at origin, and we draw frame from origin,
        # we already did this
        if !(from == origin && axis_position == :origin)
            for otherside in 1:2
                for dim in 1:N
                    if frames[N-dim+1][3-otherside]
                        p = ntuple(i-> i == dim ? limits[i][otherside] : limits[i][side], Val(N))
                        to = Point{N, Float32}(p)
                        append!(
                            linebuffer, [from, to],
                            linewidth = linewidth, color = linecolor#, linestyle = linestyle,
                        )
                    end
                end
            end
        end
    end
end

function draw_titles(
        textbuffer,
        xticks, yticks, origin, limit_widths,
        tickfont, tick_size, tick_gap, tick_title_gap,
        axis_labels,
        textcolor, textsize, rotation, align, font,
        title
    )
    tickspace_x = maximum(map(yticks) do tick
        str = last(tick)
        tick_bb = text_bb(str, to_font(tickfont[2]), tick_size[2])
        widths(tick_bb)[1]
    end)

    tickspace_y = widths(text_bb(
        last(first(xticks)), to_font(tickfont[1]), tick_size[1]
    ))[2]

    model_inv = inv(transformationmatrix(textbuffer)[])

    tickspace = transform(model_inv, (tickspace_x, tickspace_y))
    title_start = origin .- (tick_gap .+ tickspace .+ tick_title_gap)
    half_width = origin .+ (limit_widths ./ 2.0)

    posx = (half_width[1], title_start[2])
    posy = (title_start[1], half_width[2])
    positions = (posx, posy)
    for i = 1:2
        if !isempty(axis_labels[i])
            push!(
                textbuffer, axis_labels[i], positions[i],
                textsize = textsize[i], align = align[i], rotation = rotation[i],
                color = textcolor[i], font = font[i]
            )
        end
    end

    if title !== nothing
        # TODO give title own text attributes
        push!(
            textbuffer, title, (half_width[1], (origin .+ (limit_widths))[2]),
            textsize = textsize[1], align = (:center, :top), rotation = rotation[1],
            color = textcolor[1], font = font[1]
        )
    end
end

function ticks_and_labels(x)
    st, s = extrema(x)
    r = LinRange(st, s, 4)
    zip(r, string.(round.(r, 4)))
end

function transform(model::Mat4, x::T) where T
    x4d = to_ndim(Vec4f0, x, 0.0)
    to_ndim(T, model * x4d, 0.0)
end
un_transform(model::Mat4, x) = transform(inv(model), x)

to2tuple(x) = ntuple(i-> x, Val(2))
to2tuple(x::Tuple{<:Any, <: Any}) = x

function draw_axis2d(
        textbuffer,
        frame_linebuffer, grid_linebuffer,tickmarks_linebuffer,
        m, padding, limits, xyrange_labels,
        showgrid, showticks,showtickmarks,
        # grid attributes
        g_linewidth, g_linecolor, g_linestyle,

        # tick attributes
        t_linewidth, t_linecolor, t_linestyle,
        t_textcolor, t_textsize, t_rotation, t_align, t_font,
        t_gap, t_title_gap,

	#tickmark attribures
	tm_linewidth, tm_linecolor, tm_linestyle,
	tm_length,

        # frame attributes
        f_linewidth, f_linecolor, f_linestyle,
        f_axis_position, f_axis_arrow, f_arrow_size, f_frames,

        # title / axis name attributes
        ti_labels,
        ti_textcolor, ti_textsize, ti_rotation, ti_align, ti_font, ti_title
    )
    xyrange, labels = xyrange_labels
    start!(textbuffer); start!(frame_linebuffer); foreach(start!, grid_linebuffer)
    foreach(start!, tickmarks_linebuffer)
    # limits = limits Vec2f0(padding)
    # limits ((xmin, xmax), (ymin, ymax))
    limit_widths = map(x-> x[2] - x[1], limits)
    pad = (limit_widths .* to2tuple(padding))
    # pad the drawn limits
    limits = map((lim, p)-> (lim[1] - p, lim[2] + p), limits, pad)

    # recalculate widths
    limit_widths = map(x-> x[2] - x[1], limits)

    % = mean(limit_widths) / 100 # percentage
    xyticks = zip.(xyrange, labels)
    model_inv = inv(transformationmatrix(textbuffer)[])

    ti_textsize = ti_textsize .* %
    t_textsize = t_textsize .* %
    t_gap = transform(model_inv, to2tuple(t_gap .* %))
    tm_length = transform(model_inv, to2tuple(tm_length .* %))
    t_title_gap = transform(model_inv, to2tuple(t_title_gap .* %))

    origin = first.(limits)
    dirs = ((0.0, Float64(limit_widths[2])), (Float64(limit_widths[1]), 0.0))
    foreach(1:2, dirs, xyticks) do dim, dir, ticks
        if showgrid[dim]
            draw_grid(
                grid_linebuffer[dim], dim, origin, ticks, dir,
                g_linewidth, g_linecolor, g_linestyle
            )
        end
    end
    tm_offsets = ((0.0, Float64(tm_length[2])), (Float64(tm_length[1]), Float64(0.0)))
    foreach(1:2, tm_offsets, xyticks) do dim, offset, ticks
	    if showtickmarks[dim]
		    draw_tickmarks(
						   tickmarks_linebuffer[dim], dim, origin .- offset, ticks, offset,
						   t_linewidth, t_linecolor, t_linestyle,
						   t_textcolor, t_textsize, t_rotation, t_align, t_font
						   )
	    end
	end
	t_offsets = ((0.0, Float64(t_gap[2] + tm_length[2])), (Float64(t_gap[1]+tm_length[1]), Float64(0.0)))
    foreach(1:2, t_offsets, xyticks) do dim, offset, ticks
        if showticks[dim]
            draw_ticks(
                textbuffer, dim, origin .- offset, ticks,
                t_linewidth, t_linecolor, t_linestyle,
                t_textcolor, t_textsize, t_rotation, t_align, t_font
            )
	    end
    end

    draw_frame(
        frame_linebuffer, limits,
        f_linewidth, f_linecolor, f_linestyle,
        f_axis_position, f_axis_arrow, f_arrow_size,
        f_frames,
    )

    draw_titles(
        textbuffer, xyticks..., origin, limit_widths,
        t_font, t_textsize, t_gap, t_title_gap,
        ti_labels,
        ti_textcolor, ti_textsize, ti_rotation, ti_align, ti_font,
        ti_title
    )
    finish!(textbuffer); finish!(frame_linebuffer); foreach(finish!, grid_linebuffer)
	foreach(finish!, tickmarks_linebuffer)
    return
end

# for axis, we don't want to have plot!(scene, args) called on it, so we need to overload it directly
function plot!(scene::SceneLike, ::Type{<: Axis2D}, attributes::Attributes, args...)
    # create "empty" plot type - empty meaning containing no plots, just attributes + arguments
    cplot = Axis2D(scene, attributes, args)
    g_keys = (:linewidth, :linecolor, :linestyle)
    f_keys = (:linewidth, :linecolor, :linestyle, :axis_position, :axis_arrow, :arrow_size, :frames)
    t_keys = (
        :linewidth, :linecolor, :linestyle,
        :textcolor, :textsize, :rotation, :align, :font,
        :gap, :title_gap
    )
    tm_keys = (:linewidth, :linecolor, :linestyle, :length)
    ti_keys = (:axisnames, :textcolor, :textsize, :rotation, :align, :font, :title)

    g_args = getindex.(cplot.grid, g_keys)
    f_args = getindex.(cplot.frame, f_keys)
    t_args = getindex.(cplot.ticks, t_keys)
    tm_args = getindex.(cplot.tickmarks, tm_keys)
    ti_args = getindex.(cplot.names, ti_keys)

    textbuffer = TextBuffer(cplot, Point{2})

    grid_linebuffer = to_node((
        LinesegmentBuffer(
            cplot, Point{2}, transparency = true,
            linestyle = lift(first, cplot[:grid, :linestyle])
        ),
        LinesegmentBuffer(
            cplot, Point{2}, transparency = true,
            linestyle = lift(last, cplot[:grid, :linestyle])
        )
    ))
    frame_linebuffer = to_node(LinesegmentBuffer(
        cplot, Point{2}, transparency = true, linestyle = cplot.frame.linestyle
    ))
    tickmarks_linebuffer = to_node((
				    LinesegmentBuffer(
						      cplot, Point{2}, transparency = true,
						      linestyle=lift(first, cplot[:tickmarks, :linestyle])
						      ),
				    LinesegmentBuffer(
						      cplot, Point{2}, transparency = true,
						      linestyle = lift(last, cplot[:tickmarks, :linestyle])
						      )
    ))
    map_once(
        draw_axis2d,
        to_node(textbuffer),
        frame_linebuffer, grid_linebuffer,tickmarks_linebuffer,
        transformationmatrix(scene),
        cplot.padding, cplot[1], cplot.ticks.ranges_labels,
        lift.((dim2,), (cplot.showgrid, cplot.showticks, cplot.showtickmarks))...,
        g_args..., t_args..., tm_args..., f_args..., ti_args...
    )
    push!(scene, cplot)
    return cplot
end

function labelposition(ranges, dim, dir, tgap, origin::StaticVector{N}) where N
    a, b = extrema(ranges[dim])
    whalf = Float32(((b - a) / 2))
    halfaxis = unit(Point{N, Float32}, dim) .* whalf

    origin .+ (halfaxis .+ (normalize(dir) * tgap))
end

function GeometryTypes.widths(x::AbstractRange)
    mini, maxi = Float32.(extrema(x))
    maxi - mini
end

_widths(x::Tuple{<: Number, <: Number}) = x[2] - x[1]
_widths(x) = Float32(maximum(x) - minimum(x))

to3tuple(x::Tuple{Any}) = (x[1], x[1], x[1])
to3tuple(x::Tuple{Any, Any}) = (x[1], x[2], x[2])
to3tuple(x::Tuple{Any, Any, Any}) = x
to3tuple(x) = ntuple(i-> x, Val(3))

function draw_axis3d(textbuffer, linebuffer, limits, ranges_labels, args...)
    # make sure we extend all args to 3D
    ranges, labels = ranges_labels
    args3d = to3tuple.(args)
    (
        showaxis, showticks, showgrid,
        axisnames, axisnames_color, axisnames_size, axisrotation, axisalign,
        axisnames_font, titlegap,
        gridcolors, gridthickness, axiscolors,
        ttextcolor, trotation, ttextsize, talign, tfont, tgap
    ) = args3d # splat to names

    N = 3
    start!(textbuffer); start!(linebuffer)
    ranges_ticks = Pair.(ranges, labels)
    mini, maxi = first.(limits), last.(limits)

    ranges = map(i-> [mini[i]; first(ranges_ticks[i]); maxi[i]], 1:3)
    ticklabels = map(x-> [""; last(x); ""], ranges_ticks)
    origin = Point{N, Float32}(mini)
    limit_widths = maxi .- mini
    % = minimum(limit_widths) / 100 # percentage
    ttextsize = (%) .* ttextsize
    axisnames_size = (%) .* axisnames_size

    titlegap = (%) .* titlegap
    tgap = (%) .* tgap
    for i = 1:N
        axis_vec = unit(Point{N, Float32}, i)
        width = _widths(ranges[i])
        stop = origin .+ (width .* axis_vec)
        if showaxis[i]
            append!(linebuffer, [origin, stop], color = axiscolors[i], linewidth = 1.5f0)
        end
        if showticks[i]
            range = ranges[i]
            j = mod1(i + 1, N)
            tickdir = unit(Point{N, Float32}, j)
            tickdir, offset2 = if i != 2
                tickdir = unit(Vec{N, Float32}, j)
                tickdir, Float32(_widths(ranges[j]) + tgap[i]) * tickdir
            else
                tickdir = unit(Vec{N, Float32}, 1)
                tickdir, Float32(_widths(ranges[1]) + tgap[i]) * tickdir
            end
            for (j, tick) in enumerate(range)
                labels = ticklabels[i]
                if length(labels) >= j
                    str = labels[j]
                    if !isempty(str)
                        startpos = (origin .+ ((Float32(tick - range[1]) * axis_vec)) .+ offset2)
                        push!(
                            textbuffer, str, startpos,
                            color = ttextcolor[i], rotation = trotation[i],
                            textsize = ttextsize[i], align = talign[i], font = tfont[i]
                        )
                    end
                end
            end
            if !isempty(axisnames[i])
                tick_widths = if length(ticklabels[i]) >= 3
                    widths(text_bb(ticklabels[i][end-1], to_font(tfont[i]), ttextsize[i]))[1]
                else
                    0f0
                end
                pos = (labelposition(ranges, i, tickdir, titlegap[i] + tick_widths, origin) .+ offset2)
                push!(
                    textbuffer, to_latex(axisnames[i]), pos,
                    textsize = axisnames_size[i], color = axisnames_color[i],
                    rotation = axisrotation[i], align = axisalign[i], font = axisnames_font[i]
                )
            end
        end
        if showgrid[i]
            c = gridcolors[i]
            thickness = gridthickness[i]
            for _j = (i + 1):(i + N - 1)
                j = mod1(_j, N)
                dir = unit(Point{N, Float32}, j)
                range = ranges[j]
                for tick in drop(range, 1)
                    offset = Float32(tick - range[1]) * dir
                    append!(
                        linebuffer, [origin .+ offset, stop .+ offset],
                        color = c, linewidth = thickness
                    )
                end
            end
        end
        finish!(textbuffer); finish!(linebuffer)
    end
    return
end

function plot!(scene::SceneLike, ::Type{<: Axis3D}, attributes::Attributes, args...)
    axis = Axis3D(scene, attributes, args)
    textbuffer = TextBuffer(axis, Point{3}, transparency = true)
    linebuffer = LinesegmentBuffer(axis, Point{3}, transparency = true)

    tstyle, ticks, frame = to_value.(getindex.(axis, (:names, :ticks, :frame)))
    titlevals = getindex.(tstyle, (:axisnames, :textcolor, :textsize, :rotation, :align, :font, :gap))
    framevals = getindex.(frame, (:linecolor, :linewidth, :axiscolor))
    tvals = getindex.(ticks, (:textcolor, :rotation, :textsize, :align, :font, :gap))
    args = (
        getindex.(axis, (:showaxis, :showticks, :showgrid))...,
        titlevals..., framevals..., tvals...
    )
    map_once(
        draw_axis3d,
        Node(textbuffer), Node(linebuffer),
        axis[1], axis.ticks.ranges_labels, args...
    )
    push!(scene.plots, axis)
    return axis
end


text:
module TestRDA
using Test
using DataFrames
using RData

@testset "Loading RData files (version=$ver)" for ver in (2, 3)
    rdata_path = joinpath(dirname(@__FILE__), "data_v$ver")

    @testset "Reading minimal RData" begin
        df = DataFrame(num = [1.1, 2.2])
        min_rda = load(joinpath(rdata_path, "minimal.rda"), convert=false)
        rdf = min_rda["df"]
        @test rdf isa RData.RList
        @testset "class() and inherits()" begin
            # not SEXP
            @test_throws MethodError RData.class(5)
            @test_throws MethodError RData.inherits(5, "number")
            @test_throws MethodError RData.inherits(5, "number")
            @test_throws MethodError RData.inherits(5, ["number"])

            rnotobj = RData.RBuiltin("test") # not a ROBJ
            @inferred RData.class(rnotobj)
            @inferred RData.inherits(rnotobj, "dummy")
            @test RData.class(rnotobj) === RData.emptystrvec
            @test !RData.inherits(rnotobj, "dummy")

            @inferred RData.class(rdf)
            @inferred RData.inherits(rdf, "data.frame")
            @inferred RData.inherits(rdf, ["data.frame"])
            @test RData.class(rdf) == ["data.frame"]
            @test RData.inherits(rdf, "data.frame")
            @test RData.inherits(rdf, ["data.frame"])

            rnumvec = rdf.data[1]
            @test rnumvec isa RData.RNumericVector
            @test RData.class(rnumvec) != ["data.frame"]
            @test !RData.inherits(rnumvec, "data.frame")
            @test !RData.inherits(rnumvec, ["data.frame"])
        end
        @test sexp2julia(min_rda["df"]) == df
        @test load(joinpath(rdata_path, "minimal.rda"), convert=true)["df"] == df
        @test load(joinpath(rdata_path, "minimal_ascii.rda"))["df"] == df
    end

    @testset "Conversion to Julia types" begin
        df = DataFrame(num = [1.1, 2.2],
                       int = Int32[1, 2],
                       logi = [true, false],
                       chr = ["ab", "c"],
                       factor = categorical(["ab", "c"], true),
                       cplx = [1.1+0.5im, 1.0im])
        rdf = sexp2julia(load(joinpath(rdata_path, "types.rda"), convert=false)["df"])
        @test eltypes(rdf) == eltypes(df)
        @test rdf == df
        rdf_ascii = sexp2julia(load(joinpath(rdata_path, "types_ascii.rda"), convert=false)["df"])
        @test eltypes(rdf_ascii) == eltypes(df)
        @test rdf_ascii == df
    end

    @testset "NAs conversion" begin
        df = DataFrame(num = Union{Float64, Missing}[1.1, 2.2],
                       int = Union{Int32, Missing}[1, 2],
                       logi = Union{Bool, Missing}[true, false],
                       chr = Union{String, Missing}["ab", "c"],
                       factor = categorical(Union{String, Missing}["ab", "c"], true),
                       cplx = Union{ComplexF64, Missing}[1.1+0.5im, 1.0im])

        df[2, :] .= Ref(missing)
        push!(df, df[2, :]) # add another row
        df[3, :num] = NaN
        df[!, :cplx] = [missing, ComplexF64(1, NaN), NaN]
        @test isequal(sexp2julia(load(joinpath(rdata_path, "NAs.rda"), convert=false)["df"]), df)
        @test isequal(sexp2julia(load(joinpath(rdata_path, "NAs_ascii.rda"), convert=false)["df"]), df)
    end

    @testset "Column names conversion" begin
        rda_names = names(sexp2julia(load(joinpath(rdata_path, "names.rda"), convert=false)["df"]))
        expected_names = [:_end, :x!, :x1, :_B_C_, :x, :x_1]
        @test rda_names == expected_names
        rda_names = names(sexp2julia(load(joinpath(rdata_path, "names_ascii.rda"), convert=false)["df"]))
        @test rda_names == [:_end, :x!, :x1, :_B_C_, :x, :x_1]
    end

    @testset "Reading RDA with complex types (environments, closures etc)" begin
        rda_envs = load(joinpath(rdata_path, "envs.rda"), convert=false)
        rda_pairlists = load(joinpath(rdata_path, "pairlists.rda"), convert=false)
        rda_closures = load(joinpath(rdata_path, "closures.rda"), convert=false)
        rda_cmpfuns = load(joinpath(rdata_path, "cmpfun.rda"), convert=false)
    end

    @testset "Proper handling of factor and ordered" begin
        f = load(joinpath(rdata_path, "ord.rda"))
        @test !isordered(f["x"])
        @test levels(f["x"]) == ["a", "b", "c"]
        @test isordered(f["y"])
        @test levels(f["y"]) == ["b", "a", "c"]
        @test f["x"] == f["y"] == ["a", "b", "c"]
    end

end # for ver in ...

@testset "Loading AltRep-containing RData files (version=3)" begin
    altrep_rda = load(joinpath("data_v3", "altrep.rda"), convert=false)
    load(joinpath("data_v3", "altrep_ascii.rda"), convert=false)
    @test length(altrep_rda) == 3
    # test AltRep objects are recognized
    @test isa(altrep_rda["longseq"], RData.RAltRep)
    # TODO test that longseq is converted into Julia UnitRange
    @test isa(altrep_rda["wrapvec"], RData.RAltRep)
    @test sexp2julia(altrep_rda["wrapvec"]) == [1.0, 2.5, 3.0]
    @test isa(altrep_rda["nonnilpairlist"], RData.RAltRep)

    # test automatic conversion
    altrep_conv_rda = load(joinpath("data_v3", "altrep.rda"), convert=true)
    @test altrep_conv_rda["wrapvec"] == [1.0, 2.5, 3.0]
    @test isa(altrep_conv_rda["nonnilpairlist"], Matrix{Int32})
    @test size(altrep_conv_rda["nonnilpairlist"]) == (0, 10)
end

end # module TestRDA

ext: jl


text:
module SimpleWeave

using Weave
using Markdown: @md_str

export @md_str # so one can run over the md blocks
export simpleweave

abstract type Block end

struct Codeblock <: Block
    block::String
    arguments::String
end

Codeblock(block) = Codeblock(block, "")

struct Markdownblock <: Block
    block::String
end

function convert_to_blocks(inputfile)

    blocks = Block[]

    io = IOBuffer()

    in_markdown = false
    arguments = ""

    function newcode()
        push!(blocks, Codeblock(strip(String(take!(io))), arguments))
        arguments = ""
        nothing
    end
    function newmd()
        push!(blocks, Markdownblock(strip(String(take!(io)))))
        arguments = ""
        nothing
    end

    for line in readlines(inputfile)
        if !in_markdown
            if startswith(line, "##")
                newcode()
                arguments = strip(line[3:end])
            elseif startswith(line, "md\"\"\"")
                newcode()
                in_markdown = true
            else
                # try single-line markdown
                m = match(r"^md\"(.*)\"\s*$", line)
                if !isnothing(m)
                    newcode()
                    println(io, m[1])
                    newmd()
                else
                    println(io, line)
                end
            end
        else
            if startswith(line, "\"\"\"")
                in_markdown = false
                newmd()
            else
                println(io, line)
            end
        end
    end

    in_markdown ? newmd() : newcode()

    blocks = filter(!isempty ∘ (x -> x.block), blocks)
end

function blocks_to_string(blocks)
    io = IOBuffer()

    for (i, block) in enumerate(blocks)
        i > 1 && println(io)
        print_block(io, block)
    end

    String(take!(io))
end

function print_block(io, c::Codeblock)
    print(io, "```julia")
    if !isempty(c.arguments)
        println(io, " ", c.arguments)
    else
        println(io)
    end
    println(io, c.block)
    println(io, "```")
end

function print_block(io, m::Markdownblock)
    println(io, m.block)
end

function simpleweave(input, outputfolder; overwrite = true, remove_pdf_aux = true, doctype = "md2html", kwargs...)

    if isfile(outputfolder)
        error("$outputfolder is a file, not a folder")
    end
    if !isdir(outputfolder)
        @label question1
        print("Output folder $outputfolder doesn't exist. Create? y/n ")
        yn = readline()
        if yn == "n"
            return
        elseif yn != "y"
            @goto question1
        end
        mkdir(outputfolder)
    end

    blocks = convert_to_blocks(input)
    weavestring = blocks_to_string(blocks)

    jmd_filename = splitext(input)[1] * ".jmd"
    if isfile(jmd_filename)
        @label question
        print("Temporary jmd file $jmd_filename already exists. Delete? y/n ")
        yn = readline()
        if yn == "n"
            return
        elseif yn != "y"
            @goto question
        end
    end

    mktempdir() do path
        temp_output_path = joinpath(path, "weave_output")

        try
            open(jmd_filename, "w") do file
                write(file, weavestring)
            end

            Weave.weave(jmd_filename;
                doctype = doctype,
                out_path = temp_output_path,
                kwargs...
            )
        finally
            rm(jmd_filename)
        end

        # delete empty folders in the output that weave leaves there
        for (root, dirs, files) in walkdir(temp_output_path, topdown = false)
            for dir in dirs
                if isempty(readdir(joinpath(root, dir)))
                    rm(joinpath(root, dir))
                end
            end
            if remove_pdf_aux && doctype == "md2pdf"
                for file in files
                    if !endswith(file, ".pdf")
                        rm(joinpath(root, file))
                    end
                end
            end
        end

        for thing in readdir(temp_output_path)
            mv(
                joinpath(temp_output_path, thing),
                joinpath(outputfolder, thing),
                force = overwrite)
        end
    end
end

end


text:
using SACBP

function pControlUpdate(model::SimulateManipulate2D,
                        b::BelMvNormal{T},
                        UArray::Vector{<:MControl2D},
                        u_param_min::Vector{<:Real},
                        u_param_max::Vector{<:Real},
                        posGain::Real,
                        rotGain::Real) where T <: Real
    nominalPolicy = ManipulatePositionControlPolicy();
    u_val = vec(control_nominal(nominalPolicy,b,u_param_min,u_param_max,posGain,rotGain))
    UArray_new = copy(UArray)
    for ii = 1:length(UArray_new)
        if UArray_new[ii].t < round(b.t + model.dto,digits=5)
            UArray_new[ii] = MControl2D(UArray_new[ii].t,u_val);
        end
    end
    return UArray_new
end;

function pControlVal(b::BelMvNormal{T},
                     u_param_min::Vector{<:Real},
                     u_param_max::Vector{<:Real},
                     posGain::Real,
                     rotGain::Real) where T <: Real
    nominalPolicy = ManipulatePositionControlPolicy();
    u_val = vec(control_nominal(nominalPolicy,b,u_param_min,u_param_max,posGain,rotGain))
    return u_val
end;


text:
#Laguerre breaks at 16
for i in 12:20
    println("n = $i ", all(L_rec(i) .< 1/eps()))
end

# Laguerre Coefficient Calculator
# Credit to [Peter Luschny Apr. 11, 2015](https://oeis.org/A021009)
function l(n,m)
    if 0 == n == m
        1
    elseif m == -1
        0
    elseif n < m
        0
    elseif n >= 1
        (n+m+1) * l(n-1, m) - l(n-1, m-1)
    end
end

LTable(n) = [l(i,j) for i in 0:n,j in 0:n]

horner_L12(x) = @evalpoly x 1.0 -12.0 33.0 -36.666666666666664 20.625 -6.6 1.2833333333333334 -0.15714285714285714 0.012276785714285714 -0.0006062610229276896 1.8187830687830687e-5 -3.0062530062530064e-7

# Hermite Coefficient Calculator
# Credit to [Paul Barry, Aug 28, 2005.](https://oeis.org/search?q=coefficient+triangle+hermite&sort=&language=&go=Search)
function h(n,k)

    if iseven(n-k) && (n-k >= 0)
        ( (-1)^((n-k)/2) * (2^k) * factorial(n)) / (factorial(k) * factorial((n-k)/2))
    else
        0
    end
end

Htable(n) = [h(i,j) for i in 0:n, j in 0:n]

horner_H12(x) = Base.Math.@evalpoly x 665280 0 -7983360 0 13305600 0 -7096320 0 1520640 0 -135168 0 4096

# Legendre Coefficient Calculator
# Credit to Ralf Stephan, Apr. 07. 2016 https://oeis.org/A008556
PTable(n) =  [binomial(2*(i-j),i-j) * binomial(i-j,j) for i in 0:n, j in 0:n]

horner_P12(x) = Base.Math.@evalpoly x

# Tshebyshev T First Kind Coefficient Calculator
# Credit to Micahel Somos, Aug. 08, 2011 https://oeis.org/A049310
function T(n,k)
    if k < 0 || k > n || (n + k) % 2 == 1
        return 0
    end
    return (-1)^((n+k)/2 +k) * binomial(Int((n+k)/2),k)
end

TTable(n) = Int[u(i,j) for i in 0:n, j in 0:n]

# Tshebyshev U Second Kind Coefficient Calculator

function u(n,m)
    if n < m || isodd(n+m)
        0
    else
        ((-1)^((n+m)/2+m))*(2^m)*binomial(Int((n+m)/2), m)
    end
end

UTable(n) = Int[u(i,j) for i in 0:n, j in 0:n]

# Consider adding the Bell and Bernoulli polynomials

# stirling numbers of the first Kind

s1(m,n) = sum([(-1)^k*binomial(n-1+k,n-m+k)*binomial(2*n-m,n-m-k)*s2(k,n-m-k) for k in 0:n-m])
row_s1(n) = [s1(i,n) for i in 1:n]

# stirling numbers of the second kind

s2(n,k) = (1/factorial(k))*sum([(-1)^i * binomial(k,i) * (k - i)^n for i in 0:k])

row_s2(n) = [s2(n,i) for i in 1:n]

s(n,m) = sum([(-1)^(n-j) * binomial(k,j) * factorial(j*y - 1 + n) / factorial(j*y-1) * y ^ (-k) / factorial(k) for j in 0:k])

function stirlings1(n::Int, k::Int, signed::Bool=false)
    if signed == true
        return (-1)^(n - k) * stirlings1(n, k)
    end

    if n < 0
        throw(DomainError(n, "n must be nonnegative"))
    elseif n == k == 0
        return 1
    elseif n == 0 || k == 0
        return 0
    elseif n == k
        return 1
    elseif k == 1
        return factorial(n-1)
    elseif k == n - 1
        return binomial(n, 2)
    elseif k == n - 2
        return div((3 * n - 1) * binomial(n, 3), 4)
    elseif k == n - 3
        return binomial(n, 2) * binomial(n, 4)
    end

    return (n - 1) * stirlings1(n - 1, k) + stirlings1(n - 1, k - 1)
end

function stirlings2(n::Int, k::Int)
    if n < 0
        throw(DomainError(n, "n must be nonnegative"))
    elseif n == k == 0
        return 1
    elseif n == 0 || k == 0
        return 0
    elseif k == n - 1
        return binomial(n, 2)
    elseif k == 2
        return 2^(n-1) - 1
    end

    return k * stirlings2(n - 1, k) + stirlings2(n - 1, k - 1)
end

# Straight from Julia codebase
macro evalpoly(z, p...)
    a = :($(esc(p[end])))
    b = :($(esc(p[end-1])))
    as = []
    for i = length(p)-2:-1:1
        ai = Symbol("a", i)
        push!(as, :($ai = $a))
        a = :(muladd(r, $ai, $b))
        b = :($(esc(p[i])) - s * $ai) # see issue #15985 on fused mul-subtract
    end
    ai = :a0
    push!(as, :($ai = $a))
    C = Expr(:block,
             :(x = real(tt)),
             :(y = imag(tt)),
             :(r = x + x),
             :(s = muladd(x, x, y*y)),
             as...,
             :(muladd($ai, tt, $b)))
    R = Expr(:macrocall, Symbol("@horner"), (), :tt, map(esc, p)...)
    :(let tt = $(esc(z))
          isa(tt, Complex) ? $C : $R
      end)
end

#P
d(p, α,n) = binomial(n + α,n)
b(p, α, β, n, m) = (n-m+1)*(α+β+n+m)
c(p, α, m) = 2m*(α + m)
f(p, x) = 1 - x

#Ceven
d(p, α, n) = (-1)^n * (poch(α,n)/factorial(n))
b(p, α, n, m) = 2*(n-m+1)*(α+n+m-1)
c(p, m) = m*(2m-1)
f(p, x) = x^2

#Codd
d(p, α, n) = (-1)^n * (poch(α,n+1)/factorial(n))*2x
b(p, α, β, n, m) = 2*(n-m+1)*(α+β+n+m)
c(p, m) = m*(2m+1)
f(p, x) = x^2

#Teven
d(p, n) = (-1)^n
b(p, n, m) = 2(n-m+1)*(n+m-1)
c(p, m) = m*(2m-1)
f(p, x) = x^2

#Todd
d(p, n, x) = (-1^n)*(2n+1)*x
b(p, n, m) = 2*(n-m+1)*(n+m)
c(p, m) = m*(2m+1)
f(p, x) = x^2

#Ueven
d(p, n) = (-1)^n
b(p, n, m) = 2*(n-m+1)*(n+m)
c(p, m) = m*(2m+1)
f(p, x) = x^2

#Uodd
d(p, n, x) = (-1)^n * 2(n+1)*x
b(p, n, m) = 2*(n-m+1)*(n+m)
c(p, m) = m*(2m + 1)
f(p, x) = x^2

#Peven
d(p, n) = (-1/4)^n * binomial(2n,n)
b(p, n, m) = (n - m + 1) * (2n + 2m -1)
c(p, m) = m*(2m-1)
f(x) = x^2

#Podd
d(p, n, x) = (-1/4)^n * binomial(2n+1,n) * (n+1) * x
b(p, n, m) = (n - m + 1) * (2n + 2m + 1)
c(p, m) = m * (2m + 1)
f(p, x) = x^2

#L
d(p, α, n) = binomial(n + α, n)
b(p, n, m) = (n - m + 1)
c(p, m, α) = m * (α + m)
f(p, x) = x

#Heven
d(p, n) = (-1) ^ n * (factorial(2n)/factorial(n))
b(p, n, m) = 2*(n - m + 1)
c(p, m) = m*(2m - 1)
f(p, x) = x^2

#Hodd
d(p, n, x) = (-1) ^ n * (factorial(2n +1)/factorial(n))*2x
b(p, n, m) = 2*(n - m + 1)
c(p, m) = m * (2m + 1)
f(p, x) = x^2

#interpolate an expression as an expression into another expression.
ex_new = Meta.quot(ex)
quote
    still_expression = $(esc(ex_new))
end


text:
function cudnnCreate()
    handle_ref = Ref{cudnnHandle_t}()
    @check unsafe_cudnnCreate(handle_ref) CUDNN_STATUS_NOT_INITIALIZED CUDNN_STATUS_INTERNAL_ERROR
    return handle_ref[]
end

function cudnnGetProperty(property::CUDA.libraryPropertyType)
  value_ref = Ref{Cint}()
  cudnnGetProperty(property, value_ref)
  value_ref[]
end

function version()
  ver = cudnnGetVersion()
  major, ver = divrem(ver, 1000)
  minor, patch = divrem(ver, 10)

  VersionNumber(major, minor, patch)
end

function cuda_version()
  ver = cudnnGetCudartVersion()
  major, ver = divrem(ver, 1000)
  minor, patch = divrem(ver, 10)

  VersionNumber(major, minor, patch)
end


text:
import Pkg
Pkg.activate("..")

using SnoopCompile

inf_timing = @snoopi tmin=0.0 begin
    using QuantumAlgebra
    include("./snoop_statements.jl")
    include("../test/runtests.jl")
end

pc = SnoopCompile.parcel(inf_timing)

SnoopCompile.write("/tmp/precompile", pc)


text:
"""
    kernel_1d(window)
Create a kernel to perform a moving average filtering of a 1-d array. It will use window as windows size unless the given vaule is even, in which case it will add one to the given value before creating the kernel. This kernel has to be passed into imfilter or imfilter!
"""
function kernel_1d(window)
    if !isodd(window)
         window += 1
     end
     kernel_t = centered(ones(window)./(window))
 end

 """
    kernel_4d(window_h,window_t,T :: Type :: Float64)
Create a kernel to perform a moving average filtering of a 4-d array along the dimensions 1,2 and 4. It will use window_h (dimensions 1 and 2) and window_t (dimension 4)  as windows sizes unless the given vaules are even, in which case it will add one to the given value before creating the kernel. This kernel has to be passed into imfilter or imfilter!
"""
 function kernel_4d(window_h,window_t,T::Type=Float64)
    kernel_h = if (window_h == 1) | (window_h == false)
        centered([T(1.0)])
    else
        !isodd(window_h) && (window_h += 1 )
        ones(T,window_h)./(window_h)
    end
     if !isodd(window_t)
         window_t += 1
     end

    kernel_t = if (window_t == 1) | (window_t == false)
        centered([T(1.0)])
    else
        !isodd(window_t) && (window_t += 1 )
        ones(T,window_t)./(window_t)
    end

     return kernelfactors(( kernel_h, kernel_h, centered([T(1.0)]) , kernel_t ))
 end

 """
kernel_3d(window_h,window_t,T :: Type :: Float64)
Create a kernel to perform a moving average filtering of a 3-d array along the dimensions 1,2 and 3. It will use window_h (dimensions 1 and 2) and window_t (dimension 3) as windows sizes unless the given vaules are even, in which case it will add one to the given value before creating the kernel. This kernel has to be passed into imfilter or imfilter!
"""
 function kernel_3d(window_h,window_t,T::Type=Float64)
    if !isodd(window_h)
         window_h += 1
     end
     if !isodd(window_t)
         window_t += 1
     end
     kernel_h = ones(T,window_h)./(window_h)
     return  kernelfactors(( centered(kernel_h), centered(kernel_h), centered([T(1.0)]) ))
 end

 """
 kernel_2d(window_h,T :: Type :: Float64)
 Create a kernel to perform a moving average filtering of a 2-d array along the dimensions 1,2. It will use window_h as windows size unless the given vaule is even, in which case it will add one to the given value before creating the kernel. This kernel has to be passed into imfilter or imfilter!
 """
 function kernel_2d(window_h,T::Type=Float64)
    if !isodd(window_h)
         window_h += 1
     end
     kernel_h = ones(T,window_h)./(window_h)
     return  kernelfactors(( centered(kernel_h), centered(kernel_h) ))
 end

 """
 kernel_4d_t(window_t,T :: Type :: Float64)
 Create a kernel to perform a moving average filtering of a 4-d array along the dimension 4. It will use window_t as windows size unless the given vaule is even, in which case it will add one to the given value before creating the kernel. This kernel has to be passed into imfilter or imfilter!
 """
 function kernel_4d_t(window_t,T::Type=Float64)
     if !isodd(window_t)
         window_t += 1
     end
     kernel_t = ones(T,window_t)./(window_t)

     return kernelfactors(( centered([T(1.0)]), centered([T(1.0)]), centered([T(1.0)]), centered(kernel_t) ))
 end

 """
 kernel_3d_t(window_t,T :: Type :: Float64)
 Create a kernel to perform a moving average filtering of a 3-d array along the dimension 3. It will use window_t as windows size unless the given vaule is even, in which case it will add one to the given value before creating the kernel. This kernel has to be passed into imfilter or imfilter!
 """
 function kernel_3d_t(window_t,T::Type=Float64)
     if !isodd(window_t)
         window_t += 1
     end
     kernel_t = ones(T,window_t)./(window_t)
     return  kernelfactors(( centered([T(1.0)]), centered([T(1.0)]), centered(kernel_t) ))
 end

 """
     filter_array_2!(array,smooth_x,smooth_t,position)

 Filter the input array _in-place_ using a moving mean.

 In the first two dimensions the window width is smooth_x and the border is treated as circular (for doubly periodic domains). In the space dimension the width is smooth_t and the border value is replicated.

 This function calls under the hood the imfilter function of Images.jl

 """
function filter_array_2!(buf :: Array{T,4},array::Array{T,4},smooth_x,smooth_time,position) where T <: Real
    sx,sy,sz,st = size(array)
#    filtered = similar(array)
     if position == 2
         filtered = imfilter(imfilter(array, kernel4d(smooth_x,smooth_time),"circular"),kernel_4d_t(smooth_time),"inner")
     else
         ###Filtering in space

         @info typeof(array) typeof(buf)
         @inbounds for t in 1:st
             @Threads.threads for z in 1:sz
                 imfilter!(view(buf,:,:,z,t),array[:,:,z,t], kernel_2d(smooth_x),"circular")
             end
         end
         @info typeof(array) typeof(buf)
     ### filtering in time
         @inbounds for z in 1:sz, y in 1:sy
             @Threads.threads for x in 1:sx
                 array[x,y,z,:] .= imfilter(buf[x,y,z,:], kernel_1d(smooth_time),"symmetric")
             end
         end
     end
    return array
end

 """
 filter_array_2!(array,smooth_x,smooth_t,position)

Filters the input array _in-place_ using a moving mean.

In the first two dimensions the window width is smooth_x and the border is treated as circular (for doubly periodic domains). In the space dimension the width is smooth_t and the border value is replicated.

This function calls under the hood the imfilter function of Images.jl

"""
 function filter_array_2!(array::Array{T,3},smooth_x,smooth_time,position) where T <: Real
     filtered = similar(array)
     sx,sy,st = size(array)
     if position == 2
         filtered = imfilter(imfilter(array, kernel4d(smooth_x,smooth_time),"circular"),kernel4d_t(smooth_time),"inner")
     else
         ###Filtering in space
         for t in 1:st
             filtered[:,:,t] = imfilter(filtered,array[:,:,t], kernel_2d(smooth_x),"circular")
         end
         ### filtering in time
         for y in 1:sy, x in 1:sx
             array[x,y,:] = imfilter(filtered[x,y,:], kernel_1d(smooth_time),"symmetric")
         end
     end
     #return array
 end
 """
     filter_array!(buffer,array,smooth_x,smooth_t,position)

 Filters the input array _in-place_ using a moving mean. In the first two dimensions the window width is smooth_x and the border is treated as circular (for doubly periodic domains). In the space dimension the width is smooth_t and the border value is replicated.

 This function calls under the hood the imfilter function of Images.jl

 The first argument must be a buffer of the same size of array.

 """
 function filter_array!(buf::Array{T,4},array::Array{T,4},smooth_x,smooth_time,position) where T <: Real
     if position == 2
        error("Filter_array: Inner array is not implemented yet")
     else
         imfilter!(buf,array,kernel_4d(smooth_x,smooth_time,T)[1:2],"circular")
         imfilter!(array,buf,(kernel_4d_t(smooth_time,T)[4],),"symmetric")
     end
 #        return filtered
 end

 function filter_array!(buf::Array{T,3},array::Array{T,3},smooth_x,smooth_time,position) where T <: Real
     if position==2
        error("Filter_array: Inner array is not implemented yet")
     else
         imfilter!(buf,array, kernel_3d(smooth_x,smooth_time,T)[1:2],"circular")
         imfilter!(array,buf,(kernel_3d_t(smooth_time,T)[3],),"symmetric")
     end
 end

 """
 filter_array(array,smooth_x,smooth_t,position)

 Filters the input arrayusing a moving mean. In the first two dimensions the window width is smooth_x and the border is treated as circular (for doubly periodic domains). In the space dimension the width is smooth_t and the border value is replicated.

 This function calls under the hood the imfilter function of Images.jl

 The first argument must be a buffer of the same size of array.

 """
 function filter_array(array::Array{T,4},smooth_x,smooth_time,position = 1) where T <: Real
     if position == 2
        asize = size(array)
        tail_size = (asize[4]-1)÷2
        buf1 = Array{T,4}(undef,asize[1],asize[2],asize[3],asize[4] - 2*tail_size)
        buf2 = Array{T,4}(undef,asize[1],asize[2],asize[3],asize[4] - 2*tail_size)
        buf_axes = axes(buf1)
        imfilter!(OffsetArray(buf1,buf_axes[1],buf_axes[2],buf_axes[3],buf_axes[4].+tail_size),array,(kernel_4d_t(smooth_time,T)[3],),Inner())
        imfilter!(buf2,buf1, kernel_4d(smooth_x,smooth_time,T)[1:2],"circular")
     else
         buf1 = similar(array)
         buf2 = similar(array)
         imfilter!(buf1,array,kernel_4d(smooth_x,smooth_time,T)[1:2],"circular")
         imfilter!(buf2,buf1,(kernel_4d_t(smooth_time,T)[4],),"symmetric")
     end
         return buf2
 end

 """
 filter_array(array,smooth_x,smooth_t,position)

 Filters the input arrayusing a moving mean. In the first two dimensions the window width is smooth_x and the border is treated as circular (for doubly periodic domains). In the space dimension the width is smooth_t and the border value is replicated.

 This function calls under the hood the imfilter function of Images.jl

 The first argument must be a buffer of the same size of array.

 """
function filter_array(array::Array{T,3},smooth_x,smooth_time,position = 1) where T <: Real
    if !isodd(smooth_time)
         smooth_time += 1
    end
    if position==2
        asize = size(array)
        tail_size = (asize[3]-1)÷2
        buf1 = Array{T,3}(undef,asize[1],asize[2],asize[3] - 2*tail_size)
        buf2 = Array{T,3}(undef,asize[1],asize[2],asize[3] - 2*tail_size)
        buf_axes = axes(buf1)
        imfilter!(OffsetArray(buf1,buf_axes[1],buf_axes[2],buf_axes[3].+tail_size),array,(kernel_3d_t(smooth_time,T)[3],),Inner())
        imfilter!(buf2,buf1, kernel_3d(smooth_x,smooth_time,T)[1:2],"circular")

     else
        buf1 = similar(array)
        buf2 = similar(array)
        imfilter!(buf1,array, kernel_3d(smooth_x,smooth_time,T)[1:2],"circular")
        imfilter!(buf2,buf1,(kernel_3d_t(smooth_time,T)[3],),"symmetric")
     end
     return buf2
 end

 """
 filter_array_nospace(array,smooth_t,position)

 Filters the input array using a moving mean along the fourth dimension. In the 4th dimension the width is smooth_t and the border value is replicated except if position = 2, in which case it only takes the inner part of the smoothed array.

 This function calls under the hood the imfilter function of Images.jl
 """
 function filter_array_nospace(array::Array{T,4},smooth_time,position = 1) where T <: Real
     if position==2
         filtered = imfilter(array,kernel_4d_t(smooth_time),"inner")
     else
         filtered = imfilter(array,kernel_4d_t(smooth_time),"replicate")
     end
     return filtered
 end
 """
 filter_array_nospace(array,smooth_t,position)

 Filters the input array using a moving mean along the third dimension. In the 3th dimension the width is smooth_t and the border value is replicated except if position = 2, in which case it only takes the inner part of the smoothed array.

 This function calls under the hood the imfilter function of Images.jl
 """
 function filter_array_nospace(array::Array{T,3},smooth_time,position = 1) where T <: Real
     if position==2
         filtered = imfilter(array,kernel_3d_t(smooth_time),"inner")
     else

         filtered = imfilter(array,kernel_3d_t(smooth_time),"replicate")
     end
     return filtered
 end
 """
 filter_array_nospace(array,smooth_t,position=1)

 Filters the input, 1-d array. The border value is replicated except if position = 2, in which case it only takes the inner part of the smoothed array.

 This function calls under the hood the imfilter function of Images.jl
 """
 function filter_array_time(array::Array{T,1},window,position = 1) where T <: Real
     kern = centered(ones(window)./window)
     if position==2
         smooth_t = imfilter(array,kern,"inner")
     else
         smooth_t = imfilter(array,kern,"symmetric")
     end
     return smooth_t
 end


text:
convertDateTime(dt::String) = DateTime(dt[1:end-1],"y-m-dTH:M:S")

module FieldTags

using Compat

export @tag, queryTag

typealias TagDict Dict{Symbol,Dict{Symbol,Any}}

const tagRegistry = Dict{Symbol,TagDict}()

macro tag(t)
	# dump(t)
	# println()

	isa(t,Expr) || error("Tagging can only happen on a type")

	isMut = t.args[1]
	tName = t.args[2]
	# if isa(tName,Expr)
	# 	tName.head != :curly && error("Type Name is not a curly type")
	# 	tName = tName.args[1]
	# end
	fields = Any[]
	tTags = Dict{Symbol,Expr}()
	x = t.args[end].args
	for y in x
		if isa(y, LineNumberNode) || (isa(y,Expr) && y.head == :line)
			# line number node
			continue
		else isa(y,Expr)
			# println(y)
			# dump(y)
			# println()
			if y.head == :(=>)
				push!(fields,y.args[1])
				tTags[y.args[1].args[1]] = y.args[2]
			else
				push!(fields,y)
			end
		end
	end
	t.args[end].args = fields

	tagDict = Dict{Symbol,Dict{Symbol,Any}}()

	for (value,tags) in tTags
		if tags.head == :(:)
			tags = (tags,)
		end

		for tag in tags
			tagData = get(tagDict,tag.args[1],Dict{Symbol,Any}())
			tagData[value] = tag.args[2]
			tagDict[tag.args[1]] = tagData
		end
	end
	currentMod = current_module()
	if currentMod != Main
		tagRegistry[Symbol(module_name(currentMod),".",tName)] = tagDict
	else
		tagRegistry[Symbol(tName)] = tagDict
	end

	return esc(t)
end

function queryTag{T}(::Type{T}, tagName::Symbol)
	typeTags = tagRegistry[Symbol(T)]
	return get(typeTags,tagName,Dict{Symbol,Any}())
end

function queryTag(sym::Symbol, tagName::Symbol)
	typeTags = tagRegistry[sym]
	return get(typeTags,tagName,Dict{Symbol,Any}())
end

end

@generated marshalJSON{T}(data::T) = marshalJSON_impl(data)

function marshalJSON_impl(T)
	if T <: Union{Number,AbstractString,}
		return :(data)
	elseif T <: Vector
		return :(map(marshalJSON,data))
	elseif T <: Dict{String}
		return :(Dict( (k=>marshalJSON(v) for (k,v) in data)))
	else
		tagData = try
				FieldTags.queryTag(getfield(T.name.module,T.name.name),:json)
		catch
			FieldTags.TagDict()
		end
		expr = quote
		jsonData = Dict{String,Any}()
		end
		fields = fieldnames(T)

		if isempty(fields)
			error("Invalid Type")
		end

		for field in fields
			fType = fieldtype(T,field)
			fArgString = split(get(tagData,field,""),',')
			f = string(field)

			fName = isspace(fArgString[1]) ? f : fArgString[1]
			shift!(fArgString)

			fArgs = Dict{String,Nullable{String}}()

			for fArg in fArgString
				fA = split(fArg,':'; limit = 2)
				if length(fA) == 1
					fArgs[fA[1]] = Nullable{String}()
				elseif length(fA) == 2
					fArgs[fA[1]] = fA[2]
				end
			end

			if fType <: Nullable

				if fType.parameters[1] <: Dates.TimeType
					dFormat = get(get(fArgs,"format",Nullable{String}("y-m-dTH:M:S")))
					setexpr = :(Dates.format(get(getfield(data,Symbol($f))),$dFormat))
				else
					setexpr = :(marshalJSON(get(getfield(data,Symbol($f)))))
				end
				expr = quote
					$expr
					!isnull(getfield(data,Symbol($f))) && (jsonData[$(fName)] = $(setexpr))
				end
			else

				if fType <: Dates.TimeType
					dFormat = get(get(fArgs,"format",Nullable{String}("y-m-dTH:M:S")))
					setexpr = :(Dates.format(getfield(data,Symbol($f)),$dFormat))
				else
					setexpr = :(marshalJSON(getfield(data,Symbol($f))))
				end

				expr = quote
					$expr
					jsonData[$(fName)] = $setexpr
				end
			end
		end
	end
	return quote
		$expr
		return jsonData
	end
end

@generated unmarshalJSON{T}(::Type{T},data) = unmarshalJSON_impl(T,data)

function unmarshalJSON_impl{T}(::Type{T},data)
	if T <: Union{Number,AbstractString}
		return :(data)
	elseif T <: Vector
		return :(map( x -> unmarshalJSON(eltype(T),x),data))
	else
		tagData = try
			FieldTags.queryTag(getfield(T.name.module,T.name.name),:json)
		catch
			FieldTags.TagDict()
		end

		exprs = Expr[]
		fields = fieldnames(T)

		if isempty(fields)
			error("Invalid Type")
		end

		for field in fields
			fType = fieldtype(T,field)
			fArgString = split(get(tagData,field,""),',')
			f = string(field)

			fName = isspace(fArgString[1]) ? f : fArgString[1]
			shift!(fArgString)
			fArgs = Dict{String,Nullable{String}}()

			for fArg in fArgString
				fA = split(fArg,':'; limit = 2)
				if length(fA) == 1
					fArgs[fA[1]] = Nullable{String}()
				elseif length(fA) == 2
					fArgs[fA[1]] = fA[2]
				end
			end

			if fType <: Nullable

				if fType.parameters[1] <: Dates.TimeType
					dFormat = get(get(fArgs,"format",Nullable{String}("y-m-dTH:M:S")))
					getexpr = :($(fType.parameters[1])(data[$fName],$(dFormat)))
				else
					getexpr = :(unmarshalJSON($(fType.parameters[1]),data[$fName]))
				end

				expr = quote
					haskey(data,$fName) && (data[$fName] != nothing) ? $(getexpr) : Nullable{$(fType.parameters[1])}()
				end
				push!(exprs,expr);
			else

				if fType <: Dates.TimeType
					dFormat = get(get(fArgs,"format",Nullable{String}("y-m-dTH:M:S")))
					getexpr = :($(fType)(data[$fName],$(dFormat)))
				else
					getexpr = :(unmarshalJSON($(fType),data[$fName]))
				end

				expr = getexpr
				push!( exprs, expr)
			end
		end
	end
	return quote
		return $(T)($(exprs...))
	end
end


text:
using Test

@testset "arithmetic" begin
    @test 1 == 1
    @test 2 + 1 == 3
end

@test 1/1 == 1


text:
# IO routines for the Abaqus .inp file format
const abaqus_to_vtk_type = Dict{String, UInt64}(
    # 2D
    # triangle
    "CPS3"  => 5,
    "SFRU3" => 5,
    # triangle6
    "CPS6" => 22,
    # quadrilateral
    "CPS4" => 9,
    # quad8
    "CPS8" => 23,
    "M3D9" => 23
   )

function read_abaqus_2d(filepath::String)
    @info "Reading $filepath"
    # NOFE: Fhere is a crucial assumption here that elements and nodes are listed 1 to N,
    # not 8, 10, 9 or anything funky/out of order.
    name = "DefaultMeshName"
    file = open(filepath, "r")
    faces = Vector{UInt32}[]
    face_sets = Dict{String, Set{UInt32}}()
    points = Point_2D[]
    while !eof(file)
        line_split = split(readline(file))
        if length(line_split) > 0
            if occursin("**", line_split[1]) # Comment
                continue
            elseif occursin("*Heading", line_split[1])
                name = String(strip(readline(file)))
                if occursin(".inp", name)
                    name = name[1:length(name)-4]
                end
            elseif occursin("*NODE", line_split[1])
                points = read_abaqus_nodes_2d(file)
            elseif occursin("*ELEMENT", line_split[1])
                element_type = String(strip(replace(line_split[2], ("type=" => "")), ','))
                append!(faces, read_abaqus_elements(file, element_type))
            elseif occursin("*ELSET", line_split[1])
                set_name = String(replace(line_split[1], ("*ELSET,ELSET=" => "")))
                face_sets[set_name] = read_abaqus_elset(file)
            end
        end
    end
    close(file)
    # The relationship between faces, edges, and vertices is
    # Triangle
    #   F ≈ 2V
    #   E ≈ 3F/2
    # Quadrilateral
    #   F ≈ V
    #   E ≈ 2F
    # Triangle6
    #   2F ≈ V
    #   E ≈ 3F/2
    # Quadrilateral8
    #   4F ≈ V
    #   E ≈ 2F
    # We see that V ≤ 2F ≤ E, so we use 2.2F as the max number of faces, edges, or vertices
    # plus a fudge factor for small meshes, where the relationships become less accurate.
    UInt32_max = 4294967295
    nfaces = length(faces)
    if UInt32_max < 2.2*nfaces
        @error "Mesh may cause UInt32 overflow. Some poor dev now has to add UInt64 support"
    end
    face_lengths = Int64[]
    for face in faces
        l = length(face)
        if l ∉ face_lengths
            push!(face_lengths, l)
        end
    end
    if face_lengths == [3]
        return TriangleMesh_2D(name = name,
                               points = points,
                               faces = [ SVector{3, UInt32}(f) for f in faces],
                               face_sets = face_sets
                              )
    elseif face_lengths == [4]
        return QuadrilateralMesh_2D(name = name,
                                    points = points,
                                    faces = [ SVector{4, UInt32}(f) for f in faces],
                                    face_sets = face_sets
                                   )
    elseif face_lengths == [6]
        return Triangle6Mesh_2D(name = name,
                                points = points,
                                faces = [ SVector{6, UInt32}(f) for f in faces],
                                face_sets = face_sets
                               )
    else
        @error "Could not identify mesh type. Faces of lengths: $face_lengths"
        return nothing
    end
end

function read_abaqus_nodes_2d(file::IOStream)
    points = Point_2D[]
    line_split = strip.(split(readline(file)), [','])
    line_position = position(file)
    while !occursin("*", line_split[1])
        xyz = parse.(Float64, line_split[2:4])
        push!(points, Point_2D(xyz[1], xyz[2]))
        line_position = position(file)
        line_split = strip.(split(readline(file)), [','])
    end
    seek(file, line_position)
    return points
end

function read_abaqus_elements(file::IOStream, element_type::String)
    if !(element_type ∈  keys(abaqus_to_vtk_type))
        @error "$element_type is not in the abaqus to vtk type conversion dictionary"
    end
    faces = Vector{UInt32}[]
    line_split = strip.(split(readline(file)), [','])
    line_position = position(file)
    while !occursin("*", line_split[1])
        vertex_IDs = parse.(UInt32, line_split[2:length(line_split)])
        if length(vertex_IDs) == 9
            push!(faces, vertex_IDs[1:8])
        else
            push!(faces, vertex_IDs)
        end
        line_position = position(file)
        line_split = strip.(split(readline(file)), [','])
    end
    seek(file, line_position)
    return faces
end

function read_abaqus_elset(file::IOStream)
    line_split = strip.(split(readuntil(file, "*")), [','])
    seek(file, position(file)-1)
    return Set(parse.(UInt32, line_split))
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport


text:
function TimeSpanDragForce(arg0::Atmosphere, arg1::DragSensitive)
    return TimeSpanDragForce((Atmosphere, DragSensitive), arg0, arg1)
end

function TimeSpanDragForce(arg0::Atmosphere, arg1::DragSensitive, arg2::TimeScale)
    return TimeSpanDragForce((Atmosphere, DragSensitive, TimeScale), arg0, arg1, arg2)
end

function acceleration(obj::TimeSpanDragForce, arg0::FieldSpacecraftState, arg1::Vector{RealFieldElement})
    return jcall(obj, "acceleration", FieldVector3D, (FieldSpacecraftState, Vector{RealFieldElement}), arg0, arg1)
end

function acceleration(obj::TimeSpanDragForce, arg0::SpacecraftState, arg1::Vector{jdouble})
    return jcall(obj, "acceleration", Vector3D, (SpacecraftState, Vector{jdouble}), arg0, arg1)
end

function add_drag_sensitive_valid_after(obj::TimeSpanDragForce, arg0::DragSensitive, arg1::AbsoluteDate)
    return jcall(obj, "addDragSensitiveValidAfter", void, (DragSensitive, AbsoluteDate), arg0, arg1)
end

function add_drag_sensitive_valid_before(obj::TimeSpanDragForce, arg0::DragSensitive, arg1::AbsoluteDate)
    return jcall(obj, "addDragSensitiveValidBefore", void, (DragSensitive, AbsoluteDate), arg0, arg1)
end

function extract_drag_sensitive_range(obj::TimeSpanDragForce, arg0::AbsoluteDate, arg1::AbsoluteDate)
    return jcall(obj, "extractDragSensitiveRange", TimeSpanMap, (AbsoluteDate, AbsoluteDate), arg0, arg1)
end

function extract_parameters(obj::TimeSpanDragForce, arg0::Vector{RealFieldElement}, arg1::FieldAbsoluteDate)
    return jcall(obj, "extractParameters", Vector{RealFieldElement}, (Vector{RealFieldElement}, FieldAbsoluteDate), arg0, arg1)
end

function extract_parameters(obj::TimeSpanDragForce, arg0::Vector{jdouble}, arg1::AbsoluteDate)
    return jcall(obj, "extractParameters", Vector{jdouble}, (Vector{jdouble}, AbsoluteDate), arg0, arg1)
end

function get_drag_sensitive(obj::TimeSpanDragForce, arg0::AbsoluteDate)
    return jcall(obj, "getDragSensitive", DragSensitive, (AbsoluteDate,), arg0)
end

function get_drag_sensitive_span(obj::TimeSpanDragForce, arg0::AbsoluteDate)
    return jcall(obj, "getDragSensitiveSpan", TimeSpanMap_Span, (AbsoluteDate,), arg0)
end

function get_events_detectors(obj::TimeSpanDragForce)
    return jcall(obj, "getEventsDetectors", Stream, ())
end

function get_field_events_detectors(obj::TimeSpanDragForce, arg0::Field)
    return jcall(obj, "getFieldEventsDetectors", Stream, (Field,), arg0)
end

function get_parameters_drivers(obj::TimeSpanDragForce)
    return jcall(obj, "getParametersDrivers", Vector{ParameterDriver}, ())
end

function get_transitions(obj::TimeSpanDragForce)
    return jcall(obj, "getTransitions", NavigableSet, ())
end


text:
using Test

include("../src/main_model.jl")
let m = getpage("NDCs", false, false)
    run(m)

    while m[:Discontinuity,:occurdis_occurrencedummy] != [0.,0.,0.,0.,0.,0.,0.,1.,1.,1.]
        run(m)
    end

    # Climate component
    temp = m[:ClimateTemperature, :rt_g_globaltemperature]
    temp_compare = readpagedata(m, "test/validationdata/noperm/rt_g_globaltemperature.csv")
    @test temp ≈ temp_compare rtol = 1e-4

    # Abatement Costs
    abatement = m[:TotalAbatementCosts, :tct_totalcosts]
    abatement_compare = readpagedata(m, "test/validationdata/noperm/tct_totalcosts.csv")
    @test abatement ≈ abatement_compare rtol = 1e-2

    te = m[:EquityWeighting, :te_totaleffect]
    te_compare = 9.462179758012009e8
    @test te ≈ te_compare rtol = 1e4
end


text:
"""
    angleofrefraction(n1, n2, β)

Compute the angle of refraction.

# Arguments
* `n1`: refractive index of the medium in which the wave is propagating towards the interface
* `n2`: refractive index of the other medium
* `β`:  incident or reflected angle
"""
function angleofrefraction(n1, n2, β)
    asin(n1 / n2 * sin(β))
end

"""
    fresnel_x(n1, n2, β)

Compute the diagonal elements \$L_{xx}\$ of the fresnel factor.

# Arguments
* `n1`: refractive index of the medium in which the wave is propagating towards the interface
* `n2`: refractive index of the other medium
* `β`:  incident or reflected angle

https://doi.org/10.1103/PhysRevB.59.1263
"""
function fresnel_x(n1, n2, β)
    γ = angleofrefraction(n1, n2, β)
    2n1 * cos(γ) / (n1 * cos(γ) + n2 * cos(β))
end

"""
    fresnel_y(n1, n2, β)

Compute the diagonal elements \$L_{yy}\$ of the fresnel factor.

# Arguments
* `n1`: refractive index of the medium in which the wave is propagating towards the interface
* `n2`: refractive index of the other medium
* `β`:  incident or reflected angle

https://doi.org/10.1103/PhysRevB.59.1263
"""
function fresnel_y(n1, n2, β)
    γ = angleofrefraction(n1, n2, β)
    2n1 * cos(β) / (n1 * cos(β) + n2 * cos(γ))
end

"""
    fresnel_z(n1, n2, n′, β)

Compute the diagonal elements \$L_{yy}\$ of the fresnel factor.

# Arguments
* `n1`: refractive index of the medium in which the wave is propagating towards the interface
* `n2`: refractive index of the other medium
* `n′`: refractive index of the interface
* `β`:  incident or reflected angle

https://doi.org/10.1103/PhysRevB.59.1263
"""
function fresnel_z(n1, n2, n′, β)
    γ = angleofrefraction(n1, n2, β)
    2n2 * cos(β) / (n1 * cos(γ) + n2 * cos(β)) * (n1 / n′)^2
end

"""
    ninterface(n2)

Calculate the index of refraction of the interfacial layer according to

``\\left(\\frac{1}{n^\\prime}\\right)^2 = \\frac{4n_2^2 + 2}{n_2^2(n_2 ^2+5)}``.

# Arguments
* `n2`: refractive index of the medium in which the wave is **not** propagating towards the interface

https://doi.org/10.1103/PhysRevB.59.1263
"""
function ninterface(n2)
    1 / sqrt((4 * n2^2 + 2) / (n2^2 * (n2^2 + 5)))
end


text:
using Dates

function discordiandate(year::Integer, month::Integer, day::Integer)
    discordianseasons = ["Chaos", "Discord", "Confusion", "Bureaucracy", "The Aftermath"]
    holidays = Dict(
        "Chaos 5" => "Mungday",
        "Chaos 50" => "Chaoflux",
        "Discord 5" => "Mojoday",
        "Discord 50" => "Discoflux",
        "Confusion 5" => "Syaday",
        "Confusion 50" => "Confuflux",
        "Bureaucracy 5" => "Zaraday",
        "Bureaucracy 50" => "Bureflux",
        "The Aftermath 5" => "Maladay",
        "The Aftermath 50" => "Afflux")
    today = Date(year, month, day)
    isleap = isleapyear(year)
    if isleap && month == 2 && day == 29
        rst = "St. Tib's Day, YOLD " * string(year + 1166)
    else
        dy = dayofyear(today)
        if isleap && dy >= 60
            dy -= 1
        end
        rst = string(discordianseasons[div(dy, 73) + 1], " ", rem(dy, 73)) # day
        if haskey(holidays, rst)
            rst *= " ($(holidays[rst]))" # if holiday
        end
        rst *= ", YOLD $(year + 1166)" # year
    end
    return rst
end

@show discordiandate(2017, 08, 15)
@show discordiandate(1996, 02, 29)
@show discordiandate(1996, 02, 19)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#=
Type definition and methods containing the overarching backtesting object fueling the engine
=#

import Base: show

const TABWIDTH = 4
const TAB = ' ' ^ TABWIDTH

mutable struct Strategy
    universe::Universe
    indicator::Indicator
    rules::Tuple{Vararg{Rule}}
    portfolio::Portfolio
    backtest::Backtest
    function Strategy(universe::Universe,
                      indicator::Indicator,
                      rules::Tuple{Vararg{Rule}},
                      portfolio::Portfolio=Portfolio(universe))
        return new(universe, indicator, rules, portfolio, Backtest())
    end
end

function show(io::IO, strat::Strategy)
    println(io, strat.indicator.paramset)
    println(io, "# Rules:")
    for rule in strat.rules
        println(io, TAB, rule)
    end
    println()
    show(io, strat.universe)
end

function summarize_results(strat::Strategy)
    holdings = TS()
    values = TS()
    profits = TS()
    for asset in strat.universe.assets
        asset_result = strat.backtest.backtest[asset]
        holding = asset_result[:Pos]
        holdings = [holdings holding]
        values = [values cl(asset_result) * holding]
        profits = [profits asset_result[:PNL]]
    end
    # data cleaning - field assignment and missing value replacement
    holdings.fields = Symbol.(strat.universe.assets)
    profits.fields = Symbol.(strat.universe.assets)
    values.fields = Symbol.(strat.universe.assets)
    holdings.values[isnan.(holdings.values)] .= 0.0
    values.values[isnan.(values.values)] .= 0.0
    profits.values[isnan.(profits.values)] .= 0.0
    # portfolio weights, net exposure, and leverage calculations
    weights = values / apply(values, 1, fun=sum)
    weights.fields = Symbol.(strat.universe.assets)
    exposure = apply(weights, 1, fun=sum)
    leverage = apply(weights, 1, fun=x->(sum(abs.(x))))
    weights = [weights exposure leverage]
    weights.fields[end-1:end] = [:Exposure, :Leverage]
    # compute other temporal totals and return
    profits = [profits apply(profits, 1, fun=sum)]
    profits.fields[end] = :Total
    values = [values apply(values, 1, fun=sum)]
    values.fields[end] = :Total
    return weights, holdings, values, profits
end


text:
include("models/data.jl")

mordat = [ population_deaths(popdatbl; bundesland = bl, alter=a, geschlecht=g, jahr=j, deaths=(kw=kw,) )
           # or switch to interpolate deaths for kw:
           # population_deaths(alter=a, geschlecht=g, jahr=j, kw=kw )
           for bl in DeutschlandDaten.Bundesländer,
               a in DeutschlandDaten.Altersgruppen_Bundesländer,
               g in ["Männlich", "Weiblich"],
               j in 2000:2021,
               kw in 1:52 ] |> DataFrame;
filter!(x-> x.D !==missing && x.D>0, mordat)

@transform! mordat :date = endofweek.(:jahr,:kw)
dataprep!(mordat)

rowgroup=[:bundesland, :alter, :geschlecht]
@transform!(groupby(mordat, rowgroup),
            :PD_12months = gliding_mean(:PD, div(52, 2)),
            # :PD_rel2000 = :PD ./ first(skipmissing(:PD_12months)),
            :PD_relmean = :PD ./ mean(skipmissing(:PD))
#            :PD_rel_12months = :PD_12months ./ first(skipmissing(:PD_12months)),
#            :PD_relmean_12months = :PD_12months ./ mean(skipmissing(:PD_12months)),
            )

########################################
title = "Rythmen der Übersterblichkeit in %\nBundesländer Deutschland"
imagefile = "relative_Übersterblichkeit_Bundesländer"

## sortiere Bundesländer nach durchschnittlicher Übersterblichkeit 2020.6-2021.6
@transform!(groupby(mordat, :bundesland),
            :PDrelmean2021 = mean(:PD_relmean[(:date .< Date(2021,06,01)) .& (:date .> Date(2020,06,01))])
            )
sort!(mordat, [:PDrelmean2021, rowgroup..., :date])

@transform!(mordat,
            :row = map(x->tuple(x...), eachrow(mordat[:,rowgroup])),
            :col = :date,
            :value = 100 .* :PD_relmean
            )

function plot_bl_PDrel(mordat)
    rs = unique(mordat.date)
    xticks = (1:52:length(rs), string.(year(minimum(mordat.date)):1:year(maximum(mordat.date))))
    p = plot_unstacked_heatmap(
        mordat,
        title = title,
        labelf = x -> let (b, a, g) = x; (a.start == 75 && g!="Weiblich" ? (b*"") : ""); end, # * (g!="Männlich" ? alterstring(a) : "") * label_geschlecht[g]; end,
        ylabel = "Bundesländer (Alter   ♂♀)",
        xticks = xticks,
        xrotation = -90,
        xmirror=true,
        c = tot_cgrad,
        top_margin = 20px,
        bottom_margin = 20px,
        w=1024px,
        h=1500px
        ##kw...
    ) do p, m
        @show blstarts = [ i-.5 for i in 1:size(m)[1] if i == 1 || m[i,1][1] != m[i-1,1][1] ]
        #@show m[blstarts,1]
        plot!(p,blstarts, seriestype = :hline, color = :black, legend=nothing)
        plot!(p, [findfirst(e -> e > Date(2021,1,1), rs)], alpha=.9, width=.5,
              seriestype=:vline)
        p
    end
end

p = plot_bl_PDrel(mordat)

DeutschlandDaten.savefigs(p,imagefile)

p = plot_bl_PDrel(filter(x -> x.date>Date(2016,1,1), mordat ))
DeutschlandDaten.savefigs(p,imagefile*"2016_2022")

max_stars_repo_name: gkappler/DeutschlandDaten

max_issues_repo_name: gkappler/DeutschlandDaten

max_forks_repo_name: gkappler/DeutschlandDaten


text:
libcuplus=dlopen("libcuplus")

# Show cuda info
function cuda_info()
    ccall(dlsym(libcuplus, :cuda_info),
          Void, ())
end

# Allocate device memory
function cuda_malloc{T}(T::Type{T}, count::Integer)
    ptr = ccall(dlsym(libcuplus, :cuda_malloc),
                Ptr{T}, (Int32, ),
                count * sizeof(T))
    if ptr == C_NULL
        # If unable to allocate, call garbage collector
        gc()
        # Try allocating again
        ptr = ccall(dlsym(libcuplus, :cuda_malloc),
                    Ptr{T}, (Int32, ),
                    count * sizeof(T))
    end

    if ptr == C_NULL
        # If still unable to allocate, error out
        error("Can not allocate the required memory")
    end
    return ptr
end
cuda_malloc(T::Type, rows::Integer, cols::Integer) = cuda_malloc(T, rows * cols)
cuda_malloc(T::Type, dims::(Integer, Integer)) = cuda_malloc(T, dims[1] * dims[2])

# Free device memory
function cuda_free{T}(ptr::Ptr{T})
    ccall(dlsym(libcuplus, :cuda_free),
          Void, (Ptr{T}, ),
          ptr)
end

# Copy from Host to Device Memory
function mem_device{T}(dst::Ptr{T}, src::Matrix{T})
    bytes::Int32 = sizeof(T) * numel(src)
    ccall(dlsym(libcuplus, :cuda_memcpy_h2d),
          Void, (Ptr{T}, Ptr{T}, Int32),
          dst, src, bytes)
end

# Copy from Device to Host Memory
function mem_host{T}(dst::Matrix, src::Ptr{T})
    bytes::Int32 = sizeof(T) * numel(dst)
    ccall(dlsym(libcuplus, :cuda_memcpy_d2h),
          Void, (Ptr{T}, Ptr{T}, Int32),
          dst, src, bytes)
end

# Copy from Device to Device Memory
function mem_copy{T}(dst::Ptr{T}, src::Ptr{T}, bytes::Int32)
    ccall(dlsym(libcuplus, :cuda_memcpy_d2d),
          Void, (Ptr{T}, Ptr{T}, Int32),
          dst, src, bytes)
end

# Create uniform random Matrix
for (fname, elty) in ((:cudaSrand, :Float32),
                      (:cudaDrand, :Float64),
                      (:cudaCrand, :Complex64),
                      (:cudaZrand, :Complex128))
    @eval begin
        function cuda_rand(ptr::Ptr{$elty}, count::Integer)
            count = convert(Int32, count)
            ccall(dlsym(libcuplus, $string(fname)),
                  Void, (Ptr{$elty}, Int32),
                  ptr, count)
        end
    end
end

# Create normal random Matrix
for (fname, elty) in ((:cudaSrandn, :Float32),
                      (:cudaDrandn, :Float64),
                      (:cudaCrandn, :Complex64),
                      (:cudaZrandn, :Complex128))
    @eval begin
        function cuda_randn(ptr::Ptr{$elty}, count::Integer)
            count = convert(Int32, count)
            ccall(dlsym(libcuplus, $string(fname)),
                  Void, (Ptr{$elty}, Int32),
                  ptr, count)
        end
    end
end

max_stars_repo_name: stefan-k/julia-CuMatrix

max_issues_repo_name: stefan-k/julia-CuMatrix

max_forks_repo_name: stefan-k/julia-CuMatrix


text:
mutable struct GroupStatsBoolVariate <:AbstractBoolVariate
    name::String
    length::Int64
    groupstats::GroupStats
    transform::Function
end

function slice(boolvar::GroupStatsBoolVariate, fromobs::Integer, toobs::Integer, slicelength::Integer)
    factors = boolvar.groupstats.factors
    dims = Tuple(boolvar.groupstats.dims) #Tuple([length(getlevels(factor)) + 1 for factor in factors])
    k = length(factors)
    f = boolvar.transform
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    factorslices = zip([slice(factor, fromobs, toobs, slicelength) for factor in factors])
    groupstats = boolvar.groupstats
    subindexmap = covariate.groupstats.subindexmap

    subind = Vector{Int64}(undef, k)
    buffer = BitArray{1}(undef, slicelength)
    nanstats = CovariateStats(0, 0, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64, NaN64)
    map(factorslices, Bool) do slicevector
       len = length(slicevector[1])
       for i in 1:len
        for j in 1:k
            subind[j] = slicevector[j][i]
            if subind[j] > 0
                mappedindex = subindexmap[j][subind[j]]
                subind[j] = mappedindex == 0 ? 0 : mappedindex + 1
             else
                 subind[j] = 1
            end
        end
           if findfirst(subind, 0) > 0
               buffer[i] = f(nanstats)
           else
                ind = sub2ind(dims, subind...)
                mappedind = get(groupstats.linindexmap, ind, 0)
                if mappedind > 0
                    buffer[i] = f(groupstats.stats[mappedind][2])
                else
                    buffer[i] = f(nanstats)
                end
           end
       end
       view(buffer, 1:len)
    end
end

function Base.map(boolvariate::GroupStatsBoolVariate, dataframe::AbstractDataFrame)
    GroupStatsBoolVariate(boolvariate.name, length(dataframe), map(boolvariate.groupstats, dataframe), boolvariate.transform)
end

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
# Note: We can't call the method below `Trixi.include` since that is created automatically
# inside `module Trixi` to `include` source files and evaluate them within the global scope
# of `Trixi`. However, users will want to evaluate in the global scope of `Main` or something
# similar to manage dependencies on their own.
"""
    trixi_include([mod::Module=Main,] elixir::AbstractString; kwargs...)

`include` the file `elixir` and evaluate its content in the global scope of module `mod`.
You can override specific assignments in `elixir` by supplying keyword arguments.
It's basic purpose is to make it easier to modify some parameters while running Trixi from the
REPL. Additionally, this is used in tests to reduce the computational burden for CI while still
providing examples with sensible default values for users.

# Examples

```jldoctest
julia> redirect_stdout(devnull) do
         trixi_include(@__MODULE__, joinpath(examples_dir(), "1d", "elixir_advection_extended.jl"),
                       tspan=(0.0, 0.1))
         sol.t[end]
       end
0.1
```
"""
function trixi_include(mod::Module, elixir::AbstractString; kwargs...)
  Base.include(ex -> replace_assignments(ex; kwargs...), mod, elixir)
end

trixi_include(elixir::AbstractString; kwargs...) = trixi_include(Main, elixir; kwargs...)

"""
    convergence_test([mod::Module=Main,] elixir::AbstractString, iterations; kwargs...)

Run `iterations` Trixi simulations using the setup given in `elixir` and compute
the experimental order of convergence (EOC) in the ``L^2`` and ``L^\\infty`` norm.
In each iteration, the resolution of the respective mesh will be doubled.
Additional keyword arguments `kwargs...` and the optional module `mod` are passed directly
to [`trixi_include`](@ref).
"""
function convergence_test(mod::Module, elixir::AbstractString, iterations; kwargs...)
  @assert(iterations > 1, "Number of iterations must be bigger than 1 for a convergence analysis")

  # Types of errors to be calcuated
  errors = Dict(:l2 => Float64[], :linf => Float64[])

  initial_resolution = extract_initial_resolution(elixir)

  # run simulations and extract errors
  for iter in 1:iterations
    println("Running convtest iteration ", iter, "/", iterations)

    include_refined(mod, elixir, initial_resolution, iter; kwargs)

    l2_error, linf_error = mod.analysis_callback(mod.sol)

    # collect errors as one vector to reshape later
    append!(errors[:l2],   l2_error)
    append!(errors[:linf], linf_error)

    println("\n\n")
    println("#"^100)
  end

  # number of variables
  _, equations, _, _ = mesh_equations_solver_cache(mod.semi)
  variablenames = varnames(cons2cons, equations)
  nvariables = length(variablenames)

  # Reshape errors to get a matrix where the i-th row represents the i-th iteration
  # and the j-th column represents the j-th variable
  errorsmatrix = Dict(kind => transpose(reshape(error, (nvariables, iterations))) for (kind, error) in errors)

  # Calculate EOCs where the columns represent the variables
  # As dx halves in every iteration the denominator needs to be log(1/2)
  eocs = Dict(kind => log.(error[2:end, :] ./ error[1:end-1, :]) ./ log(1 / 2) for (kind, error) in errorsmatrix)

  eoc_mean_values = Dict{Symbol,Any}()
  eoc_mean_values[:variables] = variablenames

  for (kind, error) in errorsmatrix
    println(kind)

    for v in variablenames
      @printf("%-20s", v)
    end
    println("")

    for k = 1:nvariables
      @printf("%-10s", "error")
      @printf("%-10s", "EOC")
    end
    println("")

    # Print errors for the first iteration
    for k = 1:nvariables
      @printf("%-10.2e", error[1, k])
      @printf("%-10s", "-")
    end
    println("")

    # For the following iterations print errors and EOCs
    for j = 2:iterations
      for k = 1:nvariables
        @printf("%-10.2e", error[j, k])
        @printf("%-10.2f", eocs[kind][j-1, k])
      end
      println("")
    end
    println("")

    # Print mean EOCs
    mean_values = zeros(nvariables)
    for v in 1:nvariables
      mean_values[v] = sum(eocs[kind][:, v]) ./ length(eocs[kind][:, v])
      @printf("%-10s", "mean")
      @printf("%-10.2f", mean_values[v])
    end
    eoc_mean_values[kind] = mean_values
    println("")
    println("-"^100)
  end

  return eoc_mean_values
end

convergence_test(elixir::AbstractString, iterations; kwargs...) = convergence_test(Main, elixir::AbstractString, iterations; kwargs...)

# Helper methods used in the functions defined above

# Apply the function `f` to `expr` and all sub-expressions recursively.
walkexpr(f, expr::Expr) = f(Expr(expr.head, (walkexpr(f, arg) for arg in expr.args)...))
walkexpr(f, x) = f(x)

# Replace assignments to `key` in `expr` by `key = val` for all `(key,val)` in `kwargs`.
function replace_assignments(expr; kwargs...)
  # replace explicit and keyword assignemnts
  expr = walkexpr(expr) do x
    if x isa Expr
      for (key,val) in kwargs
        if (x.head === Symbol("=") || x.head === :kw) && x.args[1] === Symbol(key)
          x.args[2] = :( $val )
          # dump(x)
        end
      end
    end
    return x
  end

  return expr
end

# find a (keyword or common) assignment to `destination` in `expr`
# and return the assigned value
function find_assignment(expr, destination)
  # declare result to be able to assign to it in the closure
  local result

  # find explicit and keyword assignemnts
  walkexpr(expr) do x
    if x isa Expr
      if (x.head === Symbol("=") || x.head === :kw) && x.args[1] === Symbol(destination)
        result = x.args[2]
        # dump(x)
      end
    end
    return x
  end

  result
end

# searches the parameter that specifies the mesh reslution in the elixir
function extract_initial_resolution(elixir)
  code = read(elixir, String)
  expr = Meta.parse("begin $code end")

  try
    # get the initial_refinement_level from the elixir
    initial_refinement_level = find_assignment(expr, :initial_refinement_level)
  catch e
    if isa(e, UndefVarError)
      # get cells_per_dimension from the elixir
      cells_per_dimension = eval(find_assignment(expr, :cells_per_dimension))
    else
      throw(e)
    end
  end
end

# runs the specified elixir with a doubled resolution each time iter is increased by 1
# works for TreeMesh
function include_refined(mod, elixir, initial_refinement_level::Int, iter; kwargs)
  trixi_include(mod, elixir; kwargs..., initial_refinement_level=initial_refinement_level+iter-1)
end

# runs the specified elixir with a doubled resolution each time iter is increased by 1
# works for CurvedMesh
function include_refined(mod, elixir, cells_per_dimension::NTuple{NDIMS, Int}, iter; kwargs) where {NDIMS}
  new_cells_per_dimension = cells_per_dimension .* 2^(iter - 1)

  trixi_include(mod, elixir; kwargs..., cells_per_dimension=new_cells_per_dimension)
end


text:
# Autogenerated wrapper script for MLIR_jll for x86_64-linux-musl-cxx11
export libMLIR

using Zlib_jll
using libLLVM_jll
JLLWrappers.@generate_wrapper_header("MLIR")
JLLWrappers.@declare_library_product(libMLIR, "libMLIR.so.13jl")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, libLLVM_jll)
    JLLWrappers.@init_library_product(
        libMLIR,
        "lib/libMLIR.so",
        nothing,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    @test_exists [ntests = 100][maxtests = 1000][logto = ""][argument_data, prop]

Generates sample of values for arguments in proposition. Returns `Pass`
if proposition holds for at least one value in the sample, `Fail` if proposition
does not hold for any values, and `Error` if a part of code could not
be evaluated. `Error` is given if after `maxtests` attempts to generate
argument values, number of arguments satisfying `condition` part of `prop`
is less than `ntests`.

# Description
Similar to `@test_formany`, see it's documentation for details.

# Examples
```julia
julia> @test_exists ntests = 1000 -10<x<10, x^2>99
Test Passed
  Expression: (:(ntests = 1000),:((-10 < x < 10,x ^ 2 > 99)))

## values are generated successively in the oder x->y->z, they satisfy x<y<z automatically
julia> @test_exists 0<x<2*pi, x<y<2*pi, y<z<2*pi, sin(x)<sin(y)<sin(z)
Test Passed
  Expression: (:((0 < x < 2pi,x < y < 2pi,y < z < 2pi,sin(x) < sin(y) < sin(z))),)

## values are generated independently, the statement is checked only for those that satisfy x<y<z
julia> @test_exists 0<x<2*pi, 0<y<2*pi, 0<z<2*pi, x<y<z --> sin(x)<sin(y)<sin(z) ntests = 1000
Test Passed
  Expression: (:((0 < x < 2pi,0 < y < 2pi,0 < z < 2pi,$(Expr(:-->, :(x < y < z), :(sin(x) < sin(y) < sin(z)))))),:(ntests = 1000))

## Gives Error, because condition is too restrictive
@test_exists 0<x<1000, x>990 --> x>999.9 ntests = 100
Error During Test
  Test threw an exception of type ErrorException
  Expression: (:((0 < x < 1000,$(Expr(:-->, :(x > 990), :(x > 999.9))))),:(ntests = 100))
  Found only 14/100 values satisfying given condition.
```
"""
macro test_exists(exprs...)
inner_ex = try
	maxtests = 0; ntests = 0
	prop = :(); cond = :()
	logging  = ""
	var_data = Array(Any,1,6)
	for ex in exprs
		if isa(ex, Expr) && ex.head == :tuple
			var_data = parse_argument_data(ex)
			if ex.args[end].head == :-->
				prop = ex.args[end].args[2]
				cond = ex.args[end].args[1]
			else
				prop = ex.args[end]
				cond = true
			end
		elseif isa(ex, Expr) && ex.head == :(=)
			if ex.args[1] == :ntests
				ntests = ex.args[2]
			elseif ex.args[1] == :maxtests
				maxtests = ex.args[2]
			elseif ex.args[1]== :logto
				logging = ex.args[2]
			else
				error("Invalid macro input $ex.")
			end
		else
			error("Invalid macro input $ex.")
		end
	end
	if ntests == 0
		ntests = 100
	end
	if maxtests == 0
		maxtests = 10*ntests
	end
	num_of_vars = size(var_data,1)
	generate_values = :()
	values = :([])
	for i in 1:num_of_vars
		next_expr = :($(esc(var_data[i,1]))=
			   custom_generator($(esc(var_data[i,2])),$(esc(var_data[i,3])),
			      $(esc(var_data[i,4])),$(esc(var_data[i,5])),$(esc(var_data[i,6])))(div(n,2)+3))
		generate_values = Expr(:block,generate_values.args...,next_expr)
		values = Expr(:vect, values.args..., :($(esc(var_data[i,1]))))
	end
#####
	if logging == ""    ## no logging
		inner_ex = quote try
			for n in 1:$maxtests
				$generate_values
				if $(esc(cond))
					num_good_args += 1
					res = res || $(esc(prop))
					if res
						break
					end
					if num_good_args >= $ntests
						break
					end
				end
			end
		        # no example but insufficient amount of tests
			if !res && num_good_args<$ntests
				nt = $ntests
 				error("Found only $num_good_args/$nt values satisfying given condition.")
			end
			res ? Pass(:test,$(Expr(:quote, exprs)), nothing, nothing) :
		   	   Fail(:test,$(Expr(:quote, exprs)), fail_data, nothing)
		catch err
			Error(:test_error,$(Expr(:quote, exprs)), err, catch_backtrace())
		end end #quote #try
	else 		## logging is a path to file where the log is written to
		inner_ex = quote try
			log_file = open($logging,"a")
			writedlm(log_file,reshape([string(s) for s in $(var_data[:,1])],(1,$num_of_vars)),",")
			for n in 1:$maxtests
				$generate_values
				if $(esc(cond))
					num_good_args += 1
					res = res || $(esc(prop))
					writedlm(log_file,transpose(push!(convert(Vector{Any},$values),res)),",")
					if res
						break
					end
					if num_good_args >= $ntests
						break
					end
				end
			end
			close(log_file)
		        # no example but insufficient amount of tests
			if !res && num_good_args<$ntests
				nt = $ntests
 				error("Found only $num_good_args/$nt values satisfying given condition.")
			end
			res ? Pass(:test,$(Expr(:quote, exprs)), nothing, nothing) :
		   	   Fail(:test,$(Expr(:quote, exprs)), fail_data, nothing)
		catch err
			Error(:test_error,$(Expr(:quote, exprs)), err, catch_backtrace())
		end end #quote #try
	end
	inner_ex
#####
   catch err
	inner_ex = quote
		Error(:test_error,$(Expr(:quote, exprs)), $err, catch_backtrace())
	end
   end # try defining inner_ex
	return quote
		num_good_args = 0
		res = false
		fail_data = true
		result = $inner_ex
		record(get_testset(), result)
	end
end


text:
using RenewalProcesses
using Base.Test


text:
function log1pexp!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, res2::AbstractArray{Float64}, n::Int64=length(x))

    copy!(res2, x)
    Yeppp.exp!(res, x)
    log1p!(res, res, n)
    for i in 1:n
        @inbounds res[i]=mylog1pexp(res[i], res2[i])
    end
    # for i in 1:n
    #     @inbounds res[i]=log1pexp(x[i])
    # end
end
mylog1pexp(x::AbstractFloat, x2::AbstractFloat)=x2 < 18.0 ? x : x2 < 33.3 ? x2+exp(-x2) : x2
function log1pexp!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, n::Int64=length(x))

    res2=copy(x)
    log1pexp!(res, x, res2, n)
end
# log1pexp(x::AbstractFloat) = x < 18.0 ? log1p(exp(x)) : x < 33.3 ? x + exp(-x) : x
log1pexp!(x::AbstractArray{Float64}, n::Int64=length(x))=log1pexp!(x, x, n)

function sumexp(x::AbstractArray{T}) where T<:Real
    isempty(x) && return -Inf
    u = maximum(x)
    s = 0.
    for i in 1:length(x)
        @inbounds s += exp(x[i] - u)
    end
    s * exp(u)
end
function sumexp(x::AbstractArray{T}, coef::AbstractArray{T}) where T<:Real
    isempty(x) && return -Inf
    u = maximum(x)
    s = 0.
    for i in 1:length(x)
        @inbounds s += coef[i]*exp(x[i] - u)
    end
    s * exp(u)
end

function ratiosumexp!(x::AbstractArray{T}, coef::AbstractArray{T}, s::AbstractArray{T}, ncomponent::Int) where T<:Real
    #length(x) != length(coef) && error("Length should be the same!")
    #isempty(x) && return -Inf
    u = maximum(x)
    for i in 1:ncomponent
        @inbounds s[i] = coef[i]*exp(x[i] - u)
    end
    divide!(s, s, sum(s), ncomponent)
    #s ./ sum(s) #* exp(u)
    nothing
end
function ratiosumexp!(x::AbstractArray{T}, coef::AbstractArray{T}, s::AbstractMatrix{T}, irow::Int, ncomponent::Int) where T<:Real

    u = maximum(x)
    ssum=0.0
    @inbounds for i in 1:ncomponent
        tmp=coef[i]*exp(x[i] - u)
        s[irow, i] = tmp
        ssum += tmp
    end
    #ssum = sum(s[irow, :])
    for i in 1:ncomponent
        @inbounds s[irow, i] = s[irow, i] / ssum
    end
    nothing
end
function add!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, y::Float64, n::Int64=length(x))
   for i in 1:n
       @inbounds res[i] = x[i] + y
   end
   nothing
end

add!(res::AbstractArray{Float64}, x::Float64, y::AbstractArray{Float64}, n::Int64=length(y)) = add!(res, y, x, n)

add!(x::AbstractArray{Float64}, y::Float64, n::Int64=length(x))=add!(x, x, y, n)
add!(x::Float64, y::AbstractArray{Float64}, n::Int64=length(y))=add!(y, y, x, n)

plusone!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, n::Int64=length(x)) = add!(res, x, 1.0, n)
plusone!(x::AbstractArray{Float64}, n=length(x)) = plusone!(x, x, n)

function divide!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, y::AbstractArray{Float64}, n::Int64=length(x))
    for i in 1:n
        @inbounds res[i] = x[i] / y[i]
    end
    nothing
end

divide!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, y::Float64, n::Int64=length(x)) = multiply!(res, x, 1/y, n)
function divide!(res::AbstractArray{Float64}, x::Float64, y::AbstractArray{Float64}, n::Int64=length(y))
    for i in 1:n
        @inbounds res[i] = x / y[i]
    end
    nothing
end
divide!(x::AbstractArray{Float64}, y::AbstractArray{Float64}, n::Int64=length(x)) = divide!(x, x, y, n)
# divide!(x::Float64, y::AbstractArray{Float64}, n::Int64=length(y)) = divide!(y, x, y, n)
divide!(x::AbstractArray{Float64}, y::Float64, n::Int64=length(x)) = divide!(x, x, y, n)

rcp!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, n::Int64=length(x)) = divide!(res, 1.0, x, n)
rcp!(x::AbstractArray{Float64}, n::Int64=length(x))=rcp!(x, x, n)

function multiply!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, y::Float64, n::Int64=length(x))
    for i in 1:n
        @inbounds res[i] = x[i] * y
    end
    nothing
end

function multiply!(res::AbstractArray{Float64}, x::Float64, y::AbstractArray{Float64}, n::Int64=length(y))
    for i in 1:n
        @inbounds res[i] = x * y[i]
    end
    nothing
end

# multiply!(x::Float64, y::AbstractArray{Float64}, n::Int64=length(y)) = multiply!(y, x, y, n)
multiply!(x::AbstractArray{Float64}, y::Float64, n::Int64=length(x)) = multiply!(x, x, y, n)

function sqr!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, n::Int64=length(x))
   for i in 1:n
       @inbounds res[i] = x[i]*x[i]
   end
   nothing
end
sqr!(x::AbstractArray{Float64}, n::Int64=length(x)) = sqr!(x, x, n)

function negate!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, n::Int64=length(x))
   for i in 1:n
       @inbounds res[i] = -x[i]
   end
   nothing
end
negate!(x::AbstractArray{Float64}, n::Int64=length(x)) = negate!(x, x, n)

function negateiffalse!(x::AbstractArray{Float64}, y::AbstractArray{Bool, 1}, n::Int64=length(x))
    for i in 1:n
        @inbounds x[i] = ifelse(y[i], x[i], -x[i])
    end
    nothing
end

function negateiftrue!(x::AbstractArray{Float64}, y::AbstractArray{Bool, 1}, n::Int64=length(x))
    for i in 1:n
        @inbounds x[i] = ifelse(y[i], -x[i], x[i])
    end
    nothing
end

function log1p!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, n::Int64=length(x))
    add!(res, x, 1., n)
    Yeppp.log!(res, res)
    res
end
log1p!(x::AbstractArray{Float64}, n::Int64=length(x)) = log1p!(x, x, n)

# -log(1+exp(-xy))
function loglogistic!(x::AbstractArray{Float64}, y::AbstractArray{Bool, 1}, n::Int=length(x))
    negateiftrue!(x, y, n)
    log1pexp!(x, x, n)
    negate!(x, x, n)
    nothing
end

# 1/(1+exp(-x))
function logistic!(res::AbstractArray{Float64}, x::AbstractArray{Float64}, n::Int=length(x))

    negate!(res, x, n)
    Yeppp.exp!(res, res)
    add!(res, res, 1.0, n)
    rcp!(res, res, n)
    res
end
logistic!(x::AbstractArray{Float64}, n::Int=length(x)) = logistic!(x, x, n)

function relocate!(res::AbstractArray{Float64}, ga::Vector{Float64}, groupindex::IntegerVector, N::Int)
    for i in 1:N
        @inbounds res[i] = ga[groupindex[i]]
    end
    nothing
end

function sumby!(r::AbstractArray, y::AbstractArray, x::IntegerArray, levels::IntUnitRange)
	k = length(levels)
	length(r) == k || raise_dimerror()

	m0 = levels[1]
	m1 = levels[end]
	b = m0 - 1

	@inbounds for i in 1 : length(x)
		xi = x[i]
		if m0 <= xi <= m1
			r[xi - b] += y[i]
		end
	end
	#return r
end
function sumsqby!(r::AbstractArray, y::AbstractArray, x::IntegerArray, levels::IntUnitRange)
	k = length(levels)
	length(r) == k || raise_dimerror()

	m0 = levels[1]
	m1 = levels[end]
	b = m0 - 1

	@inbounds for i in 1 : length(x)
		xi = x[i]
		if m0 <= xi <= m1
			r[xi - b] += abs2(y[i])
		end
	end
	#return r
end

function H1(y, mu, sigmas)
    (y .- mu)./sigmas./sigmas
end
function H2(y, mu, sigmas)
    z = (y .- mu)./sigmas
    (z .^ 2 .-1) ./ sigmas^2 ./ 2
end

function H3(y, mu, sigmas)
    z = (y .- mu)./sigmas
    (z .^3 .- 3 .* z) ./ sigmas^3 ./ 6
end

function H4(y, mu, sigmas)
    z = (y .- mu)./sigmas
    (z .^4 .- 6 .* z .^ 2 .+ 3) ./ sigmas^4 ./ 24
end