text:
module jlgr

using Compat
import GR

const None = Union{}

macro _tuple(t)
    :( Tuple{$t} )
end

@compat typealias PlotArg Union{AbstractString, AbstractVector, AbstractMatrix, Function}

const gr3 = GR.gr3

const plot_kind = [:line, :scatter, :stem, :hist, :contour, :contourf, :wireframe, :surface, :plot3, :scatter3, :imshow, :isosurface, :polar, :trisurf]

const arg_fmt = [:xys, :xyac, :xyzc]

type PlotObject
    args
    kvs
end

function Figure(width=600, height=450)
    args = @_tuple(Any)
    kvs = Dict()
    kvs[:size] = (width, height)
    kvs[:ax] = false
    kvs[:subplot] = [0, 1, 0, 1]
    kvs[:clear] = true
    kvs[:update] = true
    PlotObject(args, kvs)
end

plt = Figure()

isrowvec(x::AbstractArray) = ndims(x) == 2 && size(x, 1) == 1 && size(x, 2) > 1

isvector(x::AbstractVector) = true
isvector(x::AbstractMatrix) = size(x, 1) == 1

function set_viewport(kind, subplot)
    mwidth, mheight, width, height = GR.inqdspsize()
    if haskey(plt.kvs, :figsize)
        w = 0.0254 *  width * plt.kvs[:figsize][1] / mwidth
        h = 0.0254 * height * plt.kvs[:figsize][2] / mheight
    else
        dpi = width / mwidth * 0.0254
        if dpi > 200
            w, h = [x * dpi / 100 for x in plt.kvs[:size]]
        else
            w, h = plt.kvs[:size]
        end
    end
    viewport = zeros(4)
    vp = float(subplot)
    if w > h
        ratio = float(h) / w
        msize = mwidth * w / width
        GR.setwsviewport(0, msize, 0, msize * ratio)
        GR.setwswindow(0, 1, 0, ratio)
        vp[3] *= ratio
        vp[4] *= ratio
    else
        ratio = float(w) / h
        msize = mheight * h / height
        GR.setwsviewport(0, msize * ratio, 0, msize)
        GR.setwswindow(0, ratio, 0, 1)
        vp[1] *= ratio
        vp[2] *= ratio
    end
    viewport[1] = vp[1] + 0.125 * (vp[2] - vp[1])
    viewport[2] = vp[1] + 0.925 * (vp[2] - vp[1])
    viewport[3] = vp[3] + 0.125 * (vp[4] - vp[3])
    viewport[4] = vp[3] + 0.925 * (vp[4] - vp[3])

    if w > h
        viewport[3] += (1 - (subplot[4] - subplot[3])^2) * 0.02
    end
    if kind in (:wireframe, :surface, :plot3, :scatter3, :trisurf)
        viewport[2] -= 0.0525
    end
    if kind in (:contour, :contourf, :surface, :trisurf)
        viewport[2] -= 0.1
    end
    GR.setviewport(viewport[1], viewport[2], viewport[3], viewport[4])

    plt.kvs[:viewport] = viewport
    plt.kvs[:vp] = vp
    plt.kvs[:ratio] = ratio

    if haskey(plt.kvs, :backgroundcolor)
        GR.savestate()
        GR.selntran(0)
        GR.setfillintstyle(GR.INTSTYLE_SOLID)
        GR.setfillcolorind(plt.kvs[:backgroundcolor])
        if w > h
          GR.fillrect(subplot[1], subplot[2],
                      ratio * subplot[3], ratio * subplot[4])
        else
          GR.fillrect(ratio * subplot[1], ratio * subplot[2],
                      subplot[3], subplot[4])
        end
        GR.selntran(1)
        GR.restorestate()
    end

    if kind == :polar
        xmin, xmax, ymin, ymax = viewport
        xcenter = 0.5 * (xmin + xmax)
        ycenter = 0.5 * (ymin + ymax)
        r = 0.5 * min(xmax - xmin, ymax - ymin)
        GR.setviewport(xcenter - r, xcenter + r, ycenter - r, ycenter + r)
    end
end

function fix_minmax(a, b)
    if a == b
        a -= a != 0 ? 0.1 * a : 0.1
        b += b != 0 ? 0.1 * b : 0.1
    end
    a, b
end

function minmax()
    xmin = ymin = zmin = typemax(Float64)
    xmax = ymax = zmax = typemin(Float64)
    for (x, y, z, c, spec) in plt.args
        xmin = min(minimum(x), xmin)
        xmax = max(maximum(x), xmax)
        ymin = min(minimum(y), ymin)
        ymax = max(maximum(y), ymax)
        if z != Void
            zmin = min(minimum(z), zmin)
            zmax = max(maximum(z), zmax)
        end
    end
    xmin, xmax = fix_minmax(xmin, xmax)
    ymin, ymax = fix_minmax(ymin, ymax)
    zmin, zmax = fix_minmax(zmin, zmax)
    if haskey(plt.kvs, :xlim)
        plt.kvs[:xrange] = plt.kvs[:xlim]
    else
        plt.kvs[:xrange] = xmin, xmax
    end
    if haskey(plt.kvs, :ylim)
        plt.kvs[:yrange] = plt.kvs[:ylim]
    else
        plt.kvs[:yrange] = ymin, ymax
    end
    if haskey(plt.kvs, :zlim)
        plt.kvs[:zrange] = plt.kvs[:zlim]
    else
        plt.kvs[:zrange] = zmin, zmax
    end
end

function set_window(kind)
    scale = 0
    if kind != :polar
        get(plt.kvs, :xlog, false) && (scale |= GR.OPTION_X_LOG)
        get(plt.kvs, :ylog, false) && (scale |= GR.OPTION_Y_LOG)
        get(plt.kvs, :zlog, false) && (scale |= GR.OPTION_Z_LOG)
        get(plt.kvs, :xflip, false) && (scale |= GR.OPTION_FLIP_X)
        get(plt.kvs, :yflip, false) && (scale |= GR.OPTION_FLIP_Y)
        get(plt.kvs, :zflip, false) && (scale |= GR.OPTION_FLIP_Z)
    end

    minmax()

    if kind in (:wireframe, :surface, :plot3, :scatter3, :polar, :trisurf)
        major_count = 2
    else
        major_count = 5
    end

    xmin, xmax = plt.kvs[:xrange]
    if scale & GR.OPTION_X_LOG == 0
        if !haskey(plt.kvs, :xlim)
            xmin, xmax = GR.adjustlimits(xmin, xmax)
        end
        majorx = major_count
        xtick = GR.tick(xmin, xmax) / majorx
    else
        xtick = majorx = 1
    end
    if scale & GR.OPTION_FLIP_X == 0
        xorg = (xmin, xmax)
    else
        xorg = (xmax, xmin)
    end
    plt.kvs[:xaxis] = xtick, xorg, majorx

    ymin, ymax = plt.kvs[:yrange]
    if kind in (:stem, :hist) && !haskey(plt.kvs, :ylim)
        ymin = 0
    end
    if scale & GR.OPTION_Y_LOG == 0
        if !haskey(plt.kvs, :ylim)
            ymin, ymax = GR.adjustlimits(ymin, ymax)
        end
        majory = major_count
        ytick = GR.tick(ymin, ymax) / majory
    else
        ytick = majory = 1
    end
    if scale & GR.OPTION_FLIP_Y == 0
        yorg = (ymin, ymax)
    else
        yorg = (ymax, ymin)
    end
    plt.kvs[:yaxis] = ytick, yorg, majory

    if kind in (:wireframe, :surface, :plot3, :scatter3, :trisurf)
        zmin, zmax = plt.kvs[:zrange]
        if scale & GR.OPTION_Z_LOG == 0
            if !haskey(plt.kvs, :zlim)
                zmin, zmax = GR.adjustlimits(zmin, zmax)
            end
            majorz = major_count
            ztick = GR.tick(zmin, zmax) / majorz
        else
            ztick = majorz = 1
        end
        if scale & GR.OPTION_FLIP_Z == 0
            zorg = (zmin, zmax)
        else
            zorg = (zmax, zmin)
        end
        plt.kvs[:zaxis] = ztick, zorg, majorz
    end

    plt.kvs[:window] = xmin, xmax, ymin, ymax
    if kind != :polar
        GR.setwindow(xmin, xmax, ymin, ymax)
    else
        GR.setwindow(-1, 1, -1, 1)
    end
    if kind in (:wireframe, :surface, :plot3, :scatter3, :trisurf)
        rotation = get(plt.kvs, :rotation, 40)
        tilt = get(plt.kvs, :tilt, 70)
        GR.setspace(zmin, zmax, rotation, tilt)
    end

    plt.kvs[:scale] = scale
    GR.setscale(scale)
end

function draw_axes(kind, pass=1)
    viewport = plt.kvs[:viewport]
    vp = plt.kvs[:vp]
    ratio = plt.kvs[:ratio]
    xtick, xorg, majorx = plt.kvs[:xaxis]
    ytick, yorg, majory = plt.kvs[:yaxis]

    GR.setlinecolorind(1)
    diag = sqrt((viewport[2] - viewport[1])^2 + (viewport[4] - viewport[3])^2)
    GR.setlinewidth(1)
    charheight = max(0.018 * diag, 0.012)
    GR.setcharheight(charheight)
    ticksize = 0.0075 * diag
    if kind in (:wireframe, :surface, :plot3, :scatter3, :trisurf)
        ztick, zorg, majorz = plt.kvs[:zaxis]
        if pass == 1
            GR.grid3d(xtick, 0, ztick, xorg[1], yorg[1], zorg[1], 2, 0, 2)
            GR.grid3d(0, ytick, 0, xorg[2], yorg[1], zorg[1], 0, 2, 0)
        else
            GR.axes3d(xtick, 0, ztick, xorg[1], yorg[1], zorg[1], majorx, 0, majorz, -ticksize)
            GR.axes3d(0, ytick, 0, xorg[2], yorg[1], zorg[1], 0, majory, 0, ticksize)
        end
    else
        GR.grid(xtick, ytick, 0, 0, majorx, majory)
        GR.axes(xtick, ytick, xorg[1], yorg[1], majorx, majory, ticksize)
        GR.axes(xtick, ytick, xorg[2], yorg[2], -majorx, -majory, -ticksize)
    end

    if haskey(plt.kvs, :title)
        GR.savestate()
        GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
        text(0.5 * (viewport[1] + viewport[2]), vp[4], plt.kvs[:title])
        GR.restorestate()
    end
    if kind in (:wireframe, :surface, :plot3, :scatter3, :trisurf)
        xlabel = get(plt.kvs, :xlabel, "")
        ylabel = get(plt.kvs, :ylabel, "")
        zlabel = get(plt.kvs, :zlabel, "")
        GR.titles3d(xlabel, ylabel, zlabel)
    else
        if haskey(plt.kvs, :xlabel)
            GR.savestate()
            GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_BOTTOM)
            text(0.5 * (viewport[1] + viewport[2]), vp[3] + 0.5 * charheight, plt.kvs[:xlabel])
            GR.restorestate()
        end
        if haskey(plt.kvs, :ylabel)
            GR.savestate()
            GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
            GR.setcharup(-1, 0)
            text(vp[1] + 0.5 * charheight, 0.5 * (viewport[3] + viewport[4]), plt.kvs[:ylabel])
            GR.restorestate()
        end
    end
end

function draw_polar_axes()
    viewport = plt.kvs[:viewport]
    diag = sqrt((viewport[2] - viewport[1])^2 + (viewport[4] - viewport[3])^2)
    charheight = max(0.018 * diag, 0.012)

    window = plt.kvs[:window]
    rmin, rmax = window[3], window[4]

    GR.savestate()
    GR.setcharheight(charheight)
    GR.setlinetype(GR.LINETYPE_SOLID)

    tick = 0.5 * GR.tick(rmin, rmax)
    n = round(Int, (rmax - rmin) / tick + 0.5)
    for i in 0:n
        r = float(i) / n
        if i % 2 == 0
            GR.setlinecolorind(88)
            if i > 0
                GR.drawarc(-r, r, -r, r, 0, 359)
            end
            GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
            x, y = GR.wctondc(0.05, r)
            GR.text(x, y, string(signif(rmin + i * tick, 12)))
        else
            GR.setlinecolorind(90)
            GR.drawarc(-r, r, -r, r, 0, 359)
        end
    end
    for alpha in 0:45:315
        a = alpha + 90
        sinf = sin(a * π / 180)
        cosf = cos(a * π / 180)
        GR.polyline([sinf, 0], [cosf, 0])
        GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_HALF)
        x, y = GR.wctondc(1.1 * sinf, 1.1 * cosf)
        GR.textext(x, y, string(alpha, "^o"))
    end
    GR.restorestate()
end

function inqtext(x, y, s)
    if length(s) >= 2 && s[1] == '$' && s[end] == '$'
        GR.inqtextext(x, y, s[2:end-1])
    elseif search(s, '\\') != 0 || search(s, '_') != 0 || search(s, '^') != 0
        GR.inqtextext(x, y, s)
    else
        GR.inqtext(x, y, s)
    end
end

function text(x, y, s)
    if length(s) >= 2 && s[1] == '$' && s[end] == '$'
        GR.mathtex(x, y, s[2:end-1])
    elseif search(s, '\\') != 0 || search(s, '_') != 0 || search(s, '^') != 0
        GR.textext(x, y, s)
    else
        GR.text(x, y, s)
    end
end

function draw_legend()
    viewport = plt.kvs[:viewport]
    num_labels = length(plt.kvs[:labels])
    GR.savestate()
    GR.selntran(0)
    GR.setscale(0)
    w = 0
    for label in plt.kvs[:labels]
        tbx, tby = GR.inqtextext(0, 0, label)
        w = max(w, tbx[3])
    end
    px = viewport[2] - 0.05 - w
    py = viewport[4] - 0.06
    GR.setfillintstyle(GR.INTSTYLE_SOLID)
    GR.setfillcolorind(0)
    GR.fillrect(px - 0.08, px + w + 0.02, py + 0.03, py - 0.03 * num_labels)
    GR.setlinetype(GR.LINETYPE_SOLID)
    GR.setlinecolorind(1)
    GR.setlinewidth(1)
    GR.drawrect(px - 0.08, px + w + 0.02, py + 0.03, py - 0.03 * num_labels)
    i = 0
    GR.uselinespec(" ")
    for (x, y, z, c, spec) in plt.args
        GR.savestate()
        mask = GR.uselinespec(spec)
        mask in (0, 1, 3, 4, 5) && GR.polyline([px - 0.07, px - 0.01], [py, py])
        mask & 0x02 != 0 && GR.polymarker([px - 0.06, px - 0.02], [py, py])
        GR.restorestate()
        GR.settextalign(GR.TEXT_HALIGN_LEFT, GR.TEXT_VALIGN_HALF)
        if i < num_labels
            i += 1
            text(px, py, plt.kvs[:labels][i])
        end
        py -= 0.03
    end
    GR.selntran(1)
    GR.restorestate()
end

function colorbar(off=0, colors=256)
    GR.savestate()
    viewport = plt.kvs[:viewport]
    zmin, zmax = plt.kvs[:zrange]
    GR.setwindow(0, 1, zmin, zmax)
    GR.setviewport(viewport[2] + 0.02 + off, viewport[2] + 0.05 + off,
                   viewport[3], viewport[4])
    l = zeros(Int32, 1, colors)
    l[1,:] = round(Int32, linspace(1000, 1255, colors))
    GR.cellarray(0, 1, zmax, zmin, 1, colors, l)
    diag = sqrt((viewport[2] - viewport[1])^2 + (viewport[4] - viewport[3])^2)
    charheight = max(0.016 * diag, 0.012)
    GR.setcharheight(charheight)
    if plt.kvs[:scale] & GR.OPTION_Z_LOG == 0
        ztick = 0.5 * GR.tick(zmin, zmax)
        GR.axes(0, ztick, 1, zmin, 0, 1, 0.005)
    else
        GR.setscale(GR.OPTION_Y_LOG)
        GR.axes(0, 2, 1, zmin, 0, 1, 0.005)
    end
    GR.restorestate()
end

function figure(; kv...)
    global plt
    plt = Figure()
    merge!(plt.kvs, Dict(kv))
    plt
end

function hold(flag)
    plt.kvs[:ax] = flag
    plt.kvs[:clear] = !flag
end

function subplot(nr, nc, p)
    xmin, xmax, ymin, ymax = 1, 0, 1, 0
    for i in collect(p)
        r = nr - div(i-1, nc)
        c = (i-1) % nc + 1
        xmin = min(xmin, (c-1)/nc)
        xmax = max(xmax, c/nc)
        ymin = min(ymin, (r-1)/nr)
        ymax = max(ymax, r/nr)
    end
    plt.kvs[:subplot] = [xmin, xmax, ymin, ymax]
    plt.kvs[:clear] = collect(p)[1] == 1
    plt.kvs[:update] = collect(p)[end] == nr * nc
end

function plot_img(I)
    viewport = plt.kvs[:viewport]
    vp = plt.kvs[:vp]

    if isa(I, AbstractString)
        width, height, data = GR.readimage(I)
    else
        width, height = size(I)
        data = (float(I) - minimum(I)) / (maximum(I) - minimum(I))
        data = round(Int32, 1000 + data * 255)
    end

    if width  * (viewport[4] - viewport[3]) <
       height * (viewport[2] - viewport[1])
        w = float(width) / height * (viewport[4] - viewport[3])
        xmin = max(0.5 * (viewport[1] + viewport[2] - w), viewport[1])
        xmax = min(0.5 * (viewport[1] + viewport[2] + w), viewport[2])
        ymin = viewport[3]
        ymax = viewport[4]
    else
        h = float(height) / width * (viewport[2] - viewport[1])
        xmin = viewport[1]
        xmax = viewport[2]
        ymin = max(0.5 * (viewport[4] + viewport[3] - h), viewport[3])
        ymax = min(0.5 * (viewport[4] + viewport[3] + h), viewport[4])
    end

    GR.selntran(0)
    if isa(I, AbstractString)
        GR.drawimage(xmin, xmax, ymin, ymax, width, height, data)
    else
        GR.cellarray(xmin, xmax, ymin, ymax, width, height, data)
    end

    if haskey(plt.kvs, :title)
        GR.savestate()
        GR.settextalign(GR.TEXT_HALIGN_CENTER, GR.TEXT_VALIGN_TOP)
        text(0.5 * (viewport[1] + viewport[2]), vp[4], plt.kvs[:title])
        GR.restorestate()
    end
    GR.selntran(1)
end

function plot_iso(V)
    viewport = plt.kvs[:viewport]

    if viewport[4] - viewport[3] < viewport[2] - viewport[1]
        width = viewport[4] - viewport[3]
        centerx = 0.5 * (viewport[1] + viewport[2])
        xmin = max(centerx - 0.5 * width, viewport[1])
        xmax = min(centerx + 0.5 * width, viewport[2])
        ymin = viewport[3]
        ymax = viewport[4]
    else
        height = viewport[2] - viewport[1]
        centery = 0.5 * (viewport[3] + viewport[4])
        xmin = viewport[1]
        xmax = viewport[2]
        ymin = max(centery - 0.5 * height, viewport[3])
        ymax = min(centery + 0.5 * height, viewport[4])
    end

    GR.selntran(0)
    values = round(UInt16, (V-minimum(V)) / (maximum(V)-minimum(V)) * (2^16-1))
    nx, ny, nz = size(V)
    isovalue = get(plt.kvs, :isovalue, 0.5)
    rotation = get(plt.kvs, :rotation, 40) * π / 180.0
    tilt = get(plt.kvs, :tilt, 70) * π / 180.0
    r = 2.5
    gr3.clear()
    mesh = gr3.createisosurfacemesh(values, (2/(nx-1), 2/(ny-1), 2/(nz-1)),
                                    (-1., -1., -1.),
                                    round(Int64, isovalue * (2^16-1)))
    if haskey(plt.kvs, :color)
        color = plt.kvs[:color]
    else
        color = (0.0, 0.5, 0.8)
    end
    gr3.setbackgroundcolor(1, 1, 1, 0)
    gr3.drawmesh(mesh, 1, (0, 0, 0), (0, 0, 1), (0, 1, 0), color, (1, 1, 1))
    gr3.cameralookat(r*sin(tilt)*sin(rotation), r*cos(tilt), r*sin(tilt)*cos(rotation), 0, 0, 0, 0, 1, 0)
    gr3.drawimage(xmin, xmax, ymin, ymax, 500, 500, gr3.DRAWABLE_GKS)
    gr3.deletemesh(mesh)
    GR.selntran(1)
end

function plot_polar(θ, ρ)
    window = plt.kvs[:window]
    rmin, rmax = window[3], window[4]
    ρ = (ρ - rmin) / (rmax - rmin)
    n = length(ρ)
    x, y = zeros(n), zeros(n)
    for i in 1:n
        x[i] = ρ[i] * cos(θ[i])
        y[i] = ρ[i] * sin(θ[i])
    end
    GR.polyline(x, y)
end

function plot_data(; kv...)
    merge!(plt.kvs, Dict(kv))

    if plt.args == @_tuple(Any)
        return
    end
    kind = get(plt.kvs, :kind, :line)

    plt.kvs[:clear] && GR.clearws()

    if kind in (:imshow, :isosurface)
        set_viewport(kind, plt.kvs[:subplot])
    elseif !plt.kvs[:ax]
        set_viewport(kind, plt.kvs[:subplot])
        set_window(kind)
        if kind == :polar
            draw_polar_axes()
        else
            draw_axes(kind)
        end
    end

    if haskey(plt.kvs, :colormap)
        GR.setcolormap(plt.kvs[:colormap])
    else
        GR.setcolormap(GR.COLORMAP_COOLWARM)
    end

    GR.uselinespec(" ")
    for (x, y, z, c, spec) in plt.args
        GR.savestate()
        if haskey(plt.kvs, :alpha)
            GR.settransparency(plt.kvs[:alpha])
        end
        if kind == :line
            mask = GR.uselinespec(spec)
            mask in (0, 1, 3, 4, 5) && GR.polyline(x, y)
            mask & 0x02 != 0 && GR.polymarker(x, y)
        elseif kind == :scatter
            GR.setmarkertype(GR.MARKERTYPE_SOLID_CIRCLE)
            if z != Void || c != Void
                if c != Void
                    c = (c - minimum(c)) / (maximum(c) - minimum(c))
                    cind = round(Int64, 1000 + c * 255)
                end
                for i in 1:length(x)
                    (z != Void) && GR.setmarkersize(z[i] / 100.0)
                    (c != Void )&& GR.setmarkercolorind(cind[i])
                    GR.polymarker([x[i]], [y[i]])
                end
            else
                GR.polymarker(x, y)
            end
        elseif kind == :stem
            GR.setlinecolorind(1)
            GR.polyline([plt.kvs[:window][1]; plt.kvs[:window][2]], [0; 0])
            GR.setmarkertype(GR.MARKERTYPE_SOLID_CIRCLE)
            GR.uselinespec(spec)
            for i = 1:length(y)
                GR.polyline([x[i]; x[i]], [0; y[i]])
                GR.polymarker([x[i]], [y[i]])
            end
        elseif kind == :hist
            ymin = plt.kvs[:window][3]
            for i = 2:length(y)
                GR.setfillcolorind(989)
                GR.setfillintstyle(GR.INTSTYLE_SOLID)
                GR.fillrect(x[i-1], x[i], ymin, y[i])
                GR.setfillcolorind(1)
                GR.setfillintstyle(GR.INTSTYLE_HOLLOW)
                GR.fillrect(x[i-1], x[i], ymin, y[i])
            end
        elseif kind == :contour
            zmin, zmax = plt.kvs[:zrange]
            GR.setspace(zmin, zmax, 0, 90)
            h = linspace(zmin, zmax, 20)
            if length(x) == length(y) == length(z)
                x, y, z = GR.gridit(x, y, z, 200, 200)
            end
            GR.contour(x, y, h, z, 1000)
            colorbar(0, 20)
        elseif kind == :contourf
            zmin, zmax = plt.kvs[:zrange]
            GR.setspace(zmin, zmax, 0, 90)
            if length(x) == length(y) == length(z)
                x, y, z = GR.gridit(x, y, z, 200, 200)
                z = reshape(z, 200, 200)
            end
            if plt.kvs[:scale] & GR.OPTION_Z_LOG != 0
                z = log(z)
            end
            GR.surface(x, y, z, GR.OPTION_CELL_ARRAY)
            colorbar()
        elseif kind == :wireframe
            if length(x) == length(y) == length(z)
                x, y, z = GR.gridit(x, y, z, 50, 50)
            end
            GR.setfillcolorind(0)
            GR.surface(x, y, z, GR.OPTION_FILLED_MESH)
            draw_axes(kind, 2)
        elseif kind == :surface
            if length(x) == length(y) == length(z)
                x, y, z = GR.gridit(x, y, z, 200, 200)
            end
            if get(plt.kvs, :accelerate, true)
                gr3.clear()
                GR.gr3.surface(x, y, z, GR.OPTION_COLORED_MESH)
            else
                GR.surface(x, y, z, GR.OPTION_COLORED_MESH)
            end
            draw_axes(kind, 2)
            colorbar(0.05)
        elseif kind == :plot3
            GR.polyline3d(x, y, z)
            draw_axes(kind, 2)
        elseif kind == :scatter3
            GR.polymarker3d(x, y, z)
            draw_axes(kind, 2)
        elseif kind == :imshow
            plot_img(z)
        elseif kind == :isosurface
            plot_iso(z)
        elseif kind == :polar
            GR.uselinespec(spec)
            plot_polar(x, y)
        elseif kind == :trisurf
            GR.trisurface(x, y, z)
            draw_axes(kind, 2)
            colorbar(0.05)
        end
        GR.restorestate()
    end

    if kind in (:line, :scatter, :stem) && haskey(plt.kvs, :labels)
        draw_legend()
    end

    if plt.kvs[:update]
        GR.updatews()
        if GR.isinline()
            return GR.show()
        end
    end
end

function plot_args(args; fmt=:xys)
    args = Any[args...]
    parsed_args = Any[]

    while length(args) > 0
        local x, y, z, c
        a = shift!(args)
        if isa(a, AbstractVecOrMat)
            elt = eltype(a)
            if elt <: Complex
                x = real(a)
                y = imag(a)
                z = Void
                c = Void
            elseif elt <: Real
                if fmt == :xys
                    if length(args) >= 1 &&
                       (isa(args[1], AbstractVecOrMat) && eltype(args[1]) <: Real || typeof(args[1]) == Function)
                        x = a
                        y = shift!(args)
                        z = Void
                        c = Void
                    else
                        y = a
                        n = isrowvec(y) ? size(y, 2) : size(y, 1)
                        x = linspace(1, n, n)
                        z = Void
                        c = Void
                    end
                elseif fmt == :xyac || fmt == :xyzc
                    if length(args) >= 3 &&
                        isa(args[1], AbstractVecOrMat) && eltype(args[1]) <: Real &&
                       (isa(args[2], AbstractVecOrMat) && eltype(args[2]) <: Real || typeof(args[2]) == Function) &&
                       (isa(args[3], AbstractVecOrMat) && eltype(args[3]) <: Real || typeof(args[3]) == Function)
                        x = a
                        y = shift!(args)
                        z = shift!(args)
                        c = shift!(args)
                    elseif length(args) >= 2 &&
                        isa(args[1], AbstractVecOrMat) && eltype(args[1]) <: Real &&
                       (isa(args[2], AbstractVecOrMat) && eltype(args[2]) <: Real || typeof(args[2]) == Function)
                        x = a
                        y = shift!(args)
                        z = shift!(args)
                        c = Void
                    elseif fmt == :xyac && length(args) >= 1 &&
                       (isa(args[1], AbstractVecOrMat) && eltype(args[1]) <: Real || typeof(args[1]) == Function)
                        x = a
                        y = shift!(args)
                        z = Void
                        c = Void
                    elseif fmt == :xyzc && length(args) == 0
                        z = a
                        nx, ny = size(z)
                        x = linspace(1, nx, nx)
                        y = linspace(1, ny, ny)
                        c = Void
                    end
                end
            else
                error("expected Real or Complex")
            end
        else
            error("expected array or function")
        end
        if typeof(z) == Function
            f = z
            z = Float64[f(a,b) for b in y, a in x]
        end
        spec = ""
        if fmt == :xys && length(args) > 0 && isa(args[1], AbstractString)
            spec = shift!(args)
        end
        push!(parsed_args, (x, y, z, c, spec))
    end

    pltargs = Any[]

    for arg in parsed_args
        x, y, z, c, spec = arg

        isa(x, UnitRange) && (x = collect(x))
        isa(y, UnitRange) && (y = collect(y))
        isa(z, UnitRange) && (z = collect(z))
        isa(c, UnitRange) && (c = collect(c))

        isvector(x) && (x = vec(x))
        if typeof(y) == Function
            y = [y(a) for a in x]
        else
            isvector(y) && (y = vec(y))
        end
        if z != Void
            if fmt == :xyzc && typeof(z) == Function
                z = [z(a,b) for a in x, b in y]
            else
                isvector(z) && (z = vec(z))
            end
        end
        if c != Void
            isvector(c) && (c = vec(c))
        end

        local xyzc
        if z == Void
            if isa(x, AbstractVector) && isa(y, AbstractVector)
                xyzc = [ (x, y, z, c) ]
            elseif isa(x, AbstractVector)
                xyzc = length(x) == size(y, 1) ?
                       [ (x, sub(y,:,j), z, c) for j = 1:size(y, 2) ] :
                       [ (x, sub(y,i,:), z, c) for i = 1:size(y, 1) ]
            elseif isa(y, AbstractVector)
                xyzc = size(x, 1) == length(y) ?
                       [ (sub(x,:,j), y, z, c) for j = 1:size(x, 2) ] :
                       [ (sub(x,i,:), y, z, c) for i = 1:size(x, 1) ]
            else
                @assert size(x) == size(y)
                xyzc = [ (sub(x,:,j), sub(y,:,j), z, c) for j = 1:size(y, 2) ]
            end
        elseif isa(x, AbstractVector) && isa(y, AbstractVector) &&
               (isa(z, AbstractVector) || typeof(z) == Array{Float64,2} ||
                typeof(z) == Array{Int32,2} || typeof(z) == Array{Any,2})
            xyzc = [ (x, y, z, c) ]
        else
            xyzc = [ (vec(float(x)), vec(float(y)), vec(float(z)), c) ]
        end
        for (x, y, z, c) in xyzc
            push!(pltargs, (x, y, z, c, spec))
        end
    end

    pltargs
end

function plot(args::PlotArg...; kv...)
    merge!(plt.kvs, Dict(kv))

    if plt.kvs[:ax]
        plt.args = append!(plt.args, plot_args(args))
    else
        plt.args = plot_args(args)
    end

    plot_data(kind=:line)
end

function oplot(args::PlotArg...; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = append!(plt.args, plot_args(args))

    plot_data(kind=:line)
end

function scatter(args...; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = plot_args(args, fmt=:xyac)

    plot_data(kind=:scatter)
end

function stem(args...; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = plot_args(args)

    plot_data(kind=:stem)
end

function histogram(x; kv...)
    merge!(plt.kvs, Dict(kv))

    h = Base.hist(x)
    x, y = float(collect(h[1])), float(h[2])
    plt.args = [(x, y, Void, Void, "")]

    plot_data(kind=:hist)
end

function contour(args...; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data(kind=:contour)
end

function contourf(args...; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data(kind=:contourf)
end

function wireframe(args...; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data(kind=:wireframe)
end

function surface(args...; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data(kind=:surface)
end

function plot3(args...; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data(kind=:plot3)
end

function scatter3(args...; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data(kind=:scatter3)
end

function title(s)
    plot_data(title=s)
end

function xlabel(s)
    plot_data(xlabel=s)
end

function ylabel(s)
    plot_data(ylabel=s)
end

function legend(args::AbstractString...; kv...)
    plot_data(labels=args)
end

function xlim(a)
    plot_data(xlim=a)
end

function ylim(a)
    plot_data(ylim=a)
end

function savefig(filename)
    GR.beginprint(filename)
    plot_data()
    GR.endprint()
end

function meshgrid{T}(vx::AbstractVector{T}, vy::AbstractVector{T})
    m, n = length(vy), length(vx)
    vx = reshape(vx, 1, n)
    vy = reshape(vy, m, 1)
    (repmat(vx, m, 1), repmat(vy, 1, n))
end

function meshgrid{T}(vx::AbstractVector{T}, vy::AbstractVector{T},
                     vz::AbstractVector{T})
    m, n, o = length(vy), length(vx), length(vz)
    vx = reshape(vx, 1, n, 1)
    vy = reshape(vy, m, 1, 1)
    vz = reshape(vz, 1, 1, o)
    om = ones(Int, m)
    on = ones(Int, n)
    oo = ones(Int, o)
    (vx[om, :, oo], vy[:, on, oo], vz[om, on, :])
end

function peaks(n=49)
    x = linspace(-2.5, 2.5, n)
    y = x
    x, y = meshgrid(x, y)
    3*(1-x).^2.*exp(-(x.^2) - (y+1).^2) - 10*(x/5 - x.^3 - y.^5).*exp(-x.^2-y.^2) - 1/3*exp(-(x+1).^2 - y.^2)
end

function imshow(I; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = [(Void, Void, I, Void, "")]

    plot_data(kind=:imshow)
end

function isosurface(V; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = [(Void, Void, V, Void, "")]

    plot_data(kind=:isosurface)
end

function cart2sph(x, y, z)
    azimuth = atan2(y, x)
    elevation = atan2(z, sqrt(x.^2 + y.^2))
    r = sqrt(x.^2 + y.^2 + z.^2)
    azimuth, elevation, r
end

function sph2cart(azimuth, elevation, r)
    x = r .* cos(elevation) .* cos(azimuth)
    y = r .* cos(elevation) .* sin(azimuth)
    z = r .* sin(elevation)
    x, y, z
end

function polar(args...; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = plot_args(args)

    plot_data(kind=:polar)
end

function trisurf(args...; kv...)
    merge!(plt.kvs, Dict(kv))

    plt.args = plot_args(args, fmt=:xyzc)

    plot_data(kind=:trisurf)
end

end # module


text:
mutable struct FileCachedFactor{T<:Unsigned, S<:AbstractFactor{T}} <: AbstractFactor{T}
    cachepath::String
    basefactor::S
    lockobj::Threads.TatasLock

    FileCachedFactor{T, S}(cachepath, basefactor, lockobj) where {S<:AbstractFactor{T}} where {T<:Unsigned} = (res = new(cachepath, basefactor, lockobj); finalizer(res, x -> rm(x.cachepath; force = true)); res)
end

Base.length(var::FileCachedFactor) = length(var.basefactor)

getname(var::FileCachedFactor) = getname(var.basefactor)

getlevels(var::FileCachedFactor) = getlevels(var.basefactor)

function FileCachedFactor(factor::S) where {S<:AbstractFactor{T}} where {T<:Unsigned}
    FileCachedFactor{T, S}("", factor, Threads.TatasLock())
end

function filecache(basefactor::S) where {S<:AbstractFactor{T}} where {T<:Unsigned}
    FileCachedFactor{T, S}("", basefactor, Threads.TatasLock())
end

function slice(factor::FileCachedFactor{T, S}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {T<:Unsigned} where {S<:Union{FileCachedFactor, CachedFactor, Factor, FileFactor}}
    basefactor = factor.basefactor
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    slice(basefactor, fromobs, toobs, slicelength)
end

function slice(factor::FileCachedFactor{T, S}, fromobs::Integer, toobs::Integer, slicelength::Integer) where {S<:AbstractFactor{T}} where {T<:Unsigned}
    basefactor = factor.basefactor
    slicelength = verifyslicelength(fromobs, toobs, slicelength)
    lockobj = factor.lockobj
    lock(lockobj)
    try
        if !isfile(factor.cachepath)
            v = convert(Vector{T}, basefactor)
            tmpjuml = joinpath(tempdir(), JUMLDIR)
            mkpath(tmpjuml)
            outdatpath = joinpath(tmpjuml, "$(randstring(10)).dat")
            open(outdatpath, "w") do f
                write(f, v)
            end
            factor.cachepath = outdatpath
        end
    finally
        unlock(lockobj)
    end
    sizeT = sizeof(T)
    iostream = open(factor.cachepath, "r")
    seek(iostream, sizeT * (fromobs - 1))
    buffer = Vector{T}(undef, slicelength)
    map(Seq(Vector{T}, (buffer, iostream, fromobs, toobs, slicelength), nextdatachunk), SubArray{T,1,Array{T,1},Tuple{UnitRange{Int64}},true}) do slice
        view(slice, 1:length(slice))
    end
end

function Base.map(factor::FileCachedFactor{T}, dataframe::AbstractDataFrame; permute::Bool = false) where {T<:Unsigned}
    FileCachedFactor(map(factor.basefactor, dataframe; permute = permute))
end

function isordinal(factor::FileCachedFactor{T}) where {T<:Unsigned}
    isordinal(factor.basefactor)
end

ext: jl

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct CipherSuitesEnum <: SwaggerModel

    function CipherSuitesEnum(;)
        o = new()
        o
    end
end # type CipherSuitesEnum

const _property_map_CipherSuitesEnum = Dict{Symbol,Symbol}()
const _property_types_CipherSuitesEnum = Dict{Symbol,String}()
Base.propertynames(::Type{ CipherSuitesEnum }) = collect(keys(_property_map_CipherSuitesEnum))
Swagger.property_type(::Type{ CipherSuitesEnum }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_CipherSuitesEnum[name]))}
Swagger.field_name(::Type{ CipherSuitesEnum }, property_name::Symbol) =  _property_map_CipherSuitesEnum[property_name]

function check_required(o::CipherSuitesEnum)
    true
end

function validate_property(::Type{ CipherSuitesEnum }, name::Symbol, val)
end


text:
function rootinvariant!(W2root::Vector{Float64},Q::Float64,A::Float64,Qp1::Float64,Ap1::Float64,
    c0::Float64,beta::Float64,rho::Float64,mu::Float64,diff::Float64,k::Float64,h::Float64)
    lb = Float64;
    W2 = Vector{Float64}[];
    dW2 = Float64;
    dQ = Float64;
    dA = Float64;
    Qint = Float64;
    Aint = Float64;

    # backward wave speed at proximal end
    lb = 0.5.*((Qp1./Ap1 .- sqrt.(0.5.*beta./rho).*Ap1.^0.25) .+
        (Q./A .- sqrt.(0.5.*beta./rho).*A.^0.25));
    # left-running invariant at current time
    W2 = [Q./A .- 4.*(sqrt.(0.5.*beta./rho).*A.^0.25 .- c0)];
    push!(W2,Qp1./Ap1 .- 4.*(sqrt.(0.5.*beta./rho).*Ap1.^0.25 .- c0))
    dW2 = (W2[2] .- W2[1])./k;
    # change in solution variables
    dQ = (Qp1 .- Q)./k;
    dA = (Ap1 .- A)./k;
    # interpolated state variables
    Qint = Q .- dQ.*lb.*h;
    Aint = A .- dA.*lb.*h;
    # advance proximal boundary invariant
    W2root[1] = W2[1] .- dW2.*lb.*h .- diff.*mu./rho.*h.*(Qint./Aint.^2);
end

max_stars_repo_name: dcanuto/CardioModeling

max_issues_repo_name: dcanuto/CardioModeling

max_forks_repo_name: dcanuto/CardioModeling


text:
using LinearAlgebraicRepresentation

include("./planar_arrangement.jl")
include("./dimension_travel.jl")
include("./largrid.jl")
include("./utilities.jl")
include("./interface.jl")
include("./struct.jl")
include("./mapper.jl")
include("./integr.jl")
include("./simplexn.jl")


text:
using Documenter
using WCS

makedocs(
    sitename = "WCS.jl",
    format = Documenter.HTML(),
    modules = [WCS],
    strict = true,
)

deploydocs(
    repo = "github.com/JuliaAstro/WCS.jl.git"
)


text:
module OutOfSampleBootstraps

using Docile
@docstrings

include("utils.jl")
include("recursive_ols.jl")
include("bootstrap_index.jl")
include("recursive_bootstrap.jl")

end # module


text:
module MolHandler

export Coordinate, norm, distance
export Atom, Attribute, Frame, Trajectory
export read_dcd, write_dcd, read_pdb, write_pdb, read_xyz, write_xyz
export get_frame, get_atom, clip_trajectory, geometric_center_of_mass, center_of_mass,
       pair_length_matrix, pair_length_matrix_pbc, pair_length_matrix_parallel,
       contact_bool_matrix, contact_bool_matrix_pbc, contact_bool_matrix_parallel,
       contact_probability_matrix, contact_probability_matrix_parallel,
       radius_of_gyration,
       atom_mass, residue_mass,
       distance_pbc, fix_pbc, fix_pbc!, move_pbc_center

# codes
include("variables.jl")
include("coordinate.jl")
include("component.jl")
include("fileio.jl")
include("handling.jl")

end # module


text:
struct S a end
function f(x::Int64) x.a end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
using Test

include("encounters.jl")

@testset "type abstraction" begin
    @test Dog <: Pet
    @test Cat <: Pet
end

# create some pets
buddy = Dog("Buddy")
sadie = Dog("Sadie")
minka = Cat("Minka")
felix = Cat("Felix")

@testset "names" begin
    @test name(buddy) == "Buddy"
    @test name(sadie) == "Sadie"
    @test name(minka) == "Minka"
    @test name(felix) == "Felix"
end

@testset "encounters" begin
    @test encounter(sadie, buddy) == "Sadie meets Buddy and sniffs."
    @test encounter(buddy, felix) == "Buddy meets Felix and chases."
    @test encounter(minka, sadie) == "Minka meets Sadie and hisses."
    @test encounter(felix, minka) == "Felix meets Minka and slinks."
end

# define a new type of Pet to test the fallback
# this belongs to the testset below but struct definitions within the local scope of testsets are not supported in Julia <1.1
struct Horse <: Pet
    name::String
end
name(h::Horse) = h.name

jenny = Horse("Jenny")

@testset "pet fallbacks" begin
    @test encounter(buddy, jenny) == "Buddy meets Jenny and is cautious."
    @test encounter(sadie, jenny) == "Sadie meets Jenny and is cautious."
    @test encounter(minka, jenny) == "Minka meets Jenny and is cautious."
    @test encounter(felix, jenny) == "Felix meets Jenny and is cautious."

    @test encounter(jenny, buddy) == "Jenny meets Buddy and is cautious."
    @test encounter(jenny, sadie) == "Jenny meets Sadie and is cautious."
    @test encounter(jenny, minka) == "Jenny meets Minka and is cautious."
    @test encounter(jenny, felix) == "Jenny meets Felix and is cautious."
end

struct Car
    license_plate::String
end
name(c::Car) = c.license_plate

car = Car("W-12345X")

@testset "non-pet fallback" begin
    @test encounter(buddy, car) == "Buddy meets W-12345X and runs away."
    @test encounter(minka, car) == "Minka meets W-12345X and runs away."
    @test encounter(jenny, car) == "Jenny meets W-12345X and runs away."
end

# Test generic fallback with a random type
# In practice this would be nonsense, but here we want to ensure the student doesn't implement a specific method for a given type
typename = Symbol("RandomType", rand(UInt))
@eval struct $typename end
@eval name(::$typename) = $(string(typename))
@eval t = $typename()

@testset "generic fallback" begin
    @test encounter(t, buddy) == "$(name(t)) meets Buddy and nothing happens."
    @test encounter(buddy, t) == "Buddy meets $(name(t)) and runs away."
end


text:
"""
    powermean

# Calculates the weighted powermean of a series of numbers

Calculates *order*th power mean of *values*, weighted by
*weights*. By default, *weights* are equal and *order*
is 1, so this is just the arithmetic mean.

# Arguments:

- `values`: values for which to calculate mean
- `order[s]`: order[s] of power mean
- `weights`: weights of elements, normalised to 1 inside function

# Returns:

- weighted power mean(s)
"""
function powermean end
function powermean(values::V1, order::R = 1,
                   weights::V2 = fill!(similar(values), 1)) where
    {R <: Real, FP <: AbstractFloat,
     V1 <: AbstractVector{FP}, V2 <: AbstractVector{FP}}
    length(values) == length(weights) ||
    throw(DimensionMismatch("powermean: Weight and value vectors must be the same length"))

    # Check whether all weights are zero in group.
    # In that case we want to return a NaN
    if iszero(weights)
        return convert(FP, NaN)
    end

    if isinf(order)
        if order > 0 # +Inf -> Maximum
            s = zero(FP)
            for i in eachindex(values, weights)
                @inbounds if (weights[i] > eps(FP)) & (values[i] > s)
                    s = values[i]
                end
            end
            return s
        else # -Inf -> Minimum
            s = convert(FP, Inf)
            for i in eachindex(values, weights)
                @inbounds if (weights[i] > eps(FP)) & (values[i] < s)
                    s = values[i]
                end
            end
            return s
        end
    else
        if order ≈ zero(order)
            s = zero(FP)
            for i in eachindex(values, weights)
                @inbounds if weights[i] > eps(FP)
                    s += weights[i] * log(values[i])
                end
            end
            return exp(s / sum(weights))
        else
            s = zero(FP)
            for i in eachindex(values, weights)
                @inbounds if weights[i] > eps(FP)
                    s += weights[i] * values[i] ^ order
                end
            end
            return (s / sum(weights)) ^ (one(FP) / order)
        end
    end
end

# This is the next most common case - a vector of orders
function powermean(values::V1,
                   orders::VR,
                   weights::V2 = fill!(similar(values), 1.0)) where
    {R <: Real, VR <: AbstractVector{R},
     FP <: AbstractFloat, V1 <: AbstractVector{FP}, V2 <: AbstractVector{FP}}
    return map(order -> powermean(values, order, weights), orders)
end

# This is the next most simple case - matrices with subcommunities, and an order or orders
function powermean(values::M1, orders, weights::M2 = fill!(similar(values), 1)) where
    {FP <: AbstractFloat, M1 <: AbstractMatrix{FP}, M2 <: AbstractMatrix{FP}}
    size(values) == size(weights) ||
        throw(DimensionMismatch("powermean: Weight and value matrixes " *
                                "must be the same size"))
    @views map(col -> powermean(values[:, col], orders,
                                weights[:, col]), 1:size(values, 2))
end

"""
    qD

Calculates Hill / naive-similarity diversity of order(s) *qs* of a
population with given relative proportions.

# Arguments:

- `proportions`: relative proportions of different types in population

- `qs`: single number or vector of orders of diversity measurement

# Returns:

- Diversity of order qs (single number or vector of diversities)

"""
function qD end
function qD(asm::AbstractAssemblage, qs)
    hassimilarity(asm) && error("function cannot run with $(typeof(gettypes(asm))) types as contains similarity")
    countsubcommunities(asm) == 1 ||
    throw(DimensionMismatch("Can only calculate diversity of a single community"))

    powermean(getabundance(asm), qs .- 1, getabundance(asm))[1] .^ -1
end

function qD(proportions::AbstractVector{<: Real}, qs)
    qD(Metacommunity(proportions), qs)
end

"""
    qDZ

Calculates Leinster-Cobbold / similarity-sensitive diversity of >= 1
order(s) *qs* of a population with given relative *proportions*, and
similarity matrix *Z*.

# Arguments:

- `proportions`: relative proportions of different types in a population

- `qs`: single number or vector of orders of diversity measurement

- `Z`: similarity matrix

# Returns:

- Diversity of order qs (single number or vector of diversities)

    """
function qDZ end

function qDZ(asm::AbstractAssemblage, qs)
    countsubcommunities(asm) == 1 ||
    throw(DimensionMismatch("Can only calculate diversity of a single community"))

    powermean(getordinariness!(asm), qs .- 1, getabundance(asm))[1] .^ -1
end

function qDZ(proportions::AbstractVector{<: Real}, qs,
             sim::AbstractTypes = UniqueTypes(size(proportions, 1)))
    qDZ(Metacommunity(proportions, sim), qs)
end

function qDZ(proportions::AbstractVector{<: Real}, qs, Z::AbstractMatrix{<: AbstractFloat})
    qDZ(Metacommunity(proportions, GeneralTypes(Z)), qs)
end


text:
"""
Value function for a policy on an MDP.

If `v` is a `DiscreteValueFunction`, access the value for a state with `v(s)`
"""
struct DiscreteValueFunction{M<:MDP} <: Function
    m::M
    values::Vector{Float64}
end

(v::DiscreteValueFunction)(s) = v.values[stateindex(v.m, s)]

"""
    evaluate(m::MDP, p::Policy)
    evaluate(m::MDP, p::Policy; rewardfunction=POMDPs.reward)

Calculate the value for a policy on an MDP using the approach in equation 4.2.2 of Kochenderfer, *Decision Making Under Uncertainty*, 2015.

Returns a DiscreteValueFunction, which maps states to values.

# Example
```
using POMDPModelTools, POMDPPolicies, POMDPModels
m = SimpleGridWorld()
u = evaluate(m, FunctionPolicy(x->:left))
u([1,1]) # value of always moving left starting at state [1,1]
```
"""
function evaluate(m::MDP, p::Policy; rewardfunction=POMDPs.reward)
    t = policy_transition_matrix(m, p)
    r = policy_reward_vector(m, p, rewardfunction=rewardfunction)
    u = (I-discount(m)*t)\r
    return DiscreteValueFunction(m, u)
end

"""
    policy_transition_matrix(m::Union{MDP, POMDP}, p::Policy)

Create an |S|x|S| sparse transition matrix for a given policy.

The row corresponds to the current state and column to the next state. Corresponds to ``T^π`` in equation (4.7) in Kochenderfer, *Decision Making Under Uncertainty*, 2015.
"""
function policy_transition_matrix(m::Union{MDP,POMDP}, p::Policy)
    rows = Int[]
    cols = Int[]
    probs = Float64[]
    state_space = states(m)
    ns = length(state_space)
    for s in state_space
      if !isterminal(m, s) # if terminal, the transition probabilities are all just zero
        si = stateindex(m, s)
        a = action(p, s)
        td = transition(m, s, a)
        for (sp, p) in weighted_iterator(td)
          if p > 0.0
            spi = stateindex(m, sp)
            push!(rows, si)
            push!(cols, spi)
            push!(probs, p)
          end
        end
      end
    end

    return sparse(rows, cols, probs, ns, ns)
end

function policy_reward_vector(m::Union{MDP,POMDP}, p::Policy; rewardfunction=POMDPs.reward)
    state_space = states(m)
    r = zeros(length(state_space))
    for s in state_space
        if !isterminal(m, s) # if terminal, the transition probabilities are all just zero
            si = stateindex(m, s)
            a = action(p, s)
            td = transition(m, s, a)
            for (sp, p) in weighted_iterator(td)
                if p > 0.0
                    r[si] += p*rewardfunction(m, s, a, sp)
                end
            end
        end
    end
    return r
end


text:
if Sys.iswindows()
    @info "No tests enabled for for your platform by default."
else
    include("gen.jl")

    ENV["TEST_SRVR_ASYNC"] = "true"
    include("srvr.jl")
    include("clnt.jl")

    include("memtransport_tests.jl")
    include("filetransport_tests.jl")
    include("utils_tests.jl")
end


text:
# Copyright © 2021, UChicago Argonne, LLC
# All Rights Reserved.

# This source code is licensed under the BSD-style license found in the
# LICENSE file in the root directory of this source tree.

function get_leaves(nodes)
    OrderedSet([n for n in nodes if isempty(intersect(Set([2*n, 2*n+1]), nodes))])
end

function get_nonleaves(nodes)
    leaves = get_leaves(nodes)
    setdiff(nodes, leaves)
end

function get_nodes_complete(nodes)
    """
    Return a set of nodes that its rooted subtree is a complete binary tree.
    Assume that `nodes' is a proper binary tree.
    """
    # nodes_arr = sort(collect(nodes), rev=true)
    height = Dict()
    nodes_complete = Set()
    for n in sort(collect(nodes), rev=true)
        if 2*n in nodes
            if height[2*n] >= 0 && height[2*n] == height[2*n+1]
                height[n] = height[2*n] + 1
                push!(nodes_complete, n)
            else
                height[n] = -1
            end
        else
            height[n] = 0
            push!(nodes_complete, n)
        end
    end

    OrderedSet(sort(collect(nodes_complete)))
end

function get_nodes_grandparents(nodes)
    """
    Return a set of nodes that has any of 4n, 4n+1, 4n+2, 4n+3
    """
    # nodes_arr = sort(collect(nodes), rev=true)
    nodes_grandparents = Set()
    for n in nodes
        if !isempty(intersect(Set([4*n 4*n+1 4*n+2 4*n+3]), nodes))
            push!(nodes_grandparents, n)
        end
    end

    OrderedSet(sort(collect(nodes_grandparents)))
end

function fill_single_child(nodes)
    """
    Add a missing child to make a full (proper) binary tree which every node has 0 or 2 children
    """

    nodes = Set(nodes)
    updated = true

    while updated
        updated = false
        for n in collect(nodes)
            if !isempty(intersect(nodes, Set([2*n, 2*n+1]))) && !issubset(Set([2*n, 2*n+1]), nodes)
                union!(nodes, Set([2*n, 2*n+1]))
                updated = true
                break
            end
        end
    end

    OrderedSet(sort(collect(nodes)))
end

function expand_nodes(nodes)
    """
    Attach both children at each leaf
    """

    for n in collect(nodes)
        union!(nodes, Set([n, 2*n, 2*n+1]))
    end
    OrderedSet(sort(collect(nodes)))
end

function get_nodes_by_depth(max_depth)
    OrderedSet(1:(2^(max_depth+1)-1))
end

function get_unused_nodes(ysol, nodes_ground, stepsize)
    """
    Return nodes that is large enough to serach k-neighbor of ysol given max_depth and stepsize (distance)
    """

    nodes = Set(keys(ysol))

    for i in 1:Int(floor((stepsize-1)/2))
        nodes = expand_nodes(nodes)
    end

    nodes = fill_single_child(nodes)

    nodes = setdiff(nodes_ground, nodes)

    OrderedSet(sort(collect(nodes)))
end

function get_ysol_csol(df_vars)
    df_vars_y = filter(row -> row.vartype=="y" && row.sol > 0.5, df_vars)
    ysol = Dict(row.nodeid => row.operator for row in eachrow(df_vars_y))
    csol = Dict()
    for (k,v) in ysol
        if v == 'C'
            df = filter(row -> row.vartype=="c" && row.nodeid == k, df_vars)
            csol[k] = df.sol[1]
        end
    end
    ysol, csol
end

function get_formula(ysol, csol)
    """
    ysol: Dict(nodeid -> operator)
    csol: Dict(nodeid -> constant value)
    """

    if !(1 in keys(ysol))
        return SortedDict(1=>"0")
    end

    formula = SortedDict()
    for (n,o) in ysol
        if o isa Integer
            formula[n] = "x$o"
        elseif o == 'C'
            formula[n] = "($(csol[n]))"
        elseif o == 'P'
            formula[n] = "pi"
        else
            formula[n] = nothing
        end
    end

    nodes = reverse(sort(collect(keys(formula))))
    for n in nodes
        if isnothing(formula[n])
            if ysol[n] in Set("+-*")
                formula[n] = "($(formula[2*n])$(ysol[n])$(formula[2*n+1]))"
            elseif ysol[n] == 'D'
                formula[n] = "($(formula[2*n])/$(formula[2*n+1]))"
            elseif ysol[n] == 'R'
                formula[n] = "sqrt($(formula[2*n+1]))"
            elseif ysol[n] == 'S'
                formula[n] = "sqrt(1-$(formula[2*n+1])^2)"
            elseif ysol[n] == 'E'
                formula[n] = "exp($(formula[2*n+1]))"
            elseif ysol[n] == 'L'
                formula[n] = "log($(formula[2*n+1]))"
            else
                println("tree2formula():: somrthing is wrong")
            end
        end
    end

    return formula
end

function compute_mse(y_pred, y_true; model="mse")
    err = abs.(y_pred - y_true)

    if model == "mse"
        ##### sqaured
        return Statistics.mean(err.^2)
    elseif model == "rmse"
        return sqrt(Statistics.mean(err.^2))
    elseif model == "alternative"
        ##### alternative of sqaured
        return Statistics.mean(log2.(1 .+ err.*2^30))
    else
        return 1e+06
    end
end

function compute_err_formula(obs, formula, errmodel; err2Inf=false)
    nrows, ncols = size(obs)
    nvars = ncols - 1

    y_true = Float64[]
    y_pred = Float64[]
    for i = 1:nrows
        expr = join(["x$(j) = $(obs[i,j]); " for j = 1:nvars]) * formula

        try
            push!(y_pred, eval(Meta.parse(expr)))
            push!(y_true, obs[i,end])
        catch DomainError
            nothing
            # return err2Inf ? Inf : "DomainError"
        end
    end

    if length(y_pred) <= 0.5 * nrows
        return err2Inf ? Inf : "DomainError"
    else
        return compute_mse(y_pred, y_true, model=errmodel)
    end
end

function read_parameters(file)
    params = Dict{String,Any}()

    for line in eachline(file)
        if length(line) == 0 || line[1] == '#'
            continue
        end

        chunks = split(line, limit=3, keepempty=false)
        if length(chunks) < 2
            continue
        end

        key = filter(x -> !isspace(x), chunks[1])
        value = filter(x -> !isspace(x), chunks[2])

        if length(key) >=4
            if "DBL_" == key[1:4]
                value = parse(Float64, value)
            elseif "INT_" == key[1:4]
                value = parse(Int, value)
            else
                nothing
            end
        end

        params[key] = value
    end

    params
end

function save_parameters(file, params)
    open(file, "w+") do io
        for (k,v) in params
            write(io, "$k\t$v\n")
        end
    end
end

function open_logger(logfile, loglevel)
    io = open(logfile, "a+")
    if loglevel == "Debug"
        logger = SimpleLogger(io, Logging.Debug)
    else
        logger = SimpleLogger(io, Logging.Info)
    end
    global_logger(logger)
    io
end

max_stars_repo_name: jongeunkim/STreCH

max_issues_repo_name: jongeunkim/STreCH

max_forks_repo_name: jongeunkim/STreCH


text:
@testset verbose=true "Automatic differentiation   " begin
	n = 50
	m = 5

	# create example functions
	coeffs3 = randn(n, n, n)
	coeffs1 = randn(n)
	coeffs0 = randn();
	f(x) = sum(coeffs3[i,j,k]*x[i]*x[j]*x[k] for i in 1:n, j in 1:n, k in 1:n) +
		dot(coeffs1, x) + coeffs0

	coeffs_c = randn(n-m+1, m)
	function c!(cval, x)
		for i in 1:m
			cval[i] = sum(exp.(-1.0.*view(coeffs_c, :, i).*view(x, i:n-m+i))) + x[i]^3*x[i+1]^2
		end
	end

	# ------------ test first-order derivatives -------------
	g = zeros(n)
	cval = zeros(m)
	cval2 = zeros(m)
	J = zeros(m, n)
	x0 = zeros(n)

	x = randn(n)
	xcopy = copy(x)		# make a copy to test for mutation

	actualg = copy(coeffs1)
	for l in 1:n
		actualg[l] += sum(coeffs3[l,j,k]*x[j]*x[k] for j in 1:n, k in 1:n) +
			sum(coeffs3[i,l,k]*x[i]*x[k] for i in 1:n, k in 1:n) +
			sum(coeffs3[i,j,l]*x[i]*x[j] for i in 1:n, j in 1:n)
	end

	actualJ = zeros(m, n)
	for i in 1:m
		actualJ[i, i:n-m+i] .= -view(coeffs_c, :, i) .* exp.(-1.0.*view(coeffs_c, :, i).*view(x, i:n-m+i))
		actualJ[i, i] += 3*x[i]^2*x[i+1]^2
		actualJ[i, i+1] += 2*x[i]^3*x[i+1]
	end

	# objective function gradient
	@testset "∇f" begin
		grad! = LFPSQP.generate_gradient(f, x0)
		grad!(g, x)

		@test all(g .≈ actualg)
		@test all(x .== xcopy)

		fill!(g, 0.0)
	end

	# objective function gradient with tape
	@testset "∇f (tape)" begin
		grad! = LFPSQP.generate_gradient_tape(f, x0)
		grad!(g, x)

		@test all(g .≈ actualg)
		@test all(x .== xcopy)

		fill!(g, 0.0)
	end

	# constraint Jacobian
	@testset "J_c" begin
		jac! = LFPSQP.generate_jacobian(c!, x0, m)
		jac!(J, cval, x)
		c!(cval2, x)

		@test all(J .≈ actualJ)
		@test all(cval .≈ cval2)
		@test all(x .== xcopy)

		fill!(J, 0.0)
		fill!(cval, 0.0)
	end

	# constraint Jacobian with tape
	@testset "J_c (tape)" begin
		jac! = LFPSQP.generate_jacobian_tape(c!, x0, m)
		jac!(J, cval, x)
		c!(cval2, x)

		@test all(J .≈ actualJ)
		@test all(cval .≈ cval2)
		@test all(x .== xcopy)

		fill!(J, 0.0)
		fill!(cval, 0.0)
	end

	# test Hessian vector product
	dest = zeros(n)
	dest2 = zeros(n)
	x0dual = zeros(Dual{nothing,Float64,1}, n)
	x0dual .= x0

	function actual_prod!(dest, src, x, λ)
		for l in 1:n
			dest[l] = sum(coeffs3[l,p,k]*x[k]*src[p] for k in 1:n, p in 1:n) +
				sum(coeffs3[p,l,k]*x[k]*src[p] for k in 1:n, p in 1:n) +
				sum(coeffs3[p,j,l]*x[j]*src[p] for j in 1:n, p in 1:n) +
				sum(coeffs3[l,j,p]*x[j]*src[p] for j in 1:n, p in 1:n) +
				sum(coeffs3[i,l,p]*x[i]*src[p] for i in 1:n, p in 1:n) +
				sum(coeffs3[i,p,l]*x[i]*src[p] for i in 1:n, p in 1:n)
		end

		for i in 1:m
			dest[i:n-m+i] .+= λ[i].*view(coeffs_c, :, i).^2 .* exp.(-1.0.*view(coeffs_c, :, i).*view(x, i:n-m+i)) .* view(src, i:n-m+i)

			dest[i] += λ[i]*(6*x[i]*x[i+1]^2*src[i] + 6*x[i]^2*x[i+1]*src[i+1])
			dest[i+1] += λ[i]*(6*x[i]^2*x[i+1]*src[i] + 2*x[i]^3*src[i+1])
		end
	end

	grad_dual! = LFPSQP.generate_gradient(f, x0dual)
	jac_dual! = LFPSQP.generate_jacobian(c!, x0dual, m)

	hess_lag_vec! = LFPSQP.generate_hess_lag_vec(grad_dual!, jac_dual!, x0, m)

	# try a few different products
	@testset "Hessian-vector product" begin
		for trial = 1:5
			λ = randn(m)
			λcopy = copy(λ)
			v = randn(n)
			vcopy = copy(v)

			hess_lag_vec!(dest, v, x, λ)
			actual_prod!(dest2, v, x, λ)

			@test all(dest .≈ dest2)
			@test all(λ .== λcopy)
			@test all(v .== vcopy)
		end
	end

end


text:
"""
    BitsVector{V}(data, bits_per_element, size)

Packed vector of 1, 2, or 4-bit entries.
"""
struct BitsVector <: AbstractVector{UInt8}
    data::Ptr{UInt8}
    bits_per_element::UInt8
    size::UInt
end

"""
    ScatteredBitsVector{V}(data_sectors, bits_per_element, size, n_blocks)

Scattered packed vectors of 1, 2, or 4-bit entries. Length of each block is 2^16,
except for the last one.
"""
struct ScatteredBitsVector{V} <: AbstractVector{UInt8}
    data_sectors::Vector{Base.RefValue{V}}
    bits_per_element::UInt8
    size::UInt
    n_blocks::UInt
end

"""
    ScatteredVector{T, V}(data_sectors, size, n_blocks)

Scattered vectors of type T. Length of each block is 2^16, except for the last one.
"""
struct ScatteredVector{T, V} <: AbstractVector{T}
    data_sectors::Vector{Base.RefValue{V}}
    size::Int
    n_blocks::Int
end

function Base.size(x::Union{BitsVector, ScatteredVector, ScatteredBitsVector})
    return (x.size,)
end

@inline function Base.getindex(x::BitsVector, i::Int)
    elements_per_byte = 8 ÷ x.bits_per_element
    byte_index = (i - 1) ÷ elements_per_byte + 1
    in_byte_index = (i - 1) % elements_per_byte
    byte = unsafe_load(x.data, byte_index)# x.data[][byte_index]
    (byte >> (x.bits_per_element * in_byte_index)) & mask_map[x.bits_per_element]
end

@inline function Base.getindex(x::ScatteredVector, i::Int)
    block_index = (i - 1) ÷ (2 ^ 16) + 1
    in_block_index = (i - 1) % (2 ^ 16) + 1
    x.data_sectors[block_index][][in_block_index]
end

@inline function Base.getindex(x::ScatteredBitsVector, i::Int)
    elements_per_byte = 8 ÷ x.bits_per_element
    block_index = (i - 1) ÷ (2 ^ 16) + 1
    in_block_index = ((i - 1) % (2 ^ 16)) ÷ elements_per_byte + 1
    in_byte_index = (i - 1) % elements_per_byte # 0-based
    byte = x.data_sectors[block_index][][in_block_index]
    (byte >> (x.bits_per_element * in_byte_index)) & mask_map[x.bits_per_element]
end

"""
    DiffList{V,W,X,Y}(len, sample_id_bases, last_component_sizes, has_genotypes,
    genotypes, sample_id_increments)

Data structure for difflists.
"""
mutable struct DiffList{SIBT,GT<:Union{Nothing,BitsVector}}
    len::UInt32
    sample_id_bases::SIBT
    last_component_sizes::Ptr{UInt8}
    has_genotypes::Bool
    genotypes::GT
    sample_id_increments::Ptr{UInt8}
end


text:
# The Immerse backend
# To the Cairo backend, this adds just one feature: keeping track of
# the rendered coordinates of tagged objects and svgclass("plotpanel")
# objects

export ImmerseBackend

mutable struct ImmerseBackend <: Backend
    cb::CAIROSURFACE
    coords::Dict{Symbol,Any}
    panelcoords::Vector
end

function ImmerseBackend(c::CairoSurface)
    cb = CAIROSURFACE(c)
    ImmerseBackend(cb, Dict{Symbol,Any}(), Any[])
end

function (img::ImmerseBackend)(x)
    draw(img, x)
end

root_box(backend::ImmerseBackend) = root_box(backend.cb)

push_property_frame(backend::ImmerseBackend, properties) =
        push_property_frame(backend.cb, properties)
pop_property_frame(backend::ImmerseBackend) = pop_property_frame(backend.cb)

function draw(backend::ImmerseBackend, root_container::Container)
    empty!(backend.coords)
    empty!(backend.panelcoords)
    drawpart(backend, root_container, IdentityTransform(), UnitBox(), root_box(backend))
    finish(backend)
end

function register_coords(backend::ImmerseBackend, box, units, transform, form::Form)
    if form.tag != empty_tag
        backend.coords[form.tag] = (box, units, transform)
    end
    nothing
end

function register_coords(backend::ImmerseBackend, box, units, transform, property::SVGClass)
    length(property.primitives) == 1 && property.primitives[1].value == "plotpanel" &&
            push!(backend.panelcoords, (box, units, transform))
    nothing
end

absolute_native_units(backend::ImmerseBackend, u::Float64) = absolute_native_units(backend.cb, u)
absolute_native_units(backend::ImmerseBackend, l::Length{:mm}) =
        absolute_native_units(backend.cb, l.value)
absolute_native_units(backend::ImmerseBackend, l::Tuple{Length{:mm},Length{:mm}}) =
        (absolute_native_units(backend.cb, l[1].value),
         absolute_native_units(backend.cb, l[2].value))

surface(backend::ImmerseBackend) = surface(backend.cb)

draw(backend::ImmerseBackend, form::Form) = draw(backend.cb, form)

finish(backend::ImmerseBackend) = finish(backend.cb)

iswithjs(::ImmerseBackend) = false
iswithousjs(::ImmerseBackend) = true


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5.3

# Short description of this file: Test Neural PDE low level

using Pkg
Pkg.add("ModelingToolkit")
Pkg.add("DiffEqFlux")
Pkg.add("Flux")
Pkg.add("Plots")
Pkg.add("Test")
Pkg.add("Optim")
Pkg.add("CUDA")

using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
using Plots #GalacticOptim
using PyPlot

@parameters t, x, y, θ
@variables u(..)
@derivatives Dt'~t
@derivatives Dx'~x
@derivatives Dy'~y
@derivatives Dxx''~x
@derivatives Dyy''~y

#2D PDE
eq  = Dt(u(t,x,y,θ)) + u(t,x,y,θ)*Dx(u(t,x,y,θ)) + u(t,x,y,θ)*Dy(u(t,x,y,θ)) ~ Dxx(u(t,x,y,θ)) + Dyy(u(t,x,y,θ))

# Initial and boundary conditions
bcs = [u(0,x,y,θ) ~ exp(x+y)*cos(x+y) ,
       u(2,x,y,θ) ~ exp(x+y)*cos(x+y+4*2) ,
       u(t,0,y,θ) ~ exp(y)*cos(y+4t),
       u(t,2,y,θ) ~ exp(2+y)*cos(2+y+4t) ,
       u(t,x,0,θ) ~ exp(x)*cos(x+4t),
       u(t,x,2,θ) ~ exp(x+2)*cos(x+2+4t)]

# Space and time domains
domains = [x ∈ IntervalDomain(0.0,2.0),
           y ∈ IntervalDomain(0.0,2.0),
           t ∈ IntervalDomain(0.0,2.0)]# Discretization

dx = 0.125; dy= 0.125; dt = 0.125

# Neural network
chain = FastChain(FastDense(3,16,Flux.σ),FastDense(16,16,Flux.σ),FastDense(16,1))

#strategy = GridTraining()
strategy = StochasticTraining()
discretization = PhysicsInformedNN([dt,dx,dy],chain,strategy=strategy)

indvars = [t,x,y]
depvars = [u]
dim = length(domains)

expr_pde_loss_function = build_loss_function(eq,indvars,depvars)
expr_bc_loss_functions = [build_loss_function(bc,indvars,depvars) for bc in bcs]
train_sets = generate_training_sets(domains,dx,bcs,indvars,depvars)

train_domain_set, train_bound_set, train_set = train_sets

phi = discretization.phi
autodiff = discretization.autodiff
derivative = discretization.derivative
initθ = discretization.initθ

pde_loss_function = get_loss_function(eval(expr_pde_loss_function),
                                      train_domain_set,
                                      phi,
                                      derivative,
                                      strategy)
bc_loss_function = get_loss_function(eval.(expr_bc_loss_functions),
                                     train_bound_set,
                                     phi,
                                     derivative,
                                     strategy)

function loss_function(θ,p)
    return pde_loss_function(θ) + bc_loss_function(θ)
end

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

f = NeuralPDE.OptimizationFunction(loss_function, GalacticOptim.AutoZygote())

prob = GalacticOptim.OptimizationProblem(f, initθ)

# optimizer
opt = GalacticOptim.ADAM(0.01)
res = GalacticOptim.solve(prob, opt; cb = cb, maxiters=30)
prob
phi = discretization.phi

ts,xs,ys = [domain.domain.lower:dx:domain.domain.upper for domain in domains]
#u_predict_contourf = reshape([first(phi([t,x,y],res.minimizer)) for t in ts for x in xs for y in ys] ,length(ys),length(xs),length(ts))
#Plots.plot(ts, xs, ys, u_predict_contourf, linetype=:contourf,title = "predict")

#u_predict = [[first(phi([t,x,y],res.minimizer)) for x in xs for y in ys] for t in ts ]
u_predict = [reshape([first(phi([t,x,y],res.minimizer)) for x in xs for y in ys], (length(xs),length(ys))) for t in ts ]
p1= Plots.plot(xs, u_predict[2],title = "t = 0.1");
p2= Plots.plot(xs, u_predict[6],title = "t = 0.5");
p3= Plots.plot(xs, u_predict[end],title = "t = 1");
Plots.plot(p1,p2,p3)


text:
# This script gets run once, on a developer's machine, to generate artifacts
using Pkg.Artifacts
using Downloads

function generate_artifacts(ver="2.3.0", repo="https://github.com/sglyon/PlotlyBase.jl")
    artifacts_toml = joinpath(dirname(@__DIR__), "Artifacts.toml")

    # if Artifacts.toml does not exist we also do not have to remove it
    isfile(artifacts_toml) && rm(artifacts_toml)

    plotschema_url = "https://raw.githubusercontent.com/plotly/plotly.js/v$(ver)/dist/plot-schema.json"

    plotlyartifacts_hash = create_artifact() do artifact_dir
        @show artifact_dir
        Downloads.download(plotschema_url, joinpath(artifact_dir, "plot-schema.json"))
        cp(joinpath(dirname(@__DIR__), "templates"), joinpath(artifact_dir, "templates"))
    end

    tarball_hash = archive_artifact(plotlyartifacts_hash, "plotly-base-artifacts-$ver.tar.gz")

    bind_artifact!(artifacts_toml, "plotly-base-artifacts", plotlyartifacts_hash; download_info=[
        ("$repo/releases/download/plotly-base-artifacts-$ver/plotly-base-artifacts-$ver.tar.gz", tarball_hash)
    ])
end


text:
__precompile__()
module SpikingNetworks
include("utility.jl")
include("codegen.jl")
include("group.jl")
include("unit.jl")
include("parse.jl")
include("init.jl")
include("stateMonitor.jl")
include("spike.jl")
include("reset.jl")
include("synapse.jl")
include("plasticity.jl")
end # module


text:
using Revise
using westbrick

num_objects = 5
positions = [[i,j] for (i,j) in zip([4.,3.,2.,2.5,0.3],[2.,3.,4.,2.,4.])]
objects = init_environment(num_objects, positions)
obs_dict = visualize_environment!(objects)

bobby = load_robot()
traj, obj_traj = init_visualization()

pose = [positions[1]...,π]
move!(bobby, pose, traj, obj_traj)
pick!(bobby, π, objects[1], traj, obj_traj)
move!(bobby, [-3.,-3., 0.], traj, obj_traj)
place!(bobby, -π, traj, obj_traj)
move!(bobby, [3.,-3., 0.], traj, obj_traj)

pose = [positions[5]...,π]
move!(bobby, pose, traj, obj_traj)
pick!(bobby, π, objects[5], traj, obj_traj)
move!(bobby, [-2.,-3., 0.], traj, obj_traj)
place!(bobby, -π, traj, obj_traj)
move!(bobby, [3.,-3., 0.], traj, obj_traj)

println("visualizing...")
visualize_trajectory!(bobby, traj, obj_traj, obs_dict)


text:
using DSP, Plots, LaTeXStrings, Interact

@manipulate for x₁ = -5:0.1:5, y₁ = -5:0.1:5, x₂ = -5:0.1:5, y₂ = -5:0.1:5
    O = 0
    z₁ = x₁ + y₁ * 1im
    z₂ = x₂ + y₂ * 1im

    plot( #plot z₁
        [real(O), real(z₁)], #tail point
        [imag(O), imag(z₁)], #head point
        color = :red,
        arrow = :true,
        legend = :true,
        lab = "z_1",
        alpha = 1,
        xlab = "mathrm{Re}",
        ylab = "mathrm{Im}",
        xlims = [-6, 6],
        ylims = [-6, 6],
        framestyle = :origin,
    )

    plot!( #plot z₂
        [real(O), real(z₂)], #tail point
        [imag(O), imag(z₂)], #head point
        color = :blue,
        arrow = :true,
        legend = :true,
        lab = "z_2",
        alpha = 1,
    )

    plot!(
        [real(z₁), real(z₁ + z₂)],
        [imag(z₁), imag(z₁ + z₂)],
        color = :blue,
        arrow = :true,
        legend = :true,
        lab = "z_1",
        alpha = 0.25,
    )

    plot!(
        [real(z₂), real(z₁ + z₂)],
        [imag(z₂), imag(z₁ + z₂)],
        color = :red,
        arrow = :true,
        legend = :true,
        lab = "z_2",
        alpha = 0.25,
    )

    plot!(
        [real(O), real(z₁ + z₂)],
        [imag(O), imag(z₁ + z₂)],
        color = :green,
        arrow = :true,
        legend = :true,
        lab = "z_1+z_2",
        alpha = 1,
    )
end

max_stars_repo_name: sashamae/MyTestRepo

max_issues_repo_name: sashamae/MyTestRepo

max_forks_repo_name: sashamae/MyTestRepo


text:
using Test
import VecMathOptInterface
const MOI = VecMathOptInterface
const MOIT = MOI.DeprecatedTest
const MOIU = MOI.Utilities

@testset "Basic Constraint Tests" begin
    mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
    config = MOIT.Config()
    MOIT.basic_constraint_tests(mock, config)
end

@testset "Unit Tests" begin
    # `UniversalFallback` needed for `MOI.Silent`
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
    # Optimizers attributes have to be set to default value since the mock
    # optimizer doesn't handle this
    MOI.set(mock, MOI.Silent(), true)
    MOI.set(mock, MOI.TimeLimitSec(), nothing)
    MOI.set(mock, MOI.NumberOfThreads(), nothing)
    config = MOIT.Config()
    for model in [
        mock,
        MOIU.CachingOptimizer(
            MOIU.UniversalFallback(MOIU.Model{Float64}()),
            mock,
        ),
    ]
        MOIT.unittest(
            model,
            config,
            [
                "solve_blank_obj",
                "solve_constant_obj",
                "solve_singlevariable_obj",
                "solve_with_lowerbound",
                "solve_with_upperbound",
                "solve_affine_lessthan",
                "solve_affine_greaterthan",
                "solve_affine_equalto",
                "solve_affine_interval",
                "solve_duplicate_terms_scalar_affine",
                "solve_duplicate_terms_vector_affine",
                "solve_qp_edge_cases",
                "solve_qp_zero_offdiag",
                "solve_qcp_edge_cases",
                "solve_affine_deletion_edge_cases",
                "solve_duplicate_terms_obj",
                "solve_integer_edge_cases",
                "solve_objbound_edge_cases",
                "raw_status_string",
                "solve_time",
                "solve_zero_one_with_bounds_1",
                "solve_zero_one_with_bounds_2",
                "solve_zero_one_with_bounds_3",
                "solve_start_soc",
                "solve_unbounded_model",
                "solve_single_variable_dual_min",
                "solve_single_variable_dual_max",
                "solve_result_index",
                "solve_farkas_equalto_lower",
                "solve_farkas_equalto_upper",
                "solve_farkas_lessthan",
                "solve_farkas_greaterthan",
                "solve_farkas_interval_lower",
                "solve_farkas_interval_upper",
                "solve_farkas_variable_lessthan",
                "solve_farkas_variable_lessthan_max",
                "solve_twice",
            ],
        )
        MOI.empty!(model)
    end

    @testset "solve_blank_obj" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1]),
                MOI.FEASIBLE_POINT,
            ),
        )
        MOIT.solve_blank_obj(mock, config)
        # The objective is blank so any primal value ≥ 1 is correct
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [2]),
                MOI.FEASIBLE_POINT,
            ),
        )
        MOIT.solve_blank_obj(mock, config)
    end
    @testset "solve_constant_obj" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1]),
                MOI.FEASIBLE_POINT,
            ),
        )
        MOIT.solve_constant_obj(mock, config)
    end
    @testset "solve_singlevariable_obj" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1]),
                MOI.FEASIBLE_POINT,
            ),
        )
        MOIT.solve_singlevariable_obj(mock, config)
    end
    @testset "solve_with_lowerbound" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1]),
                MOI.FEASIBLE_POINT,
                (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [2.0],
                (MOI.VariableIndex, MOI.LessThan{Float64}) => [0.0],
            ),
        )
        # x has two variable constraints
        mock.eval_variable_constraint_dual = false
        MOIT.solve_with_lowerbound(mock, config)
        mock.eval_variable_constraint_dual = true
    end
    @testset "solve_with_upperbound" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1]),
                MOI.FEASIBLE_POINT,
                (MOI.VariableIndex, MOI.LessThan{Float64}) => [-2.0],
                (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [0.0],
            ),
        )
        # x has two variable constraints
        mock.eval_variable_constraint_dual = false
        MOIT.solve_with_upperbound(mock, config)
        mock.eval_variable_constraint_dual = true
    end
    @testset "solve_affine_lessthan" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.5]),
                MOI.FEASIBLE_POINT,
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-0.5],
            ),
        )
        MOIT.solve_affine_lessthan(mock, config)
    end
    @testset "solve_affine_greaterthan" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.5]),
                MOI.FEASIBLE_POINT,
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [0.5],
            ),
        )
        MOIT.solve_affine_greaterthan(mock, config)
    end
    @testset "solve_affine_equalto" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.5]),
                MOI.FEASIBLE_POINT,
                (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                    [0.5],
            ),
        )
        MOIT.solve_affine_equalto(mock, config)
    end
    @testset "solve_affine_interval" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [2.0]),
                MOI.FEASIBLE_POINT,
                (MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}) => [-1.5],
            ),
        )
        MOIT.solve_affine_interval(mock, config)
    end
    @testset "solve_duplicate_terms_scalar_affine" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.5]),
                MOI.FEASIBLE_POINT,
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-0.5],
            ),
        )
        MOIT.solve_duplicate_terms_scalar_affine(mock, config)
    end
    @testset "solve_duplicate_terms_vector_affine" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.5]),
                MOI.FEASIBLE_POINT,
                (MOI.VectorAffineFunction{Float64}, MOI.Nonpositives) =>
                    [[-0.5]],
            ),
        )
        MOIT.solve_duplicate_terms_vector_affine(mock, config)
    end

    @testset "solve_qcp_edge_cases" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.5, 0.5]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.5, (√13 - 1) / 4]),
            ),
        )
        MOIT.solve_qcp_edge_cases(mock, config)
    end

    @testset "solve_qp_edge_cases" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0, 2.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0, 2.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0, 2.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0, 2.0]),
            ),
        )
        MOIT.solve_qp_edge_cases(mock, config)
    end
    @testset "solve_qp_zero_offdiag" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0, 2.0]),
            ),
        )
        MOIT.solve_qp_zero_offdiag(mock, config)
    end
    @testset "solve_duplicate_terms_obj" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1]),
                MOI.FEASIBLE_POINT,
            ),
        )
        MOIT.solve_duplicate_terms_obj(mock, config)
    end
    @testset "solve_affine_deletion_edge_cases" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [2.0]),
            ),
        )
        MOIT.solve_affine_deletion_edge_cases(mock, config)
    end
    @testset "solve_integer_edge_cases" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [2.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.0]),
            ),
        )
        MOIT.solve_integer_edge_cases(mock, config)
    end
    @testset "solve_objbound_edge_cases" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> begin
                MOI.set(mock, MOI.ObjectiveBound(), 3.0)
                MOIU.mock_optimize!(
                    mock,
                    MOI.OPTIMAL,
                    (MOI.FEASIBLE_POINT, [2.0]),
                )
            end,
            (mock::MOIU.MockOptimizer) -> begin
                MOI.set(mock, MOI.ObjectiveBound(), 3.0)
                MOIU.mock_optimize!(
                    mock,
                    MOI.OPTIMAL,
                    (MOI.FEASIBLE_POINT, [1.0]),
                )
            end,
            (mock::MOIU.MockOptimizer) -> begin
                MOI.set(mock, MOI.ObjectiveBound(), 2.0)
                MOIU.mock_optimize!(
                    mock,
                    MOI.OPTIMAL,
                    (MOI.FEASIBLE_POINT, [1.5]),
                )
            end,
            (mock::MOIU.MockOptimizer) -> begin
                MOI.set(mock, MOI.ObjectiveBound(), 4.0)
                MOIU.mock_optimize!(
                    mock,
                    MOI.OPTIMAL,
                    (MOI.FEASIBLE_POINT, [1.5]),
                )
            end,
        )
        MOIT.solve_objbound_edge_cases(mock, config)
    end

    @testset "raw_status_string" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> begin
                MOI.set(
                    mock,
                    MOI.RawStatusString(),
                    "Mock solution set by `mock_optimize!`.",
                )
                MOIU.mock_optimize!(
                    mock,
                    MOI.OPTIMAL,
                    (MOI.FEASIBLE_POINT, [0.0]),
                )
            end,
        )
        MOIT.raw_status_string(mock, config)
    end

    @testset "solve_time" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> begin
                MOI.set(mock, MOI.SolveTimeSec(), 0.0)
                MOIU.mock_optimize!(
                    mock,
                    MOI.OPTIMAL,
                    (MOI.FEASIBLE_POINT, [0.0]),
                )
            end,
        )
        MOIT.solve_time(mock, config)
    end

    @testset "solve_zero_one_with_bounds" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> begin
                MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [1.0]))
            end,
            (mock::MOIU.MockOptimizer) -> begin
                MOIU.mock_optimize!(mock, MOI.OPTIMAL, (MOI.FEASIBLE_POINT, [0.0]))
            end,
            (mock::MOIU.MockOptimizer) -> begin
                MOIU.mock_optimize!(mock, MOI.INFEASIBLE)
            end,
        )
        MOIT.solve_zero_one_with_bounds_1(mock, config)
        MOIT.solve_zero_one_with_bounds_2(mock, config)
        MOIT.solve_zero_one_with_bounds_3(mock, config)
    end

    @testset "solve_start_soc" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
                MOI.FEASIBLE_POINT,
                (MOI.VectorAffineFunction{Float64}, MOI.SecondOrderCone) =>
                    [[1.0, -1.0]],
            ),
        )
        MOIT.solve_start_soc(mock, config)
    end

    @testset "solve_unbounded_model" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> begin
                MOIU.mock_optimize!(mock, MOI.DUAL_INFEASIBLE)
            end,
        )
        MOIT.solve_unbounded_model(mock, config)
    end

    @testset "solve_single_variable_dual_min" begin
        flag = mock.eval_variable_constraint_dual
        mock.eval_variable_constraint_dual = false
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
                MOI.FEASIBLE_POINT,
                (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [1.0],
                (MOI.VariableIndex, MOI.LessThan{Float64}) => [0.0],
            ),
        )
        MOIT.solve_single_variable_dual_min(mock, config)
        mock.eval_variable_constraint_dual = flag
    end

    @testset "solve_single_variable_dual_max" begin
        flag = mock.eval_variable_constraint_dual
        mock.eval_variable_constraint_dual = false
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
                MOI.FEASIBLE_POINT,
                (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [0.0],
                (MOI.VariableIndex, MOI.LessThan{Float64}) => [-1.0],
            ),
        )
        MOIT.solve_single_variable_dual_max(mock, config)
        mock.eval_variable_constraint_dual = flag
    end

    @testset "solve_result_index" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
                MOI.FEASIBLE_POINT,
                (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [1.0],
            ),
        )
        MOIT.solve_result_index(mock, config)
    end

    @testset "solve_farkas_equalto_upper" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.INFEASIBLE,
                (MOI.NO_SOLUTION, [NaN, NaN]),
                MOI.INFEASIBILITY_CERTIFICATE,
                (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [2.0, 1.0],
                (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                    [-1.0],
            ),
        )
        MOIT.solve_farkas_equalto_upper(mock, config)
    end

    @testset "solve_farkas_equalto_lower" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.INFEASIBLE,
                (MOI.NO_SOLUTION, [NaN, NaN]),
                MOI.INFEASIBILITY_CERTIFICATE,
                (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [2.0, 1.0],
                (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) =>
                    [1.0],
            ),
        )
        MOIT.solve_farkas_equalto_lower(mock, config)
    end

    @testset "solve_farkas_lessthan" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.INFEASIBLE,
                (MOI.NO_SOLUTION, [NaN, NaN]),
                MOI.INFEASIBILITY_CERTIFICATE,
                (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [2.0, 1.0],
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-1.0],
            ),
        )
        MOIT.solve_farkas_lessthan(mock, config)
    end

    @testset "solve_farkas_greaterthan" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.INFEASIBLE,
                (MOI.NO_SOLUTION, [NaN, NaN]),
                MOI.INFEASIBILITY_CERTIFICATE,
                (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [2.0, 1.0],
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [1.0],
            ),
        )
        MOIT.solve_farkas_greaterthan(mock, config)
    end

    @testset "solve_farkas_interval_upper" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.INFEASIBLE,
                (MOI.NO_SOLUTION, [NaN, NaN]),
                MOI.INFEASIBILITY_CERTIFICATE,
                (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [2.0, 1.0],
                (MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}) => [-1.0],
            ),
        )
        MOIT.solve_farkas_interval_upper(mock, config)
    end

    @testset "solve_farkas_interval_lower" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.INFEASIBLE,
                (MOI.NO_SOLUTION, [NaN, NaN]),
                MOI.INFEASIBILITY_CERTIFICATE,
                (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [2.0, 1.0],
                (MOI.ScalarAffineFunction{Float64}, MOI.Interval{Float64}) => [1.0],
            ),
        )
        MOIT.solve_farkas_interval_lower(mock, config)
    end

    @testset "solve_farkas_variable_lessthan" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.INFEASIBLE,
                (MOI.NO_SOLUTION, [NaN, NaN]),
                MOI.INFEASIBILITY_CERTIFICATE,
                (MOI.VariableIndex, MOI.LessThan{Float64}) => [-2.0, -1.0],
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [1.0],
            ),
        )
        MOIT.solve_farkas_variable_lessthan(mock, config)
    end

    @testset "solve_farkas_variable_lessthan_max" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.INFEASIBLE,
                (MOI.NO_SOLUTION, [NaN, NaN]),
                MOI.INFEASIBILITY_CERTIFICATE,
                (MOI.VariableIndex, MOI.LessThan{Float64}) => [-2.0, -1.0],
                (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => [1.0],
            ),
        )
        MOIT.solve_farkas_variable_lessthan_max(mock, config)
    end

    @testset "solve_twice" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
            ),
        )
        MOIT.solve_twice(mock, config)
    end
end

@testset "modifications" begin
    # `UniversalFallback` needed for `MOI.Silent`
    mock = MOIU.MockOptimizer(MOIU.UniversalFallback(MOIU.Model{Float64}()))
    config = MOIT.Config()
    @testset "delete_variables_in_a_batch" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0, 1.0, 1.0]),
                MOI.FEASIBLE_POINT,
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
                MOI.FEASIBLE_POINT,
            ),
        )
        MOIT.delete_variables_in_a_batch(mock, config)
    end
    @testset "solve_set_singlevariable_lessthan" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
                MOI.FEASIBLE_POINT,
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [2.0]),
                MOI.FEASIBLE_POINT,
            ),
        )
        MOIT.solve_set_singlevariable_lessthan(mock, config)
    end
    @testset "solve_transform_singlevariable_lessthan" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
                MOI.FEASIBLE_POINT,
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [2.0]),
                MOI.FEASIBLE_POINT,
            ),
        )
        MOIT.solve_transform_singlevariable_lessthan(mock, config)
    end
    @testset "solve_set_scalaraffine_lessthan" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
                MOI.FEASIBLE_POINT,
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-1.0],
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [2.0]),
                MOI.FEASIBLE_POINT,
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-1.0],
            ),
        )
        MOIT.solve_set_scalaraffine_lessthan(mock, config)
    end
    @testset "solve_coef_scalaraffine_lessthan" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
                MOI.FEASIBLE_POINT,
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-1.0],
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.5]),
                MOI.FEASIBLE_POINT,
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-0.5],
            ),
        )
        MOIT.solve_coef_scalaraffine_lessthan(mock, config)
    end
    @testset "solve_func_scalaraffine_lessthan" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
                MOI.FEASIBLE_POINT,
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-1.0],
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.5]),
                MOI.FEASIBLE_POINT,
                (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}) => [-0.5],
            ),
        )
        MOIT.solve_func_scalaraffine_lessthan(mock, config)
    end
    @testset "solve_const_vectoraffine_nonpos" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.0, 0.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0, 0.75]),
            ),
        )
        MOIT.solve_const_vectoraffine_nonpos(mock, config)
    end
    @testset "solve_multirow_vectoraffine_nonpos" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.5]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [0.25]),
            ),
        )
        MOIT.solve_multirow_vectoraffine_nonpos(mock, config)
    end
    @testset "solve_const_scalar_objective" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
            ),
        )
        MOIT.solve_const_scalar_objective(mock, config)
    end
    @testset "solve_coef_scalar_objective" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
            ),
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
            ),
        )
        MOIT.solve_coef_scalar_objective(mock, config)
    end
    @testset "delete_variable_with_single_variable_obj" begin
        MOIU.set_mock_optimize!(
            mock,
            (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
                mock,
                MOI.OPTIMAL,
                (MOI.FEASIBLE_POINT, [1.0]),
                MOI.FEASIBLE_POINT,
                (MOI.VariableIndex, MOI.GreaterThan{Float64}) => [1.0],
            ),
        )
        MOIT.delete_variable_with_single_variable_obj(mock, config)
    end
end