text:
import J4Org:
    create_documented_item,
    tokenized,
    find_tag,
    org_string_comment,
    org_string_code,
    tag_idx

@testset "documented_item" begin
    filename="$(dirname(@__FILE__))/code_examples/basic.jl"

    code=readcode(filename)
    t=tokenized(code)

    r = find_tag(t,1)
end;

# sigsegv... with 0.7 beta ???
if VERSION<v"0.7.0-alpha"
    @testset "di_real_code_1" begin
        filename="$(dirname(@__FILE__))/code_examples/real_code_1.jl"
        code=readcode(filename)
        t=tokenized(code)

        r=find_tag(t,1)
        di=create_documented_item(t,tag_idx(r),filename=filename)
        @test org_string_comment(di,[di],[di],"","BoxingModule") == "#+BEGIN_QUOTE\nA *central* structure containing documented item\n#+END_QUOTE\n"
        @test org_string_code(di) == "#+BEGIN_SRC julia :eval never :exports code\nstruct Documented_Item\n#+END_SRC\n"

        r=find_tag(t,8)
        di=create_documented_item(t,tag_idx(r),filename=filename)
        @test org_string_comment(di,[di],[di],"","BoxingModule") == ""
        @test org_string_code(di) == "#+BEGIN_SRC julia :eval never :exports code\nfilename(di::Documented_Item)::String\n#+END_SRC\n"
    end;
end

@testset "di_real_code_problematic_1" begin
    filename="$(dirname(@__FILE__))/code_examples/real_code_problematic_1.jl"

    tok=tokenized(String(read(filename)))

    idx=1
    extracted_tag = extract_tag(tok,idx)
    @assert extracted_tag != nothing
    idx=idx+1
    doc=extract_comment(tok,idx)
end;

@testset "enum" begin
    filename="$(dirname(@__FILE__))/code_examples/enum.jl"
    code=readcode(filename)
    t=tokenized(code)

    r=find_tag(t,1)
    di=create_documented_item(t,tag_idx(r),filename=filename)
    @test org_string_comment(di,[di],[di],"","BoxingModule") == "#+BEGIN_QUOTE\nThis is an enum example [[target][@enum Alpha]]\n#+END_QUOTE\n"
    @test org_string_code(di) == "#+BEGIN_SRC julia :eval never :exports code\n@enum Alpha A  B=1 C\n#+END_SRC\n"
end;

@testset "variable" begin
    filename="$(dirname(@__FILE__))/code_examples/variable.jl"
    code=readcode(filename)
    t=tokenized(code)

    r=find_tag(t,1)
    di=create_documented_item(t,tag_idx(r),filename=filename)
    @test org_string_comment(di,[di],[di],"","BoxingModule") == "#+BEGIN_QUOTE\nThis is an variable example [[target][variable A]]\n#+END_QUOTE\n"
    @test org_string_code(di) == "#+BEGIN_SRC julia :eval never :exports code\nA\n#+END_SRC\n"
end;

@testset "macro" begin
    filename="$(dirname(@__FILE__))/code_examples/macro.jl"
    code=readcode(filename)
    t=tokenized(code)

    r=find_tag(t,1)
    di=create_documented_item(t,tag_idx(r),filename=filename)
    @test org_string_comment(di,[di],[di],"","BoxingModule") == "#+BEGIN_QUOTE\nThis is an macro example [[target][@swap]]\n#+END_QUOTE\n"
    @test org_string_code(di) == "#+BEGIN_SRC julia :eval never :exports code\nmacro swap(x,y)\n#+END_SRC\n"
end;


text:
using Plots
using KernelAbstractions: CPU
using ClimateMachine.MPIStateArrays: array_device
using ClimateMachine.BalanceLaws: Prognostic, Auxiliary, GradientFlux

"""
    plot_friendly_name(ϕ)

Get plot-friendly string, since many Unicode
characters do not render in plot labels.
"""
function plot_friendly_name(ϕ)
    s = ϕ
    s = replace(s, "ρ" => "rho")
    s = replace(s, "α" => "alpha")
    s = replace(s, "∂" => "partial")
    s = replace(s, "∇" => "nabla")
    return s
end

"""
    export_plot(z, all_data, ϕ_all, filename;
                xlabel, ylabel, time_data, round_digits, horiz_layout)

Export plot of all variables, or all
available time-steps in `all_data`.
"""
function export_plot(
    z,
    all_data::Array,
    ϕ_all,
    filename;
    xlabel,
    ylabel,
    time_data,
    round_digits = 2,
    horiz_layout = false,
)
    ϕ_all isa Tuple || (ϕ_all = (ϕ_all,))
    single_var = ϕ_all[1] == xlabel && length(ϕ_all) == 1
    p = plot()
    for (t, data) in zip(time_data, all_data)
        for ϕ in ϕ_all
            ϕ_string = String(ϕ)
            ϕ_name = plot_friendly_name(ϕ_string)
            ϕ_data = data[ϕ_string][:]
            label = single_var ? "t=$(round(t, digits=round_digits))" :
                "$(ϕ_string), t=$(round(t, digits=round_digits))"
            if !horiz_layout
                plot!(
                    ϕ_data,
                    z;
                    xlabel = xlabel,
                    ylabel = ylabel,
                    label = label,
                )
            else
                plot!(
                    z,
                    ϕ_data;
                    xlabel = xlabel,
                    ylabel = ylabel,
                    label = label,
                )
            end
        end
    end
    savefig(filename)
end

function save_binned_surface_plots(
    x,
    y,
    z,
    title,
    filename,
    n_plots = (3, 3),
    z_label_prefix = "z",
    n_digits = 5,
)
    n_z_partitions = prod(n_plots)
    z_min_global = min(z...)
    z_max_global = max(z...)
    Δz = (z_max_global - z_min_global) / n_z_partitions
    z_min = ntuple(i -> z_min_global + (i - 1) * Δz, n_z_partitions)
    z_max = ntuple(i -> z_min_global + (i) * Δz, n_z_partitions)
    p = []
    for i in 1:n_z_partitions
        mask = z_min[i] .<= z .<= z_max[i]
        x_i = x[mask]
        y_i = y[mask]
        sz_min = string(z_min[i])[1:min(n_digits, length(string(z_min[i])))]
        sz_max = string(z_max[i])[1:min(n_digits, length(string(z_max[i])))]
        p_i = plot(
            x_i,
            y_i,
            title = "$(title), in ($sz_min, $sz_max)",
            seriestype = :scatter,
            markersize = 5,
        )
        push!(p, p_i)
    end
    plot(p..., layout = n_plots, legend = false)
    savefig(filename)
end;

state_prefix(::Prognostic) = "prog_"
state_prefix(::Auxiliary) = "aux_"
state_prefix(::GradientFlux) = "grad_flux_"

"""
    plot_results(solver_config, all_data, time_data, output_dir)

Exports plots of states given
 - `solver_config` a `SolverConfiguration`
 - `all_data` an array of dictionaries, returned from `dict_of_nodal_states`
 - `time_data` an array of time values
 - `output_dir` output directory
"""
function export_state_plots(
    solver_config,
    all_data,
    time_data,
    output_dir;
    state_types = (Prognostic(), Auxiliary()),
)
    FT = eltype(solver_config.Q)
    z = get_z(solver_config.dg.grid)
    z = array_device(solver_config.Q) isa CPU ? z : Array(z)
    mkpath(output_dir)
    for st in state_types
        vs = vars_state(solver_config.dg.balance_law, st, FT)
        for fn in flattenednames(vs)
            file_name = state_prefix(st) * replace(fn, "." => "_")
            export_plot(
                z,
                all_data,
                (fn,),
                joinpath(output_dir, "$(file_name).png");
                xlabel = fn,
                ylabel = "z [m]",
                time_data = time_data,
                round_digits = 5,
            )
        end
    end
end

ext: jl

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
# TODO move this back to GraphLayout

function _tree_rectangle(x, y, rect)
    SimpleRectangle(rect.x + x - rect.w/2, rect.y + y - rect.w/2, rect.w, rect.h)
end

"""
    Creates an arrow between two rectangles in the tree

    Arguments:
    o_x, o_y, o_h   Origin x, y, and height
    d_x, d_y, d_h   Destination x, y, and height
"""
function _arrow_tree{T}(origin::Point{2,T}, o_h::T, destination::Point{2,T}, d_h::T)
    x1, y1 = origin[1], origin[2] + o_h/2
    x2, y2 = destination[1], destination[2] - d_h/2
    Δx, Δy = x2 - x1, y2 - y1
    θ = atan2(Δy, Δx)
    # Put an arrow head only if destination isn't dummy
    head = d_h != 0 ? _arrow_heads(θ, x2, y2, T(2)) : []
    start = Point{2,T}(x1, y1)
    endpoint = Point{2,T}(x2, y2)
    LineSegment{Point{2,T}}[LineSegment{Point{2,T}}(start, endpoint), head...]
end
"""
    Creates an arrow head given the angle of the arrow and its destination.

    Arguments:
    θ               Angle of arrow (radians)
    dest_x, dest_y  End of arrow
    λ               Length of arrow head tips
    ϕ               Angle of arrow head tips relative to angle of arrow
"""
function _arrow_heads{T}(θ::T, dest_x::T, dest_y::T, λ::T, ϕ=T(0.125π))
    left = Point{2,T}(dest_x - λ*cos(θ+ϕ), dest_y - λ*sin(θ+ϕ))
    pointy = Point{2,T}(dest_x, dest_y)
    right = Point{2,T}(dest_x - λ*cos(θ-ϕ), dest_y - λ*sin(θ-ϕ))
    LineSegment{Point{2,T}}[
        (left, pointy)
        (pointy, left)
    ]
end

function layout_tree{T<:Integer, R<:SimpleRectangle}(
        adj_list::AdjList{T},
        labels::Vector{R};
        ordering    = :optimal,
        coord       = :optimal,
        xsep        = 50,
        ysep        = 120,
        scale       = 0.05,
        labelpad    = 1.2,
    )
    # Calculate the original number of vertices
    n = length(adj_list)
    TV = fieldtype(R, :x)

    # 2     Layering
    # 2.1   Assign a layer to each vertex
    layers = GraphLayout._layer_assmt_longestpath(adj_list)
    num_layers = maximum(layers)
    # 2.2  Create dummy vertices for long edges
    adj_list, layers = GraphLayout._layer_assmt_dummy(adj_list, layers)
    orig_n, n = n, length(adj_list)

    # 3     Vertex ordering [to reduce crossings]
    # 3.1   Build initial permutation vectors
    layer_verts = [L => Int[] for L in 1:num_layers]
    for i in 1:n
        push!(layer_verts[layers[i]], i)
    end
    # 3.2  Reorder permutations to reduce crossings
    if ordering == :barycentric
        layer_verts = GraphLayout._ordering_barycentric(adj_list, layers, layer_verts)
    elseif ordering == :optimal
        layer_verts = GraphLayout._ordering_ip(adj_list, layers, layer_verts)
    end

    # 4     Vertex coordinates [to straighten edges]
    # 4.1   Place y coordinates in layers
    locs_y = zeros(TV, n)
    for L in 1:num_layers
        for (x,v) in enumerate(layer_verts[L])
            locs_y[v] = (L-1)*ysep
        end
    end
    # 4.2   Get widths of each label, if there are any
    widths_  = ones(TV, n); widths_[orig_n+1:n]  = 0
    heights = ones(TV, n); heights[orig_n+1:n] = 0
    # Note that we will convert these sizes into "absolute" units
    # and then work in these same units throughout. The font size used
    # here is just arbitrary, and unchanging. This hack arises because it
    # is meaningless to ask for the size of the font in "relative" units
    # but we don't want to collapse to absolute units until the end.
    if length(labels) == orig_n
        @inbounds for (i, rect) in enumerate(labels)
            w, h = Vec{2, TV}(widths(rect))
            widths_[i]  = w
            heights[i] = h
        end
    end
    locs_x = convert(Vector{TV}, GraphLayout._coord_ip(adj_list, layers, layer_verts, orig_n, widths_, xsep))
    # 4.3   Summarize vertex info
    max_x, max_y = maximum(locs_x), maximum(locs_y)
    max_w, max_h = maximum(widths_), maximum(heights)

    # 5     Draw the tree
    # 5.1   Create the vertices
    positioned_rectangles = [_tree_rectangle(locs_x[i], locs_y[i], labels[i]) for i in 1:orig_n]
    # 5.2   Create the arrows
    lines = LineSegment{Point{2,TV}}[]
    for L in 1:num_layers, i in layer_verts[L], j in adj_list[i]
        node_lines = _arrow_tree(
            Point{2,TV}(locs_x[i], locs_y[i]), i<=orig_n ? max_h : TV(0),
            Point{2,TV}(locs_x[j], locs_y[j]), j<=orig_n ? max_h : TV(0)
        )
        append!(lines, node_lines)
    end

    return positioned_rectangles, lines
end

"""
Takes a signal and creates a dictionary of unique nodes and
their neighbours from the signal graph.
"""
function to_adjency_dict(s::Signal, adjecency_dict = Dict{Signal, Vector{Any}}(), neighbours = [])
    adjecency_dict[s] = neighbours
    for action in s.actions
        neighour_s = action.recipient.value
        if !haskey(adjecency_dict, neighour_s)
            push!(neighbours, neighour_s)
            to_adjency_dict(neighour_s, adjecency_dict)
        end
    end
    for parent_s in s.parents
        if !haskey(adjecency_dict, parent_s)
            to_adjency_dict(parent_s, adjecency_dict, Any[s])
        end
    end
    adjecency_dict
end

"""
Takes the signal `s` and turns it into a adjecency list plus a list of the unique
nodes in the correct order for the adjecency list.
"""
function to_adjency_list(s::Signal)
    dict = to_adjency_dict(s)
    adjecency_list = Vector{Int}[]
    labels = Signal[]
    pos_in_adjecency = Dict{Signal, Int}()
    # Initialize neighbours in adjecency_list
    for (signal, v) in dict
        push!(adjecency_list, Int[])
        push!(labels, signal)
        pos_in_adjecency[signal] = length(adjecency_list)
    end
    # fill neighbours
    for (signal, neighbours) in dict
        neighbourlist = adjecency_list[pos_in_adjecency[signal]]
        for neighbour in neighbours
            push!(neighbourlist, pos_in_adjecency[neighbour])
        end
    end
    adjecency_list, labels
end


text:
include("physical_constants.jl")
include("tools.jl")
include("gap_functions.jl")
include("kernel.jl")

"""Return Z at several temperatures or frequencies.

Parameters
----------
input_vector : list-like
    List of either reduced temperatures (T/Tc) or reduced frequencies
    (h*f/delta0), where Tc is the superconducting critical temperature, h is
    Planck"s constant and delta0 is the zero-temperature superconducting
    energy gap.

tc : float
    Critical temperature in Kelvin.

vf : float
    Fermi velocity in m/s.

london0 : float
    London penetration depth at zero temperature in m.

axis : "string" (optional)
    Acceptable values are ``"temperature" or "frequency" or "mean free
    path"``. Specifies what the units of the ``input_vector`` parameter are.
    Default is ``"temperature"``.

Keyword Arguments
-----------------
mfp, tr, fr : float (required)
    If ``axis == "temperature"`` must specify reduced frequency value fr =
    h*f/delta0 and mean-free-path in meters. If ``axis == "frequency"`` must specify a reduced
    temperature value tr = T/Tc and mean-free-path in meters. If ``axis == "mean free path"`` must
    specify both fr and tr.

bcs : float (optional)
    The constant that gives the zero-temperature superconducting energy gap
    delta0 according to the equation delta0 = bcs*kB*Tc, where kB is
    Boltzmann"s constant and Tc is the superconducting critical temperature.
    Default value is the Bardeen-Cooper-Schrieffer value of 1.76.

gap : python function or string (optional)
    Python function that gives a value for the reduced superconducting
    energy gap deltar(T) = delta(T)/delta0, where delta is the
    superconducting energy gap and delta0 = delta(T=0). Function signature
    is float(float) with return value between zero and one. Default is
    tabulated values from Muhlschlegel (1959) via the deltar_bcs function.
    Optionally, one may pass the string ``"cos"`` to use the built-in cosine
    approximation of the gap.

output_depths : bool (optional)
    Sets the output units. False returns complex impedance in Ohms. True
    converts the complex impedance to a skin-depth (real part) and a
    superconducting penetration depth (imaginary part), both in meters.
    Default is False.

boundary : string (optional)
    Options are ``"diffuse"``/``"d"`` or ``"specular"``/``"s"``. Determines
    whether the impedance calculation assumes diffuse or specular scattering
    at the boundaries of the superconductor. Default is ``"diffuse"``.

verbose : int
    Whether to print out some debugging information. Very much a work in
    progress. Default is 0.

    * 0: Don"t print out any debugging information
    * 1: Print out minimal debugging information
    * 2: Print out full output from quad routines

Returns
-------
zs : numpy array
    The complex impedance calculated at each value of input_vector.

Note
----
fr = 0 may return slightly wrong number for specular boundary conditions.
Not sure why."""
function get_Zvec(input_vector, tc, vf, london0, axis="temperature"; kwargs...)

    allowed_axes = ["temperature", "t", "frequency", "f", "mean free path", "mfp", "l"]

    @assert axis in allowed_axes "Invalid axis."

    kwargs = Dict(kwargs)

    #This is the value that specifies delta0 = bcs*kB*Tc
    bcs = pop!(kwargs, :bcs, 1.76)
    delta0 = bcs*kb_eV*tc

    #Kwargs to pass to the cmplx_impedance function
    zs_kwargs = Dict()

    #Allow for passing in a custom gap function
    gap = pop!(kwargs, :gap, nothing)
    if gap != nothing
        if gap == "cos"
            zs_kwargs[:gap] = deltar_cos
        else
            zs_kwargs[:gap] = gap
        end
    end

    verbose = pop!(kwargs, :verbose, 0)

    zs_kwargs[:verbose] = verbose
    #Optionally can output penetration/skin depth in meters instead of Ohms
    output_depths = pop!(kwargs, :output_depths, false)
    zs_kwargs[:output_depths] = output_depths

    #See if specular reflection is requested
    boundary = pop!(kwargs, :boundary, "diffuse")
    @assert boundary in ["diffuse", "d", "specular", "s"] "Invalid boundary type."
    zs_kwargs[:boundary] = boundary

    if axis in ["mfp", "mean free path", "l"]
        @assert haskey(kwargs, :fr) "Must supply reduced frequency"
        fr = kwargs[:fr]

        @assert haskey(kwargs, :tr) "Must supply reduced temperature"
        tr = kwargs[:tr]

        mfps = input_vector

        zs = Array{Complex128}(length(input_vector))

        for (ix, mfp) in enumerate(mfps)
            #Convert physical data to params
            params_dict = init_from_physical_data(tc, vf, london0, mfp, bcs)

            #Add the parameters from physical data into the kwargs dict
            merge!(zs_kwargs, params_dict)

            xk = pop!(zs_kwargs, :xk, nothing)
            xm = pop!(zs_kwargs, :xm, nothing)
            vf = pop!(zs_kwargs, :vf, nothing)

            #Calculate the next impedance
            zs[ix] = cmplx_impedance(tr, fr, tc, xk, xm, vf; zs_kwargs...)
        end

    else
        @assert haskey(kwargs, :mfp) "Must supply mean free path"
        mfp = kwargs[:mfp]

        #Convert physical data to params
        params_dict = init_from_physical_data(tc, vf, london0, mfp, bcs)

        #Add the parameters from physical data into the kwargs dict
        merge!(zs_kwargs, params_dict)

        xk = pop!(zs_kwargs, :xk, nothing)
        xm = pop!(zs_kwargs, :xm, nothing)
        vf = pop!(zs_kwargs, :vf, nothing)

        if axis in ["temperature", "t"]
            @assert haskey(kwargs, :fr) "Must supply reduced frequency"
            fr = kwargs[:fr]

            trs = input_vector
            zs = Array{Complex128}(length(input_vector))

            for (ix, tr) in enumerate(trs)
              zs[ix] = cmplx_impedance(tr, fr, tc, xk, xm, vf; zs_kwargs...)
            end

        elseif axis in ["frequency", "f"]
            @assert haskey(kwargs, :tr) "Must supply reduced temperature"
            tr = kwargs[:tr]

            frs = input_vector
            zs = Array{Complex128}(length(input_vector))

            for (ix, fr) in enumerate(frs)
              zs[ix] = cmplx_impedance(tr, fr, tc, xk, xm, vf; zs_kwargs...)
            end
        end
    end

    return zs
end

"""Calculate the complex surface impedance (Z) of a superconductor at a
given temperature and frequency.

Parameters
----------
tr : float
    Reduced temperature tr = T/Tc, where T is temperature and Tc is the
    superconducting critical temperature.

fr : float
    Reduced frequency fr = h*f/delta0 where h is Planck"s constant, f is the
    frequency in Hz, and delta0 is the zero-temperature superconducting
    energy gap.

tc : float
    Superconducting critical temperature in Kelvin.

xk : float
    BCS coherence length (ksi0) divided by mean free path (mfp). ksi0 =
    hbar*vf/(pi*delta0), where hbar is Planck"s constant divided by 2*pi, vf
    is the Fermi velocity, and delta0 is the zero-temprature superconducting
    energy gap.

xm : float
    Mean free path (mfp) divided by the zero-temperature London penetration
    depth (london0).

vf : float
    Fermi velocity in m/s.

Keyword Arguments
-----------------
bcs : float (optional)
    The constant that gives the zero-temperature superconducting energy gap
    delta0 according to the equation delta0 = bcs*kB*Tc, where kB is
    Boltzmann"s constant and Tc is the superconducting critical temperature.
    Default value is the Bardeen-Cooper-Schrieffer value of 1.76.

boundary : string (optional)
    Options are ``"diffuse"``/``"d"`` or ``"specular"``/``"s"``. Determines
    whether the impedance calculation assumes diffuse or specular scattering
    at the boundaries of the superconductor. Default is ``"diffuse"``.

gap : python function or string (optional)
    Python function that gives a value for the reduced superconducting
    energy gap deltar(T) = delta(T)/delta0, where delta is the
    superconducting energy gap and delta0 = delta(T=0). Function signature
    is float(float) with return value between zero and one. Default is
    tabulated values from Muhlschlegel (1959) via the deltar_bcs function.
    Optionally, one may pass the string ``"cos"`` to use the built-in cosine
    approximation of the gap.

output_depths : bool (optional)
    Sets the output units. False returns complex impedance in Ohms. True
    converts the complex impedance to a skin-depth (real part) and a
    superconducting penetration depth (imaginary part), both in meters.
    Default is False.

verbose : int
    Whether to print out some debugging information. Very much a work in
    progress. Default is 0.

    * 0: Don"t print out any debugging information
    * 1: Print out minimal debugging information
    * 2: Print out full output from quad routines

Returns
-------
Z: The complex surface impedance in Ohms.

Note
----
fr = 0 may return slightly wrong number for specular boundary conditions.
Not sure why."""
function cmplx_impedance(tr, fr, tc, xk, xm, vf; kwargs...)

    kwargs = Dict(kwargs)
    #Optionally can output penetration/skin depth in meters instead of Ohms
    output_depths = pop!(kwargs, :output_depths, false)

    units = "Ohms"
    if output_depths
        units = "meters"
    end

    verbose = pop!(kwargs, :verbose, false)
    boundary = pop!(kwargs, :boundary, "diffuse")

    gap = pop!(kwargs, :gap, nothing)
    if gap == nothing
        gap = deltar_bcs
    elseif gap == "cos"
        gap = deltar_cos
    end

    dr = gap(tr)
    bcs = pop!(kwargs, :bcs, 1.76)
    delta0 = bcs*kb_eV*tc

    #Calculate the prefactor. Mostly this doesn"t matter since we care about ratios.
    if output_depths
        #Units here are meters
        prefactor = vf*hbar_eV/(xk*delta0)
    else
        #Units here are Ohms
        prefactor = fr*mu_0*vf/xk
    end

    k(x) = cmplx_kernel(tr, fr, x, xk, xm, dr, bcs, verbose)

    if (boundary == "diffuse") || (boundary == "d")
        #Need two integrals to help with singularity removal

        invZint_a(x) = log(x^2+k(x))-2*log(x)
        invZint_b(x) = log(x^2+k(x))

        invZ_a, invZerr = quadgk(invZint_a, 1, Inf)
        invZ_b, invZerr = quadgk(invZint_b, 0, 1)
        #quadgk(x->-2*log(x), 0, 1) = 2

        invZ = invZ_a + invZ_b + 2

        Z = 1.0/invZ

    elseif (boundary == "specular") || (boundary == "s")
        zint(x) = 1.0/(x^2+k(x))

        Z, Zerr = quadgk(zint, -Inf, Inf)

        Z = Z/π^2
    end

    Z *= prefactor*1im

    return Z
end


text:
# unit test he SVRG solver
include("../src/DMD_Util.jl")
include("../src/DMD_Type.jl")
include("../src/DMD_Func.jl")
include("../src/DMD_SVRG.jl")

# Generate DMD Synthetic Data
#------------------------------------------------------------------------------
# dimensions
m = 1000;    # temporal dimension
n = 100;    # spatial dimension
k = 8;      # number of modes
T = Float64;
sigma = 1e-2;
# generate data
X, t, at, Bt = simDMD(m, n, k, T; seed=123456);
Xdat = X + sigma*Random.randn!(copy(X))
# specify the loss
function lossFunc(r)
	return sum(abs2, r)
end
#
function lossGrad(r)
	conj!(r)
end
#
# Create object
#------------------------------------------------------------------------------
params = DMDParams(k, Xdat, t, lossFunc, lossGrad);
#
# Initial Guess
#------------------------------------------------------------------------------
# create initial guess using trap. rule with exact dmd
a0, B0 = dmdexactestimate(m, n, k, Xdat, t);
#
# Create Object
#------------------------------------------------------------------------------
params = DMDParams(k, Xdat, t, lossFunc, lossGrad);
copyto!(params.a, a0);
 copyto!(params.B, B0);
#Random.randn!(params.ar);
#
# Apply Solver
#------------------------------------------------------------------------------
τ = 20;
η = 1e-1;
itm = 1000;
tol = 1e-9;
ptf = 100;

function prox_stab(ar)
    # project for numerical stability
    k = length(ar) >> 1;
    @inbounds for i = 1:k
        ar[i<<1-1] = min(0.0, ar[i<<1-1]);
    end
end

opts = DMD_SVRG_Options(τ, η, params, itm=itm, tol=tol, ptf=ptf,
                        prox=prox_stab);

obj_his, err_his = solveDMD_withSVRG(params, opts);

println(besterrperm(params.a,at))

println(aBFunc(params)/lossFunc(X))

params2 = DMDParams(k,Xdat,t,lossFunc,lossGrad);
copyto!(params2.a,at)
println(aBFunc(params2)/lossFunc(X))
copyto!(params2.a,a0)
println(aBFunc(params2)/lossFunc(X))


text:
module EventShapes
using StaticArrays
using LinearAlgebra

include("thrust.jl")
include("foxWolfram.jl")

end # module


text:
"""Runge Kutta time stepping scheme diagnostic cariables collected in a struct."""
@with_kw struct RungeKuttaVars{T<:Real}

    # to be specified
    nx::Int
    ny::Int
    bc::String
    halo::Int
    haloη::Int

    nux::Int = if (bc == "periodic") nx else nx-1 end   # u-grid in x-direction
    nuy::Int = ny                                       # u-grid in y-direction
    nvx::Int = nx                                       # v-grid in x-direction
    nvy::Int = ny-1                                     # v-grid in y-direction
    nqx::Int = if (bc == "periodic") nx else nx+1 end   # q-grid in x-direction
    nqy::Int = ny+1                                     # q-grid in y-direction

    # EDGE POINT (1 = yes, 0 = no)
    ep::Int = if bc == "periodic" 1 else 0 end  # is there a u-point on the left edge?

    u0::Array{T,2} = zeros(T,nux+2*halo,nuy+2*halo)     # u-velocities for RK updates
    u1::Array{T,2} = zeros(T,nux+2*halo,nuy+2*halo)
    v0::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo)     # v-velocities for RK updates
    v1::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo)
    η0::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)     # sea surface height for RK updates
    η1::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)
end

"""Generator function for RungeKutta VarCollection."""
function RungeKuttaVars{T}(G::Grid) where {T<:Real}

    @unpack nx,ny,bc = G
    @unpack halo,haloη = G

    return RungeKuttaVars{T}(nx=nx,ny=ny,bc=bc,halo=halo,haloη=haloη)
end

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

"""Tendencies collected in a struct."""
@with_kw struct TendencyVars{T<:Real}

    # to be specified
    nx::Int
    ny::Int
    bc::String
    halo::Int
    haloη::Int

    nux::Int = if (bc == "periodic") nx else nx-1 end   # u-grid in x-direction
    nuy::Int = ny                                       # u-grid in y-direction
    nvx::Int = nx                                       # v-grid in x-direction
    nvy::Int = ny-1                                     # v-grid in y-direction
    nqx::Int = if (bc == "periodic") nx else nx+1 end   # q-grid in x-direction
    nqy::Int = ny+1                                     # q-grid in y-direction

    # EDGE POINT (1 = yes, 0 = no)
    ep::Int = if bc == "periodic" 1 else 0 end  # is there a u-point on the left edge?

    du::Array{T,2} = zeros(T,nux+2*halo,nuy+2*halo)     # tendency of u without time step
    dv::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo)     # tendency of v without time step
    dη::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)     # tendency of η without time step

    # sum of tendencies (incl time step) over all sub-steps
    du_sum::Array{T,2} = zeros(T,nux+2*halo,nuy+2*halo)
    dv_sum::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo)
    dη_sum::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)

    # compensation for tendencies (variant of Kahan summation)
    du_comp::Array{T,2} = zeros(T,nux+2*halo,nuy+2*halo)
    dv_comp::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo)
    dη_comp::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)
end

"""Generator function for Tendencies VarCollection."""
function TendencyVars{T}(G::Grid) where {T<:Real}

    @unpack nx,ny,bc = G
    @unpack halo,haloη = G

    return TendencyVars{T}(nx=nx,ny=ny,bc=bc,halo=halo,haloη=haloη)
end

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

"""VolumeFluxes collected in a struct."""
@with_kw struct VolumeFluxVars{T<:Real}

    # to be specified
    nx::Int
    ny::Int
    bc::String
    halo::Int
    haloη::Int

    nux::Int = if (bc == "periodic") nx else nx-1 end   # u-grid in x-direction
    nuy::Int = ny                                       # u-grid in y-direction
    nvx::Int = nx                                       # v-grid in x-direction
    nvy::Int = ny-1                                     # v-grid in y-direction
    nqx::Int = if (bc == "periodic") nx else nx+1 end   # q-grid in x-direction
    nqy::Int = ny+1                                     # q-grid in y-direction

    # EDGE POINT (1 = yes, 0 = no)
    ep::Int = if bc == "periodic" 1 else 0 end      # is there a u-point on the left edge?

    h::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)         # layer thickness
    h_u::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη)     # layer thickness on u-grid
    U::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη)       # U=uh volume flux

    h_v::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη-1)     # layer thickness on v-grid
    V::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη-1)       # V=vh volume flux

    dUdx::Array{T,2} = zeros(T,nx+2*haloη-2,ny+2*haloη)    # gradients thereof
    dVdy::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη-2)
end

"""Generator function for VolumeFluxes VarCollection."""
function VolumeFluxVars{T}(G::Grid) where {T<:Real}

    @unpack nx,ny,bc = G
    @unpack halo,haloη = G

    return VolumeFluxVars{T}(nx=nx,ny=ny,bc=bc,halo=halo,haloη=haloη)
end

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

"""Vorticity variables collected in a struct."""
@with_kw struct VorticityVars{T<:Real}

    # to be specified
    nx::Int
    ny::Int
    bc::String
    halo::Int
    haloη::Int

    nux::Int = if (bc == "periodic") nx else nx-1 end   # u-grid in x-direction
    nuy::Int = ny                                       # u-grid in y-direction
    nvx::Int = nx                                       # v-grid in x-direction
    nvy::Int = ny-1                                     # v-grid in y-direction
    nqx::Int = if (bc == "periodic") nx else nx+1 end   # q-grid in x-direction
    nqy::Int = ny+1                                     # q-grid in y-direction

    # EDGE POINT (1 = yes, 0 = no)
    ep::Int = if bc == "periodic" 1 else 0 end      # is there a u-point on the left edge?

    h_q::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη-1)  # layer thickness h interpolated on q-grid
    q::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη-1)    # potential vorticity

    q_v::Array{T,2} = zeros(T,nx+2*haloη-2,ny+2*haloη-1)  # q interpolated on v-grid
    U_v::Array{T,2} = zeros(T,nx+2*haloη-2,ny+2*haloη-1)  # mass flux U=uh on v-grid

    q_u::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη-2)  # q interpolated on u-grid
    V_u::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη-2)  # mass flux V=vh on v-grid

    qhu::Array{T,2} = zeros(T,nvx,nvy)            # potential vorticity advection term u-component
    qhv::Array{T,2} = zeros(T,nux,nuy)            # potential vorticity advection term v-component

    u_v::Array{T,2} = zeros(T,nux+2*halo-1,nuy+2*halo-1)  # u-velocity on v-grid
    v_u::Array{T,2} = zeros(T,nvx+2*halo-1,nvy+2*halo-1)  # v-velocity on u-grid

    dudx::Array{T,2} = zeros(T,nux+2*halo-1,nuy+2*halo)   # ∂u/∂x
    dudy::Array{T,2} = zeros(T,nux+2*halo,nuy+2*halo-1)   # ∂u/∂y

    dvdx::Array{T,2} = zeros(T,nvx+2*halo-1,nvy+2*halo)   # ∂v/∂x
    dvdy::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo-1)   # ∂v/∂y
end

"""Generator function for Vorticity VarCollection."""
function VorticityVars{T}(G::Grid) where {T<:Real}

    @unpack nx,ny,bc = G
    @unpack halo,haloη = G

    return VorticityVars{T}(nx=nx,ny=ny,bc=bc,halo=halo,haloη=haloη)
end

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

"""Bernoulli variables collected in a struct."""
@with_kw struct BernoulliVars{T<:Real}

    # to be specified
    nx::Int
    ny::Int
    bc::String
    halo::Int
    haloη::Int

    nux::Int = if (bc == "periodic") nx else nx-1 end   # u-grid in x-direction
    nuy::Int = ny                                       # u-grid in y-direction
    nvx::Int = nx                                       # v-grid in x-direction
    nvy::Int = ny-1                                     # v-grid in y-direction
    nqx::Int = if (bc == "periodic") nx else nx+1 end   # q-grid in x-direction
    nqy::Int = ny+1                                     # q-grid in y-direction

    # EDGE POINT (1 = yes, 0 = no)
    ep::Int = if bc == "periodic" 1 else 0 end      # is there a u-point on the left edge?

    u²::Array{T,2} = zeros(T,nux+2*halo,nuy+2*halo)         # u-velocity squared
    v²::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo)         # v-velocity squared

    KEu::Array{T,2} = zeros(T,nux+2*halo-1,nuy+2*halo)      # u-velocity squared on T-grid
    KEv::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo-1)      # v-velocity squared on T-grid

    p::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)          # Bernoulli potential
    dpdx::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη)     # ∂p/∂x
    dpdy::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη-1)     # ∂p/∂y
end

"""Generator function for Bernoulli VarCollection."""
function BernoulliVars{T}(G::Grid) where {T<:Real}

    @unpack nx,ny,bc = G
    @unpack halo,haloη = G

    return BernoulliVars{T}(nx=nx,ny=ny,bc=bc,halo=halo,haloη=haloη)
end

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

"""Bottomdrag variables collected in a struct."""
@with_kw struct BottomdragVars{T<:Real}

    # to be specified
    nx::Int
    ny::Int
    bc::String
    halo::Int
    haloη::Int

    nux::Int = if (bc == "periodic") nx else nx-1 end   # u-grid in x-direction
    nuy::Int = ny                                       # u-grid in y-direction
    nvx::Int = nx                                       # v-grid in x-direction
    nvy::Int = ny-1                                     # v-grid in y-direction
    nqx::Int = if (bc == "periodic") nx else nx+1 end   # q-grid in x-direction
    nqy::Int = ny+1                                     # q-grid in y-direction

    # EDGE POINT (1 = yes, 0 = no)
    ep::Int = if bc == "periodic" 1 else 0 end      # is there a u-point on the left edge?

    sqrtKE::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)       # sqrt of kinetic energy
    sqrtKE_u::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη)   # interpolated on u-grid
    sqrtKE_v::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη-1)   # interpolated on v-grid

    Bu::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη)         # bottom friction term u-component
    Bv::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη-1)         # bottom friction term v-component
end

"""Generator function for Bottomdrag VarCollection."""
function BottomdragVars{T}(G::Grid) where {T<:Real}

    @unpack nx,ny,bc = G
    @unpack halo,haloη = G

    return BottomdragVars{T}(nx=nx,ny=ny,bc=bc,halo=halo,haloη=haloη)
end

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

"""ArakawaHsu variables collected in a struct."""
@with_kw struct ArakawaHsuVars{T<:Real}

    # to be specified
    nx::Int
    ny::Int
    bc::String
    halo::Int
    haloη::Int

    nux::Int = if (bc == "periodic") nx else nx-1 end   # u-grid in x-direction
    nuy::Int = ny                                       # u-grid in y-direction
    nvx::Int = nx                                       # v-grid in x-direction
    nvy::Int = ny-1                                     # v-grid in y-direction
    nqx::Int = if (bc == "periodic") nx else nx+1 end   # q-grid in x-direction
    nqy::Int = ny+1                                     # q-grid in y-direction

    # EDGE POINT (1 = yes, 0 = no)
    ep::Int = if bc == "periodic" 1 else 0 end      # is there a u-point on the left edge?

    # Linear combination of potential vorticity
    qα::Array{T,2} = zeros(T,nx+2*haloη-2,ny+2*haloη-2)
    qβ::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη-2)
    qγ::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη-2)
    qδ::Array{T,2} = zeros(T,nx+2*haloη-2,ny+2*haloη-2)
end

"""Generator function for ArakawaHsu VarCollection."""
function ArakawaHsuVars{T}(G::Grid) where {T<:Real}

    @unpack nx,ny,bc = G
    @unpack halo,haloη = G

    return ArakawaHsuVars{T}(nx=nx,ny=ny,bc=bc,halo=halo,haloη=haloη)
end

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

"""Laplace variables collected in a struct."""
@with_kw struct LaplaceVars{T<:Real}

    # to be specified
    nx::Int
    ny::Int
    bc::String
    halo::Int
    haloη::Int

    nux::Int = if (bc == "periodic") nx else nx-1 end   # u-grid in x-direction
    nuy::Int = ny                                       # u-grid in y-direction
    nvx::Int = nx                                       # v-grid in x-direction
    nvy::Int = ny-1                                     # v-grid in y-direction
    nqx::Int = if (bc == "periodic") nx else nx+1 end   # q-grid in x-direction
    nqy::Int = ny+1                                     # q-grid in y-direction

    # EDGE POINT (1 = yes, 0 = no)
    ep::Int = if bc == "periodic" 1 else 0 end      # is there a u-point on the left edge?

    Lu::Array{T,2} = zeros(T,nux+2*halo-2,nuy+2*halo-2)         # ∇²u
    Lv::Array{T,2} = zeros(T,nvx+2*halo-2,nvy+2*halo-2)         # ∇²v

    # Derivatives of Lu,Lv
    dLudx::Array{T,2} = zeros(T,nux+2*halo-3,nuy+2*halo-2)
    dLudy::Array{T,2} = zeros(T,nux+2*halo-2,nuy+2*halo-3)
    dLvdx::Array{T,2} = zeros(T,nvx+2*halo-3,nvy+2*halo-2)
    dLvdy::Array{T,2} = zeros(T,nvx+2*halo-2,nvy+2*halo-3)
end

"""Generator function for Laplace VarCollection."""
function LaplaceVars{T}(G::Grid) where {T<:Real}

    @unpack nx,ny,bc = G
    @unpack halo,haloη = G

    return LaplaceVars{T}(nx=nx,ny=ny,bc=bc,halo=halo,haloη=haloη)
end

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

"""Smagorinsky variables collected in a struct."""
@with_kw struct SmagorinskyVars{T<:Real}

    # to be specified
    nx::Int
    ny::Int
    bc::String
    halo::Int
    haloη::Int

    nux::Int = if (bc == "periodic") nx else nx-1 end   # u-grid in x-direction
    nuy::Int = ny                                       # u-grid in y-direction
    nvx::Int = nx                                       # v-grid in x-direction
    nvy::Int = ny-1                                     # v-grid in y-direction
    nqx::Int = if (bc == "periodic") nx else nx+1 end   # q-grid in x-direction
    nqy::Int = ny+1                                     # q-grid in y-direction

    # EDGE POINT (1 = yes, 0 = no)
    ep::Int = if bc == "periodic" 1 else 0 end      # is there a u-point on the left edge?

    DT::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)       # Tension squared (on the T-grid)
    DS::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)       # Shearing strain squared (on the T-grid)
    νSmag::Array{T,2} = zeros(T,nx+2*haloη,ny+2*haloη)    # Viscosity coefficient

    # Tension squared on the q-grid
    DS_q::Array{T,2} = zeros(T,nvx+2*halo-1,nvy+2*halo)

    # Smagorinsky viscosity coefficient on the q-grid
    νSmag_q::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη-1)

    # Entries of the Smagorinsky viscous tensor
    S12::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη-1)
    S21::Array{T,2} = zeros(T,nx+2*haloη-1,ny+2*haloη-1)

    S11::Array{T,2} = zeros(T,nux+2*halo-3,nuy+2*halo-2)
    S22::Array{T,2} = zeros(T,nvx+2*halo-2,nvy+2*halo-3)

    # u- and v-components 1 and 2 of the biharmonic diffusion tendencies
    LLu1::Array{T,2} = zeros(T,nux+2*halo-4,nuy+2*halo-2)
    LLu2::Array{T,2} = zeros(T,nx+1,ny)

    LLv1::Array{T,2} = zeros(T,nx,ny+1)
    LLv2::Array{T,2} = zeros(T,nvx+2*halo-2,nvy+2*halo-4)
end

"""Generator function for Smagorinsky VarCollection."""
function SmagorinskyVars{T}(G::Grid) where {T<:Real}

    @unpack nx,ny,bc = G
    @unpack halo,haloη = G

    return SmagorinskyVars{T}(nx=nx,ny=ny,bc=bc,halo=halo,haloη=haloη)
end

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

"""SemiLagrange variables collected in a struct."""
@with_kw struct SemiLagrangeVars{T<:Real}

    # to be specified
    nx::Int
    ny::Int
    bc::String
    halo::Int
    haloη::Int
    halosstx::Int
    halossty::Int

    nux::Int = if (bc == "periodic") nx else nx-1 end   # u-grid in x-direction
    nuy::Int = ny                                       # u-grid in y-direction
    nvx::Int = nx                                       # v-grid in x-direction
    nvy::Int = ny-1                                     # v-grid in y-direction
    nqx::Int = if (bc == "periodic") nx else nx+1 end   # q-grid in x-direction
    nqy::Int = ny+1                                     # q-grid in y-direction

    # EDGE POINT (1 = yes, 0 = no)
    ep::Int = if bc == "periodic" 1 else 0 end      # is there a u-point on the left edge?

    xd::Array{T,2} = zeros(T,nx,ny)                         # departure points x-coord
    yd::Array{T,2} = zeros(T,nx,ny)                         # departure points y-coord

    um::Array{T,2} = zeros(T,nux+2*halo,nuy+2*halo)         # u-velocity temporal mid-point
    vm::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo)         # v-velocity temporal mid-point

    u_T::Array{T,2} = zeros(T,nux+2*halo-1,nuy+2*halo)      # u-velocity interpolated on T-grid
    um_T::Array{T,2} = zeros(T,nux+2*halo-1,nuy+2*halo)     # um interpolated on T-grid
    v_T::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo-1)      # v-velocity interpolated on T-grid
    vm_T::Array{T,2} = zeros(T,nvx+2*halo,nvy+2*halo-1)     # vm interpolated on T-grid

    uinterp::Array{T,2} = zeros(T,nx,ny)                    # u interpolated on mid-point xd,yd
    vinterp::Array{T,2} = zeros(T,nx,ny)                    # v interpolated on mid-point xd,yd

    ssti::Array{T,2} = zeros(T,nx+2*halosstx,ny+2*halossty) # sst interpolated on departure points
    sst_ref::Array{T,2} = zeros(T,nx+2*halosstx,ny+2*halossty) # sst initial conditions for relaxation

    # compensated summation
    dsst_comp::Array{T,2} = zeros(T,nx+2*halosstx,ny+2*halossty)
end

"""Generator function for SemiLagrange VarCollection."""
function SemiLagrangeVars{T}(G::Grid) where {T<:Real}

    @unpack nx,ny,bc = G
    @unpack halo,haloη = G
    @unpack halosstx,halossty = G

    return SemiLagrangeVars{T}(nx=nx,ny=ny,bc=bc,halo=halo,haloη=haloη,
                            halosstx=halosstx,halossty=halossty)
end

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

struct DiagnosticVars{T,Tprog}
    RungeKutta::RungeKuttaVars{Tprog}
    Tendencies::TendencyVars{Tprog}
    VolumeFluxes::VolumeFluxVars{T}
    Vorticity::VorticityVars{T}
    Bernoulli::BernoulliVars{T}
    Bottomdrag::BottomdragVars{T}
    ArakawaHsu::ArakawaHsuVars{T}
    Laplace::LaplaceVars{T}
    Smagorinsky::SmagorinskyVars{T}
    SemiLagrange::SemiLagrangeVars{T}
    PrognosticVarsRHS::PrognosticVars{T}        # low precision version
end

"""Preallocate the diagnostic variables and return them as matrices in structs."""
function preallocate(   ::Type{T},
                        ::Type{Tprog},
                        G::Grid) where {T<:Real,Tprog<:Real}

    RK = RungeKuttaVars{Tprog}(G)
    TD = TendencyVars{Tprog}(G)
    VF = VolumeFluxVars{T}(G)
    VT = VorticityVars{T}(G)
    BN = BernoulliVars{T}(G)
    BD = BottomdragVars{T}(G)
    AH = ArakawaHsuVars{T}(G)
    LP = LaplaceVars{T}(G)
    SM = SmagorinskyVars{T}(G)
    SL = SemiLagrangeVars{T}(G)
    PV = PrognosticVars{T}(G)

    return DiagnosticVars{T,Tprog}(RK,TD,VF,VT,BN,BD,AH,LP,SM,SL,PV)
end


text:
using Images
using OCReract
using Test

pkg_path = abspath(joinpath(dirname(pathof(OCReract)), ".."))

TEST_ITEMS = Dict(
    "simple" => Dict(
        "path" => "$(pkg_path)/test/files/testocr.png",
        "text" => "This is a lot of 12 point text to test the\n"*
                  "ocr code and see if it works on all types\n"*
                  "of file format.\n\n"*
                  "The quick brown dog jumped over the\n"*
                  "lazy fox. The quick brown dog jumped\n"*
                  "over the lazy fox. The quick brown dog\n"*
                  "jumped over the lazy fox. The quick\n"*
                  "brown dog jumped over the lazy fox."
    ),
    "noisy" => Dict(
        "path" => "$(pkg_path)/test/files/noisy.png",
        "text" => "Noisy image\nto test\nOCReract.jl"
    )
)

function test_run_tesseract()
    tmp_path = tempname()
    mkdir(tmp_path)
    path_to_res = tmp_path*"/res.txt"
    res = run_tesseract(TEST_ITEMS["simple"]["path"], path_to_res)
    @test res == true

    result_txt = ""
    open(path_to_res, "r") do f
        result_txt = read(f, String)
    end
    @test strip(result_txt) == strip(TEST_ITEMS["simple"]["text"])

    # Again to test warning about already existing output file
    res = run_tesseract(TEST_ITEMS["simple"]["path"], path_to_res)
    @test res == true
end

function test_run_and_get_output()
    result_txt = run_tesseract(Images.load(TEST_ITEMS["simple"]["path"]))
    @test strip(TEST_ITEMS["simple"]["text"]) == strip(result_txt)
end

function test_run_with_kwargs()
    # Set an option that will make OCR bad for this image
    result_txt = run_tesseract(Images.load(TEST_ITEMS["noisy"]["path"]), tessedit_pageseg_mode=7)
    @test strip(TEST_ITEMS["noisy"]["text"]) != strip(result_txt)
end

function test_path_exist()
    res = run_tesseract("/tmp/not_existing_image.png", "/tmp/res")
    @test res == false
end

function test_bad_arguments()
    # Bad Language
    tmp_path = tempname()
    mkdir(tmp_path)
    path_to_res = tmp_path*"/res.txt"
    res = run_tesseract(TEST_ITEMS["simple"]["path"], path_to_res, lang="not-existing-lang")
    @test res == false

    # Despite bad parameters oem and psm, this will work due to default changes
    res = run_tesseract(TEST_ITEMS["simple"]["path"], path_to_res, oem=100, psm=100)
    @test res == true
end

@testset "RunTesseract" begin
    # Normal run
    test_run_tesseract()
    # Run and get
    test_run_and_get_output()
    # Test with non existent image
    test_path_exist()
    # Test with bad arguments
    test_bad_arguments()
    # Test with kwargs
    test_run_with_kwargs()
end


text:
# some model
mgrid = [range(0.0, stop=10.,step=0.05), range(0.0, stop=10.,step=0.05)];
model = GeoPhyInv.Models.Seismic_zeros(mgrid);
vp0=[2100.,2200.];vs0=[-1., -1.]; ρ0=[2100., 2300.]
GeoPhyInv.Models.adjust_bounds!(model, vp0,vs0,ρ0);

nznx=prod(length.(model.mgrid))

for parameterization in [[:χKI, :χρI, :null],[:χKI, :null, :null],[:χvp, :χρ, :null],[:χvp, :null, :null]]
	δx1=randn(count(parameterization .≠ :null)*nznx)
	δxout1=zero(δx1)

	J.Models.pert_reparameterize!(δxout1, δx1, model, parameterization)

	δxout2=randn(size(δxout1));
	δx2=zero(δx1)
	J.Models.pert_gradient_chainrule!(δx2, δxout2, model, parameterization)

	@test LinearAlgebra.dot(δxout2,δxout1) ≈ LinearAlgebra.dot(δx1,δx2)
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "hsakmt_roct"
version = v"3.8.0"

# Collection of sources required to build ROCT-Thunk-Interface
sources = [
    ArchiveSource("https://github.com/RadeonOpenCompute/ROCT-Thunk-Interface/archive/rocm-$(version).tar.gz",
                  "cd5440f31f592737b5d05448704bd01f91f73cfcab8a7829922e81332575cfe8"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd ${WORKSPACE}/srcdir/ROCT-Thunk-Interface*/

# fix for musl (but works on glibc too)
atomic_patch -p1 "${WORKSPACE}/srcdir/patches/0001-Build-correctly-on-musl.patch"

mkdir build && cd build
cmake -DCMAKE_BUILD_TYPE=Release \
      -DCMAKE_INSTALL_PREFIX=${prefix} \
      -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
      -DCMAKE_C_FLAGS="-Dstatic_assert=_Static_assert" \
      ..
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
# ROCT-Thunk-Interface only supports Linux, seemingly only 64bit
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "linux"; libc="musl"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct(["libhsakmt"], :libhsakmt),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("NUMA_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies,
               preferred_gcc_version=v"8")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil