text:
function viewVort2D(tev::Matrix{Float64}, lev::Matrix{Float64}, bv::Matrix{Float64})
    scatter(tev[:,2], tev[:,3], s=5, c=tev[:,1], edgecolors="none")
    sc = scatter(lev[:,2], lev[:,3], s=5, c=lev[:,1], edgecolors="none")
    sc2 = scatter(bv[:,2], bv[:,3], s=5, c=bv[:,1], edgecolors="none")
    #plot(bv[:,2], bv[:,3], color = "black", linewidth=1.0)
end

#For quick views while debugging
function viewVort2D(tev::Vector{TwoDVort}, lev::Vector{TwoDVort}, bv::Vector{TwoDVort})
    tev = hcat(map(q->q.s, tev), map(q->q.x, tev), map(q->q.z, tev))
    lev = hcat(map(q->q.s, lev), map(q->q.x, lev), map(q->q.z, lev))
    bv = hcat(map(q->q.s, bv), map(q->q.x, bv), map(q->q.z, bv))

    scatter(tev[:,2], tev[:,3], s=5, c=tev[:,1], edgecolors="none")
    sc = scatter(lev[:,2], lev[:,3], s=5, c=lev[:,1], edgecolors="none")
    sc2 = scatter(bv[:,2], bv[:,3], s=5, c=bv[:,1], edgecolors="none")
    #plot(bv[:,2], bv[:,3], color = "black", linewidth=1.0)
end

function viewVortConnect2D(tev::Matrix{Float64}, lev::Matrix{Float64}, bv::Matrix{Float64})
    scatter(tev[:,2], tev[:,3], s=5, c=tev[:,1], edgecolors="none")
    sc = scatter(lev[:,2], lev[:,3], s=5, c=lev[:,1], edgecolors="none")
    sc2 = scatter(bv[:,2], bv[:,3], s=5, c=bv[:,1], edgecolors="none")

    plot(tev[:,2], tev[:,3], color = "red", linewidth=1.0)
    plot(lev[:,2], lev[:,3], color = "blue", linewidth=1.0)
    plot(bv[:,2], bv[:,3], color = "black", linewidth=1.0)
end

function makeVortPlots2D(tt = "all")
 # tt is a target time to make the plot for. It will choose the closest time available
    try
        cd("Step Files")
    catch
        println(" Error: No 'Step Files' directory. Make sure 'wflag' is enabled before running simulation.")
        return
    end

    dirvec = readdir()
    dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
    if tt != "all"
        figure()
        ttidx = findmin(abs.(dirresults.-tt))[2]
    end
    #Determine axis limits
    dirmax = maximum(dirresults)
    cd("$(dirmax)")

    multsurfflag = 0
    if isfile("boundv-1") == true  #only 1 surface
        multsurfflag = 1
    end

    if multsurfflag == 0 #single surface

        tev = DelimitedFiles.readdlm("tev", '\t', Float64,  skipstart = 1)
        lev =   try
            DelimitedFiles.readdlm("lev", '\t', Float64,  skipstart = 1)
        catch
            zeros(0,3)
        end
        bv = DelimitedFiles.readdlm("boundv", '\t', Float64,  skipstart = 1)
        #tupletest = minimum([UnsteadyFlowSolvers.tev[:,2];UnsteadyFlowSolvers.lev[:,2];UnsteadyFlowSolvers.bv[:,2]])
        xmin = minimum([tev[:,2];lev[:,2];bv[:,2]])
        zmin = minimum([tev[:,3];lev[:,3];bv[:,3]])
        xmax = maximum([tev[:,2];lev[:,2]])
        zmax = maximum([bv[:,2];tev[:,3];lev[:,3];bv[:,3]])

        cd("../..")

        if tt == "all"
            dirvec = readdir()
            if "vortPlots" in dirvec
                rm("vortPlots", recursive=true)
            end
            mkdir("vortPlots")
        end
        cd("Step Files")

        for i=1:length(dirresults)
            if tt != "all" && ttidx != i
                continue
            end
            if dirresults[i] != 0
                dirstr="$(dirresults[i])"
                cd(dirstr)
                tev= DelimitedFiles.readdlm("tev", '\t', Float64,  skipstart = 1)
                lev =   try
                    DelimitedFiles.readdlm("lev", '\t', Float64,  skipstart = 1)
                catch
                    zeros(0,3)
                end
                bv= DelimitedFiles.readdlm("boundv", '\t', Float64,  skipstart = 1)

                viewVort2D(tev, lev, bv)
                axis([xmin-1, xmax+1, zmin-1, zmax+1])

                if tt == "all"
                    savefig("../../vortPlots/$(dirresults[i]).png")
                    close()
                end
                cd("..")
            end
        end
        if tt == "all"
            println(" All plots saved under the 'vortPlots' directory.")
        else
            pause(0.0001)
            show(block = false)
        end

    else #multiple surfaces

        tev= DelimitedFiles.readdlm("tev", '\t', Float64,  skipstart = 1)
        lev =   try
            DelimitedFiles.readdlm("lev", '\t', Float64,  skipstart = 1)
        catch
            zeros(0,3)
        end

        nsurf = 0
        for i = 1:1000
            if isfile("boundv-$i") == true
                nsurf += 1

            end
        end
        bv= DelimitedFiles.readdlm("boundv-1", '\t', Float64,  skipstart = 1)
        for i = 2:nsurf
            bvt = DelimitedFiles.readdlm("boundv-$i", '\t', Float64,  skipstart = 1)
            bv = [bv;bvt;]
        end

        xmin = minimum([tev[:,2];lev[:,2];bv[:,2]])
        zmin = minimum([tev[:,3];lev[:,3];bv[:,3]])
        xmax = maximum([tev[:,2];lev[:,2]])
        zmax = maximum([bv[:,2];tev[:,3];lev[:,3];bv[:,3]])

        cd("..")
    if tt == "all"
            dirvec = readdir()
        if "vortPlots" in dirvec
            rm("vortPlots", recursive=true)
        end

        cd("..")
        mkpath("vortPlots")
        cd("Step Files")
        for i=1:length(dirresults)
            if dirresults[i] != 0
                dirstr="$(dirresults[i])"
                cd(dirstr)
                tev = DelimitedFiles.readdlm("tev", '\t', Float64,  skipstart = 1)
                lev =   try
                    DelimitedFiles.readdlm("lev", '\t', Float64,  skipstart = 1)
                catch
                    zeros(0,3)
                end
                bv= DelimitedFiles.readdlm("boundv-1", '\t', Float64,  skipstart = 1)
                for i = 2:nsurf
                    bvt = DelimitedFiles.readdlm("boundv-$i", '\t', Float64,  skipstart = 1)
                    bv = [bv;bvt;]
                end
                viewVort2D(tev, lev, bv)
                axis([xmin-1, xmax+1, zmin-1, zmax+1])
                savefig("../../vortPlots/$(dirresults[i]).png")
                close()
                cd("..")
            end
        end
    end
end
    cd("..")

function makeForcePlots2D(mat = "")

    dirvec = readdir()
    if "forcePlots" in dirvec
        rm("forcePlots", recursive=true)
    end
    mkdir("forcePlots")

    if mat == ""
       mat = DelimitedFiles.readdlm("resultsSummary", Float64,  skipstart = 1)
    end

    if length(mat[1,:])==8 #only 1 surface
        t = mat[:,1]
        alpha = mat[:,2]*180/pi
        h = mat[:,3]
        cl = mat[:,6]
        cd = mat[:,7]
        cm = mat[:,8]

        plot(t, alpha)
        len = length(t)
        xmin = t[1]
        xmax = t[end]
        zmin = minimum(alpha) - 0.1*abs(minimum(alpha))
        zmax = maximum(alpha) + 0.1*abs(maximum(alpha))
        axis([xmin, xmax, zmin, zmax])
        xlabel("t*")
        ylabel(L"$\alpha$ (deg)")
        savefig("forcePlots/pitch.png")
        close()

        plot(t, h)
        xmin = t[1]
        xmax = t[end]
        zmin = minimum(h) - 0.1*abs(minimum(h))
        zmax = maximum(h) + 0.1*abs(maximum(h))
        axis([xmin, xmax, zmin, zmax])
        xlabel("t*")
        ylabel(L"$h/c$")
        savefig("forcePlots/plunge.png")
        close()

        plot(t, cl)
        range = Int(round(0.05*len)):Int(round(0.95*len))
        xmin = t[1]
        xmax = t[end]
        zmin = minimum(cl[range]) - 0.1*abs(minimum(cl[range]))
        zmax = maximum(cl[range]) + 0.1*abs(maximum(cl[range]))
        axis([xmin, xmax, zmin, zmax])
        xlabel("t*")
        ylabel(L"$C_l$")
        savefig("forcePlots/cl.png")
        close()

        plot(t, cd)
        range = Int(round(0.05*len)):Int(round(0.95*len))
        xmin = t[1]
        xmax = t[end]
        zmin = minimum(cd[range]) - 0.1*abs(minimum(cd[range]))
        zmax = maximum(cd[range]) + 0.1*abs(maximum(cd[range]))
        axis([xmin, xmax, zmin, zmax])
        xlabel("t*")
        ylabel(L"$C_d$")
        savefig("forcePlots/cd.png")
        close()

        plot(t, cm)
        range = Int(round(0.05*len)):Int(round(0.95*len))
        xmin = t[1]
        xmax = t[end]
        zmin = minimum(cm[range]) - 0.1*abs(minimum(cm[range]))
        zmax = maximum(cm[range]) + 0.1*abs(maximum(cm[range]))
        axis([xmin, xmax, zmin, zmax])
        xlabel("t*")
        ylabel(L"$C_m$")
        savefig("forcePlots/cm.png")
        close()

    else #multiple surfaces

        nsurf = (length(mat[1,:]) - 1)/7
        t = mat[:,1]

        for i = 1:nsurf
            alpha = mat[:,Int((i-1)*7+2)]*180/pi
            h = mat[:,Int((i-1)*7+3)]
            cl = mat[:,Int((i-1)*7+6)]
            cd = mat[:,Int((i-1)*7+7)]
            cm = mat[:,Int((i-1)*7+8)]

            plot(t, alpha)
            len = length(t)
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(alpha) - 0.1*abs(minimum(alpha))
            zmax = maximum(alpha) + 0.1*abs(maximum(alpha))
            axis([xmin, xmax, zmin, zmax])
            xlabel("t*")
            ylabel(L"$\alpha$ (deg)")
            savefig("forcePlots/pitch-$i.png")
            close()

            plot(t, h)
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(h) - 0.1*abs(minimum(h))
            zmax = maximum(h) + 0.1*abs(maximum(h))
            axis([xmin, xmax, zmin, zmax])
            xlabel("t*")
            ylabel(L"$h/c$")
            savefig("forcePlots/plunge-$i.png")
            close()

            plot(t, cl)
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(cl[range]) - 0.1*abs(minimum(cl[range]))
            zmax = maximum(cl[range]) + 0.1*abs(maximum(cl[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel("t*")
            ylabel(L"$C_l$")
            savefig("forcePlots/cl-$i.png")
            close()

            plot(t, cd)
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(cd[range]) - 0.1*abs(minimum(cd[range]))
            zmax = maximum(cd[range]) + 0.1*abs(maximum(cd[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel("t*")
            ylabel(L"$C_d$")
            savefig("forcePlots/cd-$i.png")
            close()

            plot(t, cm)
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(cm[range]) - 0.1*abs(minimum(cm[range]))
            zmax = maximum(cm[range]) + 0.1*abs(maximum(cm[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel("t*")
            ylabel(L"$C_m$")
            savefig("forcePlots/cm-$i.png")
            close()
        end

    end

end

function checkConverge(k::Float64)

    mat = DelimitedFiles.readdlm("resultsSummary", '\t', Float64,  skipstart = 1)

    T = pi/k

    end_cycle = mat[end,1]

    #Find number of cycles
    ncyc = 0
    for i = 1:1000
        t_l = real(i)*T
        if t_l > end_cycle
            break
        end
        ncyc = ncyc + 1
    end

    #Lift convergence
    figure
    for i = 1:ncyc
        start_t = real(i-1)*T
        end_t = real(i)*T
        start_ind = argmin(abs.(mat[:,1] .- start_t))
        end_ind = argmin(abs.(mat[:,1] .- end_t))
        plot((mat[start_ind:end_ind,1] .- start_t)/T, mat[start_ind:end_ind,6])
    end
    xmin = 0.
    xmax = 1.
    zmin = minimum(mat[:,6])
    zmax = maximum(mat[:,6])
    axis([xmin, xmax, zmin, zmax])

    xlabel("t*")
    ylabel(L"$C_l$")
    savefig("forcePlots/cl-convergence.png")
    close()

    #Drag convergence
    figure
    for i = 1:ncyc
        start_t = real(i-1)*T
        end_t = real(i)*T
        start_ind = argmin(abs.(mat[:,1] .- start_t))
        end_ind = argmin(abs.(mat[:,1] .- end_t))
        plot((mat[start_ind:end_ind,1] .- start_t)/T, mat[start_ind:end_ind,7])
    end
    xmin = 0.
    xmax = 1.
    zmin = minimum(mat[:,7])
    zmax = maximum(mat[:,7])
    axis([xmin, xmax, zmin, zmax])

    xlabel("t*")
    ylabel(L"$C_d$")
    savefig("forcePlots/cd-convergence.png")
    close()

    #Pitching moment convergence
    figure
    for i = 1:ncyc
        start_t = real(i-1)*T
        end_t = real(i)*T
        start_ind = argmin(abs.(mat[:,1] .- start_t))
        end_ind = argmin(abs.(mat[:,1] .- end_t))
        plot((mat[start_ind:end_ind,1] .- start_t)/T, mat[start_ind:end_ind,8])
    end
    xmin = 0.
    xmax = 1.
    zmin = minimum(mat[:,8])
    zmax = maximum(mat[:,8])
    axis([xmin, xmax, zmin, zmax])

    xlabel("t*")
    ylabel(L"$C_m$")
    savefig("forcePlots/cm-convergence.png")
    close()
end

function makeKinemClVortPlots2D(mat = "", tt = "all")
# tt is a target time to make the plot for. It will choose the closest time available

    if mat == ""
        mat = DelimitedFiles.readdlm("resultsSummary", '\t', Float64,  skipstart = 1)
    end

    t = mat[:,1]
    len = length(t)
    cl = mat[:,6]
    alpha = mat[:,2]*180/pi
    h = mat[:,3]
    u = mat[:,4]

    try
        cd("Step Files")
    catch
        println(" Error: No 'Step Files' directory. Make sure wflag is enabled.")
        return
    end

    dirvec = readdir()
    dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
    if tt != "all"
        ttidx = findmin(abs.(dirresults.-tt))[2]
    end

    #Determine axis limits
    dirmax = maximum(dirresults)

    cd("$(dirmax)")

    multsurfflag = 0
    if isfile("boundv-1") == true  #only 1 surface
        error("this plot function is only written for single surface")
    end

    tev = DelimitedFiles.readdlm("tev", '\t', Float64,  skipstart = 1)
    lev =   try
        DelimitedFiles.readdlm("lev", '\t', Float64,  skipstart = 1)[1]
    catch
        zeros(0,3)
    end
    bv = DelimitedFiles.readdlm("boundv", '\t', Float64,  skipstart = 1)

    xminv = minimum([tev[:,2];lev[:,2];bv[:,2]])
    zminv = minimum([tev[:,3];lev[:,3];bv[:,3]])
    xmaxv = maximum([tev[:,2];lev[:,2]])
    zmaxv = maximum([bv[:,2];tev[:,3];lev[:,3];bv[:,3]])

    cd("../..")
    if tt == "all"
        dirvec = readdir()
        if "infoPlots" in dirvec
            rm("infoPlots", recursive=true)
        end
        mkdir("infoPlots")
    end

    cd("Step Files")
    for i=1:length(dirresults)
        if tt != "all" && ttidx != i
            continue
        end
        if dirresults[i] != 0

            dirstr="$(dirresults[i])"
            cd(dirstr)

            t_cur = dirresults[i]

            figure(figsize=(8,8))

            subplot2grid((6, 2), (0, 0))
            plot(t, alpha)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(alpha[range]) - 0.1*abs(minimum(alpha[range]))
            zmax = maximum(alpha[range]) + 0.1*abs(maximum(alpha[range]))
            axis([xmin, xmax, zmin, zmax])
            #xlabel("t*")
            xticks([])
            ylabel(L"$\alpha (deg)$")

            subplot2grid((6, 2), (1, 0))
            plot(t, h)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(h[range]) - 0.1*abs(minimum(h[range]))
            zmax = maximum(h[range]) + 0.1*abs(maximum(h[range]))
            axis([xmin, xmax, zmin, zmax])
            #xlabel("t*")
            xticks([])
            ylabel(L"$h$")

            subplot2grid((6, 2), (2, 0))
            plot(t, u)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(u[range]) - 0.1*abs(minimum(u[range]))
            zmax = maximum(u[range]) + 0.1*abs(maximum(u[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel("t*")
            ylabel(L"$u$")

            subplot2grid((6, 2), (0, 1), rowspan=3)
            plot(t, cl)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(cl[range]) - 0.1*abs(minimum(cl[range]))
            zmax = maximum(cl[range]) + 0.1*abs(maximum(cl[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel("t*")
            ylabel(L"$C_l$")

            subplot2grid((6, 2), (3, 0), rowspan=3, colspan=2)
            tev= DelimitedFiles.readdlm("tev", '\t', Float64,  skipstart = 1)
            lev =   try
                DelimitedFiles.readdlm("lev", '\t', Float64,  skipstart = 1)
            catch
                zeros(0,3)
            end
            bv = DelimitedFiles.readdlm("boundv", '\t', Float64,  skipstart = 1)
            viewVort2D(tev, lev, bv)
            axis([xminv-1, xmaxv+1, zminv-1, zmaxv+1])

            tight_layout()

            if tt == "all"
                savefig("../../infoPlots/$(dirresults[i]).png")
                close()
            end
            cd("..")
        end
    end
    cd("..")
    if tt == "all"
        println(" All plots saved under the 'infoPlots' directory.")
    else
        pause(0.0001)
        show(block = false)
    end
end

function makeKinemVelVortPlots2D(mat = "",tt = "all")

    if mat == ""
        mat = DelimitedFiles.readdlm("resultsSummary", '\t', Float64,  skipstart = 1)
    end

    t = mat[:,1]
    len = length(t)
    cl = mat[:,6]
    alpha = mat[:,2]*180/pi
    h = mat[:,3]
    u = mat[:,4]

    try
        cd("Step Files")
    catch
        println(" Error: No 'Step File' directory. Make sure wflag is enabled.")
        return
    end

    dirvec = readdir()
    dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
    if tt != "all"
        ttidx = findmin(abs(dirresults-tt))[2]
    end

    #Determine axis limits
    dirmax = maximum(dirresults)

    cd("$(dirmax)")

    multsurfflag = 0
    if isfile("boundv-1") == true  #only 1 surface
        error("this plot function is only written for single surface")
    end

    tev = DelimitedFiles.readdlm("tev", '\t', Float64,  skipstart = 1)
    lev =   try
        DelimitedFiles.readdlm("lev", '\t', Float64,  skipstart = 1)
    catch
        zeros(0,3)
    end
    bv = DelimitedFiles.readdlm("boundv", '\t', Float64,  skipstart = 1)

    xminv = minimum([tev[:,2];lev[:,2];bv[:,2]])
    zminv = minimum([tev[:,3];lev[:,3];bv[:,3]])
    xmaxv = maximum([tev[:,2];lev[:,2]])
    zmaxv = maximum([bv[:,2];tev[:,3];lev[:,3];bv[:,3]])

    cd("../..")

    dirvec = readdir()
    if "infoplots" in dirvec
        rm("infoplots", recursive=true)
    end
    mkdir("infoplots")

    cd("Step Files")
    for i=1:length(dirresults)
        if tt != "all" && ttidx != i
            continue
        end
        if dirresults[i] != 0
            dirstr="$(dirresults[i])"
            cd(dirstr)

            t_cur = dirresults[i]

            figure(figsize=(8,8))

            subplot2grid((6, 2), (0, 0))
            plot(t, alpha)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(alpha[range]) - 0.1*abs(minimum(alpha[range]))
            zmax = maximum(alpha[range]) + 0.1*abs(maximum(alpha[range]))
            axis([xmin, xmax, zmin, zmax])
            #xlabel("t*")
            xticks([])
            ylabel(L"$\alpha (deg)$")

            subplot2grid((6, 2), (1, 0))
            plot(t, h)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(h[range]) - 0.1*abs(minimum(h[range]))
            zmax = maximum(h[range]) + 0.1*abs(maximum(h[range]))
            axis([xmin, xmax, zmin, zmax])
            #xlabel("t*")
            xticks([])
            ylabel(L"$h$")

            subplot2grid((6, 2), (2, 0))
            plot(t, u)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(u[range]) - 0.1*abs(minimum(u[range]))
            zmax = maximum(u[range]) + 0.1*abs(maximum(u[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel("t*")
            ylabel(L"$u$")

            subplot2grid((6, 2), (0, 1), rowspan=3)
            plot(t, cl)
            plot([t_cur; t_cur], [-10000; 10000], "k-")
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(cl[range]) - 0.1*abs(minimum(cl[range]))
            zmax = maximum(cl[range]) + 0.1*abs(maximum(cl[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel("t*")
            ylabel(L"$C_l$")

            subplot2grid((6, 2), (3, 0), rowspan=3, colspan=2)
            tev = DelimitedFiles.readdlm("tev", '\t', Float64,  skipstart = 1)
            lev =   try
                DelimitedFiles.readdlm("lev", '\t', Float64,  skipstart = 1)
            catch
                zeros(0,3)
            end
            bv = DelimitedFiles.readdlm("boundv", '\t', Float64,  skipstart = 1)
            viewVort2D(tev, lev, bv)
            axis([xminv-1, xmaxv+1, zminv-1, zmaxv+1])

            tight_layout()

            savefig("../../infoplots/$(dirresults[i]).png")
            close()
            cd("..")
        end
    end
    cd("..")
end

function subPlotForces(mat, xVar :: Int64, xVarLabel :: String, saveName :: String)
    # Create Cl,Cd,Cm,LESP vs xVar plot
    # mat = [t , alpha , h , u , LESP , cl , cd , cm] for each timestep
    xVar = mat[:,xVar] ;
    cl = mat[:,6] ;
    cd = mat[:,7] ;
    cm = mat[:,8] ;
    LESP = mat[:,5] ;

    subplot(221)
    plot(xVar,cl) ;
    xlabel(xVarLabel)
    ylabel("Cl")

    subplot(222)
    plot(xVar,cd) ;
    xlabel(xVarLabel)
    ylabel("Cd")

    subplot(223)
    plot(xVar,cm) ;
    xlabel(xVarLabel)
    ylabel("Cm")

    subplot(224)
    plot(xVar,LESP) ;
    xlabel(xVarLabel)
    ylabel("LESP")

    tight_layout()
    savefig(saveName)
    close()
end

## UI Plots
function initPlot(A,C,t,a,h,u)
    # Create intitilization plots that shows airfoil and camber and motion parameters

    # Camber Plot
    if A != 0 # airfoil defined
        subplot(121,aspect = 1)
        plot(A[:,1],A[:,2],color = "black") # full airfoil
        plot(C[:,1],C[:,2], color = "blue", marker = "|") # Camber
        ylim(-.5,.5)
        xlabel("x")
        ylabel("z")
        title("Geometry")
    end

    # Motion plots
    if sum(a) != 0 # alphadef defined
        subplot(122)
        plot(t,h, color = "black") # height
        plot(t,u, color = "blue") # velocity
        plot(t,a, color = "red") # alpha
        xlabel("t*")
        legend(["Height","Velocity","Alpha [deg]"])
        title("Motion")
    end
    pause(0.0001)
    show(block = false)
end

function plotMat(mat,pltCnt,plotx,ploty,matIdx)
    # plot subplots (number determined by plotCnt) with an x-axis of plotx
    #   and a y axis of ploty. plotx/y are indices of the columns of mat.

    figure()
    # find subplot input dimensions
    if pltCnt <= 2 # plot on one row
        subRow = 1
        subCol = pltCnt
    else
        subRow = 2
        subCol = ceil(pltCnt/2)
    end
    subDim = subRow*100 + subCol*10

    for i = 1:pltCnt
        subplot(subDim + i)
        x = matIdx[lowercase(plotx)]
        y = matIdx[lowercase(ploty[i])]
        plot(mat[:,x],mat[:,y],"b-")
        xlabel(plotx)
        ylabel(ploty[i])
        title("$(ploty[i]) vs $plotx")
    end
    tight_layout()
    pause(0.0001)
    show(block = false)
end


text:
export cellaverage, cellcenters, coordinates
export visualize

"""
function cellaverage(Q; M = nothing)

# Description
Compute the cell-average of Q given the mass matrix M.
Assumes that Q and M are the same size

# Arguments
`Q`: MPIStateArrays (array)

# Keyword Arguments
`M`: Mass Matrix (array)

# Return
The cell-average of Q
"""
function cellaverage(Q; M = nothing)
    if M==nothing
        return nothing
    end
    return (sum(M .* Q, dims = 1) ./ sum(M , dims = 1))[:]
end

"""
function coordinates(grid::DiscontinuousSpectralElementGrid)

# Description
Gets the (x,y,z) coordinates corresponding to the grid

# Arguments
- `grid`: DiscontinuousSpectralElementGrid

# Return
- `x, y, z`: views of x, y, z coordinates
"""
function coordinates(grid::DiscontinuousSpectralElementGrid)
    x = view(grid.vgeo, :, grid.x1id, :)   # x-direction
    y = view(grid.vgeo, :, grid.x2id, :)   # y-direction
    z = view(grid.vgeo, :, grid.x3id, :)   # z-direction
    return x, y, z
end

"""
function cellcenters(Q; M = nothing)

# Description
Get the cell-centers of every element in the grid

# Arguments
- `grid`: DiscontinuousSpectralElementGrid

# Return
- Tuple of cell-centers
"""
function cellcenters(grid::DiscontinuousSpectralElementGrid)
    x, y, z = coordinates(grid)
    M = view(grid.vgeo, :, grid.Mid, :)  # mass matrix
    xC = cellaverage(x, M = M)
    yC = cellaverage(y, M = M)
    zC = cellaverage(z, M = M)
    return xC[:], yC[:], zC[:]
end

"""
function visualize(g::DiscontinuousSpectralElementGrid)

# Description
- Use Makie to visualize a 3D grid

# Arguments
- `g`: A DiscontinuousSpectralElementGrid
"""
function visualize(g::DiscontinuousSpectralElementGrid)
    x, y, z = coordinates(g)

    e = collect(1:size(x)[2])
    nfaces = 6
    faces = collect(1:nfaces)
    opacities = collect(range(0,1, length = 10))

    scene, layout = layoutscene()

    element = Node{Any}(e[1])
    face = Node{Any}(faces[6])
    opacity = Node{Any}(opacities[1])

    tmpx = @lift(x[:, $element])
    tmpy = @lift(y[:, $element])
    tmpz = @lift(z[:, $element])

    tmpxf = @lift(x[g.vmap⁻[:, $face, $element]])
    tmpyf = @lift(y[g.vmap⁻[:, $face, $element]])
    tmpzf = @lift(z[g.vmap⁻[:, $face, $element]])

    total_color = @lift(RGBA(0,0,0, $opacity))

    lscene = layout[1:4, 2:4] = LScene(scene)
    Makie.scatter!(lscene, x[:], y[:], z[:], color = total_color, markersize = 100.0, strokewidth = 0)
    Makie.scatter!(lscene, tmpx, tmpy, tmpz, color = RGBA(0,0,0,0.5), markersize = 100.0, strokewidth = 0, camera = cam3d!)
    Makie.scatter!(lscene, tmpxf, tmpyf, tmpzf, color = RGBA(1,0,0,1.0), markersize = 100.0, strokewidth = 0, camera = cam3d!)
    supertitle = layout[1,2] = LText(scene, " "^10 * " Gauss-Lobatto Points " * " "^10, textsize = 50, color = :black)

    menu  = LMenu(scene, options = zip(e,e))
    menu2 = LMenu(scene, options = zip(faces,faces))
    menu3 = LMenu(scene, options = zip(opacities,opacities))
    layout[1, 1] = vgrid!(
        LText(scene, "Element", width = nothing),
        menu,
        LText(scene, "Face", width = nothing),
        menu2,
        LText(scene, "Opacity", width = nothing),
        menu3,
    )
    on(menu.selection) do s
        element[] = s
    end
    on(menu2.selection) do s
        face[] = s
    end
    on(menu3.selection) do s
        opacity[] = s
    end
    display(scene)
    return scene
end

ext: jl

max_stars_repo_name: CliMA/ImperoHooks

max_issues_repo_name: CliMA/ImperoHooks

max_forks_repo_name: CliMA/ImperoHooks


text:
export MongoCursor,
       query, fields, limit, skip

type MongoCursor
    client::MongoClient
    namespace::String
    query::BSONObject
    fields::BSONObject
    limit::Int
    skip::Int
    _cursor::Ptr{Void}

    function MongoCursor(client::MongoClient, namespace::String, q, f, l, s)
        p = ccall((:mongo_cursor_alloc, MONGO_LIB), Ptr{Void}, ())
        if p == C_NULL
            error("Unable to create mongo cursor – mongo_cursor_alloc() failed")
        end

        ccall((:mongo_cursor_init, MONGO_LIB),
            Void,
            (Ptr{Void}, Ptr{Void}, Ptr{Uint8}),
            p, client._mongo, bytestring(namespace))

        cursor = new(client, namespace, q, f, l, s, p)
        finalizer(cursor, destroy)

        query(cursor, q)
        fields(cursor, f)
        limit(cursor, l)
        skip(cursor, s)

        return cursor
    end
end

MongoCursor(client::MongoClient, namespace::String) = MongoCursor(client, namespace, BSONObject(), BSONObject(), 0, 0)

query(cursor::MongoCursor, query::Dict) = query(cursor, BSONObject(query))
function query(cursor::MongoCursor, query::BSONObject)
    cursor.query = query
    ccall((:mongo_cursor_set_query, MONGO_LIB), Void, (Ptr{Void}, Ptr{Void}), cursor._cursor, query._bson)
end

fields(cursor::MongoCursor, fields::Dict) = query(cursor, BSONObject(fields))
function fields(cursor::MongoCursor, fields::BSONObject)
    cursor.fields = fields
    ccall((:mongo_cursor_set_fields, MONGO_LIB), Void, (Ptr{Void}, Ptr{Void}), cursor._cursor, fields._bson)
end

function limit(cursor::MongoCursor, limit::Int)
    cursor.limit = limit
    ccall((:mongo_cursor_set_limit, MONGO_LIB), Void, (Ptr{Void}, Int32), cursor._cursor, limit)
end

function skip(cursor::MongoCursor, skip::Int)
    cursor.skip = skip
    ccall((:mongo_cursor_set_skip, MONGO_LIB), Void, (Ptr{Void}, Int32), cursor._cursor, skip)
end

## Iterator ##
import Base.start, Base.next, Base.done

start(c::MongoCursor) = ccall((:mongo_cursor_next, MONGO_LIB), Int32, (Ptr{Void},), c._cursor)
done(c::MongoCursor, errno::Int32) = begin
    if errno != MONGO_OK
        reset(c)
        return true
    end
    return false
end

next(c::MongoCursor, errno::Int32) = begin
    _current_bson = ccall((:mongo_cursor_bson, MONGO_LIB), Ptr{Void}, (Ptr{Void},), c._cursor)

    # Create a copy of the bson
    _bson = ccall((:bson_alloc, MONGO_LIB), Ptr{Void}, ())
    errno = ccall((:bson_copy, MONGO_LIB), Int32, (Ptr{Void}, Ptr{Void}), _bson, _current_bson)
    if errno == BSON_ERROR
      error("Could not copy BSON object")
    end

    errno = ccall((:mongo_cursor_next, MONGO_LIB), Int32, (Ptr{Void},), c._cursor)

    (dict(BSONObject(_bson)), errno)
end

## Private Methods ##

function destroy(cursor::MongoCursor)
    errno = ccall((:mongo_cursor_destroy, MONGO_LIB), Int32, (Ptr{Void},), cursor._cursor)
    if errno == MONGO_ERROR
        error("Unable to destroy mongo cursor – mongo_cursor_destroy() failed")
    end
end

function reset(cursor::MongoCursor)
    ccall((:mongo_cursor_init, MONGO_LIB), Void,
      (Ptr{Void}, Ptr{Void}, Ptr{Uint8}),
      cursor._cursor, cursor.client._mongo, bytestring(cursor.namespace))
end


text:
import DielectricUpdateTechique as DUTCH
using Test

a = rand(9);
b = rand(9);
c = zeros(9);

scalar = rand(1)[1]

DUTCH.caxpy!(c, scalar, a, b);
@test c == scalar.*a + b


text:
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/ida/ida.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function IDACreate()
    ccall((:IDACreate, libsundials_ida), IDAMemPtr, ())
end

function IDASetErrHandlerFn(ida_mem, ehfun::IDAErrHandlerFn, eh_data)
    ccall((:IDASetErrHandlerFn, libsundials_ida), Cint, (IDAMemPtr, IDAErrHandlerFn, Ptr{Void}), ida_mem, ehfun, eh_data)
end

function IDASetErrFile(ida_mem, errfp)
    ccall((:IDASetErrFile, libsundials_ida), Cint, (IDAMemPtr, Ptr{FILE}), ida_mem, errfp)
end

function IDASetUserData(ida_mem, user_data)
    ccall((:IDASetUserData, libsundials_ida), Cint, (IDAMemPtr, Any), ida_mem, user_data)
end

function __IDASetMaxOrd(ida_mem, maxord::Cint)
    ccall((:IDASetMaxOrd, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxord)
end

function IDASetMaxOrd(ida_mem, maxord)
    __IDASetMaxOrd(ida_mem, convert(Cint, maxord))
end

function __IDASetMaxNumSteps(ida_mem, mxsteps::Clong)
    ccall((:IDASetMaxNumSteps, libsundials_ida), Cint, (IDAMemPtr, Clong), ida_mem, mxsteps)
end

function IDASetMaxNumSteps(ida_mem, mxsteps)
    __IDASetMaxNumSteps(ida_mem, convert(Clong, mxsteps))
end

function IDASetInitStep(ida_mem, hin::realtype)
    ccall((:IDASetInitStep, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, hin)
end

function IDASetMaxStep(ida_mem, hmax::realtype)
    ccall((:IDASetMaxStep, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, hmax)
end

function IDASetStopTime(ida_mem, tstop::realtype)
    ccall((:IDASetStopTime, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, tstop)
end

function IDASetNonlinConvCoef(ida_mem, epcon::realtype)
    ccall((:IDASetNonlinConvCoef, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, epcon)
end

function __IDASetMaxErrTestFails(ida_mem, maxnef::Cint)
    ccall((:IDASetMaxErrTestFails, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxnef)
end

function IDASetMaxErrTestFails(ida_mem, maxnef)
    __IDASetMaxErrTestFails(ida_mem, convert(Cint, maxnef))
end

function __IDASetMaxNonlinIters(ida_mem, maxcor::Cint)
    ccall((:IDASetMaxNonlinIters, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxcor)
end

function IDASetMaxNonlinIters(ida_mem, maxcor)
    __IDASetMaxNonlinIters(ida_mem, convert(Cint, maxcor))
end

function __IDASetMaxConvFails(ida_mem, maxncf::Cint)
    ccall((:IDASetMaxConvFails, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxncf)
end

function IDASetMaxConvFails(ida_mem, maxncf)
    __IDASetMaxConvFails(ida_mem, convert(Cint, maxncf))
end

function __IDASetSuppressAlg(ida_mem, suppressalg::Cint)
    ccall((:IDASetSuppressAlg, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, suppressalg)
end

function IDASetSuppressAlg(ida_mem, suppressalg)
    __IDASetSuppressAlg(ida_mem, convert(Cint, suppressalg))
end

function __IDASetId(ida_mem, id::N_Vector)
    ccall((:IDASetId, libsundials_ida), Cint, (IDAMemPtr, N_Vector), ida_mem, id)
end

function IDASetId(ida_mem, id)
    __id = convert(NVector, id)
    __IDASetId(ida_mem, convert(N_Vector, __id))
end

function __IDASetConstraints(ida_mem, constraints::N_Vector)
    ccall((:IDASetConstraints, libsundials_ida), Cint, (IDAMemPtr, N_Vector), ida_mem, constraints)
end

function IDASetConstraints(ida_mem, constraints)
    __constraints = convert(NVector, constraints)
    __IDASetConstraints(ida_mem, convert(N_Vector, __constraints))
end

function IDASetRootDirection(ida_mem, rootdir)
    ccall((:IDASetRootDirection, libsundials_ida), Cint, (IDAMemPtr, Ptr{Cint}), ida_mem, rootdir)
end

function IDASetNoInactiveRootWarn(ida_mem)
    ccall((:IDASetNoInactiveRootWarn, libsundials_ida), Cint, (IDAMemPtr,), ida_mem)
end

function __IDAInit(ida_mem, res::IDAResFn, t0::realtype, yy0::N_Vector, yp0::N_Vector)
    ccall((:IDAInit, libsundials_ida), Cint, (IDAMemPtr, IDAResFn, realtype, N_Vector, N_Vector), ida_mem, res, t0, yy0, yp0)
end

function IDAInit(ida_mem, res, t0, yy0, yp0)
    __yy0 = convert(NVector, yy0)
    __yp0 = convert(NVector, yp0)
    __IDAInit(ida_mem, IDAResFn_wrapper(res), t0, convert(N_Vector, __yy0), convert(N_Vector, __yp0))
end

function __IDAReInit(ida_mem, t0::realtype, yy0::N_Vector, yp0::N_Vector)
    ccall((:IDAReInit, libsundials_ida), Cint, (IDAMemPtr, realtype, N_Vector, N_Vector), ida_mem, t0, yy0, yp0)
end

function IDAReInit(ida_mem, t0, yy0, yp0)
    __yy0 = convert(NVector, yy0)
    __yp0 = convert(NVector, yp0)
    __IDAReInit(ida_mem, t0, convert(N_Vector, __yy0), convert(N_Vector, __yp0))
end

function IDASStolerances(ida_mem, reltol::realtype, abstol::realtype)
    ccall((:IDASStolerances, libsundials_ida), Cint, (IDAMemPtr, realtype, realtype), ida_mem, reltol, abstol)
end

function __IDASVtolerances(ida_mem, reltol::realtype, abstol::N_Vector)
    ccall((:IDASVtolerances, libsundials_ida), Cint, (IDAMemPtr, realtype, N_Vector), ida_mem, reltol, abstol)
end

function IDASVtolerances(ida_mem, reltol, abstol)
    __abstol = convert(NVector, abstol)
    __IDASVtolerances(ida_mem, reltol, convert(N_Vector, __abstol))
end

function IDAWFtolerances(ida_mem, efun::IDAEwtFn)
    ccall((:IDAWFtolerances, libsundials_ida), Cint, (IDAMemPtr, IDAEwtFn), ida_mem, efun)
end

function IDASetNonlinConvCoefIC(ida_mem, epiccon::realtype)
    ccall((:IDASetNonlinConvCoefIC, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, epiccon)
end

function __IDASetMaxNumStepsIC(ida_mem, maxnh::Cint)
    ccall((:IDASetMaxNumStepsIC, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxnh)
end

function IDASetMaxNumStepsIC(ida_mem, maxnh)
    __IDASetMaxNumStepsIC(ida_mem, convert(Cint, maxnh))
end

function __IDASetMaxNumJacsIC(ida_mem, maxnj::Cint)
    ccall((:IDASetMaxNumJacsIC, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxnj)
end

function IDASetMaxNumJacsIC(ida_mem, maxnj)
    __IDASetMaxNumJacsIC(ida_mem, convert(Cint, maxnj))
end

function __IDASetMaxNumItersIC(ida_mem, maxnit::Cint)
    ccall((:IDASetMaxNumItersIC, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxnit)
end

function IDASetMaxNumItersIC(ida_mem, maxnit)
    __IDASetMaxNumItersIC(ida_mem, convert(Cint, maxnit))
end

function __IDASetLineSearchOffIC(ida_mem, lsoff::Cint)
    ccall((:IDASetLineSearchOffIC, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, lsoff)
end

function IDASetLineSearchOffIC(ida_mem, lsoff)
    __IDASetLineSearchOffIC(ida_mem, convert(Cint, lsoff))
end

function IDASetStepToleranceIC(ida_mem, steptol::realtype)
    ccall((:IDASetStepToleranceIC, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, steptol)
end

function __IDASetMaxBacksIC(ida_mem, maxbacks::Cint)
    ccall((:IDASetMaxBacksIC, libsundials_ida), Cint, (IDAMemPtr, Cint), ida_mem, maxbacks)
end

function IDASetMaxBacksIC(ida_mem, maxbacks)
    __IDASetMaxBacksIC(ida_mem, convert(Cint, maxbacks))
end

function __IDARootInit(ida_mem, nrtfn::Cint, g::IDARootFn)
    ccall((:IDARootInit, libsundials_ida), Cint, (IDAMemPtr, Cint, IDARootFn), ida_mem, nrtfn, g)
end

function IDARootInit(ida_mem, nrtfn, g)
    __IDARootInit(ida_mem, convert(Cint, nrtfn), IDARootFn_wrapper(g))
end

function __IDACalcIC(ida_mem, icopt::Cint, tout1::realtype)
    ccall((:IDACalcIC, libsundials_ida), Cint, (IDAMemPtr, Cint, realtype), ida_mem, icopt, tout1)
end

function IDACalcIC(ida_mem, icopt, tout1)
    __IDACalcIC(ida_mem, convert(Cint, icopt), tout1)
end

function __IDASolve(ida_mem, tout::realtype, tret, yret::N_Vector, ypret::N_Vector, itask::Cint)
    ccall((:IDASolve, libsundials_ida), Cint, (IDAMemPtr, realtype, Ptr{realtype}, N_Vector, N_Vector, Cint), ida_mem, tout, tret, yret, ypret, itask)
end

function IDASolve(ida_mem, tout, tret, yret, ypret, itask)
    __yret = convert(NVector, yret)
    __ypret = convert(NVector, ypret)
    __IDASolve(ida_mem, tout, tret, convert(N_Vector, __yret), convert(N_Vector, __ypret), convert(Cint, itask))
end

function __IDAGetDky(ida_mem, t::realtype, k::Cint, dky::N_Vector)
    ccall((:IDAGetDky, libsundials_ida), Cint, (IDAMemPtr, realtype, Cint, N_Vector), ida_mem, t, k, dky)
end

function IDAGetDky(ida_mem, t, k, dky)
    __dky = convert(NVector, dky)
    __IDAGetDky(ida_mem, t, convert(Cint, k), convert(N_Vector, __dky))
end

function IDAGetWorkSpace(ida_mem, lenrw, leniw)
    ccall((:IDAGetWorkSpace, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}), ida_mem, lenrw, leniw)
end

function IDAGetNumSteps(ida_mem, nsteps)
    ccall((:IDAGetNumSteps, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nsteps)
end

function IDAGetNumResEvals(ida_mem, nrevals)
    ccall((:IDAGetNumResEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nrevals)
end

function IDAGetNumLinSolvSetups(ida_mem, nlinsetups)
    ccall((:IDAGetNumLinSolvSetups, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nlinsetups)
end

function IDAGetNumErrTestFails(ida_mem, netfails)
    ccall((:IDAGetNumErrTestFails, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, netfails)
end

function IDAGetNumBacktrackOps(ida_mem, nbacktr)
    ccall((:IDAGetNumBacktrackOps, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nbacktr)
end

function __IDAGetConsistentIC(ida_mem, yy0_mod::N_Vector, yp0_mod::N_Vector)
    ccall((:IDAGetConsistentIC, libsundials_ida), Cint, (IDAMemPtr, N_Vector, N_Vector), ida_mem, yy0_mod, yp0_mod)
end

function IDAGetConsistentIC(ida_mem, yy0_mod, yp0_mod)
    __yy0_mod = convert(NVector, yy0_mod)
    __yp0_mod = convert(NVector, yp0_mod)
    __IDAGetConsistentIC(ida_mem, convert(N_Vector, __yy0_mod), convert(N_Vector, __yp0_mod))
end

function IDAGetLastOrder(ida_mem, klast)
    ccall((:IDAGetLastOrder, libsundials_ida), Cint, (IDAMemPtr, Ptr{Cint}), ida_mem, klast)
end

function IDAGetCurrentOrder(ida_mem, kcur)
    ccall((:IDAGetCurrentOrder, libsundials_ida), Cint, (IDAMemPtr, Ptr{Cint}), ida_mem, kcur)
end

function IDAGetActualInitStep(ida_mem, hinused)
    ccall((:IDAGetActualInitStep, libsundials_ida), Cint, (IDAMemPtr, Ptr{realtype}), ida_mem, hinused)
end

function IDAGetLastStep(ida_mem, hlast)
    ccall((:IDAGetLastStep, libsundials_ida), Cint, (IDAMemPtr, Ptr{realtype}), ida_mem, hlast)
end

function IDAGetCurrentStep(ida_mem, hcur)
    ccall((:IDAGetCurrentStep, libsundials_ida), Cint, (IDAMemPtr, Ptr{realtype}), ida_mem, hcur)
end

function IDAGetCurrentTime(ida_mem, tcur)
    ccall((:IDAGetCurrentTime, libsundials_ida), Cint, (IDAMemPtr, Ptr{realtype}), ida_mem, tcur)
end

function IDAGetTolScaleFactor(ida_mem, tolsfact)
    ccall((:IDAGetTolScaleFactor, libsundials_ida), Cint, (IDAMemPtr, Ptr{realtype}), ida_mem, tolsfact)
end

function __IDAGetErrWeights(ida_mem, eweight::N_Vector)
    ccall((:IDAGetErrWeights, libsundials_ida), Cint, (IDAMemPtr, N_Vector), ida_mem, eweight)
end

function IDAGetErrWeights(ida_mem, eweight)
    __eweight = convert(NVector, eweight)
    __IDAGetErrWeights(ida_mem, convert(N_Vector, __eweight))
end

function __IDAGetEstLocalErrors(ida_mem, ele::N_Vector)
    ccall((:IDAGetEstLocalErrors, libsundials_ida), Cint, (IDAMemPtr, N_Vector), ida_mem, ele)
end

function IDAGetEstLocalErrors(ida_mem, ele)
    __ele = convert(NVector, ele)
    __IDAGetEstLocalErrors(ida_mem, convert(N_Vector, __ele))
end

function IDAGetNumGEvals(ida_mem, ngevals)
    ccall((:IDAGetNumGEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, ngevals)
end

function IDAGetRootInfo(ida_mem, rootsfound)
    ccall((:IDAGetRootInfo, libsundials_ida), Cint, (IDAMemPtr, Ptr{Cint}), ida_mem, rootsfound)
end

function IDAGetIntegratorStats(ida_mem, nsteps, nrevals, nlinsetups, netfails, qlast, qcur, hinused, hlast, hcur, tcur)
    ccall((:IDAGetIntegratorStats, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Clong}, Ptr{Cint}, Ptr{Cint}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}, Ptr{realtype}), ida_mem, nsteps, nrevals, nlinsetups, netfails, qlast, qcur, hinused, hlast, hcur, tcur)
end

function IDAGetNumNonlinSolvIters(ida_mem, nniters)
    ccall((:IDAGetNumNonlinSolvIters, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nniters)
end

function IDAGetNumNonlinSolvConvFails(ida_mem, nncfails)
    ccall((:IDAGetNumNonlinSolvConvFails, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nncfails)
end

function IDAGetNonlinSolvStats(ida_mem, nniters, nncfails)
    ccall((:IDAGetNonlinSolvStats, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}), ida_mem, nniters, nncfails)
end

function __IDAGetReturnFlagName(flag::Clong)
    ccall((:IDAGetReturnFlagName, libsundials_ida), Cstring, (Clong,), flag)
end

function IDAGetReturnFlagName(flag)
    __IDAGetReturnFlagName(convert(Clong, flag))
end

function IDAFree(ida_mem)
    ccall((:IDAFree, libsundials_ida), Void, (Ref{IDAMemPtr},), ida_mem)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/ida/ida_bbdpre.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function IDABBDPrecInit(ida_mem, Nlocal::sunindextype, mudq::sunindextype, mldq::sunindextype, mukeep::sunindextype, mlkeep::sunindextype, dq_rel_yy::realtype, Gres::IDABBDLocalFn, Gcomm::IDABBDCommFn)
    ccall((:IDABBDPrecInit, libsundials_ida), Cint, (IDAMemPtr, sunindextype, sunindextype, sunindextype, sunindextype, sunindextype, realtype, IDABBDLocalFn, IDABBDCommFn), ida_mem, Nlocal, mudq, mldq, mukeep, mlkeep, dq_rel_yy, Gres, Gcomm)
end

function IDABBDPrecReInit(ida_mem, mudq::sunindextype, mldq::sunindextype, dq_rel_yy::realtype)
    ccall((:IDABBDPrecReInit, libsundials_ida), Cint, (IDAMemPtr, sunindextype, sunindextype, realtype), ida_mem, mudq, mldq, dq_rel_yy)
end

function IDABBDPrecGetWorkSpace(ida_mem, lenrwBBDP, leniwBBDP)
    ccall((:IDABBDPrecGetWorkSpace, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}), ida_mem, lenrwBBDP, leniwBBDP)
end

function IDABBDPrecGetNumGfnEvals(ida_mem, ngevalsBBDP)
    ccall((:IDABBDPrecGetNumGfnEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, ngevalsBBDP)
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/ida/ida_direct.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function IDADlsSetLinearSolver(ida_mem, LS::SUNLinearSolver, A::SUNMatrix)
    ccall((:IDADlsSetLinearSolver, libsundials_ida), Cint, (IDAMemPtr, SUNLinearSolver, SUNMatrix), ida_mem, LS, A)
end

function IDADlsSetJacFn(ida_mem, jac::IDADlsJacFn)
    ccall((:IDADlsSetJacFn, libsundials_ida), Cint, (IDAMemPtr, IDADlsJacFn), ida_mem, jac)
end

function IDADlsGetWorkSpace(ida_mem, lenrwLS, leniwLS)
    ccall((:IDADlsGetWorkSpace, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}), ida_mem, lenrwLS, leniwLS)
end

function IDADlsGetNumJacEvals(ida_mem, njevals)
    ccall((:IDADlsGetNumJacEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, njevals)
end

function IDADlsGetNumResEvals(ida_mem, nfevalsLS)
    ccall((:IDADlsGetNumResEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nfevalsLS)
end

function IDADlsGetLastFlag(ida_mem, flag)
    ccall((:IDADlsGetLastFlag, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, flag)
end

function __IDADlsGetReturnFlagName(flag::Clong)
    ccall((:IDADlsGetReturnFlagName, libsundials_ida), Cstring, (Clong,), flag)
end

function IDADlsGetReturnFlagName(flag)
    __IDADlsGetReturnFlagName(convert(Clong, flag))
end
# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/ida/ida_impl.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /home/osboxes/.julia/v0.6/Sundials/deps/usr/include/ida/ida_spils.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function IDASpilsSetLinearSolver(ida_mem, LS::SUNLinearSolver)
    ccall((:IDASpilsSetLinearSolver, libsundials_ida), Cint, (IDAMemPtr, SUNLinearSolver), ida_mem, LS)
end

function IDASpilsSetPreconditioner(ida_mem, pset::IDASpilsPrecSetupFn, psolve::IDASpilsPrecSolveFn)
    ccall((:IDASpilsSetPreconditioner, libsundials_ida), Cint, (IDAMemPtr, IDASpilsPrecSetupFn, IDASpilsPrecSolveFn), ida_mem, pset, psolve)
end

function IDASpilsSetJacTimes(ida_mem, jtsetup::IDASpilsJacTimesSetupFn, jtimes::IDASpilsJacTimesVecFn)
    ccall((:IDASpilsSetJacTimes, libsundials_ida), Cint, (IDAMemPtr, IDASpilsJacTimesSetupFn, IDASpilsJacTimesVecFn), ida_mem, jtsetup, jtimes)
end

function IDASpilsSetEpsLin(ida_mem, eplifac::realtype)
    ccall((:IDASpilsSetEpsLin, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, eplifac)
end

function IDASpilsSetIncrementFactor(ida_mem, dqincfac::realtype)
    ccall((:IDASpilsSetIncrementFactor, libsundials_ida), Cint, (IDAMemPtr, realtype), ida_mem, dqincfac)
end

function IDASpilsGetWorkSpace(ida_mem, lenrwLS, leniwLS)
    ccall((:IDASpilsGetWorkSpace, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}, Ptr{Clong}), ida_mem, lenrwLS, leniwLS)
end

function IDASpilsGetNumPrecEvals(ida_mem, npevals)
    ccall((:IDASpilsGetNumPrecEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, npevals)
end

function IDASpilsGetNumPrecSolves(ida_mem, npsolves)
    ccall((:IDASpilsGetNumPrecSolves, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, npsolves)
end

function IDASpilsGetNumLinIters(ida_mem, nliters)
    ccall((:IDASpilsGetNumLinIters, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nliters)
end

function IDASpilsGetNumConvFails(ida_mem, nlcfails)
    ccall((:IDASpilsGetNumConvFails, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nlcfails)
end

function IDASpilsGetNumJtimesEvals(ida_mem, njvevals)
    ccall((:IDASpilsGetNumJtimesEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, njvevals)
end

function IDASpilsGetNumResEvals(ida_mem, nrevalsLS)
    ccall((:IDASpilsGetNumResEvals, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, nrevalsLS)
end

function IDASpilsGetLastFlag(ida_mem, flag)
    ccall((:IDASpilsGetLastFlag, libsundials_ida), Cint, (IDAMemPtr, Ptr{Clong}), ida_mem, flag)
end

function __IDASpilsGetReturnFlagName(flag::Clong)
    ccall((:IDASpilsGetReturnFlagName, libsundials_ida), Cstring, (Clong,), flag)
end

function IDASpilsGetReturnFlagName(flag)
    __IDASpilsGetReturnFlagName(convert(Clong, flag))
end


text:
###
# PseudoBlockArray apply
###

arguments(LAY::MemoryLayout, A::PseudoBlockArray) = arguments(LAY, A.blocks)

##########
# BlockVcat
##########

struct BlockVcat{T, N, Arrays} <: AbstractBlockArray{T,N}
    arrays::Arrays
    function BlockVcat{T,N,Arrays}(arrays::Arrays) where {T,N,Arrays}
        length(arrays) == 1 || blockisequal(axes.(arrays,2)...) || throw(ArgumentError("Blocks must match"))
        new{T,N,Arrays}(arrays)
    end
end

BlockVcat{T,N}(arrays::AbstractArray...) where {T,N} = BlockVcat{T,N,typeof(arrays)}(arrays)
BlockVcat{T}(arrays::AbstractArray{<:Any,N}...) where {T,N} = BlockVcat{T,N}(arrays...)
BlockVcat(arrays::AbstractArray...) = BlockVcat{mapreduce(eltype, promote_type, arrays)}(arrays...)
blockvcat(a) = a
blockvcat(a, b...) = BlockVcat(a, b...)

# all
_vcat_axes(ax::OneTo{Int}...) = blockedrange(SVector(map(length,ax)...))
_vcat_axes(ax...) = blockedrange(vcat(map(blocklengths,ax)...))
axes(b::BlockVcat{<:Any,1}) = (_vcat_axes(axes.(b.arrays,1)...),)
axes(b::BlockVcat{<:Any,2}) = (_vcat_axes(axes.(b.arrays,1)...),axes(b.arrays[1],2))

copy(b::BlockVcat{T,N}) where {T,N} = BlockVcat{T,N}(map(copy, b.arrays)...)
copy(b::AdjOrTrans{<:Any,<:BlockVcat}) = copy(parent(b))'
AbstractArray{T}(B::BlockVcat{<:Any,N}) where {T,N} = BlockVcat{T,N}(map(AbstractArray{T}, B.arrays)...)
AbstractArray{T,N}(B::BlockVcat{<:Any,N}) where {T,N} = BlockVcat{T,N}(map(AbstractArray{T}, B.arrays)...)
convert(::Type{AbstractArray{T}}, B::BlockVcat{<:Any,N}) where {T,N} = BlockVcat{T,N}(convert.(AbstractArray{T}, B.arrays)...)
convert(::Type{AbstractArray{T,N}}, B::BlockVcat{<:Any,N}) where {T,N} = BlockVcat{T,N}(convert.(AbstractArray{T}, B.arrays)...)
convert(::Type{AbstractArray{T}}, B::BlockVcat{T,N}) where {T,N} = B
convert(::Type{AbstractArray{T,N}}, B::BlockVcat{T,N}) where {T,N} = B

# avoid making naive view
_viewifblocked(::Tuple{Vararg{OneTo}}, a, kj::Block{1}) = a
_viewifblocked(::Tuple{Vararg{OneTo}}, a, kj::Block{2}) = a
_viewifblocked(::Tuple{OneTo}, a::AbstractVector, kj::Block{2}) = a
_viewifblocked(_, a, kj::Block) = view(a, kj)
function _viewifblocked(_, a::AbstractVector, kj::Block{2})
    k,j = kj.n
    @assert j == 1
    view(a, Block(k))
end
_viewifblocked(a, kj) = _viewifblocked(axes(a), a, kj)

_findvcatblock(k) = throw(BoundsError())

function _findvcatblock(k::Block{1}, a, b...)
    n = Block(blocklength(a))
    k ≤ n && return _viewifblocked(a, k)
    _findvcatblock(k - n, b...)
end

function _findvcatblock(kj::Block{2}, a, b...)
    k,j= kj.n
    n = blocksize(a,1)
    k ≤ n && return _viewifblocked(a, kj)
    _findvcatblock(Block(k-n, j), b...)
end

viewblock(b::BlockVcat, k::Block) = _findvcatblock(k, b.arrays...)
getindex(b::BlockVcat, Kk::BlockIndex{1}) = view(b,block(Kk))[Kk.α...]
getindex(b::BlockVcat, Kk::BlockIndex{2}) = view(b,block(Kk))[Kk.α...]
getindex(b::BlockVcat{<:Any,1}, k::Integer) = b[findblockindex(axes(b,1), k)]
getindex(b::BlockVcat{<:Any,2}, k::Integer, j::Integer) = b[findblockindex(axes(b,1),k), findblockindex(axes(b,2),j)]

MemoryLayout(::Type{<:BlockVcat}) = ApplyLayout{typeof(vcat)}()
arguments(::ApplyLayout{typeof(vcat)}, b::BlockVcat) = b.arrays

sub_materialize(lay::ApplyLayout{typeof(vcat)}, V::AbstractVector, ::Tuple{<:BlockedUnitRange}) = blockvcat(arguments(lay, V)...)
sub_materialize(lay::ApplyLayout{typeof(vcat)}, V::AbstractMatrix, ::Tuple{<:BlockedUnitRange,<:BlockedUnitRange}) = blockvcat(arguments(lay, V)...)
sub_materialize(lay::ApplyLayout{typeof(vcat)}, V::AbstractMatrix, ::Tuple{<:BlockedUnitRange,<:AbstractUnitRange}) = blockvcat(arguments(lay, V)...)

LazyArrays._vcat_sub_arguments(lay::ApplyLayout{typeof(vcat)}, A, V, kr::BlockSlice{<:BlockRange{1}}) =
    arguments(lay, A)[Int.(kr.block)]

_split2blocks(KR) = ()
function _split2blocks(KR, ax::OneTo, C...)
    if isempty(KR)
        (Base.OneTo(0), _split2blocks(Block.(1:0), C...)...)
    elseif first(KR) ≠ Block(1)
        (Base.OneTo(0), _split2blocks((KR[1] - Block(1)):(KR[end] - Block(1)), C...)...)
    elseif length(KR) == 1
        (ax, _split2blocks(Block.(1:0), C...)...)
    else
        (ax, _split2blocks((KR[2]- Block(1)):(KR[end]-Block(1)), C...)...)
    end
end
function _split2blocks(KR, A, C...)
    M = blocklength(A)
    if Int(last(KR)) ≤ M
        (KR, _split2blocks(Block.(1:0), C...)...)
    else
        (KR[1]:Block(M), _split2blocks(Block(1):(last(KR)-Block(M)), C...)...)
    end
end

function arguments(lay::ApplyLayout{typeof(vcat)}, V::SubArray{<:Any,1,<:Any,<:Tuple{BlockSlice{<:BlockRange1}}})
    kr, = parentindices(V)
    P = parent(V)
    a = arguments(lay, P)
    KR = _split2blocks(kr.block, axes.(a,1)...)
    filter(!isempty,getindex.(a, KR))
end

function arguments(lay::ApplyLayout{typeof(vcat)}, V::SubArray{<:Any,2,<:Any,<:Tuple{BlockSlice{<:BlockRange1},BlockSlice}})
    kr, jr = parentindices(V)
    P = parent(V)
    a = arguments(lay, P)
    KR = _split2blocks(kr.block, axes.(a,1)...)
    filter(!isempty,getindex.(a, KR, Ref(jr.block)))
end

##########
# BlockHcat
##########

struct BlockHcat{T, Arrays} <: AbstractBlockMatrix{T}
    arrays::Arrays
    function BlockHcat{T,Arrays}(arrays::Arrays) where {T,Arrays}
        length(arrays) == 1 || blockisequal(axes.(arrays,1)...) || throw(ArgumentError("Blocks must match"))
        new{T,Arrays}(arrays)
    end
end

BlockHcat{T}(arrays::AbstractArray...) where T = BlockHcat{T,typeof(arrays)}(arrays)
BlockHcat(arrays::AbstractArray...) = BlockHcat{mapreduce(eltype, promote_type, arrays)}(arrays...)
blockhcat(a) = a
blockhcat(a, b...) = BlockHcat(a, b...)

axes(b::BlockHcat) = (axes(b.arrays[1],1),_vcat_axes(axes.(b.arrays,2)...))
axes(b::BlockHcat{<:Any, <:Tuple{Vararg{AbstractVector}}}) = (axes(b.arrays[1],1),blockedrange(Ones{Int}(length(b.arrays))))

copy(b::BlockHcat{T}) where T = BlockHcat{T}(map(copy, b.arrays)...)
copy(b::AdjOrTrans{<:Any,<:BlockHcat}) = copy(parent(b))'
AbstractArray{T}(B::BlockHcat) where T = BlockHcat{T}(map(AbstractArray{T}, B.arrays)...)
AbstractMatrix{T}(B::BlockHcat) where T = BlockHcat{T}(map(AbstractArray{T}, B.arrays)...)
convert(::Type{AbstractArray{T}}, B::BlockHcat) where T = BlockHcat{T}(convert.(AbstractArray{T}, B.arrays)...)
convert(::Type{AbstractMatrix{T}}, B::BlockHcat) where T = BlockHcat{T}(convert.(AbstractArray{T}, B.arrays)...)
convert(::Type{AbstractArray{T}}, B::BlockHcat{T}) where T = B
convert(::Type{AbstractMatrix{T}}, B::BlockHcat{T}) where T = B

_blocksize2(a::AbstractVector) = 1
_blocksize2(a) = blocksize(a,2)

function _findhcatblock(kj::Block{2}, a, b...)
    k,j = kj.n
    n = _blocksize2(a)
    j ≤ n && return _viewifblocked(a, kj)
    _findhcatblock(Block(k, j-n), b...)
end

viewblock(b::BlockHcat, kj::Block{2}) = _findhcatblock(kj, b.arrays...)
getindex(b::BlockHcat, Kk::BlockIndex{1}, Jj::BlockIndex{1}) = view(b,block(Kk), block(Jj))[blockindex(Kk), blockindex(Jj)]
getindex(b::BlockHcat, k::Integer, j::Integer) = b[findblockindex(axes(b,1),k), findblockindex(axes(b,2),j)]

arguments(::ApplyLayout{typeof(hcat)}, b::BlockHcat) = b.arrays

sub_materialize(lay::ApplyLayout{typeof(hcat)}, V::AbstractMatrix, ::Tuple{<:BlockedUnitRange,<:BlockedUnitRange}) = blockhcat(arguments(lay, V)...)
sub_materialize(lay::ApplyLayout{typeof(hcat)}, V::AbstractMatrix, ::Tuple{<:BlockedUnitRange,<:AbstractUnitRange}) = blockhcat(arguments(lay, V)...)
sub_materialize(lay::ApplyLayout{typeof(hcat)}, V::AbstractMatrix, ::Tuple{<:AbstractUnitRange,<:BlockedUnitRange}) = blockhcat(arguments(lay, V)...)

function arguments(lay::ApplyLayout{typeof(hcat)}, V::SubArray{<:Any,2,<:Any,<:Tuple{BlockSlice,BlockSlice{<:BlockRange1}}})
    kr, jr = parentindices(V)
    P = parent(V)
    a = arguments(lay, P)
    JR = _split2blocks(jr.block, axes.(a,2)...)
    filter(!isempty,getindex.(a, Ref(kr.block), JR))
end

function arguments(lay::ApplyLayout{typeof(hcat)}, V::SubArray{<:Any,2,<:Any,<:Tuple{BlockSlice{<:BlockRange1},BlockSlice{<:Block1}}})
    kr, jr = parentindices(V)
    J = jr.block
    P = parent(V)
    arguments(lay, view(P,kr,J:J))
end

for adj in (:adjoint, :transpose)
    @eval begin
        $adj(A::BlockHcat{T}) where T = BlockVcat{T}(map($adj,A.arrays)...)
        $adj(A::BlockVcat{T,2}) where T = BlockHcat{T}(map($adj,A.arrays)...)
    end
end

##########
# BlockHvcat
##########

struct BlockHvcat{T, Args} <: AbstractBlockMatrix{T}
    n::Int
    args::Args
    function BlockHvcat{T,Args}(n::Int, args::Args) where {T,Args}
        new{T,Args}(n, args)
    end
end

BlockHvcat{T}(n::Int, args...) where T = BlockHvcat{T,typeof(args)}(n, args)
BlockHvcat(n::Int, args...) = BlockHvcat{mapreduce(eltype, promote_type, args)}(n, args...)

axes(b::BlockHvcat) = (_vcat_axes(axes.(b.args[1:b.n:end],1)...),_vcat_axes(axes.(b.args[1:b.n],2)...))

copy(b::BlockHvcat{T}) where T = BlockHvcat{T}(b.n, map(copy, b.args)...)
copy(b::AdjOrTrans{<:Any,<:BlockHvcat}) = copy(parent(b))'

# _hvcat_viewifblocked(::OneTo, a::AbstractMatrix, k) = a
# _hvcat_viewifblocked(::OneTo, a::AbstractVector, k) = a
# _hvcat_viewifblocked(_, a::AbstractMatrix, k) = view(a, Block(k,1))
# _hvcat_viewifblocked(_, a::AbstractVector, k) = view(a, Block(k))
# _hvcat_viewifblocked(a, k) = _hvcat_viewifblocked(axes(a,1), a, k)

function viewblock(b::BlockHvcat, kj::Block{2})
    k,j = kj.n
    b.args[b.n*(k-1) + j] # TODO: blocked
end
getindex(b::BlockHvcat, Kk::BlockIndex{1}, Jj::BlockIndex{1}) = view(b,block(Kk), block(Jj))[blockindex(Kk), blockindex(Jj)]
getindex(b::BlockHvcat, k::Integer, j::Integer) = b[findblockindex(axes(b,1),k), findblockindex(axes(b,2),j)]

# MemoryLayout(::Type{<:BlockHvcat}) = ApplyLayout{typeof(hvcat)}()
# arguments(::ApplyLayout{typeof(hvcat)}, b::BlockHvcat) = b.args

# sub_materialize(lay::ApplyLayout{typeof(hvcat)}, V::AbstractMatrix, ::Tuple{<:BlockedUnitRange,<:BlockedUnitRange}) = BlockHvcat(arguments(lay, V)...)
# sub_materialize(lay::ApplyLayout{typeof(hvcat)}, V::AbstractMatrix, ::Tuple{<:AbstractUnitRange,<:BlockedUnitRange}) = BlockHvcat(arguments(lay, V)...)

# function arguments(lay::ApplyLayout{typeof(hvcat)}, V::SubArray{<:Any,2,<:Any,<:Tuple{BlockSlice{<:BlockRange1},BlockSlice{<:Block1}}})
#     kr, jr = parentindices(V)
#     @assert kr.block ≡ Block(1)
#     arguments(lay, parent(V))[Int.(jr.block)]
# end

#############
# BlockApplyArray
#############

"""
   BlockBroadcastArray(f, A, B, C...)

is a block array corresponding to `f.(blocks(A), blocks(B), ...)`,
except if `A` is scalar.
"""

struct BlockBroadcastArray{T, N, FF, Args} <: AbstractBlockArray{T, N}
    f::FF
    args::Args
end

const BlockBroadcastVector{T, FF, Args} = BlockBroadcastArray{T, 1, FF, Args}
const BlockBroadcastMatrix{T, FF, Args} = BlockBroadcastArray{T, 2, FF, Args}

_blocks(A::AbstractArray) = blocks(A)
_blocks(A::Number) = A

BlockBroadcastArray{T,N}(f, args...) where {T,N} = BlockBroadcastArray{T,N,typeof(f),typeof(args)}(f, args)
BlockBroadcastArray{T}(bc::Broadcasted) where T = BlockBroadcastArray{T}(bc.f, bc.args...)
BlockBroadcastArray(bc::Broadcasted) = BlockBroadcastArray{eltype(Base.Broadcast.combine_eltypes(bc.f, bc.args))}(bc)
BlockBroadcastArray(f, args...) = BlockBroadcastArray(instantiate(broadcasted(f, args...)))

BlockBroadcastArray{T}(::typeof(hcat), args...) where T = BlockBroadcastMatrix{T}(hcat, args...)
BlockBroadcastArray{T}(::typeof(vcat), args::AbstractVector...) where T = BlockBroadcastVector{T}(vcat, args...)
BlockBroadcastArray{T}(::typeof(vcat), args...) where T = BlockBroadcastMatrix{T}(vcat, args...)
BlockBroadcastArray{T}(::typeof(hvcat), args...) where T = BlockBroadcastMatrix{T}(hvcat, args...)
BlockBroadcastArray{T}(::typeof(Diagonal), args...) where T = BlockBroadcastMatrix{T}(Diagonal, args...)

_block_vcat_axes(ax...) = BlockArrays._BlockedUnitRange(1,+(map(blocklasts,ax)...))

_block_interlace_axes(::Int, ax::Tuple{BlockedUnitRange{OneTo{Int}}}...) = _block_vcat_axes(ax...)

function _block_interlace_axes(nbc::Int, ax::NTuple{2,BlockedUnitRange{OneTo{Int}}}...)
    n,m = max(length.(first.(ax))...),max(length.(last.(ax))...)
    (blockedrange(Fill(length(ax) ÷ nbc, n)),blockedrange(Fill(mod1(length(ax),nbc), m)))
end

axes(A::BlockBroadcastVector{<:Any,typeof(vcat)}) = (_block_vcat_axes(axes.(A.args,1)...),)
axes(A::BlockBroadcastMatrix{<:Any,typeof(hcat)}) = (axes(A.args[1],1), _block_vcat_axes(axes.(A.args,2)...))

axes(A::BlockBroadcastMatrix{<:Any,typeof(hvcat)}) = _block_interlace_axes(A.args[1], map(axes,A.args[2:end])...)
axes(A::BlockBroadcastMatrix{<:Any,typeof(Diagonal)}) = (_block_vcat_axes(axes.(A.args,1)...), _block_vcat_axes(axes.(A.args,2)...))
# size(A::BlockBroadcastArray) = map(length, axes(A))

copy(b::BlockBroadcastArray{T,N}) where {T,N} = BlockBroadcastArray{T,N}(b.f, map(copy, b.args)...)
copy(b::AdjOrTrans{<:Any,<:BlockBroadcastArray}) = copy(parent(b))'

AbstractArray{T}(B::BlockBroadcastArray{<:Any,N}) where {T,N} = BlockBroadcastArray{T,N}(B.f, map(AbstractArray{T}, B.args)...)
AbstractArray{T,N}(B::BlockBroadcastArray{<:Any,N}) where {T,N} = BlockBroadcastArray{T,N}(B.f, map(AbstractArray{T}, B.args)...)
convert(::Type{AbstractArray{T}}, B::BlockBroadcastArray{<:Any,N}) where {T,N} = BlockBroadcastArray{T,N}(B.f, convert.(AbstractArray{T}, B.args)...)
convert(::Type{AbstractArray{T,N}}, B::BlockBroadcastArray{<:Any,N}) where {T,N} = BlockBroadcastArray{T,N}(B.f, convert.(AbstractArray{T}, B.args)...)
convert(::Type{AbstractArray{T}}, B::BlockBroadcastArray{T,N}) where {T,N} = B
convert(::Type{AbstractArray{T,N}}, B::BlockBroadcastArray{T,N}) where {T,N} = B

AbstractArray{T}(B::BlockBroadcastArray{<:Any,N,typeof(hvcat)}) where {T,N} = BlockBroadcastArray{T,N}(B.f, first(B.args), map(AbstractArray{T}, tail(B.args))...)
AbstractArray{T,N}(B::BlockBroadcastArray{<:Any,N,typeof(hvcat)}) where {T,N} = BlockBroadcastArray{T,N}(B.f, first(B.args), map(AbstractArray{T}, tail(B.args))...)
convert(::Type{AbstractArray{T}}, B::BlockBroadcastArray{<:Any,N,typeof(hvcat)}) where {T,N} = BlockBroadcastArray{T,N}(B.f, first(B.args), convert.(AbstractArray{T}, tail(B.args))...)
convert(::Type{AbstractArray{T,N}}, B::BlockBroadcastArray{<:Any,N,typeof(hvcat)}) where {T,N} = BlockBroadcastArray{T,N}(B.f, first(B.args), convert.(AbstractArray{T}, tail(B.args))...)
convert(::Type{AbstractArray{T}}, B::BlockBroadcastArray{T,N,typeof(hvcat)}) where {T,N} = B
convert(::Type{AbstractArray{T,N}}, B::BlockBroadcastArray{T,N,typeof(hvcat)}) where {T,N} = B

Base.BroadcastStyle(::Type{<:BlockBroadcastArray{T,N}}) where {T,N} = LazyArrayStyle{N}()

function getindex(A::BlockBroadcastVector{<:Any,typeof(vcat)}, k::Int)
    K = findblockindex(axes(A,1), k)
    A.args[blockindex(K)][Int(block(K))]
end

function getindex(A::BlockBroadcastMatrix{<:Any,typeof(hcat)}, k::Int, j::Int)
    J = findblockindex(axes(A,2), j)
    A.args[blockindex(J)][k,Int(block(J))]
end

function getindex(A::BlockBroadcastMatrix{<:Any,typeof(hvcat)}, k::Int, j::Int)
    K,J = findblockindex.(axes(A), (k,j))
    A.args[(blockindex(K)-1)*A.args[1] + blockindex(J)+1][Int(block(K)), Int(block(J))]
end

function getindex(A::BlockBroadcastMatrix{T,typeof(Diagonal)}, k::Int, j::Int) where T
    K,J = findblockindex.(axes(A), (k,j))
    blockindex(K) == blockindex(J) || return zero(T)
    A.args[blockindex(K)][Int(block(K)),Int(block(J))]
end

viewblock(A::BlockBroadcastVector{<:Any,typeof(vcat)}, k::Block{1}) = Vcat(getindex.(A.args, Int(k))...)
viewblock(A::BlockBroadcastMatrix{<:Any,typeof(hcat)}, kj::Block{2}) = Hcat(getindex.(A.args, kj.n...)...)
viewblock(A::BlockBroadcastMatrix{<:Any,typeof(hvcat)}, kj::Block{2}) = ApplyArray(hvcat, first(A.args), getindex.(tail(A.args), kj.n...)...)

###
# blockbandwidths
###
blockbandwidths(A::BlockBroadcastMatrix{<:Any,typeof(hvcat)}) = max.(map(blockbandwidths,Base.tail(A.args))...)
blockbandwidths(A::BlockBroadcastMatrix{<:Any,typeof(Diagonal)}) = max.(map(blockbandwidths,A.args)...)
subblockbandwidths(A::BlockBroadcastMatrix{<:Any,typeof(Diagonal)}) = (0,0)

function blockbandwidths(M::BlockVcat{<:Any,2})
    cs = tuple(0, _cumsum(blocksize.(M.arrays[1:end-1],1)...)...) # cumsum of sizes
    (maximum(cs .+ blockbandwidth.(M.arrays,1)), maximum(blockbandwidth.(M.arrays,2) .- cs))
end
isblockbanded(M::BlockVcat{<:Any,2}) = all(isblockbanded, M.arrays)

function blockbandwidths(M::BlockHcat)
    cs = tuple(0, _cumsum(blocksize.(M.arrays[1:end-1],2)...)...) # cumsum of sizes
    (maximum(blockbandwidth.(M.arrays,1) .- cs), maximum(blockbandwidth.(M.arrays,2) .+ cs))
end
isblockbanded(M::BlockHcat) = all(isblockbanded, M.arrays)

function subblockbandwidths(B::BlockBroadcastMatrix{<:Any,typeof(hvcat)})
    p = B.args[1]
    ret = p .* subblockbandwidths(B.args[2]) # initialise with first
    shft = 0
    rw = 0
    for a in Base.tail(Base.tail(B.args))
        shft += 1
        if shft == p
            # next row
            rw += 1
            shft = -rw
        end
        ret = max.(ret, p .* subblockbandwidths(a) .+ (-shft,shft))
    end
    ret
end

###
# MemoryLayout
#
# sometimes we get block/banded when concatenting block/banded matrices
###

blockhcatlayout(_...) = ApplyLayout{typeof(hcat)}()
# at the moment we just support hcat for a special case of a subview of Eye concatenated with a block banded.
# This can be generalised later as needed
blockhcatlayout(::AbstractBandedLayout, ::AbstractBlockBandedLayout) = ApplyBlockBandedLayout{typeof(hcat)}()
MemoryLayout(::Type{<:BlockHcat{<:Any,Arrays}}) where Arrays = blockhcatlayout(LazyArrays.tuple_type_memorylayouts(Arrays)...)
sublayout(::ApplyBlockBandedLayout{typeof(hcat)}, ::Type{<:Tuple{<:BlockSlice{<:BlockRange1}, <:BlockSlice{<:BlockRange1}}}) = ApplyBlockBandedLayout{typeof(hcat)}()

_copyto!(_, LAY::ApplyBlockBandedLayout{typeof(hcat)}, dest::AbstractMatrix, H::AbstractMatrix) =
    block_hcat_copyto!(dest, arguments(LAY,H)...)
function block_hcat_copyto!(dest::AbstractMatrix, arrays...)
    nrows = blocksize(dest, 1)
    ncols = 0
    dense = true
    for a in arrays
        dense &= isa(a,Array)
        nd = ndims(a)
        ncols += (nd==2 ? blocksize(a,2) : 1)
    end

    nrows == blocksize(first(arrays),1) || throw(DimensionMismatch("Destination rows must match"))
    ncols == blocksize(dest,2) || throw(DimensionMismatch("Destination columns must match"))

    pos = 1
    for a in arrays
        p1 = pos+(isa(a,AbstractMatrix) ? blocksize(a, 2) : 1)-1
        copyto!(view(dest,:, Block.(pos:p1)), a)
        pos = p1+1
    end
    return dest
end

struct BlockBandedInterlaceLayout <: AbstractLazyBlockBandedLayout end

sublayout(::BlockBandedInterlaceLayout, ::Type{<:NTuple{2,BlockSlices}}) = BlockBandedInterlaceLayout()

arguments(::BlockBandedInterlaceLayout, A::BlockBroadcastMatrix{<:Any,typeof(vcat)}) = A.args

# avoid extra types, since we are using int indexing for now...
# TODO: rewrite when other block sizes are allowed
deblock(A::PseudoBlockArray) = A.blocks
deblock(A::Zeros{T}) where T = Zeros{T}(size(A)...)
function arguments(::BlockBandedInterlaceLayout, A::SubArray)
    P = parent(A)
    args = arguments(BlockBandedInterlaceLayout(), P)
    KR,JR = parentindices(A)
    kr,jr = Int.(KR.block),Int.(JR.block)
    tuple(first(args), view.(deblock.(tail(args)), Ref(kr), Ref(jr))...)
end

# function _copyto!(::BlockBandedColumnMajor, ::BlockBandedInterlaceLayout, dest::AbstractMatrix, src::AbstractMatrix)
#     args = arguments(BlockBandedInterlaceLayout(), src)
#     N,M = blocksize(dest)

# end

blockinterlacelayout(_...) = LazyLayout()
blockinterlacelayout(::Union{ZerosLayout,PaddedLayout,AbstractBandedLayout}...) = BlockBandedInterlaceLayout()

MemoryLayout(::Type{<:BlockBroadcastMatrix{<:Any,typeof(hvcat),Arrays}}) where Arrays = blockinterlacelayout(Base.tail(LazyArrays.tuple_type_memorylayouts(Arrays))...)

# temporary hack, need to think of how to flag as lazy for infinite case.
MemoryLayout(::Type{<:BlockBroadcastMatrix{<:Any,typeof(hcat),Arrays}}) where Arrays = LazyLayout()

MemoryLayout(::Type{<:BlockBroadcastMatrix{<:Any,typeof(Diagonal),Arrays}}) where Arrays = LazyBandedBlockBandedLayout()

##
# special for unitblocks
blockbandwidths(A::PseudoBlockMatrix{<:Any,<:Any,<:NTuple{2,BlockedUnitRange{<:AbstractUnitRange{Int}}}}) = bandwidths(A.blocks)
subblockbandwidths(A::PseudoBlockMatrix{<:Any,<:Any,<:NTuple{2,BlockedUnitRange{<:AbstractUnitRange{Int}}}}) = (0,0)

###
# work around bug in dat
###

LazyArrays._lazy_getindex(dat::PseudoBlockArray, kr::UnitRange) = view(dat.blocks,kr)
LazyArrays._lazy_getindex(dat::PseudoBlockArray, kr::OneTo) = view(dat.blocks,kr)

###
# block col/rowsupport
###
blockcolsupport(M::BlockVcat, j) = first(blockcolsupport(first(M.arrays),j)):(Block(blocksize(BlockVcat(most(M.arrays)...),1))+last(blockcolsupport(last(M.arrays),j)))
blockrowsupport(M::BlockHcat, k) = first(blockrowsupport(first(M.arrays),k)):(Block(blocksize(BlockHcat(most(M.arrays)...),1))+last(blockrowsupport(last(M.arrays),k)))
function blockcolsupport(H::BlockHcat, j::Integer)
    ξ = j
    for A in arguments(H)
        n = blocksize(A,2)
        ξ ≤ n && return blockcolsupport(A, ξ)
        ξ -= n
    end
    return Block.(1:0)
end

function blockrowsupport(H::BlockVcat, k::Integer)
    ξ = k
    for A in arguments(H)
        n = blocksize(A,1)
        ξ ≤ n && return blockrowsupport(A, ξ)
        ξ -= n
    end
    return Block.(1:0)
end

blockcolsupport(H::BlockHcat, J::Block{1}) = blockcolsupport(H, Int(J))
blockrowsupport(H::BlockVcat, K::Block{1}) = blockrowsupport(H, Int(K))

blockcolsupport(A::BlockBroadcastMatrix{<:Any,typeof(hvcat)}, j) = Block.(convexunion(colsupport.(tail(A.args), Ref(Int.(j)))...))
blockrowsupport(A::BlockBroadcastMatrix{<:Any,typeof(hvcat)}, k) = Block.(convexunion(rowsupport.(tail(A.args), Ref(Int.(k)))...))

blockcolsupport(A::BlockBroadcastVector{<:Any,typeof(vcat)}, j) = Block.(convexunion(colsupport.(tail(A.args), Ref(Int.(j)))...))

blockbroadcastlayout(FF, args...) = UnknownLayout()
blockbroadcastlayout(::Type{typeof(vcat)}, ::PaddedLayout...) = PaddedLayout{UnknownLayout}()

function paddeddata(B::BlockBroadcastVector{T,typeof(vcat)}) where T
    dats = map(paddeddata,B.args)
    N = max(map(length,dats)...)
    all(length.(dats) .== N) || error("differening padded lengths not supported")
    BlockBroadcastVector{T}(vcat, dats...)
end

MemoryLayout(::Type{BlockBroadcastArray{T,N,FF,Args}}) where {T,N,FF,Args} = blockbroadcastlayout(FF, tuple_type_memorylayouts(Args)...)

resize!(c::BlockBroadcastVector{T,typeof(vcat)}, N::Block{1}) where T = BlockBroadcastVector{T}(vcat, resize!.(c.args, N)...)


text:
module PlutoVista
using UUIDs
using Colors
using ColorSchemes
using DocStringExtensions
using GridVisualize

abstract type AbstractPlutoVistaBackend end

include("plutovtkplot.jl")
export PlutoVTKPlot
export tricontour!,contour!
export trimesh!
export tetcontour!
export tetmesh!
export quiver2d!

include("plutoplotlyplot.jl")
export plot!,plot
export PlutoPlotlyPlot

# Experimental:
export triupdate!, triplot!
export axis3d!, axis2d!

include("plutocanvasplot.jl")
export polygon!,linecolor!, fillcolor!
export textcolor!,textsize!,text!
export polyline!,linecolor!
export polygon!,fillcolor!
export axis!
export PlutoCanvasPlot

# API
include("plutovistaplot.jl")
export PlutoVistaPlot
export  tricontour,trimesh,contour,tetcontour,tetmesh,quiver2d

# Tools
include("plutoutil.jl")
export ScreenWidthGrabber, PlutoCellWidener

end # module


text:
p_1 = createmodel(:CapitalAdjustModel; nK = 40, nz = 11, γ = 2., intdim = :Separable)
sol_1_prebuild = solve(p_1; rewardcall=:pre)
sol_1_prebuild_partial = solve(p_1; rewardcall=:pre_partial)
sol_1_nobuild = solve(p_1; rewardcall=:jit)

p_2 = createmodel(:CapitalAdjustModel2; nK = 15, nN = 10, nz = 3, intdim=:Separable_ExogStates,)
diopt = Dict(:concavity => [true,true], :monotonicity => [true,true])
sol_2_prebuild = solve(p_2; diopt..., rewardcall=:pre)
sol_2_prebuild_partial = solve(p_2; diopt..., rewardcall=:pre_partial)
sol_2_nobuild = solve(p_2; diopt..., rewardcall=:jit)
sol_1_nobuild

@testset "reward matrix options" begin

	@testset "SingleChoiceVar" begin
		@test isapprox(sol_1_nobuild, sol_1_prebuild_partial, rtol=mytol)
		@test isapprox(sol_1_nobuild, sol_1_prebuild, rtol=mytol)
		@test isapprox(sol_1_prebuild, sol_1_prebuild_partial, rtol=mytol)
	end

	@testset "TwoChoiceVar" begin
		@test isapprox(sol_2_nobuild, sol_2_prebuild_partial, rtol=mytol)
		@test isapprox(sol_2_nobuild, sol_2_prebuild, rtol=mytol)
		@test isapprox(sol_2_prebuild, sol_2_prebuild_partial, rtol=mytol)
	end
end


text:
include("testIdealGas.jl")
include("testHDSim.jl")

max_stars_repo_name: bolverk/Sirion

max_issues_repo_name: bolverk/Sirion

max_forks_repo_name: bolverk/Sirion


text:
function IMPLEMENT_ONE_ARG_FUNC(meth, symnm; lib=:basic_)
     @eval begin
        function ($meth)(b::SymbolicType)
            a = Basic()
            ccall(($(string(lib,symnm)), libsymengine), Void, (Ptr{Basic}, Ptr{Basic}), &a, &b)
            return a
        end
    end
end

function IMPLEMENT_TWO_ARG_FUNC(meth, symnm; lib=:basic_)
     @eval begin
        function ($meth)(b1::SymbolicType, b2::Number)
            a = Basic()
            b1, b2 = promote(b1, b2)
            ccall(($(string(lib,symnm)), libsymengine), Void, (Ptr{Basic}, Ptr{Basic}, Ptr{Basic}), &a, &b1, &b2)
            return a
        end
    end
end

## import from base one argument functions
## these are from cwrapper.cpp, one arg func
for (meth, libnm) in [
                      (:abs,:abs),
                      (:sin,:sin),
                      (:cos,:cos),
                      (:tan,:tan),
                      (:csc,:csc),
                      (:sec,:sec),
                      (:cot,:cot),
                      (:asin,:asin),
                      (:acos,:acos),
                      (:asec,:asec),
                      (:acsc,:acsc),
                      (:atan,:atan),
                      (:acot,:acot),
                      (:sinh,:sinh),
                      (:cosh,:cosh),
                      (:tanh,:tanh),
                      (:csch,:csch),
                      (:sech,:sech),
                      (:coth,:coth),
                      (:asinh,:asinh),
                      (:acosh,:acosh),
                      (:asech,:asech),
                      (:acsch,:acsch),
                      (:atanh,:atanh),
                      (:acoth,:acoth),
                      (:zeta,:zeta),
                      (:gamma,:gamma),
                      (:eta,:dirichlet_eta),
                      (:log,:log),
                      (:sqrt,:sqrt),
                      (:exp,:exp),
                      ]
    eval(Expr(:import, :Base, meth))
    IMPLEMENT_ONE_ARG_FUNC(meth, libnm)
end

# export not import
for  (meth, libnm) in [
                       (:lambertw,:lambertw)   # in add-on packages, not base
                       ]
    IMPLEMENT_ONE_ARG_FUNC(meth, libnm)
    eval(Expr(:export, meth))
end

## add these in until they are wrapped
Base.cbrt(a::SymbolicType) = a^(1//3)

for (meth, fn) in [(:sind, :sin), (:cosd, :cos), (:tand, :tan), (:secd, :sec), (:cscd, :csc), (:cotd, :cot)]
    eval(Expr(:import, :Base, meth))
    @eval begin
        $(meth)(a::SymbolicType) = $(fn)(a*PI/180)
    end
end

## Number theory module from cppwrapper
for (meth, libnm) in [(:gcd, :gcd),
                      (:lcm, :lcm),
                      (:div, :quotient),
                      (:mod, :mod_f),
                      ]
    eval(Expr(:import, :Base, meth))
    IMPLEMENT_TWO_ARG_FUNC(meth, libnm, lib=:ntheory_)
end

Base.binomial(n::Basic, k::Number) = binomial(N(n), N(k))  #ntheory_binomial seems wrong
Base.rem(a::SymbolicType, b::SymbolicType) = a - (a ÷ b) * b
Base.factorial(n::SymbolicType, k) = factorial(N(n), N(k))

## but not (:fibonacci,:fibonacci), (:lucas, :lucas) (Basic type is not the signature)
for (meth, libnm) in [(:nextprime,:nextprime)
                      ]
    IMPLEMENT_ONE_ARG_FUNC(meth, libnm, lib=:ntheory_)
    eval(Expr(:export, meth))
end

function Base.convert{T}(::Type{CVecBasic}, x::Vector{T})
    vec = CVecBasic()
    for i in x
       b::Basic = Basic(i)
       ccall((:vecbasic_push_back, libsymengine), Void, (Ptr{Void}, Ptr{Basic}), vec.ptr, &b)
    end
    return vec
end

Base.convert(::Type{CVecBasic}, x...) = Base.convert(CVecBasic, collect(promote(x...)))

type SymFunction
    name::String
end

SymFunction(s::Symbol) = SymFunction(string(s))

@compat function (f::SymFunction)(x::CVecBasic)
    a = Basic()
    ccall((:function_symbol_set, libsymengine), Void, (Ptr{Basic}, Ptr{Int8}, Ptr{Void}), &a, f.name, x.ptr)
    return a
end

@compat (f::SymFunction){T}(x::Vector{T}) = (f::SymFunction)(convert(CVecBasic, x))
@compat (f::SymFunction)(x...) = (f::SymFunction)(convert(CVecBasic, x...))

macro funs(x...)
    q=Expr(:block)
    if length(x) == 1 && isa(x[1],Expr)
        @assert x[1].head === :tuple "@funs expected a list of symbols"
        x = x[1].args
    end
    for s in x
        @assert isa(s,Symbol) "@funs expected a list of symbols"
        push!(q.args, Expr(:(=), esc(s), Expr(:call, :(SymEngine.SymFunction), Expr(:quote, s))))
    end
    push!(q.args, Expr(:tuple, map(esc, x)...))
    q
end


text:
@testcase "issue 87" begin
    @test 1 == @match 1 begin
        ::Union{Int, String} => 1
        ::Int => 2
    end
end


text:
using CSV
using CairoMakie
using Colors

color_lst = [colorant"#332288" colorant"#555555" colorant"#117733" colorant"#882255"]

f = CSV.File("DimensionScaling.csv")

crlb_trace = f.crlb_trace
random_crlb_trace = f.random_crlb_trace

noto_sans = "../resources/NotoSans-Regular.ttf"

tickfontsize = 28
labelfontsize = 32

f = Figure(resolution = (850,400), font=noto_sans)

ax = Axis(  f[1,1],
            xticklabelsize=tickfontsize,
            yticklabelsize=tickfontsize,
            yticklabelpad=2,
            xlabel="Ambient Dimension", ylabel = "Tr(Σ)",
            xlabelsize=labelfontsize,
            ylabelsize=labelfontsize,
#            yscale = log10,
#            xscale = log10,
            yminorticksvisible = true,
            yminorgridvisible = true,
            yminorticks = IntervalsBetween(10),
            xminorticksvisible = true,
            xminorgridvisible = true,
            xminorticks = IntervalsBetween(10))

scaling = sum(crlb_trace)/length(crlb_trace)
l1 = lines!(ax, 2 .+ (1:length(crlb_trace)), crlb_trace/scaling,        color = color_lst[4], linewidth=2)
l2 = lines!(ax, 2 .+ (1:length(crlb_trace)), random_crlb_trace/scaling, color = color_lst[1], linewidth=2)

xlims!(ax,(0,52))

l1elem = LineElement(color=color_lst[4], linewidth=10)
l2elem = LineElement(color=color_lst[1], linewidth=10)

Legend(f[1,2], [l2elem,l1elem], ["Random Sampling","1D Limit Set"], labelsize=labelfontsize)

colgap!(f.layout,60)

save("ScalingPlot.pdf",f)


text:
###############
# API methods #
###############

"""
    ForwardDiff.hessian(f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, x), check=Val{true}())

Return `H(f)` (i.e. `J(∇(f))`) evaluated at `x`, assuming `f` is called as `f(x)`.

This method assumes that `isa(f(x), Real)`.

Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
"""
function hessian(f, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {T,CHK}
    CHK && checktag(T, f, x)
    ∇f = y -> gradient(f, y, cfg.gradient_config, Val{false}())
    return jacobian(∇f, x, cfg.jacobian_config, Val{false}())
end

"""
    ForwardDiff.hessian!(result::AbstractArray, f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, x), check=Val{true}())

Compute `H(f)` (i.e. `J(∇(f))`) evaluated at `x` and store the result(s) in `result`,
assuming `f` is called as `f(x)`.

This method assumes that `isa(f(x), Real)`.

Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
"""
function hessian!(result::AbstractArray, f, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, x), ::Val{CHK}=Val{true}()) where {T,CHK}
    CHK && checktag(T, f, x)
    ∇f = y -> gradient(f, y, cfg.gradient_config, Val{false}())
    jacobian!(result, ∇f, x, cfg.jacobian_config, Val{false}())
    return result
end

# We use this struct below instead of an
# equivalent closure in order to avoid
# JuliaLang/julia#15276-related performance
# issues. See #316.
mutable struct InnerGradientForHess{R,C,F}
    result::R
    cfg::C
    f::F
end

function (g::InnerGradientForHess)(y, z)
    inner_result = DiffResult(zero(eltype(y)), y)
    gradient!(inner_result, g.f, z, g.cfg.gradient_config, Val{false}())
    g.result = DiffResults.value!(g.result, value(DiffResults.value(inner_result)))
    return y
end

"""
    ForwardDiff.hessian!(result::DiffResult, f, x::AbstractArray, cfg::HessianConfig = HessianConfig(f, result, x), check=Val{true}())

Exactly like `ForwardDiff.hessian!(result::AbstractArray, f, x::AbstractArray, cfg::HessianConfig)`, but
because `isa(result, DiffResult)`, `cfg` is constructed as `HessianConfig(f, result, x)` instead of
`HessianConfig(f, x)`.

Set `check` to `Val{false}()` to disable tag checking. This can lead to perturbation confusion, so should be used with care.
"""
function hessian!(result::DiffResult, f, x::AbstractArray, cfg::HessianConfig{T} = HessianConfig(f, result, x), ::Val{CHK}=Val{true}()) where {T,CHK}
    CHK && checktag(T, f, x)
    ∇f! = InnerGradientForHess(result, cfg, f)
    jacobian!(DiffResults.hessian(result), ∇f!, DiffResults.gradient(result), x, cfg.jacobian_config, Val{false}())
    return ∇f!.result
end

hessian(f, x::StaticArray) = jacobian(y -> gradient(f, y), x)
hessian(f, x::StaticArray, cfg::HessianConfig) = hessian(f, x)
hessian(f, x::StaticArray, cfg::HessianConfig, ::Val) = hessian(f, x)

hessian!(result::AbstractArray, f, x::StaticArray) = jacobian!(result, y -> gradient(f, y), x)

hessian!(result::MutableDiffResult, f, x::StaticArray) = hessian!(result, f, x, HessianConfig(f, result, x))

hessian!(result::ImmutableDiffResult, f, x::StaticArray, cfg::HessianConfig) = hessian!(result, f, x)
hessian!(result::ImmutableDiffResult, f, x::StaticArray, cfg::HessianConfig, ::Val) = hessian!(result, f, x)

function hessian!(result::ImmutableDiffResult, f, x::StaticArray)
    T = typeof(Tag(f, eltype(x)))
    d1 = dualize(T, x)
    d2 = dualize(T, d1)
    fd2 = f(d2)
    val = value(T,value(T,fd2))
    grad = extract_gradient(T,value(T,fd2), x)
    hess = extract_jacobian(T,partials(T,fd2), x)
    result = DiffResults.hessian!(result, hess)
    result = DiffResults.gradient!(result, grad)
    result = DiffResults.value!(result, val)
    return result
end


text:
using DynamicalSystemsBase
using LinearAlgebra:norm
using Test

println("\nTesting continuous system evolution...")

@testset "Lorenz63 oop" begin

    lo11 = Systems.lorenz() #with Jac
    ts1 = trajectory(lo11, 1.0, Δt = 0.01)
    @test size(ts1) == (101, 3)
    @test ts1[1, :] == SVector{3}(lo11.u0)
    @test !any(x -> abs(x) > 1e3, ts1[end])

    ts1 = trajectory(lo11, 1.0; Δt=0.1, diffeq = (abstol=1e-8, reltol=1e-8))
    @test size(ts1) == (11, 3)
    @test ts1[1, :] == SVector{3}(lo11.u0)

    data2 = trajectory(lo11, 100; save_idxs = 1:2)
    @test size(data2)[2] == 2
end

@testset "Lorenz96 iin" begin
    u = ones(5)
    lo11 = Systems.lorenz96(5, u)
    ts1 = trajectory(lo11, 2.0, Δt = 0.01)
    @test size(ts1) == (201, 5)
    @test ts1[1, :] == SVector{5}(u)

    ts1 = trajectory(lo11, 2.0; Δt=0.1, diffeq = (abstol=1e-8, reltol=1e-8))
    @test size(ts1) == (21, 5)
    @test ts1[1, :] == SVector{5}(u)
    @test !any(x -> abs(x) > 1e3, ts1[end])
end

println("\nTesting stroboscopic map...")

@testset "Duffing stroboscopic map" begin
    F = 0.27; ω = 0.1;  # smooth boundary
    ds = Systems.duffing(zeros(2), ω = ω, f = F, d = 0.15, β = -1)
    smap = stroboscopicmap(ds, 2*pi/ω)
    reinit!(smap,[1., 1.])
    for j in 1:100
      step!(smap)
    end
    u = get_state(smap)
    @test abs(sum(u - [1.11, 0])) < 0.01
end


text:
# extend to vectors
function CommonSolve.solve(eqns::Vector{T}, args...; kwargs...) where {T <: SymbolicObject}
    ↑((Py(sympy).solve)(↓(eqns), ↓(args)...; kwargs...))
end

function dsolve(eqns::Vector{T}, args...; kwargs...) where {T <: SymbolicObject}
    ↑((Py(sympy).dsolve)(↓(eqns), ↓(args)...; kwargs...))
end

function nsolve(eqns::Vector{T}, args...; kwargs...) where {T <: SymbolicObject}
    ↑((Py(sympy).nsolve)(↓(eqns), ↓(args)...; kwargs...))
end

function linsolve(eqns::AbstractArray{T}, args...; kwargs...) where {T <: SymbolicObject}
    ↑((Py(sympy).linsolve)(↓(eqns), ↓(args)...; kwargs...))
end
function linsolve(eqns::NTuple{N, T}, args...; kwargs...) where {N, T <: SymbolicObject}
    ↑((Py(sympy).linsolve)(↓(eqns), ↓(args)...; kwargs...))
end

function nonlinsolve(eqns::Vector{T}, args...; kwargs...) where {T <: SymbolicObject}
    ↑((Py(sympy).nonlinsolve)(↓(eqns), ↓(args)...; kwargs...))
end

function solveset(eqns::Vector{T}, args...; kwargs...) where {T <: SymbolicObject}
    ↑((Py(sympy).solveset)(↓(eqns), ↓(args)...; kwargs...))
end

# disambiguate
SpecialFunctions.beta(x::Sym, y::Sym) = sympy.beta(x, y)
Base.log(b::Sym, x::Sym) = log(x)/log(b)

# special exports
export limit, subs,
    lhs, rhs, doit

export free_symbols
export elements

function limit(ex::Sym, xc::Pair; kwargs...)
    ↑((Py(sympy).limit)(↓(ex), ↓(first(xc)), ↓(last(xc)); kwargs...))
end

function subs(ex::T, xs...; kwargs...) where {T <: SymbolicObject}
    fs = free_symbols(ex)
    length(fs) == length(xs) || throw(ArgumentError("free symbols mismatch"))
    subs(ex, (k=>v for (k,v) ∈ zip(fs, xs))...)
end

function subs(ex::T, d::Pair...; kwargs...)  where {T <: SymbolicObject}
    isnothing(pygetattr(getpy(ex), "subs", nothing)) && return ex
    ↑(ex.py.subs(↓([(first(p), last(p)) for p in d])))
end

function subs(ex::T, d::Dict; kwargs...)  where {T <: SymbolicObject}
    for (k,v) ∈ d
        ex = ex.subs(k, v)
    end
    ex
end

# this one is in a test, but ...
subs(ex::T, y::Tuple{Any, Any}; kwargs...)          where {T <: SymbolicObject} = ex.subs(y[1], Sym(y[2]); kwargs...)
subs(ex::T, y::Tuple{Any, Any}, args...; kwargs...) where {T <: SymbolicObject} = subs(subs(ex, y), args...; kwargs...)

## curried versions to use with |>
#subs(x::SymbolicObject, y; kwargs...) = ex -> subs(ex, x, y; kwargs...)
subs(;kwargs...)                      = ex -> subs(ex; kwargs...)
subs(dict::Dict; kwargs...)           = ex -> subs(ex, dict...; kwargs...)
subs(d::Pair...; kwargs...)           = ex -> subs(ex, [(p.first, p.second) for p in d]...; kwargs...)

lhs(ex::SymbolicObject) = ↑(pygetattr(ex, "lhs", ↓(ex)))
rhs(ex::SymbolicObject) = ↑(pygetattr(ex, "rhs", ↓(ex)))
doit(ex::T; deep::Bool=false) where {T<:SymbolicObject} = ex.doit(deep=deep)
doit(; deep::Bool=false)      where {T<:SymbolicObject} = ex -> doit(ex, deep=deep)

rewrite(x::Sym, args...; kwargs...) = ↑(getpy(x).rewrite(↓(args)...; kwargs...))

# extend to non-sym
simplify(x, args...; kwargs...) = x

# N is tricky
function N(x::Sym)
    length(free_symbols(x)) > 0 && return x # really

    u = getpy(x)
    if is_(:real, x) || is_(:Float, x)
        is_(:infinite, x) && return is_(:negative, x) ? -Inf : Inf

        if is_(:integer, x)
            return pyconvert(Integer, u)
#            is_(:zero, x) && return 0
#            T = Bool(getpy(abs(x)) < typemax(Int)) ? Int : BigInt
#            return pyconvert(T, getpy(x))

        elseif is_(:rational, x)
            return N(numerator(x)) // N(denominator(x))
        else
            # Floating point
            is_(:zero, x) && return 0.0

            _prec = pygetattr(x, "_prec", nothing)
            if !isnothing(_prec)
                prec = pyconvert(Int, _prec)
                prec <= 64 && return pyconvert(Float64, x)
                return setprecision(BigFloat, prec) do
                    pyconvert(BigFloat, x)
                end
            else
                x == PI && return π
                x == sympy.E && return ℯ
                # _from_mpmath   ??
                return N(u.evalf(16))
                #error("not _prec?")
            end
        end
    elseif is_(:imaginary, x)
        return N(real(x)) * im
    elseif is_(:complex, x)

        is_(:infinite, x) && return complex(Inf)
        return complex(N(real(x)), N(imag(x)))

    else
        is_(:boolean, x) && return pyconvert(Bool, getpy(x))
        # ...
        # check for Python object int, float, complex, mpfm mpc ,...
        u = getpy(x)
        pyisinstance(u, pybuiltins.int) &&
            return pyconvert(Bool(u.bit_length() > 64) ? BigInt : Int, u)
        if pyisinstance(u, pybuiltins.float)
            prec = pyconvert(Int, pygetattr(u, "_prec", 64))
            prec <= 64 && return pyconvert(Float64, u)
            return setprecision(BigFloat, prec) do
                pyconvert(BigFloat, u)
            end
        end
        pyisinstance(u, pybuiltins.complex) &&
            return pyconvert(Complex, u)
        if pygetattr(u, "_numerator", nothing) !== nothing
            return N(↑(pygetattr(u, "_numerator"))) // N(↑(pygetattr(u, "_denominator")))
        end

    end

    _evalf = pygetattr(u, "evalf", nothing)
    !isnothing(_evalf) && return N(u.evalf(16))

    throw(ArgumentError("Could not convert $x"))
end
N(x) = N(Sym(x))
N(x::Sym, prec) = N(x.evalf(prec))

function Base.Bool(x::Sym)
    Bool(getpy(x) == pybuiltins.True) && return true
    Bool(getpy(x) == pybuiltins.False) && return false
    throw(ArgumentError("x is not boolean"))
end

##
function free_symbols(x::Sym)
    u = pygetattr(getpy(x), "free_symbols", nothing)
    u == nothing && return Vector{Sym}[]

    fs = [Sym(uᵢ) for uᵢ ∈ PythonCall.pyconvert(PySet, u)]
    fs[sortperm(string.(fs))]   # some sorting order to rely on
end

free_symbols(x::Vector{T}) where {T <: Sym} = reduce(union, free_symbols.(x))

## setes
Base.contains(s::Sym, val) = Bool(↓(s.contains(val)))

# # Matrices
# # call SymMatrix method on Matrix{Sym}
# ## Eg. A.norm() where A = [x 1; 1 x], say
# function Base.getproperty(A::AbstractArray{T}, k::Symbol) where {T <: SymbolicObject}
#     if k in fieldnames(typeof(A))
#         return getfield(A,k)
#     else
#         u = ↓(A)
#         val = pygetattr(u, string(k), nothing)
#         if pycallable(val)
#             return SymbolicCallable(val)
#         end
#         val
#     end
# end

# Base.getindex(M::Sym, i::Int, j::Int) = M.py.__getitem__((i,j))
# Base.getindex(V::Sym, i::Int) = V.py.__getitem__((i,0))
# Base.getindex(M::Sym) = M

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall