text:
@testset "Jacobian" begin
    @testset "LineSegment" begin
        for T ∈ Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(1, 0, 0)
            l = LineSegment(P₁, P₂)
            @test jacobian(l, 1) == [1,0,0]
            @test jacobian(l, 0) == [1,0,0]
        end
    end

    @testset "QuadraticSegment" begin
        for T ∈ Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(1, 0, 0)
            P₃ = Point{3,T}(1//2, 0, 0)
            q = QuadraticSegment(P₁, P₂, P₃)
            @test jacobian(q, 0)    == [1,0,0]
            @test jacobian(q, 1//2) == [1,0,0]
            @test jacobian(q, 1)    == [1,0,0]
        end
    end

    @testset "Triangle" begin
        for T ∈ Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(1, 0, 0)
            P₃ = Point{3,T}(0, 1, 0)
            tri = Triangle(P₁, P₂, P₃)
            @test jacobian(tri, 0, 0) == [1 0; 0 1; 0 0]
            @test jacobian(tri, 1, 0) == [1 0; 0 1; 0 0]
            @test jacobian(tri, 0, 1) == [1 0; 0 1; 0 0]
        end
    end

    @testset "Quadrilateral" begin
        for T ∈ Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(1, 0, 0)
            P₃ = Point{3,T}(1, 1, 0)
            P₄ = Point{3,T}(0, 1, 0)
            quad = Quadrilateral(P₁, P₂, P₃, P₄)
            @test jacobian(quad, 0, 0) == [1 0; 0 1; 0 0]
            @test jacobian(quad, 1, 0) == [1 0; 0 1; 0 0]
            @test jacobian(quad, 1, 1) == [1 0; 0 1; 0 0]
            @test jacobian(quad, 0, 1) == [1 0; 0 1; 0 0]
        end
    end

    @testset "Quadratic Triangle" begin
        for T ∈ Floats
            P₁ = Point{3,T}(   0,    0, 0)
            P₂ = Point{3,T}(   1,    0, 0)
            P₃ = Point{3,T}(   0,    1, 0)
            P₄ = Point{3,T}(1//2,    0, 0)
            P₅ = Point{3,T}(1//2, 1//2, 0)
            P₆ = Point{3,T}(   0, 1//2, 0)
            tri6 = QuadraticTriangle(P₁, P₂, P₃, P₄, P₅, P₆)
            @test jacobian(tri6,    0,    0) == [1 0; 0 1; 0 0]
            @test jacobian(tri6,    1,    0) == [1 0; 0 1; 0 0]
            @test jacobian(tri6,    0,    1) == [1 0; 0 1; 0 0]
            @test jacobian(tri6, 1//2,    0) == [1 0; 0 1; 0 0]
            @test jacobian(tri6,    0, 1//2) == [1 0; 0 1; 0 0]
            @test jacobian(tri6, 1//2, 1//2) == [1 0; 0 1; 0 0]
        end
    end

    @testset "Quadratic Quadrilateral" begin
        for T ∈ Floats
            P₁ = Point{3,T}(   0,    0, 0)
            P₂ = Point{3,T}(   1,    0, 0)
            P₃ = Point{3,T}(   1,    1, 0)
            P₄ = Point{3,T}(   0,    1, 0)
            P₅ = Point{3,T}(1//2,    0, 0)
            P₆ = Point{3,T}(   1, 1//2, 0)
            P₇ = Point{3,T}(1//2,    1, 0)
            P₈ = Point{3,T}(   0, 1//2, 0)
            quad8 = QuadraticQuadrilateral(P₁, P₂, P₃, P₄, P₅, P₆, P₇, P₈)
            @test jacobian(quad8,    0,    0) == [1 0; 0 1; 0 0]
            @test jacobian(quad8,    1,    0) == [1 0; 0 1; 0 0]
            @test jacobian(quad8,    1,    1) == [1 0; 0 1; 0 0]
            @test jacobian(quad8,    0,    1) == [1 0; 0 1; 0 0]
            @test jacobian(quad8, 1//2,    0) == [1 0; 0 1; 0 0]
            @test jacobian(quad8,    1, 1//2) == [1 0; 0 1; 0 0]
            @test jacobian(quad8, 1//2,    1) == [1 0; 0 1; 0 0]
            @test jacobian(quad8,    0, 1//2) == [1 0; 0 1; 0 0]
        end
    end

    @testset "Tetrahedron" begin
        for T ∈ Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(1, 0, 0)
            P₃ = Point{3,T}(0, 1, 0)
            P₄ = Point{3,T}(0, 0, 1)
            tet = Tetrahedron(P₁, P₂, P₃, P₄)
            @test jacobian(tet, 0, 0, 0) == I
            @test jacobian(tet, 1, 0, 0) == I
            @test jacobian(tet, 0, 1, 0) == I
            @test jacobian(tet, 0, 0, 1) == I
        end
    end

    @testset "Hexahedron" begin
        for T ∈ Floats
            P₁ = Point{3,T}(0, 0, 0)
            P₂ = Point{3,T}(1, 0, 0)
            P₃ = Point{3,T}(1, 1, 0)
            P₄ = Point{3,T}(0, 1, 0)
            P₅ = Point{3,T}(0, 0, 1)
            P₆ = Point{3,T}(1, 0, 1)
            P₇ = Point{3,T}(1, 1, 1)
            P₈ = Point{3,T}(0, 1, 1)
            hex = Hexahedron(P₁, P₂, P₃, P₄, P₅, P₆, P₇, P₈)
            @test jacobian(hex, 0, 0, 0) == I
            @test jacobian(hex, 1, 0, 0) == I
            @test jacobian(hex, 1, 1, 0) == I
            @test jacobian(hex, 0, 1, 0) == I
            @test jacobian(hex, 0, 0, 1) == I
            @test jacobian(hex, 1, 0, 1) == I
            @test jacobian(hex, 1, 1, 1) == I
            @test jacobian(hex, 0, 1, 1) == I
        end
    end

    @testset "QuadraticTetrahedron" begin
        for T ∈ Floats
            P₁  = Point{3,T}(0, 0, 0)
            P₂  = Point{3,T}(1, 0, 0)
            P₃  = Point{3,T}(0, 1, 0)
            P₄  = Point{3,T}(0, 0, 1)
            P₅  = Point{3,T}(1//2,    0,    0)
            P₆  = Point{3,T}(1//2, 1//2,    0)
            P₇  = Point{3,T}(0,    1//2,    0)
            P₈  = Point{3,T}(0,       0, 1//2)
            P₉  = Point{3,T}(1//2,    0, 1//2)
            P₁₀ = Point{3,T}(0,    1//2, 1//2)
            tet = QuadraticTetrahedron(P₁, P₂, P₃, P₄, P₅, P₆, P₇, P₈, P₉, P₁₀)
            @test jacobian(tet, 0, 0, 0) == I
            @test jacobian(tet, 1, 0, 0) == I
            @test jacobian(tet, 0, 1, 0) == I
            @test jacobian(tet, 0, 0, 1) == I
            @test jacobian(tet, 1//2,    0,    0) == I
            @test jacobian(tet, 1//2, 1//2,    0) == I
            @test jacobian(tet, 0,    1//2,    0) == I
            @test jacobian(tet, 0,       0, 1//2) == I
            @test jacobian(tet, 1//2,    0, 1//2) == I
            @test jacobian(tet, 0,    1//2, 1//2) == I
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
@testset "code_transformation/util" begin

    # Test Nabla.unionise_type.
    @test Nabla.unionise_type(:Real) == :(Union{Real, Node{<:Real}})
    @test Nabla.unionise_type(:(T{V})) == :(Union{T{V}, Node{<:T{V}}})
    @test Nabla.unionise_type(:(Vararg)) == :(Vararg{Union{Any, Node{<:Any}}})
    @test Nabla.unionise_type(:(Vararg{T})) == :(Vararg{Union{T, Node{<:T}}})
    @test Nabla.unionise_type(:(T{V} where V)) ==
        :(Union{T{V} where V, Node{<:(T{V} where V)}})
    @test Nabla.unionise_type(:(Vararg{T} where T)) ==
        :(Vararg{Union{T where T, Node{<:(T where T)}}})
    @test Nabla.unionise_type(:(Vararg{<:T})) == :(Vararg{Union{<:T, Node{<:T}}})

    # Test Nabla.replace_body.
    @test Nabla.replace_body(:Real, :Float64) == :Float64
    @test Nabla.replace_body(:Real, :(Union{Float64, foo})) == :(Union{Float64, foo})
    @test Nabla.replace_body(:(T where T), :(Node{T})) == :(Node{T} where T)
    @test Nabla.replace_body(:(T where T <: A{V} where V), :(Node{T})) ==
        :(Node{T} where T <: A{V} where V)

    # Test Nabla.get_body.
    @test Nabla.get_body(:Real) == :Real
    @test Nabla.get_body(:(Node{T})) == :(Node{T})
    @test Nabla.get_body(:(Node{T} where T)) == :(Node{T})
    @test Nabla.get_body(:(Node{T} where T<:A{Q} where Q)) == :(Node{T})
    @test Nabla.get_body(:(Node{A{T} where T})) == :(Node{A{T} where T})
    @test Nabla.get_body(:(Node{A{T} where T} where A)) == :(Node{A{T} where T})

    # Test Nabla.get_types.
    @test Nabla.get_types(:(Tuple{T})) == [:T]
    @test Nabla.get_types(:(Tuple{T, T})) == [:T, :T]
    @test Nabla.get_types(:(Tuple{T{V}})) == [:(T{V})]
    @test Nabla.get_types(:(Tuple{T{<:V}})) == [:(T{<:V})]

    # Test Nabla.isa_vararg.
    @test Nabla.isa_vararg(:foo) == false
    @test Nabla.isa_vararg(:Vararg) == true
    @test Nabla.isa_vararg(:(Node{T} where T)) == false
    @test Nabla.isa_vararg(:(Vararg{T} where T)) == true
    @test Nabla.isa_vararg(:(Vararg{T, N} where T where N)) == true

    # Test Nabla.remove_vararg.
    @test Nabla.remove_vararg(:Vararg) == (:Any, :Vararg)
    @test Nabla.remove_vararg(:(Node{T})) == (:(Node{T}), :nothing)
    @test Nabla.remove_vararg.([:Vararg]) == [(:Any, :Vararg)]
    @test Nabla.remove_vararg.(:Vararg) == (:Any, :Vararg)
    @test Nabla.remove_vararg.([:Vararg, :Vararg]) == [(:Any, :Vararg), (:Any, :Vararg)]
    @test Nabla.remove_vararg(:Real) == (:Real, :nothing)
    @test Nabla.remove_vararg(:(Vararg{T} where T)) == (:(T where T), :Vararg)
    @test Nabla.remove_vararg(:(Vararg{T, N} where T<:Real)) == (:(T where T<:Real), :N)
    # Redundant local `where N` (rather than leaving the `N` to be outside in the sig)
    @test Nabla.remove_vararg(:(Vararg{Real, N} where N)) == (:Real, :Vararg)
    @test Nabla.remove_vararg(:(Vararg{T, N} where {N,T})) == (:(T where T), :Vararg)
    @test Nabla.remove_vararg(:(Vararg{T, N} where {T<:Real, N})) == (:(T where T<:Real), :Vararg)
    # This case doesn't work but never occurs in practice
    @test_broken Nabla.remove_vararg(:(Vararg{T, N} where T where N)) == (:(T where T), :Vararg)

    # Test Nabla.replace_vararg.
    @test Nabla.replace_vararg(:(U{T, N{T}}), (:V, :nothing)) == :(U{T, N{T}})
    @test Nabla.replace_vararg(:(Real), (:Any, :Vararg)) == :(Vararg{Real})
    @test Nabla.replace_vararg(:(T where T), (:T, :Vararg)) == :(Vararg{T} where T)
    @test Nabla.replace_vararg(:(T where T), (:T, :2)) == :(Vararg{T, 2} where T)
    @test Nabla.replace_vararg(:(U{T} where T), (:T, :N)) == :(Vararg{U{T}, N} where T)
    @test Nabla.replace_vararg(:(U{T} where N), (:T, :N)) == :(Vararg{U{T}, N} where N)

    @testset "parse_kwargs" begin
        @test Nabla.parse_kwargs(:(())) == NamedTuple()
        @test Nabla.parse_kwargs(:(NamedTuple())) == NamedTuple()
        @test Nabla.parse_kwargs(:((a = 1, b = 2))) == NamedTuple{(:a, :b)}((1, 2))
        @test Nabla.parse_kwargs(:((; a = 1, b = 2))) == NamedTuple{(:a, :b)}((1, 2))
        @test Nabla.parse_kwargs(:((a = sum(1:10), b = :c))) ==
            NamedTuple{(:a, :b)}((:(sum(1:10)), :(:c)))
        @test Nabla.parse_kwargs(:((; a = sum(1:10), b = :c))) ==
            NamedTuple{(:a, :b)}((:(sum(1:10)), :(:c)))
        @test_throws ArgumentError Nabla.parse_kwargs(:([a => 2]))
    end

    @testset "parse_is_node" begin
        @test Nabla.parse_is_node(:([])) == Bool[]
        @test Nabla.parse_is_node(:([true])) == [true]
        @test Nabla.parse_is_node(:([true, false])) == [true, false]
        @test_throws ArgumentError Nabla.parse_is_node(:((true, false)))
    end

    @testset "node_type" begin
        # special case for a redudant local where N in a Vararg
        @test Nabla.node_type(:(Vararg{Int64, N} where N)) == :(Vararg{Node{<:Int64}})

        @test Nabla.node_type(:Float32) == :(Node{<:Float32})
    end
end

ext: jl


text:
import Pkg
Pkg.activate(normpath(@__DIR__))

using ProjectMoana
import Images: save
import Base: write
import ColorTypes: RGB

using ArgParse

function parse_commandline()
    s = ArgParseSettings(description = "This program generates an image of 10 spheres. Try me!",
                               usage = "usage: [--help] [--width W] [--height H] [--camera C] [--angle α] [--distance D] [--file_out FILENAME] [--render_alg ALG] [--a A]",
                              epilog = "Let's try again!")

    @add_arg_table s begin
        "--width"
            help = "width of the image"
            required = true
            arg_type = Int
        "--height"
            help = "height of the image"
            required = true
            arg_type = Int
        "--camera"
            help = "type of the camera [O, for Orthogonal, or P, for Perspective]"
            required = false
            default = "O"
            arg_type = String
        "--angle"
            help = "angle of the z-axis rotation applied to camera (in degrees)"
            required = false
            default = 0.
            arg_type = Float64
        "--distance"
            help = "distance of a Perspective Camera"
            required = false
            default = 1.
            arg_type = Float64
        "--file_out"
            help = "name of the output file (without extension)"
            required = false
            default = "demo"
            arg_type = String
        "--render_alg"
            help = "type of algortihm to use for render [on_off, flat, path_tracer]"
            required = false
            default = "flat"
            arg_type = String
        "--a"
            help = "a_factor for nomralize image luminosity during the convertion"
            required = false
            default = 1.
            arg_type = Float64
    end

    return parse_args(s)
end

function main()
    params = parse_commandline()

    w = params["width"]
    h = params["height"]
    a = w/h
    d = params["distance"]
    camera_tr = rotation_z(params["angle"]*π/180.0) * translation(Vec(-1.0,0.,0.))
    image = HdrImage(w, h)
    file_out_pfm = "$(params["file_out"]).pfm"
    file_out_png = "$(params["file_out"]).png"
    algorithm = params["render_alg"]

# Creating WORLD with 10 spheres
    material1 = Material(DiffuseBRDF(UniformPigment(RGB(0.7, 0.3, 0.2))))
    material2 = Material(DiffuseBRDF(CheckeredPigment(RGB(0.2, 0.7, 0.3), RGB(0.3, 0.2, 0.7), 8)))

    sphere_texture = HdrImage(2, 2)
    set_pixel(sphere_texture, 1, 1, RGB(0.1, 0.2, 0.3))
    set_pixel(sphere_texture, 1, 2, RGB(0.2, 0.1, 0.3))
    set_pixel(sphere_texture, 2, 1, RGB(0.3, 0.2, 0.1))
    set_pixel(sphere_texture, 2, 2, RGB(0.3, 0.1, 0.2))

    material3 = Material(DiffuseBRDF(ImagePigment(sphere_texture)))
    world = World()

    for x in [-0.5, 0.5]
        for y in [-0.5, 0.5]
            for z in [-0.5, 0.5]
               add_shape(world, Sphere(translation(Vec(x, y, z)) * scaling(Vec(0.1, 0.1, 0.1)), material1 ))
            end
        end
    end

   add_shape(world, Sphere(translation(Vec( 0.0, 0.5, 0.0)) * scaling(Vec(0.1, 0.1, 0.1)), material2 ))
   add_shape(world, Sphere(translation(Vec( -0.5, 0.0,-0.5)) * scaling(Vec(.1,.1,.1)), material3 ))

   println("World objects created.")

# Creating a Perspective of Orthogonal CAMERA
    if params["camera"] == "O"
        camera = OrthogonalCamera(a, camera_tr)
    elseif params["camera"] == "P"
        camera = PerspectiveCamera(a, camera_tr, d)
    end

# Creating an ImageTracer object
    tracer = ImageTracer(image, camera)

    println("Observer initialized.")

    print("Computing ray intersection ")
    if params["render_alg"] == "flat"
        println("using Flat renderer")
        renderer = Flat_Renderer(world, RGB(0.4,0.4,0.4))
        fire_all_rays(tracer, Flat, renderer)
    else
        println("using On/Off renderer")
        renderer = OnOff_Renderer(world)
        fire_all_rays(tracer, OnOff, renderer)
    end

    println("Ray intersections evaluated.")

# Saving the PFM FILE
    write(file_out_pfm, tracer.image)

    println("$(file_out_pfm) has been written to disk.")

# Automatic CONVERSION TO JPEG FILE
    # normalize_image(tracer.image, params["a"])
    # clamp_image(tracer.image)

    # matrix_pixels = reshape(tracer.image.pixels, (tracer.image.width, tracer.image.height))

    # save(file_out_png, matrix_pixels')

    # println("$(file_out_png) has been automatically written to disk.")

end

main()

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


text:
using FFTW
using Printf
using HDF5

struct Pressure{TF <: Union{Float32, Float64}}
    fft_forward
    fft_backward
    bmati::Vector{TF}
    bmatj::Vector{TF}
    a::Vector{TF}
    c::Vector{TF}
    work3d::Array{TF, 3}
    work2d::Array{TF, 2}
    b::Array{TF, 3}
    p_nogc::Array{TF, 3}
end

function Pressure(g::Grid, TF)
    nthreads = (Threads.nthreads() == 1) ? 1 : 2*Threads.nthreads()
    FFTW.set_num_threads(nthreads)

    # We set the type of rand here explictly to trigger the right precision in FFTW
    tmp = rand(TF, g.itot, g.jtot, g.ktot)
    fft_plan_f = FFTW.plan_r2r(tmp, FFTW.R2HC, [1, 2], flags=FFTW.MEASURE)
    fft_plan_b = FFTW.plan_r2r(tmp, FFTW.HC2R, [1, 2], flags=FFTW.MEASURE)

    bmati = zeros(g.itot)
    bmatj = zeros(g.jtot)
    a = zeros(g.ktot)
    c = zeros(g.ktot)

    dxidxi = g.dxi^2
    dyidyi = g.dyi^2

    for j in 0:g.jtot÷2
        bmatj[j+1] = 2 * (cos(2pi*j/g.jtot) - 1) * dyidyi;
    end

    for j in g.jtot÷2+1:g.jtot-1
        bmatj[j+1] = bmatj[g.jtot-j+1];
    end

    for i in 0:g.itot÷2
        bmati[i+1] = 2 * (cos(2pi*i/g.itot) - 1) * dxidxi;
    end

    for i in g.itot÷2+1:g.itot-1
        bmati[i+1] = bmati[g.itot-i+1];
    end

    for k in 1:g.ktot
        a[k] = g.dz[k+g.kgc] * g.dzhi[k+g.kgc  ];
        c[k] = g.dz[k+g.kgc] * g.dzhi[k+g.kgc+1];
    end

    work3d = zeros(g.itot, g.jtot, g.ktot)
    work2d = zeros(g.itot, g.jtot)
    b = zeros(g.itot, g.jtot, g.ktot)
    p_nogc = zeros(g.itot, g.jtot, g.ktot)

    Pressure{TF}(fft_plan_f, fft_plan_b, bmati, bmatj, a, c, work3d, work2d, b, p_nogc)
end

function input_kernel!(
    p,
    u, v, w,
    ut, vt, wt,
    dxi, dyi, dzi, dti,
    is, ie, js, je, ks, ke)

    @fast3d begin
        p[i-is+1, j-js+1, k-ks+1] = @fd (u, v, w, ut, vt, wt) gradx(ut + u*dti) + grady(vt + v*dti) + gradz(wt + w*dti)
    end
end

function solve_pre_kernel!(
    p, b,
    dz, bmati, bmatj, a, c,
    itot, jtot, ktot, kgc)

    @tturbo unroll=8 for k in 1:ktot
        for j in 1:jtot
            for i in 1:itot
                b[i, j, k] = dz[k+kgc]^2 * (bmati[i]+bmatj[j]) - (a[k]+c[k]);
                p[i, j, k] *= dz[k+kgc]^2
            end
        end
    end

    # Set the BCs of all wave numbers to Neumann = 0.
    @tturbo for j in 1:jtot
        for i in 1:itot
            b[i, j, 1] += a[1]
            b[i, j, ktot] += c[ktot]
        end
    end

    # Set the wave number 0 to Dirichlet = 0.
    b[1, 1, ktot] -= 2c[ktot]
end

function solve_tdma_kernel!(
    p, work3d, work2d,
    a, b, c,
    itot, jtot, ktot)

    @tturbo for j in 1:jtot
        for i in 1:itot
            work2d[i, j] = b[i, j, 1]
            p[i, j, 1] /= work2d[i, j]
        end
    end

    for k in 2:ktot
        @tturbo for j in 1:jtot
            for i in 1:itot
                work3d[i, j, k] = c[k-1] / work2d[i, j]
                work2d[i, j] = b[i, j, k] - a[k]*work3d[i, j, k]
                p[i, j, k] -= a[k] * p[i, j, k-1]
                p[i, j, k] /= work2d[i, j]
            end
        end
    end

    for k in ktot-1:-1:1
        @tturbo for j in 1:jtot
            for i in 1:itot
                p[i, j, k] -= work3d[i, j, k+1] * p[i, j, k+1]
            end
        end
    end
end

function solve_post_kernel!(
    p, p_nogc,
    is, ie, js, je, ks, ke,
    igc, jgc, kgc)

    @fast3d begin
        p[i, j, k] = p_nogc[i-igc, j-jgc, k-kgc]
    end
end

function output_kernel!(
    ut, vt, wt,
    p,
    dxi, dyi, dzhi,
    is, ie, js, je, ks, ke)

    @fast3d begin
        @fd (ut, p) ut -= gradx(p)
        @fd (vt, p) vt -= grady(p)
        @fd (wt, p) wt -= gradz(p)
    end
end

function calc_pressure_tend!(f::Fields, g::Grid, t::Timeloop, p::Pressure)
    # Set the cyclic boundaries for the tendencies.
    boundary_cyclic_kernel!(
        f.u_tend, g.is, g.ie, g.js, g.je, g.igc, g.jgc)
    boundary_cyclic_kernel!(
        f.v_tend, g.is, g.ie, g.js, g.je, g.igc, g.jgc)

    # Set the boundaries of wtend to zero.
    # CvH: Fix this later.
    @tturbo for j in 1:g.jcells
        for i in 1:g.icells
            f.w_tend[i, j, g.ks ] = 0
            f.w_tend[i, j, g.keh] = 0
        end
    end

    # Convert the dti to the type of the arrays to prevent casting in kernel.
    dti_sub = convert(typeof(f.u_tend[1]), 1/get_sub_dt(t))

    input_kernel!(
        p.p_nogc,
        f.u, f.v, f.w,
        f.u_tend, f.v_tend, f.w_tend,
        g.dxi, g.dyi, g.dzi, dti_sub,
        g.is, g.ie, g.js, g.je, g.ks, g.ke)

    p_fft = p.fft_forward * p.p_nogc

    solve_pre_kernel!(
        p_fft, p.b,
        g.dz, p.bmati, p.bmatj, p.a, p.c,
        g.itot, g.jtot, g.ktot, g.kgc)

    solve_tdma_kernel!(
        p_fft, p.work3d, p.work2d,
        p.a, p.b, p.c,
        g.itot, g.jtot, g.ktot)

    p.p_nogc[:, :, :] = (p.fft_backward * p_fft) ./ (g.itot * g.jtot)

    solve_post_kernel!(
        f.p, p.p_nogc,
        g.is, g.ie, g.js, g.je, g.ks, g.ke,
        g.igc, g.jgc, g.kgc)

    # Set the bot and top bc's
    @tturbo for j in 1:g.jcells
        for i in 1:g.icells
            f.p[i, j, g.ks-1] = f.p[i, j, g.ks]
            f.p[i, j, g.ke+1] = f.p[i, j, g.ke]
        end
    end

    boundary_cyclic_kernel!(
        f.p, g.is, g.ie, g.js, g.je, g.igc, g.jgc)

    output_kernel!(
        f.u_tend, f.v_tend, f.w_tend,
        f.p,
        g.dxi, g.dyi, g.dzhi,
        g.is, g.ie, g.js, g.je, g.ks, g.ke)
end


text:
using Merge
using Base.Test

a = [1, 2, 5, 7, 8, 11]
b = [1, 4, 5, 6, 9, 10, 12, 16]
expected = [1, 1, 2, 4, 5, 5, 6, 7, 8, 9, 10, 11, 12, 16]
@test collect(merge_sorted([a; b])) == expected
@test collect(merge_sorted([a; b], x -> x * -1)) == reverse(expected)

a = ["abcd", "acbd", "bcad", "db", "dbac", "dcab"]
b = ["abcde", "acbdf", "baad", "cabd", "dba", "dc"]
c = ["bbcde", "cfbdf", "daad", "dabd", "dba", "ec"]
expected = sort([a; b; c])
@test collect(merge_sorted(Vector{ASCIIString}[a, b, c])) == expected

a = [1, 2, 7]
b = []
@test collect(merge_sorted(Array[a, b])) == [1, 2, 7]

a = []
b = [2, 5, 8]
@test collect(merge_sorted(Array[a, b])) == [2, 5, 8]

a = []
b = []
@test collect(merge_sorted(Array[a, b])) == []

a = [1, 2, 7]
b = []
c = []
d = [2, 5, 8]
e = [100]
f = []
@test collect(merge_sorted(Array[a, b, c, d, e, f])) == [1, 2, 2, 5, 7, 8, 100]


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/InterfaceMechanics.jl/blob/master/LICENSE

@with_kw mutable struct CouloumbDriverState <: AbstractInterfaceState
    time :: Float64 = zero(Float64)
    displacements :: SVector{3,Float64} = zero(SVector{3,Float64})
end

@with_kw mutable struct CouloumbParameterState <: AbstractInterfaceState
    mu :: Float64 = 0.0
    elastic_slip :: Float64 = 0.000
    p0 :: Float64 = 0.0
    c0 :: Float64 = 0.0
end

@with_kw mutable struct CouloumbVariableState <: AbstractInterfaceState
    traction :: SVector{3,Float64} = zero(SVector{3,Float64})
    slip :: SVector{2,Float64} = zero(SVector{2,Float64})
    frictional_energy :: Float64 = 0.0
    jacobian :: SMatrix{3,3,Float64} = zero(SMatrix{3,3,Float64})
end

@with_kw mutable struct Couloumb <: AbstractInterface
    drivers :: CouloumbDriverState = CouloumbDriverState()
    ddrivers :: CouloumbDriverState = CouloumbDriverState()
    variables :: CouloumbVariableState = CouloumbVariableState()
    variables_new :: CouloumbVariableState = CouloumbVariableState()
    parameters :: CouloumbParameterState = CouloumbParameterState()
    dparameters :: CouloumbParameterState = CouloumbParameterState()
end

function integrate_interface!(interface::Couloumb)
    par = interface.parameters
    var = interface.variables
    ddri = interface.ddrivers
    dri = interface.drivers
    @unpack mu, elastic_slip, p0, c0 = par
    @unpack displacements, time = dri
    h0, u10, u20 = displacements
    dtime = ddri.time

    # Drivers at start state
    @unpack traction, slip, frictional_energy, jacobian = var

    # Begin integration
    function dvars(ddisplacement)
        dh, du1, du2 = ddisplacement
        h = h0 + dh
        # Define exponential pressure-overclosure
        p_func(h) = p0/(exp(1.0)-1.0)*((h/c0 + 1.0)*(exp(h/c0 + 1.0) - 1.0))
        dp_func(h) = ForwardDiff.derivative(p_func, h)
        if h <= -c0 # No contact
            p = 0.0
            tau1 = 0.0
            tau2 = 0.0
        elseif h<=6*c0 # Exponential pressure-overclosure
            p = p_func(h)
        else # Linearized continuation after h>6*c0
            p = p_func(6.0*c0) + dp_func(6.0*c0)*(h - 6.0*c0)
        end
        dp = p - traction[1]
        # Trial values for interface tangential movement
        u1 = u10 + du1
        u2 = u20 + du2
        tau_cr = mu*p
        gamma1_tr = u1 - slip[1]
        gamma2_tr = u2 - slip[2]
        gamma_eq_tr = sqrt(gamma1_tr^2 + gamma2_tr^2)
        k = tau_cr/elastic_slip
        if gamma_eq_tr == 0 # Avoid division by zero
            n1 = 1.0/sqrt(2.0)
            n2 = 1.0/sqrt(2.0)
        else
            n1 = gamma1_tr/gamma_eq_tr
            n2 = gamma2_tr/gamma_eq_tr
        end
        jacobian_ = zeros(3,3)
        if gamma_eq_tr>elastic_slip # Slipping
            dslip1 = n1*(gamma_eq_tr - elastic_slip)
            dslip2 = n2*(gamma_eq_tr - elastic_slip)
            tau1 = n1*tau_cr
            tau2 = n2*tau_cr
            dtau1 = tau1 - traction[2]
            dtau2 = tau2 - traction[3]
            dfrictional_energy_ = tau_cr*(gamma_eq_tr - elastic_slip)
            jacobian_[1,1] = dp_func(h)
            jacobian_[2,1] = n1*mu*jacobian_[1,1]
            jacobian_[2,2] = (1.0-n1*n1)*tau_cr/gamma_eq_tr
            jacobian_[2,3] = -n1*n2*tau_cr/gamma_eq_tr
            jacobian_[3,1] = n2*mu*jacobian_[1,1]
            jacobian_[3,2] = -n1*n2*tau_cr/gamma_eq_tr
            jacobian_[3,3] = (1.0-n2*n2)*tau_cr/gamma_eq_tr
        else # Stick/reversal
            dtau1 = k*du1
            dtau2 = k*du2
            dslip1 = 0.0
            dslip2 = 0.0
            dfrictional_energy_ = 0.0
            jacobian_[1,1] = dp_func(h)
            jacobian_[2,1] = (traction[2]+dtau1)/tau_cr*mu
            jacobian_[2,2] = k*jacobian_[1,1]
            jacobian_[3,1] = (traction[3]+dtau2)/tau_cr*mu
            jacobian_[3,3] = k*jacobian_[1,1]
        end
        return [dp, dtau1, dtau2], [dslip1, dslip2], dfrictional_energy_, jacobian_
    end
    dtraction, dslip, dfrictional_energy, jacobian = dvars(ddri.displacements)
    @info("ddisplacements = $(ddri.displacements)")
    @info("traction = $traction, slip = $slip, frictional_energy = $frictional_energy")
    @info("dtraction = $dtraction, dslip = $dslip, dfrictional_energy = $dfrictional_energy")
    #jacobian = ForwardDiff.jacobian(x -> dvars(x)[1], ddri.displacements)
    @info("jacobian = $jacobian")
    frictional_energy += dfrictional_energy
    traction += SVector{3, Float64}(dtraction)
    slip += SVector{2, Float64}(dslip)
    variables_new = CouloumbVariableState(traction=SVector{3, Float64}(traction),
                slip = SVector{2, Float64}(slip),
                frictional_energy = frictional_energy,
                jacobian = SMatrix{3,3,Float64}(jacobian))
    interface.variables_new = variables_new
end


text:
using Test, EndowmentsLH, Random

function fixed_percentiles_test(addNoise)
    rng = MersenneTwister(45);
    @testset "Draw fixed percentiles" begin
        n = 15;
        ed = EndowmentsLH.make_test_endowment_draws(n);
        nameV = names(ed);

        n2 = 5;
        pct = 0.7;
        draws = draw_fixed_percentiles(ed, n2, pct; addNoise = addNoise);
        @test names(draws) == nameV;
        for eName ∈ nameV
            v = get_draws(draws, eName);
            # Check that element type was not changed
            @test eltype(v) == eltype(ed, eName)
            if eltype(v) <: Integer
                @test all(v .== v[1])
            elseif eltype(v) <: Real
                @test all(isapprox.(v, v[1], atol = addNoise + 1e-5))
            end
        end
    end
end

@testset "Fixed percentiles" begin
    for addNoise ∈ (0.0, 0.1)
        fixed_percentiles_test(addNoise);
    end
end

# -------------


text:
using Test
import LanguageServer: Document, Range, Position,
                       get_text, get_offset, get_line_offsets, get_position_at,
                       get_open_in_editor, set_open_in_editor,
                       is_workspace_file,
                       applytextdocumentchanges, TextDocumentContentChangeEvent

@testset "Document changes" begin
    doc = Document("file:///example/path/example.jl", "function foo()", false)
    c1 = TextDocumentContentChangeEvent(Range(Position(0,14), Position(0,14)),
                                        0, "\n")
    c2 = TextDocumentContentChangeEvent(Range(Position(1,0), Position(1,0)),
                                        0, "    ")
    c3 = TextDocumentContentChangeEvent(missing, missing, "println(\"Hello World\")")

    applytextdocumentchanges(doc, c1)
    @test get_text(doc) == "function foo()\n"
    # Implicitly test for issue #403
    applytextdocumentchanges(doc, c2)
    @test get_text(doc) == "function foo()\n    "
    applytextdocumentchanges(doc, c3)
    @test get_text(doc) == "println(\"Hello World\")"
    # doc currently has only one line, applying change to 2nd line should throw
    @test_throws BoundsError applytextdocumentchanges(doc, c2)
end


text:
################################################# FILE DESCRIPTION #########################################################

# This file contains graph algorithms.

################################################# IMPORT/EXPORT ############################################################
export
# External interfaces
hoplist, hoptree, hopgraph, mutualcount, mutual
################################################# BFSLIST ##################################################################

""" Standard BFS implementation that returns a parent vector """
function bfs(g::Graph, seed::Vector{Int}, hopend::Number=Inf)
   N = nv(g)

   parvec = fill(-1, N)
   parvec[seed] = 0

   Q = sizehint!(copy(seed), N)
   Q_ = size
   adj = sizehint!(Vector{Int}(N), N)

   nhops = 1

   lsize = length(seed)
   lcount = 0

   while !isempty(Q) && nhops <= hopend
      u = shift!(Q)

      for v in fadj!(g, u, adj)
         if parvec[v] == -1
            parvec[v] = u
            push!(Q, v)
         end
      end

      lcount += 1

      if lcount == lsize
         lsize = length(Q)
         lcount = 0
         nhops += 1
      end
   end

   return parvec
end

bfs(g::Graph, seed::AbstractVector, nhops::Number=Inf) = bfs(g, collect(seed), nhops)
bfs(g::Graph, seed::Int, nhops::Number=Inf) = bfs(g, Int[seed], nhops)

################################################# BFSLIST ##################################################################

""" Get the list of vertices at a certain distance from the seed """
function bfs_list(g::Graph, seed::Vector{Int}, hopstart::Int=1, hopend::Number=Inf)
   N = nv(g)
   vs = sizehint!(Vector{Int}(0), N)

   visited = falses(N)
   visited[seed] = true

   Q = sizehint!(copy(seed), N)
   adj = sizehint!(Vector{Int}(0), N)

   nhops = 1

   lsize = length(seed)
   lcount = 0

   # Visit one entire level of level_size elements
   while !isempty(Q) && nhops <= hopend
      u = shift!(Q)

      # Visit u's adjacencies
      for v in fadj!(g, u, adj)
         # If v hasn't been visited, visit it
         if !visited[v]
            visited[v] = true

            # Return v if it's in the hoprange
            if nhops >= hopstart
               push!(vs, v)
            end
            push!(Q, v)
         end
      end

      lcount += 1

      if lcount == lsize
         lsize = length(Q)
         lcount = 0
         nhops += 1
      end
   end

   return vs
end

bfs_list(g::Graph, seed::AbstractVector, hopstart::Int=1, hopend::Number=Inf) = bfs_list(g, collect(seed), hopstart, hopend)
bfs_list(g::Graph, seed::Int, hopstart::Int=1, hopend::Number=Inf) = bfs_list(g, [seed], hopstart, hopend)

"""
Get a list of vertices at a certain hop distance from a labelled vertex
"""
function hoplist(g::Graph, x, hopstart::Int, hopend::Number=Inf)
   encode(g, bfs_list(g, decode(g, x), hopstart, hopend))
end

################################################# BFSTREE ###############################################################

"""
Returns a BFS tree, containing explored vertices and only tree edges.
"""
function bfs_tree(g::Graph, seed::Int, hopend::Number=Inf)
   parvec = bfs(g, seed, hopend)

   us = sizehint!(Int[], nv(g))
   vs = sizehint!(Int[], nv(g))

   for v in eachindex(parvec)
      u = parvec[v]
      if u > 0
         push!(us, u)
         push!(vs, v)
      end
   end
   vlist = vcat(seed, vs)
   eit = EdgeIter(length(us), us, vs)

   subgraph(g, vlist, eit)
end

"""
Get the bfs tree starting from the input labelled vertex, and ending at a certain hop distance
"""
function hoptree(g::Graph, x, hopend::Number=Inf)
   bfs_tree(g, decode(g, x), hopend)
end

################################################# BFSSUBGRAPH #############################################################

"""
Returns a BFS subgraph, containing explored vertices and all edges between
them
"""
function bfs_subgraph(g::Graph, seed, hopend::Number=Inf)
   vs = append!(bfs_list(g, seed, 1, hopend), seed)
   subgraph(g, vs)
end

"""
Get a subgraph containing vertices and edges within a certain hop distance from the input labelled
vertex
"""
function hopgraph(g::Graph, x, hopend::Number=Inf)
   bfs_subgraph(g, decode(g, x), hopend)
end

################################################# MUTUAL FRIENDS ###########################################################

# Hacky to ensure performance
function count_mutual_adj(g::Graph, u::VertexID, v::VertexID)
   x = indxs(g)

   i = x.colptr[u]
   j = x.colptr[v]

   ei = x.colptr[u+1] - 1
   ej = x.colptr[v+1] - 1

   count = 0

   while i <= ei && j <= ej
      @inbounds vi = x.rowval[i]
      @inbounds vj = x.rowval[j]

      if vi < vj
         i += 1
      elseif vi == vj
         count += 1
         i += 1
         j += 1
      else
         j += 1
      end
   end

   return count
end

mutualcount(g::Graph, ul, vl) = count_mutual_adj(g, decode(g, ul), decode(g, vl))

function mutual_adj(g::Graph, u::VertexID, v::VertexID)
   x = indxs(g)

   aru = nzrange(x, u)
   arv = nzrange(x, v)

   i = start(aru)
   j = start(arv)

   res = sizehint!(Int[], count_mutual_adj(g, u, v))

   while i in aru && j in arv
      if x.rowval[i] < x.rowval[j]
         i += 1
      elseif x.rowval[i] == x.rowval[j]
         push!(res, x.rowval[i])
         i += 1
         j += 1
      else
         j += 1
      end
   end

   return res
end

mutual(g::Graph, ul, vl) = encode(g, mutual_adj(g, decode(g, ul), decode(g, vl)))


text:
using MUMPS, MATLAB

#### TO DO:
#
# - clean up normal equation stuff

######################################################################
# Each instance of abstract_linear_system_solver should contain
#	initialize::Function
# 	ls_factor::Function
# 	ls_solve!::Function
# i.e. see class_mumps_linear_system_solver
#######################################################################

# standard function for solving either the whole system or normal equations

# solve whole system
function construct_K(qp::class_quadratic_program, my_settings::class_settings)
	return sparse([[qp._H qp.A']; [qp.A -my_settings.diagonal_modification*speye(qp.m,qp.m)]])
end

function update_K!(K::SparseMatrixCSC{Float64,Int64}, vars::class_variables, qp::class_quadratic_program)
	x = vars.x();
	s = vars.s();

	H = qp.H(vars);
	for i = 1:(qp.n)
		K[i,i] = H[i,i] + s[i]/x[i];
	end
end

# solve normal equations A*A'

function form_normal_equations(vars::class_variables, qp::class_quadratic_program)
	H = qp.H(vars);
	diH = diag(H);
	@assert( sum(diH .>  0) == length(nonzeros(H)) ) # check that nothing is on the diagonals of H

	invH = spdiagm( 1./(diH + vars.s() ./ vars.x() ) );

	GLOBAL_timer.start("Factor/Schur")
	normal_matrix = qp.A * invH * qp.A' + settings.diagonal_modification*speye(qp.m,qp.m);
	GLOBAL_timer.stop("Factor/Schur")

	return normal_matrix, invH;
end

function standard_ls_initialize(qp, my_settings, is_normal)
	if is_normal
		return speye(0);
	else
		return construct_K(qp, my_settings);
	end
end

function standard_ls_before_factor!(Sparse_Matrix, vars, qp, is_normal)
	#this.SparseMatrix, this.invH = before_factor(this.SparseMatrix, vars, qp, this.options.normal);
	if is_normal
		return form_normal_equations(vars, qp);
	else
		update_K!(Sparse_Matrix, vars,  qp);
		return Sparse_Matrix, speye(0)
	end
end

function standard_ls_before_solve(my_sol::Array{Float64}, my_rhs::Array{Float64}, invH::SparseMatrixCSC{Float64,Int64}, qp::class_quadratic_program, is_normal::Bool)
	if is_normal
		n = qp.n
		m = qp.m
		temp_rhs = A * invH * my_rhs[1:n] - my_rhs[(1 + n):(m + n)];
		temp_sol = my_sol[(n + 1):(n + m)];

		return temp_sol, temp_rhs
	else
		return my_sol, my_rhs
	end
end

function standard_ls_after_solve!(my_sol::Array{Float64}, temp_sol::Array{Float64}, my_rhs::Array{Float64}, invH::SparseMatrixCSC{Float64,Int64}, qp::class_quadratic_program, is_normal::Bool)
	if is_normal
		my_sol[1:qp.n] = invH*(my_rhs[1:qp.n] - qp.A' * temp_sol);
		my_sol[(qp.n + 1):(qp.n + qp.m)] = temp_sol;
	end
	# don't need to do anything if equations are not normal
	return my_sol
end

#########################
# compute the inertia
#########################

function compute_inertia(diag_vec::Array{Float64}, num_vars::Int64, num_constraints::Int64)
	try
		@assert(length(diag_vec) == num_vars + num_constraints)
		pos_eigs = sum(diag_vec .> 0.0);
		neg_eigs = sum(diag_vec .< 0.0);

		if pos_eigs == num_vars && neg_eigs == num_constraints
			return true
		elseif pos_eigs > num_vars || neg_eigs + pos_eigs != num_vars + num_constraints
			println("Warning: numerical instability in LDL factorization")
			if pos_eigs > num_vars
				println("more positive eigenvalues than variables")
				return false
			else
				println("zero eigenvalues")
				return false
			end
		else
			return false
		end
	catch e
		println("ERROR in compute_inertia")
		throw(e)
	end
end

# norm_rhs = normal_rhs(rhs,qp)
# dy = factor \ norm_rhs
# dx = solve_normal(dy, b)
# dx = H\ A'*

#########################################
# Linear solver options
#########################################

abstract abstract_linear_solver_options

type standard_linear_solver_options <: abstract_linear_solver_options
	sym::Int64
	normal::Bool

	function standard_linear_solver_options()
		this = new();

		this.sym = 0; # use LU factorization by default
		this.normal = false;

		return this;
	end
end

##########################################################################
# THE DIFFERENT LINEAR SOLVERS
##########################################################################

####
#### add any new linear solver here ...
####

abstract abstract_linear_system_solver

type linear_solver_MUMPS <: abstract_linear_system_solver
	initialize::Function
	ls_factor::Function
	ls_solve!::Function

	SparseMatrix::SparseMatrixCSC{Float64,Int64}
	factor::MUMPSfactorization{Float64}

	# options
	options::standard_linear_solver_options

	qp::class_quadratic_program
	invH::SparseMatrixCSC{Float64,Int64}

	function linear_solver_MUMPS()
		this = new();
		this.options = standard_linear_solver_options();

		this.initialize = function(qp::class_quadratic_program, my_settings::class_settings)
			this.qp = qp;
			this.invH = speye(0);
			this.SparseMatrix = standard_ls_initialize(qp, my_settings, this.options.normal);
		end

		this.ls_factor = function(vars::class_variables, qp::class_quadratic_program)
			#if isdefined(this,:factor)
			#	destroyMUMPS(this.factor);
			#end
			this.SparseMatrix, this.invH = standard_ls_before_factor!(this.SparseMatrix, vars, qp, this.options.normal);

			GLOBAL_timer::class_algorithm_timer
			GLOBAL_timer.start("Factor/MUMPs")
			this.factor = factorMUMPS(this.SparseMatrix, this.options.sym);
			GLOBAL_timer.stop("Factor/MUMPs")

			return 1; # inertia
		end

		this.ls_solve! = function(my_rhs::Array{Float64}, my_sol::Array{Float64}, qp::class_quadratic_program)
			is_normal = this.options.normal;

			temp_sol, temp_rhs = standard_ls_before_solve(my_sol, my_rhs, this.invH, this.qp, is_normal)

			GLOBAL_timer::class_algorithm_timer
			GLOBAL_timer.start("Solve/MUMPs")
			applyMUMPS!(this.factor, temp_rhs, temp_sol);
			GLOBAL_timer.stop("Solve/MUMPs")

			standard_ls_after_solve!(my_sol, temp_sol, my_rhs, this.invH, this.qp, is_normal)
		end

		return this;
	end
end

type linear_solver_julia <: abstract_linear_system_solver
	initialize::Function
	ls_factor::Function
	ls_solve!::Function

	SparseMatrix::SparseMatrixCSC{Float64,Int64}
	factor

	# options
	options::standard_linear_solver_options

	qp::class_quadratic_program
	invH::SparseMatrixCSC{Float64,Int64}

	function linear_solver_julia()
		this = new();
		this.options = standard_linear_solver_options()
		this.options.sym = 0; # use LU factorization by default

		this.initialize = function(qp::class_quadratic_program, my_settings::class_settings)
			this.invH = speye(0);
			this.SparseMatrix = standard_ls_initialize(qp, my_settings, this.options.normal);
		end

		this.ls_factor = function(vars::class_variables,qp::class_quadratic_program)
			this.SparseMatrix, this.invH = standard_ls_before_factor!(this.SparseMatrix, vars, qp, this.options.normal);

			if this.options.sym == 0
				this.factor = lufact(this.SparseMatrix)
			elseif this.options.sym == 1
				this.factor = cholfact(this.SparseMatrix)
			else
				error("this.options.sym = " * string(this.options.sym) * " not supported")
			end

			return 1; # inertia
		end

		this.ls_solve! = function(my_rhs::Array{Float64}, my_sol::Array{Float64}, qp::class_quadratic_program)
			is_normal = this.options.normal;

			temp_sol, temp_rhs = standard_ls_before_solve(my_sol, my_rhs, this.invH, qp, is_normal)

			GLOBAL_timer::class_algorithm_timer
			GLOBAL_timer.start("Solve/julia")
			temp_sol[1:length(temp_sol)] = this.factor \ temp_rhs; #::UmfpackLU{Float64,Int64}
			GLOBAL_timer.stop("Solve/julia")

			standard_ls_after_solve!(my_sol, temp_sol, my_rhs, this.invH, qp, is_normal)
		end

		return this;
	end
end

type linear_solver_matlab <: abstract_linear_system_solver
	initialize::Function
	ls_factor::Function
	ls_solve!::Function

	SparseMatrix::SparseMatrixCSC{Float64,Int64}

	# options
	options::abstract_linear_solver_options

	function linear_solver_matlab()
		this = new();
		this.options = standard_linear_solver_options()
		this.options.sym = 0; # use LU factorization by default

		this.initialize = function(qp::class_quadratic_program, my_settings::class_settings)
			this.SparseMatrix = construct_K(qp, my_settings);
		end

		this.ls_factor = function(vars::class_variables, qp::class_quadratic_program)
			update_K!(this.SparseMatrix, vars, qp);

			sparse_matrix = this.SparseMatrix
			@mput sparse_matrix

			@matlab begin
				mat_L, mat_D, mat_P, mat_S = ldl(sparse_matrix,1e-6);
			end

			@mget mat_D

			return compute_inertia(diag(mat_D), qp.n, qp.m)
		end

		this.ls_solve! = function(my_rhs::Array{Float64}, my_sol::Array{Float64}, qp::class_quadratic_program)
			@mput my_rhs

			@matlab begin
				mat_solution = mat_S * ( mat_P* (mat_L'\ (mat_D\ (mat_L\ (mat_P' * ( mat_S * my_rhs ) ) ) ) ) );
			end

			@mget mat_solution

			my_sol[1:length(my_sol)] =  mat_solution;
		end

		return this;
	end
end

max_forks_repo_name: ohinder/Everest


text:
using Distributed

@everywhere using Ipopt

include("code-2-lib/parsers.jl")
include("code-2-lib/lib.jl")

function compute_solution2(con_file::String, inl_file::String, raw_file::String, rop_file::String, time_limit::Int, scoring_method::Int, network_model::String; output_dir::String="", scenario_id::String="none")
    time_data_start = time()
    PowerModels.silence()
    goc_data = parse_goc_files(con_file, inl_file, raw_file, rop_file, scenario_id=scenario_id)
    network = build_pm_model(goc_data)

    sol = read_solution1(network, output_dir=output_dir)
    PowerModels.update_data!(network, sol)

    check_network_solution(network)

    network_tmp = deepcopy(network)
    balance = compute_power_balance_deltas!(network_tmp)

    if balance.p_delta_abs_max > 0.01 || balance.q_delta_abs_max > 0.01
        error(LOGGER, "solution1 power balance requirements not satified (all power balance values should be below 0.01). $(balance)")
    end
    load_time = time() - time_data_start

    ###### Prepare Solution 2 ######

    time_contingencies_start = time()

    gen_cont_total = length(network["gen_contingencies"])
    branch_cont_total = length(network["branch_contingencies"])
    cont_total = gen_cont_total + branch_cont_total

    cont_order = contingency_order(network)

    workers = Distributed.workers()

    process_data = []

    cont_per_proc = cont_total/length(workers)

    for p in 1:length(workers)
        cont_start = trunc(Int, ceil(1+(p-1)*cont_per_proc))
        cont_end = min(cont_total, trunc(Int,ceil(p*cont_per_proc)))
        pd = (
            pid = p,
            processes = length(workers),
            con_file = con_file,
            inl_file = inl_file,
            raw_file = raw_file,
            rop_file = rop_file,
            scenario_id = scenario_id,
            output_dir = output_dir,
            cont_range = cont_start:cont_end,
        )
        push!(process_data, pd)
    end

    for (i,pd) in enumerate(process_data)
        info(LOGGER, "worker task $(pd.pid): $(length(pd.cont_range)) / $(pd.cont_range)")
    end

    solution2_files = pmap(solution2_solver, process_data, retry_delays = zeros(3))
    sort!(solution2_files)
    #println("pmap result: $(solution2_files)")

    time_contingencies = time() - time_contingencies_start
    info(LOGGER, "contingency eval time: $(time_contingencies)")
    info(LOGGER, "time per contingency: $(time_contingencies/cont_total)")

    info(LOGGER, "combine $(length(solution2_files)) solution2 files")
    combine_files(solution2_files, "solution2.txt"; output_dir=output_dir)
    remove_files(solution2_files)

    println("")

    data = [
        "----",
        "scenario id",
        "bus",
        "branch",
        "gen_cont",
        "branch_cont",
        "runtime (sec.)",
    ]
    println(join(data, ", "))

    data = [
        "DATA_SSS",
        goc_data.scenario,
        length(network["bus"]),
        length(network["branch"]),
        length(network["gen_contingencies"]),
        length(network["branch_contingencies"]),
        time_contingencies,
    ]
    println(join(data, ", "))
end

@everywhere function solution2_solver(process_data)
    #println(process_data)
    time_data_start = time()
    PowerModels.silence()
    goc_data = parse_goc_files(
        process_data.con_file, process_data.inl_file, process_data.raw_file,
        process_data.rop_file, scenario_id=process_data.scenario_id)
    network = build_pm_model(goc_data)

    sol = read_solution1(network, output_dir=process_data.output_dir)
    PowerModels.update_data!(network, sol)
    correct_voltage_angles!(network)

    time_data = time() - time_data_start

    for (i,bus) in network["bus"]
        if haskey(bus, "evhi")
            bus["vmax"] = bus["evhi"]
        end
        if haskey(bus, "evlo")
            bus["vmin"] = bus["evlo"]
        end
    end

    for (i,branch) in network["branch"]
        if haskey(branch, "rate_c")
            branch["rate_a"] = branch["rate_c"]
        end
    end

    contingencies = contingency_order(network)[process_data.cont_range]

    for (i,branch) in network["branch"]
        g, b = PowerModels.calc_branch_y(branch)
        tr, ti = PowerModels.calc_branch_t(branch)
        branch["g"] = g
        branch["b"] = b
        branch["tr"] = tr
        branch["ti"] = ti
    end

    bus_gens = gens_by_bus(network)

    network["delta"] = 0
    for (i,bus) in network["bus"]
        bus["vm_base"] = bus["vm"]
        bus["vm_start"] = bus["vm"]
        bus["va_start"] = bus["va"]
        bus["vm_fixed"] = length(bus_gens[i]) != 0

    end

    for (i,gen) in network["gen"]
        gen["pg_base"] = gen["pg"]
        gen["pg_start"] = gen["pg"]
        gen["qg_start"] = gen["qg"]
        gen["pg_fixed"] = false
        gen["qg_fixed"] = false
    end

    #nlp_solver = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6, mu_init=1e-6, hessian_approximation="limited-memory", print_level=0)
    nlp_solver = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6, print_level=0)
    #nlp_solver = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6)
    #nlp_solver = JuMP.with_optimizer(Ipopt.Optimizer, tol=1e-6, hessian_approximation="limited-memory")

    pad_size = trunc(Int, ceil(log(10,process_data.processes)))
    padded_pid = lpad(string(process_data.pid), pad_size, "0")
    solution_filename = "solution2-$(padded_pid).txt"

    if length(process_data.output_dir) > 0
        solution_path = joinpath(process_data.output_dir, solution_filename)
    else
        solution_path = solution_filename
    end
    if isfile(solution_path)
        warn(LOGGER, "removing existing solution2 file $(solution_path)")
        rm(solution_path)
    end
    open(solution_path, "w") do sol_file
        # creates an empty file in the case of workers without contingencies
    end

    #network_tmp = deepcopy(network)
    for cont in contingencies
        if cont.type == "gen"
            info(LOGGER, "working on: $(cont.label)")
            time_start = time()
            network_tmp = deepcopy(network)
            debug(LOGGER, "contingency copy time: $(time() - time_start)")
            network_tmp["cont_label"] = cont.label

            cont_gen = network_tmp["gen"]["$(cont.idx)"]
            cont_gen["contingency"] = true
            cont_gen["gen_status"] = 0
            pg_lost = cont_gen["pg"]

            time_start = time()
            result = run_fixpoint_pf_v2_2!(network_tmp, pg_lost, ACRPowerModel, nlp_solver, iteration_limit=5)
            debug(LOGGER, "second-stage contingency solve time: $(time() - time_start)")

            result["solution"]["label"] = cont.label
            result["solution"]["feasible"] = (result["termination_status"] == LOCALLY_SOLVED)
            result["solution"]["cont_type"] = "gen"
            result["solution"]["cont_comp_id"] = cont.idx

            result["solution"]["gen"]["$(cont.idx)"]["pg"] = 0.0
            result["solution"]["gen"]["$(cont.idx)"]["qg"] = 0.0

            correct_contingency_solution!(network, result["solution"])
            open(solution_path, "a") do sol_file
                sol2 = write_solution2_contingency(sol_file, network, result["solution"])
            end

            network_tmp["gen"]["$(cont.idx)"]["gen_status"] = 1
        elseif cont.type == "branch"
            info(LOGGER, "working on: $(cont.label)")
            time_start = time()
            network_tmp = deepcopy(network)
            debug(LOGGER, "contingency copy time: $(time() - time_start)")
            network_tmp["cont_label"] = cont.label
            network_tmp["branch"]["$(cont.idx)"]["br_status"] = 0

            time_start = time()
            result = run_fixpoint_pf_v2_2!(network_tmp, 0.0, ACRPowerModel, nlp_solver, iteration_limit=5)
            debug(LOGGER, "second-stage contingency solve time: $(time() - time_start)")

            result["solution"]["label"] = cont.label
            result["solution"]["feasible"] = (result["termination_status"] == LOCALLY_SOLVED)
            result["solution"]["cont_type"] = "branch"
            result["solution"]["cont_comp_id"] = cont.idx

            correct_contingency_solution!(network, result["solution"])
            open(solution_path, "a") do sol_file
                sol2 = write_solution2_contingency(sol_file, network, result["solution"])
            end

            network_tmp["branch"]["$(cont.idx)"]["br_status"] = 1
        else
            @assert("contingency type $(cont.type) not known")
        end
    end

    return solution_path
end


text:
export RPLoadBalancer, parallel_lb_rp

"""
Load balancing through random polling.
"""
immutable RPLoadBalancer <: AbstractLoadBalancer
    msg_chls::Array{RemoteChannel{Channel{Message}}}
    stat_chl::RemoteChannel{Channel{Message}}
    res_chl::RemoteChannel{Channel{Message}}
    statuses::Array{Symbol}
end
RPLoadBalancer(cap::Integer=64) = RPLoadBalancer(create_msg_chls(cap),
                                                 RemoteChannel(()->Channel{Message}(cap), 1),
                                                 RemoteChannel(()->Channel{Message}(cap), 1),
                                                 fill!(Array{Symbol}(nworkers()), :unstarted))

"""
    parallel_lb(balancer::RPLoadBalancer, work::WorkUnit)

Compute the task associated with `work` using random polling load balancing.
"""
function parallel_lb(balancer::RPLoadBalancer, work::WorkUnit)
    @sync begin
        # start the worker processes
        for wid in workers()
            @spawnat wid worker(balancer)
        end

        @async recv_results(balancer)

        @sync begin
            balancer.statuses[1] = :started

            # send initial work
            @async put!(balancer.msg_chls[1], Message(:work, myid(), work))

            @async status_manager(balancer)

        end

        for i = 1:nworkers()
            put!(balancer.msg_chls[i], Message(:end, -1))
        end
    end
end

"""
    parallel_lb_rp(work::WorkUnit)

Compute the task associated with `work` using the default RPLoadBalancer
"""
parallel_lb_rp(work::WorkUnit) = parallel_lb(RPLoadBalancer(), work)

"""
    status_manager(balancer::RPLoadBalancer)

Track worker status through update messages.
"""
function status_manager(balancer::RPLoadBalancer)
    # while there are any started, nonidle nodes
    if nworkers() > 1
        info("status_manager is waking up idle workers...")
        for w_idx in 2:nworkers()
            info("\t... ", workers()[w_idx])
            put!(balancer.msg_chls[w_idx], Message(:_idle, myid()))
        end
    end
    while any((balancer.statuses .!= :unstarted) .& (balancer.statuses .!= :idle))
        status_msg = take!(balancer.stat_chl)

        #w_idx = nprocs() > 1 ? status_msg.data - 1 : status_msg.data
        wid = status_msg.data

        if status_msg.kind == :idle
            # mark this worker as idle
            balancer.statuses[w_idx(wid)] = :idle

        elseif status_msg.kind == :nonidle
            balancer.statuses[w_idx(wid)] = :nonidle

        else
            error("Invalid message received by status_manager")

        end
    end
end

"""
Receive messages on its remote channel. Depending on the message,
different actions will be taken:

External Messages: (Could come from anywhere)
- :work - Pass from the remote channel to local_chl
- :jlance - Start a new task that will attempt to send work in local_chl
              to the remote worker specified by this message
- :nowork - Another worker failed to send work to this worker; request
              work from a random worker.
- :end  - Pass to local_chl and exit.

Internal Messages: (Expect to receive these only from other tasks on this process)
- :_idle - Send an :idle message to the controller via stat_chl;
             Request work from a random worker.
- :_nonidle - Send a :nonidle message to the controller via stat_chl
"""
function _msg_handler(balancer::RPLoadBalancer,
                      local_chl::Channel{Message})

    msg_chl = get_msg_chl(balancer)

    while true
        let msg = take!(msg_chl)
        # start soft local scope
        if msg.kind == :end
            put!(local_chl, msg)
            break

        elseif msg.kind == :work
            put!(local_chl, msg)

        elseif msg.kind == :nowork
            put!(balancer.stat_chl, Message(:idle, myid()))
            other_wid = rand(workers())
            other_msg_chl = get_msg_chl(balancer, other_wid)
            put!(other_msg_chl, Message(:jlance, myid()))

        elseif msg.kind == :jlance && msg.data > 0
            # attempt to load balance
            @schedule send_work(balancer, local_chl, msg)

        elseif msg.kind == :_idle
            put!(balancer.stat_chl, Message(:idle, myid()))
            # if there is only 1 worker, there is no one else to get work from
            #    and, if the only worker is idle, that means we should finish
            info("worker ", myid(), " idle")
            if nworkers() > 1
                other_wid = rand(workers())
                other_msg_chl = get_msg_chl(balancer, other_wid)
                put!(other_msg_chl, Message(:jlance, myid()))
            end

        elseif msg.kind == :_nonidle
            # put! on remote chl blocks, so schedule in different task
            @schedule put!(balancer.stat_chl, Message(:nonidle, myid()))
        end
        # end soft local scope
        end
    end
end

max_stars_repo_name: dillondaudert/jlance

max_issues_repo_name: dillondaudert/jlance

max_forks_repo_name: dillondaudert/jlance


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## number-theoretic functions ##

"""
    gcd(x,y)

Greatest common (positive) divisor (or zero if `x` and `y` are both zero).

# Examples
```jldoctest
julia> gcd(6,9)
3

julia> gcd(6,-9)
3
```
"""
function gcd(a::T, b::T) where T<:Integer
    while b != 0
        t = b
        b = rem(a, b)
        a = t
    end
    checked_abs(a)
end

# binary GCD (aka Stein's) algorithm
# about 1.7x (2.1x) faster for random Int64s (Int128s)
function gcd(a::T, b::T) where T<:Union{Int8,UInt8,Int16,UInt16,Int32,UInt32,Int64,UInt64,Int128,UInt128}
    a == 0 && return abs(b)
    b == 0 && return abs(a)
    za = trailing_zeros(a)
    zb = trailing_zeros(b)
    k = min(za, zb)
    u = unsigned(abs(a >> za))
    v = unsigned(abs(b >> zb))
    while u != v
        if u > v
            u, v = v, u
        end
        v -= u
        v >>= trailing_zeros(v)
    end
    r = u << k
    # T(r) would throw InexactError; we want OverflowError instead
    r > typemax(T) && __throw_gcd_overflow(a, b)
    r % T
end
@noinline __throw_gcd_overflow(a, b) = throw(OverflowError("gcd($a, $b) overflows"))

"""
    lcm(x,y)

Least common (non-negative) multiple.

# Examples
```jldoctest
julia> lcm(2,3)
6

julia> lcm(-2,3)
6
```
"""
function lcm(a::T, b::T) where T<:Integer
    # explicit a==0 test is to handle case of lcm(0,0) correctly
    if a == 0
        return a
    else
        return checked_abs(a * div(b, gcd(b,a)))
    end
end

gcd(a::Integer) = a
lcm(a::Integer) = a
gcd(a::Integer, b::Integer) = gcd(promote(a,b)...)
lcm(a::Integer, b::Integer) = lcm(promote(a,b)...)
gcd(a::Integer, b::Integer...) = gcd(a, gcd(b...))
lcm(a::Integer, b::Integer...) = lcm(a, lcm(b...))

lcm(abc::AbstractArray{<:Integer}) = reduce(lcm, abc; init=one(eltype(abc)))

function gcd(abc::AbstractArray{<:Integer})
    a = zero(eltype(abc))
    for b in abc
        a = gcd(a,b)
        if a == 1
            return a
        end
    end
    return a
end

# return (gcd(a,b),x,y) such that ax+by == gcd(a,b)
"""
    gcdx(x,y)

Computes the greatest common (positive) divisor of `x` and `y` and their Bézout
coefficients, i.e. the integer coefficients `u` and `v` that satisfy
``ux+vy = d = gcd(x,y)``. ``gcdx(x,y)`` returns ``(d,u,v)``.

# Examples
```jldoctest
julia> gcdx(12, 42)
(6, -3, 1)

julia> gcdx(240, 46)
(2, -9, 47)
```

!!! note
    Bézout coefficients are *not* uniquely defined. `gcdx` returns the minimal
    Bézout coefficients that are computed by the extended Euclidean algorithm.
    (Ref: D. Knuth, TAoCP, 2/e, p. 325, Algorithm X.)
    For signed integers, these coefficients `u` and `v` are minimal in
    the sense that ``|u| < |y/d|`` and ``|v| < |x/d|``. Furthermore,
    the signs of `u` and `v` are chosen so that `d` is positive.
    For unsigned integers, the coefficients `u` and `v` might be near
    their `typemax`, and the identity then holds only via the unsigned
    integers' modulo arithmetic.
"""
function gcdx(a::T, b::T) where T<:Integer
    # a0, b0 = a, b
    s0, s1 = oneunit(T), zero(T)
    t0, t1 = s1, s0
    # The loop invariant is: s0*a0 + t0*b0 == a
    while b != 0
        q = div(a, b)
        a, b = b, rem(a, b)
        s0, s1 = s1, s0 - q*s1
        t0, t1 = t1, t0 - q*t1
    end
    a < 0 ? (-a, -s0, -t0) : (a, s0, t0)
end
gcdx(a::Integer, b::Integer) = gcdx(promote(a,b)...)

# multiplicative inverse of n mod m, error if none

"""
    invmod(x,m)

Take the inverse of `x` modulo `m`: `y` such that ``x y = 1 \\pmod m``,
with ``div(x,y) = 0``. This is undefined for ``m = 0``, or if
``gcd(x,m) \\neq 1``.

# Examples
```jldoctest
julia> invmod(2,5)
3

julia> invmod(2,3)
2

julia> invmod(5,6)
5
```
"""
function invmod(n::T, m::T) where T<:Integer
    g, x, y = gcdx(n, m)
    g != 1 && throw(DomainError((n, m), "Greatest common divisor is $g."))
    m == 0 && throw(DomainError(m, "`m` must not be 0."))
    # Note that m might be negative here.
    # For unsigned T, x might be close to typemax; add m to force a wrap-around.
    r = mod(x + m, m)
    # The postcondition is: mod(r * n, m) == mod(T(1), m) && div(r, m) == 0
    r
end
invmod(n::Integer, m::Integer) = invmod(promote(n,m)...)

# ^ for any x supporting *
to_power_type(x) = convert(Base._return_type(*, Tuple{typeof(x), typeof(x)}), x)
@noinline throw_domerr_powbysq(::Any, p) = throw(DomainError(p,
    string("Cannot raise an integer x to a negative power ", p, '.',
           "\nConvert input to float.")))
@noinline throw_domerr_powbysq(::Integer, p) = throw(DomainError(p,
   string("Cannot raise an integer x to a negative power ", p, '.',
          "\nMake x a float by adding a zero decimal (e.g., 2.0^$p instead ",
          "of 2^$p), or write 1/x^$(-p), float(x)^$p, or (x//1)^$p")))
@noinline throw_domerr_powbysq(::AbstractMatrix, p) = throw(DomainError(p,
   string("Cannot raise an integer matrix x to a negative power ", p, '.',
          "\nMake x a float matrix by adding a zero decimal ",
          "(e.g., [2.0 1.0;1.0 0.0]^$p instead ",
          "of [2 1;1 0]^$p), or write float(x)^$p or Rational.(x)^$p")))
function power_by_squaring(x_, p::Integer)
    x = to_power_type(x_)
    if p == 1
        return copy(x)
    elseif p == 0
        return one(x)
    elseif p == 2
        return x*x
    elseif p < 0
        isone(x) && return copy(x)
        isone(-x) && return iseven(p) ? one(x) : copy(x)
        throw_domerr_powbysq(x, p)
    end
    t = trailing_zeros(p) + 1
    p >>= t
    while (t -= 1) > 0
        x *= x
    end
    y = x
    while p > 0
        t = trailing_zeros(p) + 1
        p >>= t
        while (t -= 1) >= 0
            x *= x
        end
        y *= x
    end
    return y
end
power_by_squaring(x::Bool, p::Unsigned) = ((p==0) | x)
function power_by_squaring(x::Bool, p::Integer)
    p < 0 && !x && throw_domerr_powbysq(x, p)
    return (p==0) | x
end

^(x::T, p::T) where {T<:Integer} = power_by_squaring(x,p)
^(x::Number, p::Integer)  = power_by_squaring(x,p)

# x^p for any literal integer p is lowered to Base.literal_pow(^, x, Val(p))
# to enable compile-time optimizations specialized to p.
# However, we still need a fallback that calls the function ^ which may either
# mean Base.^ or something else, depending on context.
# We mark these @inline since if the target is marked @inline,
# we want to make sure that gets propagated,
# even if it is over the inlining threshold.
@inline literal_pow(f, x, ::Val{p}) where {p} = f(x,p)

# Restrict inlining to hardware-supported arithmetic types, which
# are fast enough to benefit from inlining.
const HWReal = Union{Int8,Int16,Int32,Int64,UInt8,UInt16,UInt32,UInt64,Float32,Float64}
const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}}

# Core.Compiler has complicated logic to inline x^2 and x^3 for
# numeric types.  In terms of Val we can do it much more simply.
# (The first argument prevents unexpected behavior if a function ^
# is defined that is not equal to Base.^)
@inline literal_pow(::typeof(^), x::HWNumber, ::Val{0}) = one(x)
@inline literal_pow(::typeof(^), x::HWNumber, ::Val{1}) = x
@inline literal_pow(::typeof(^), x::HWNumber, ::Val{2}) = x*x
@inline literal_pow(::typeof(^), x::HWNumber, ::Val{3}) = x*x*x

# don't use the inv(x) transformation here since float^p is slightly more accurate
@inline literal_pow(::typeof(^), x::AbstractFloat, ::Val{p}) where {p} = x^p
@inline literal_pow(::typeof(^), x::AbstractFloat, ::Val{-1}) = inv(x)

# for other types, define x^-n as inv(x)^n so that negative literal powers can
# be computed in a type-stable way even for e.g. integers.
@inline @generated function literal_pow(f::typeof(^), x, ::Val{p}) where {p}
    if p < 0
        :(literal_pow(^, inv(x), $(Val{-p}())))
    else
        :(f(x,$p))
    end
end

# note: it is tempting to add optimized literal_pow(::typeof(^), x, ::Val{n})
#       methods here for various n, but this easily leads to method ambiguities
#       if anyone has defined literal_pow(::typeof(^), x::T, ::Val).

# b^p mod m

"""
    powermod(x::Integer, p::Integer, m)

Compute ``x^p \\pmod m``.

# Examples
```jldoctest
julia> powermod(2, 6, 5)
4

julia> mod(2^6, 5)
4

julia> powermod(5, 2, 20)
5

julia> powermod(5, 2, 19)
6

julia> powermod(5, 3, 19)
11
```
"""
function powermod(x::Integer, p::Integer, m::T) where T<:Integer
    p < 0 && return powermod(invmod(x, m), -p, m)
    p == 0 && return mod(one(m),m)
    (m == 1 || m == -1) && return zero(m)
    b = oftype(m,mod(x,m))  # this also checks for divide by zero

    t = prevpow(2, p)
    r::T = 1
    while true
        if p >= t
            r = mod(widemul(r,b),m)
            p -= t
        end
        t >>>= 1
        t <= 0 && break
        r = mod(widemul(r,r),m)
    end
    return r
end

# optimization: promote the modulus m to BigInt only once (cf. widemul in generic powermod above)
powermod(x::Integer, p::Integer, m::Union{Int128,UInt128}) = oftype(m, powermod(x, p, big(m)))

_nextpow2(x::Unsigned) = oneunit(x)<<((sizeof(x)<<3)-leading_zeros(x-oneunit(x)))
_nextpow2(x::Integer) = reinterpret(typeof(x),x < 0 ? -_nextpow2(unsigned(-x)) : _nextpow2(unsigned(x)))
_prevpow2(x::Unsigned) = one(x) << unsigned((sizeof(x)<<3)-leading_zeros(x)-1)
_prevpow2(x::Integer) = reinterpret(typeof(x),x < 0 ? -_prevpow2(unsigned(-x)) : _prevpow2(unsigned(x)))

"""
    ispow2(n::Integer) -> Bool

Test whether `n` is a power of two.

# Examples
```jldoctest
julia> ispow2(4)
true

julia> ispow2(5)
false
```
"""
ispow2(x::Integer) = x > 0 && count_ones(x) == 1

"""
    nextpow(a, x)

The smallest `a^n` not less than `x`, where `n` is a non-negative integer. `a` must be
greater than 1, and `x` must be greater than 0.

# Examples
```jldoctest
julia> nextpow(2, 7)
8

julia> nextpow(2, 9)
16

julia> nextpow(5, 20)
25

julia> nextpow(4, 16)
16
```

See also [`prevpow`](@ref).
"""
function nextpow(a::Real, x::Real)
    x <= 0 && throw(DomainError(x, "`x` must be positive."))
    # Special case fast path for x::Integer, a == 2.
    # This is a very common case. Constant prop will make sure that a call site
    # specified as `nextpow(2, x)` will get this special case inlined.
    a == 2 && isa(x, Integer) && return _nextpow2(x)
    a <= 1 && throw(DomainError(a, "`a` must be greater than 1."))
    x <= 1 && return one(a)
    n = ceil(Integer,log(a, x))
    p = a^(n-1)
    # guard against roundoff error, e.g., with a=5 and x=125
    p >= x ? p : a^n
end

"""
    prevpow(a, x)

The largest `a^n` not greater than `x`, where `n` is a non-negative integer.
`a` must be greater than 1, and `x` must not be less than 1.

# Examples
```jldoctest
julia> prevpow(2, 7)
4

julia> prevpow(2, 9)
8

julia> prevpow(5, 20)
5

julia> prevpow(4, 16)
16
```
See also [`nextpow`](@ref).
"""
function prevpow(a::Real, x::Real)
    x < 1 && throw(DomainError(x, "`x` must be ≥ 1."))
    # See comment in nextpos() for a == special case.
    a == 2 && isa(x, Integer) && return _prevpow2(x)
    a <= 1 && throw(DomainError(a, "`a` must be greater than 1."))
    n = floor(Integer,log(a, x))
    p = a^(n+1)
    p <= x ? p : a^n
end

## ndigits (number of digits) in base 10 ##

# decimal digits in an unsigned integer
const powers_of_ten = [
    0x0000000000000001, 0x000000000000000a, 0x0000000000000064, 0x00000000000003e8,
    0x0000000000002710, 0x00000000000186a0, 0x00000000000f4240, 0x0000000000989680,
    0x0000000005f5e100, 0x000000003b9aca00, 0x00000002540be400, 0x000000174876e800,
    0x000000e8d4a51000, 0x000009184e72a000, 0x00005af3107a4000, 0x00038d7ea4c68000,
    0x002386f26fc10000, 0x016345785d8a0000, 0x0de0b6b3a7640000, 0x8ac7230489e80000,
]
function bit_ndigits0z(x::Base.BitUnsigned64)
    lz = (sizeof(x)<<3)-leading_zeros(x)
    nd = (1233*lz)>>12+1
    nd -= x < powers_of_ten[nd]
end
function bit_ndigits0z(x::UInt128)
    n = 0
    while x > 0x8ac7230489e80000
        x = div(x,0x8ac7230489e80000)
        n += 19
    end
    return n + ndigits0z(UInt64(x))
end

ndigits0z(x::BitSigned) = bit_ndigits0z(unsigned(abs(x)))
ndigits0z(x::BitUnsigned) = bit_ndigits0z(x)
ndigits0z(x::Integer) = ndigits0zpb(x, 10)

## ndigits with specified base ##

# The suffix "nb" stands for "negative base"
function ndigits0znb(x::Integer, b::Integer)
    d = 0
    if x isa Unsigned
        d += (x != 0)::Bool
        x = -signed(fld(x, -b))
    end
    # precondition: b < -1 && !(typeof(x) <: Unsigned)
    while x != 0
        x = cld(x,b)
        d += 1
    end
    return d
end

# do first division before conversion with signed here, which can otherwise overflow
ndigits0znb(x::Bool, b::Integer) = x % Int

# The suffix "pb" stands for "positive base"
function ndigits0zpb(x::Integer, b::Integer)
    # precondition: b > 1
    x == 0 && return 0
    b = Int(b)
    x = abs(x)
    if x isa Base.BitInteger
        x = unsigned(x)::Unsigned
        b == 2  && return sizeof(x)<<3 - leading_zeros(x)
        b == 8  && return (sizeof(x)<<3 - leading_zeros(x) + 2) ÷ 3
        b == 16 && return sizeof(x)<<1 - leading_zeros(x)>>2
        b == 10 && return bit_ndigits0z(x)
    end

    d = 0
    while x > typemax(Int)
        x = div(x,b)
        d += 1
    end
    x = div(x,b)
    d += 1

    m = 1
    while m <= x
        m *= b
        d += 1
    end
    return d
end

ndigits0zpb(x::Bool, b::Integer) = x % Int

# The suffix "0z" means that the output is 0 on input zero (cf. #16841)
"""
    ndigits0z(n::Integer, b::Integer=10)

Return 0 if `n == 0`, otherwise compute the number of digits in
integer `n` written in base `b` (i.e. equal to `ndigits(n, base=b)`
in this case).
The base `b` must not be in `[-1, 0, 1]`.

# Examples
```jldoctest
julia> Base.ndigits0z(0, 16)
0

julia> Base.ndigits(0, base=16)
1

julia> Base.ndigits0z(0)
0

julia> Base.ndigits0z(10, 2)
4

julia> Base.ndigits0z(10)
2
```

See also [`ndigits`](@ref).
"""
function ndigits0z(x::Integer, b::Integer)
    if b < -1
        ndigits0znb(x, b)
    elseif b > 1
        ndigits0zpb(x, b)
    else
        throw(DomainError(b, "The base must not be in `[-1, 0, 1]`."))
    end
end

"""
    ndigits(n::Integer; base::Integer=10, pad::Integer=1)

Compute the number of digits in integer `n` written in base `base`
(`base` must not be in `[-1, 0, 1]`), optionally padded with zeros
to a specified size (the result will never be less than `pad`).

# Examples
```jldoctest
julia> ndigits(12345)
5

julia> ndigits(1022, base=16)
3

julia> string(1022, base=16)
"3fe"

julia> ndigits(123, pad=5)
5
```
"""
ndigits(x::Integer; base::Integer=10, pad::Integer=1) = max(pad, ndigits0z(x, base))

## integer to string functions ##

function bin(x::Unsigned, pad::Integer, neg::Bool)
    i = neg + max(pad,sizeof(x)<<3-leading_zeros(x))
    a = StringVector(i)
    while i > neg
        @inbounds a[i] = 48+(x&0x1)
        x >>= 1
        i -= 1
    end
    if neg; @inbounds a[1]=0x2d; end
    String(a)
end

function oct(x::Unsigned, pad::Integer, neg::Bool)
    i = neg + max(pad,div((sizeof(x)<<3)-leading_zeros(x)+2,3))
    a = StringVector(i)
    while i > neg
        @inbounds a[i] = 48+(x&0x7)
        x >>= 3
        i -= 1
    end
    if neg; @inbounds a[1]=0x2d; end
    String(a)
end

function dec(x::Unsigned, pad::Integer, neg::Bool)
    i = neg + ndigits(x, base=10, pad=pad)
    a = StringVector(i)
    while i > neg
        @inbounds a[i] = 48+rem(x,10)
        x = oftype(x,div(x,10))
        i -= 1
    end
    if neg; @inbounds a[1]=0x2d; end
    String(a)
end

function hex(x::Unsigned, pad::Integer, neg::Bool)
    i = neg + max(pad,(sizeof(x)<<1)-(leading_zeros(x)>>2))
    a = StringVector(i)
    while i > neg
        d = x & 0xf
        @inbounds a[i] = 48+d+39*(d>9)
        x >>= 4
        i -= 1
    end
    if neg; @inbounds a[1]=0x2d; end
    String(a)
end

const base36digits = ['0':'9';'a':'z']
const base62digits = ['0':'9';'A':'Z';'a':'z']

function _base(b::Integer, x::Integer, pad::Integer, neg::Bool)
    (x >= 0) | (b < 0) || throw(DomainError(x, "For negative `x`, `b` must be negative."))
    2 <= abs(b) <= 62 || throw(DomainError(b, "base must satisfy 2 ≤ abs(base) ≤ 62"))
    digits = abs(b) <= 36 ? base36digits : base62digits
    i = neg + ndigits(x, base=b, pad=pad)
    a = StringVector(i)
    @inbounds while i > neg
        if b > 0
            a[i] = digits[1+rem(x,b)]
            x = div(x,b)
        else
            a[i] = digits[1+mod(x,-b)]
            x = cld(x,b)
        end
        i -= 1
    end
    if neg; a[1]='-'; end
    String(a)
end

split_sign(n::Integer) = unsigned(abs(n)), n < 0
split_sign(n::Unsigned) = n, false

"""
    string(n::Integer; base::Integer = 10, pad::Integer = 1)

Convert an integer `n` to a string in the given `base`,
optionally specifying a number of digits to pad to.

```jldoctest
julia> string(5, base = 13, pad = 4)
"0005"

julia> string(13, base = 5, pad = 4)
"0023"
```
"""
function string(n::Integer; base::Integer = 10, pad::Integer = 1)
    if base == 2
        (n_positive, neg) = split_sign(n)
        bin(n_positive, pad, neg)
    elseif base == 8
        (n_positive, neg) = split_sign(n)
        oct(n_positive, pad, neg)
    elseif base == 10
        (n_positive, neg) = split_sign(n)
        dec(n_positive, pad, neg)
    elseif base == 16
        (n_positive, neg) = split_sign(n)
        hex(n_positive, pad, neg)
    else
        _base(base, base > 0 ? unsigned(abs(n)) : convert(Signed, n), pad, (base>0) & (n<0))
    end
end

string(b::Bool) = b ? "true" : "false"

"""
    bitstring(n)

A string giving the literal bit representation of a number.

# Examples
```jldoctest
julia> bitstring(4)
"0000000000000000000000000000000000000000000000000000000000000100"

julia> bitstring(2.2)
"0100000000000001100110011001100110011001100110011001100110011010"
```
"""
function bitstring end

bitstring(x::Union{Bool,Int8,UInt8})           = string(reinterpret(UInt8,x), pad = 8, base = 2)
bitstring(x::Union{Int16,UInt16,Float16})      = string(reinterpret(UInt16,x), pad = 16, base = 2)
bitstring(x::Union{Char,Int32,UInt32,Float32}) = string(reinterpret(UInt32,x), pad = 32, base = 2)
bitstring(x::Union{Int64,UInt64,Float64})      = string(reinterpret(UInt64,x), pad = 64, base = 2)
bitstring(x::Union{Int128,UInt128})            = string(reinterpret(UInt128,x), pad = 128, base = 2)

"""
    digits([T<:Integer], n::Integer; base::T = 10, pad::Integer = 1)

Return an array with element type `T` (default `Int`) of the digits of `n` in the given
base, optionally padded with zeros to a specified size. More significant digits are at
higher indices, such that `n == sum([digits[k]*base^(k-1) for k=1:length(digits)])`.

# Examples
```jldoctest
julia> digits(10, base = 10)
2-element Array{Int64,1}:
 0
 1

julia> digits(10, base = 2)
4-element Array{Int64,1}:
 0
 1
 0
 1

julia> digits(10, base = 2, pad = 6)
6-element Array{Int64,1}:
 0
 1
 0
 1
 0
 0
```
"""
digits(n::Integer; base::Integer = 10, pad::Integer = 1) =
    digits(typeof(base), n, base = base, pad = pad)

function digits(T::Type{<:Integer}, n::Integer; base::Integer = 10, pad::Integer = 1)
    digits!(zeros(T, ndigits(n, base=base, pad=pad)), n, base=base)
end

"""
    hastypemax(T::Type) -> Bool

Return `true` if and only if `typemax(T)` is defined.
"""
hastypemax(::Base.BitIntegerType) = true
hastypemax(::Type{T}) where {T} = applicable(typemax, T)

"""
    digits!(array, n::Integer; base::Integer = 10)

Fills an array of the digits of `n` in the given base. More significant digits are at higher
indices. If the array length is insufficient, the least significant digits are filled up to
the array length. If the array length is excessive, the excess portion is filled with zeros.

# Examples
```jldoctest
julia> digits!([2,2,2,2], 10, base = 2)
4-element Array{Int64,1}:
 0
 1
 0
 1

julia> digits!([2,2,2,2,2,2], 10, base = 2)
6-element Array{Int64,1}:
 0
 1
 0
 1
 0
 0
```
"""
function digits!(a::AbstractVector{T}, n::Integer; base::Integer = 10) where T<:Integer
    2 <= abs(base) || throw(DomainError(base, "base must be ≥ 2 or ≤ -2"))
    hastypemax(T) && abs(base) - 1 > typemax(T) &&
        throw(ArgumentError("type $T too small for base $base"))
    isempty(a) && return a

    if base > 0
        for i in eachindex(a)
            n, d = divrem(n, base)
            a[i] = d
        end
    else
        # manually peel one loop iteration for type stability
        n, d = fldmod(n, -base)
        a[firstindex(a)] = d
        n = -signed(n)
        for i in firstindex(a)+1:lastindex(a)
            n, d = fldmod(n, -base)
            a[i] = d
            n = -n
        end
    end
    return a
end

"""
    isqrt(n::Integer)

Integer square root: the largest integer `m` such that `m*m <= n`.

```jldoctest
julia> isqrt(5)
2
```
"""
isqrt(x::Integer) = oftype(x, trunc(sqrt(x)))

function isqrt(x::Union{Int64,UInt64,Int128,UInt128})
    x==0 && return x
    s = oftype(x, trunc(sqrt(x)))
    # fix with a Newton iteration, since conversion to float discards
    # too many bits.
    s = (s + div(x,s)) >> 1
    s*s > x ? s-1 : s
end

"""
    factorial(n::Integer)

Factorial of `n`. If `n` is an [`Integer`](@ref), the factorial is computed as an
integer (promoted to at least 64 bits). Note that this may overflow if `n` is not small,
but you can use `factorial(big(n))` to compute the result exactly in arbitrary precision.

# Examples
```jldoctest
julia> factorial(6)
720

julia> factorial(21)
ERROR: OverflowError: 21 is too large to look up in the table; consider using `factorial(big(21))` instead
Stacktrace:
[...]

julia> factorial(big(21))
51090942171709440000
```

# See also
* [`binomial`](@ref)

# External links
* [Factorial](https://en.wikipedia.org/wiki/Factorial) on Wikipedia.
"""
function factorial(n::Integer)
    n < 0 && throw(DomainError(n, "`n` must be nonnegative."))
    f::typeof(n*n) = 1
    for i::typeof(n*n) = 2:n
        f *= i
    end
    return f
end

"""
    binomial(n::Integer, k::Integer)

The _binomial coefficient_ ``\\binom{n}{k}``, being the coefficient of the ``k``th term in
the polynomial expansion of ``(1+x)^n``.

If ``n`` is non-negative, then it is the number of ways to choose `k` out of `n` items:
```math
\\binom{n}{k} = \\frac{n!}{k! (n-k)!}
```
where ``n!`` is the [`factorial`](@ref) function.

If ``n`` is negative, then it is defined in terms of the identity
```math
\\binom{n}{k} = (-1)^k \\binom{k-n-1}{k}
```

# Examples
```jldoctest
julia> binomial(5, 3)
10

julia> factorial(5) ÷ (factorial(5-3) * factorial(3))
10

julia> binomial(-5, 3)
-35
```

# See also
* [`factorial`](@ref)

# External links
* [Binomial coeffient](https://en.wikipedia.org/wiki/Binomial_coefficient) on Wikipedia.
"""
function binomial(n::T, k::T) where T<:Integer
    n0, k0 = n, k
    k < 0 && return zero(T)
    sgn = one(T)
    if n < 0
        n = -n + k -1
        if isodd(k)
            sgn = -sgn
        end
    end
    k > n && return zero(T)
    (k == 0 || k == n) && return sgn
    k == 1 && return sgn*n
    if k > (n>>1)
        k = (n - k)
    end
    x::T = nn = n - k + 1
    nn += 1
    rr = 2
    while rr <= k
        xt = div(widemul(x, nn), rr)
        x = xt % T
        x == xt || throw(OverflowError("binomial($n0, $k0) overflows"))
        rr += 1
        nn += 1
    end
    convert(T, copysign(x, sgn))
end

max_stars_repo_name: EricForgy/julia

max_issues_repo_name: EricForgy/julia

max_forks_repo_name: EricForgy/julia


text:
module Raytracing
  using CG, Camera, SceneObjects;
  export Ray, generateRay, intersect;

  type Ray
    origin::Vec4f;
    direction::Vec4f;
  end

  function generateRay(camera::PinholeCamera, i::Int, j::Int)
    u = -camera.w / 2 + camera.w * (i - 0.5) / camera.nx
  	v = +camera.h / 2 - camera.h * (j - 0.5) / camera.ny
  	# pixel position
  	p = Vec4f(u, v, 0, 1)
  	# eye position
  	o = Vec4f(0, 0, camera.d, 1)

  	# origin in world coordinates
  	origin = camera.camToWorld * p
  	# unit vector for ray direction
  	direction = camera.camToWorld * (unitize(p-o))
    direction.e4 = 0;
  	return Ray(origin, direction)
  end

  function intersect(ray::Ray, sphere::Sphere)
    d = ray.direction;
    o = ray.origin;
    c = sphere.center;
    r = sphere.radius;

    A = d * d;
    B = 2 * d * (o - c);
    C = (o - c) * (o - c) - r^2;

    radical = B * B - 4 * A * C;
    if radical < 0
      return false, 0.0f0;
    else
      t1 = -B / (2 * A) + sqrt(radical);
      t2 = -B / (2 * A) - sqrt(radical);

      return true, C > 0 ? t2 : t1;
    end
  end

  function intersect(ray::Ray, aabb::AABB)
    tmin::Float32 = -Inf32
  	tmax::Float32 = Inf32
  	p = aabb.center-ray.origin
  	# reduce hiting problem to 1D
  	# for each axis aligned band calculate the intersection points
  	for (f,g,h) in [(p.e1,ray.direction.e1,aabb.hx),(p.e2,ray.direction.e2,aabb.hy),(p.e3,ray.direction.e3,aabb.hz)]
  		t₁ = (f+h)/g
  		t₂ = (f-h)/g
  		if t₁ > t₂
  			tmp = t₁
  			t₁ = t₂
  			t₂ = tmp
  		end
  		t₁ > tmin && (tmin=t₁)
  		t₂ < tmax && (tmax=t₂)
  		# if cube hit this wont hold true
  		if tmin > tmax
  			return false, 0.0f0
  		# if tmax < 0 cube is behind camera
  		elseif tmax < 0
  			return false, 0.0f0
  		end
  	end
  	# cube is in front of us
  	if tmin > 0
  		return true, tmin
  	# cam is inside the cube
  	else
  		return true, tmax
  	end
  end

  function intersect(ray::Ray, scene::Scene)
    tmin = Inf;
    hitObject = nothing;
    for obj in scene.sceneObjects
      hit, t = intersect(ray, obj);
      if (hit && t < tmin)
        tmin = t;
        hitObject = obj;
      end
    end

    return hitObject != nothing, tmin, hitObject;
  end

  function surfaceNormal(ray::Ray, t::Float32, sphere::Sphere)
    p = ray.origin + t * ray.direction;
    n = unitize(p - sphere.center);
    n.e4 = 0;
    return n;
  end

  function surfaceNormal(ray::Ray, t::Float32, aabb::AABB)
    # center of cube
    c = aabb.center
    #point ray hits on aabb
    h = ray.origin+t*ray.direction
    ch = h-c
    x = abs(ch.e1/aabb.hx)
    y = abs(ch.e2/aabb.hy)
    z = abs(ch.e3/aabb.hz)
    if x > y
      if x > z
        return sign(ch.e1/aabb.hx)*Vec4f(1,0,0,0)
      else
        return sign(ch.e3/aabb.hz)*Vec4f(0,0,1,0)
      end
    else
      if y > z
        return sign(ch.e2/aabb.hy)*Vec4f(0,1,0,0)
      else
        return sign(ch.e3/aabb.hz)*Vec4f(0,0,1,0)
      end
    end
  end
end

max_stars_repo_name: fpoppinga/cg

max_issues_repo_name: fpoppinga/cg

max_forks_repo_name: fpoppinga/cg


text:
function makeForcePlots3D()

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

    mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true)

    nspan = (length(mat[1,:]) - 4)/8

    t = mat[:,1]
    len = length(t)
    cl = mat[:,2]
    cd = mat[:,3]
    cm = mat[:,4]

    #Plot wing coefficients
    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(L"$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(L"$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(L"$t^*$")
    ylabel(L"$C_M$")
    savefig("forcePlots/cm.png")
    close()

    #Co-plot kinematics for the various sections
    alpha_all = Float64[]
    h_all = Float64[]
    for i = 1:nspan
        alpha = mat[:,Int((i-1)*8+5)]*180/pi
        alpha_all = [alpha_all; alpha]
        h = mat[:,Int((i-1)*8+6)]
        h_all = [h_all; h]
    end
    for i = 1:nspan
        alpha = mat[:,Int((i-1)*8+5)]*180/pi
        plot(t, alpha)
    end
    xmin = t[1]
    xmax = t[end]
    zmin = minimum(alpha_all) - 0.1*abs(minimum(alpha_all))
    zmax = maximum(alpha_all) + 0.1*abs(maximum(alpha_all))
    axis([xmin, xmax, zmin, zmax])
    xlabel(L"$t^*$")
    ylabel(L"$\alpha$ (deg)")
    savefig("forcePlots/pitch-spanwise.png")
    close()

    for i = 1:nspan
        h = mat[:,Int((i-1)*8+6)]
        plot(t, h)
    end
    xmin = t[1]
    xmax = t[end]
    zmin = minimum(h_all) - 0.1*abs(minimum(h_all))
    zmax = maximum(h_all) + 0.1*abs(maximum(h_all))
    axis([xmin, xmax, zmin, zmax])
    xlabel(L"$t^*$")
    ylabel(L"$h/c$")
    savefig("forcePlots/plunge-spanwise.png")
    close()

    #Coplot spanwise force coefficients
    cl_all = Float64[]
    cd_all = Float64[]
    cm_all = Float64[]
    for i = 1:nspan
        cl = mat[range,Int((i-1)*8+9)]./mat[range,2]
        cl_all = [cl_all; cl]
        cd = mat[range,Int((i-1)*8+10)]./mat[range,3]
        cd_all = [cd_all; cd]
        cm = mat[range,Int((i-1)*8+11)]./mat[range,4]
        cm_all = [cm_all; cm]
    end

    for i = 1:nspan
        cl = mat[:,Int((i-1)*8+9)]./mat[:,2]
        plot(t, cl)
    end
    xmin = t[1]
    xmax = t[end]
    zmin = minimum(cl_all[:]) - 0.1*abs(minimum(cl_all[:]))
    zmax = maximum(cl_all[:]) + 0.1*abs(maximum(cl_all[:]))
    axis([xmin, xmax, zmin, zmax])
    xlabel(L"$t^*$")
    ylabel(L"$C_l c/C_L c_{ref}$")
    savefig("forcePlots/cl-spanwise.png")
    close()

    for i = 1:nspan
        cd = mat[:,Int((i-1)*8+10)]./mat[:,3]
        plot(t, cd)
    end
    xmin = t[1]
    xmax = t[end]
    zmin = minimum(cd_all[:]) - 0.1*abs(minimum(cd_all[:]))
    zmax = maximum(cd_all[:]) + 0.1*abs(maximum(cd_all[:]))
    axis([xmin, xmax, zmin, zmax])
    xlabel(L"$t^*$")
    ylabel(L"$C_d c/C_D c_{ref}$")
    savefig("forcePlots/cd-spanwise.png")
    close()

    for i = 1:nspan
        cm = mat[:,Int((i-1)*8+11)]./mat[:,4]
        plot(t, cm)
    end
    xmin = t[1]
    xmax = t[end]
    zmin = minimum(cm_all[:]) - 0.1*abs(minimum(cm_all[:]))
    zmax = maximum(cm_all[:]) + 0.1*abs(maximum(cm_all[:]))
    axis([xmin, xmax, zmin, zmax])
    xlabel(L"$t^*$")
    ylabel(L"$C_m c^2/C_M c_{ref}^2$")
    savefig("forcePlots/cm-spanwise.png")
    close()

    #Coplot spanwise a0 variation
    a0_all = Float64[]
    a03d_all = Float64[]
    for i = 1:nspan
        a0 = mat[:,Int((i-1)*8+8)]
        a0_all = [a0_all; a0]
        a03d = mat[:,Int((i-1)*8+12)]
        a03d_all = [a03d_all; a03d]
    end
    for i = 1:nspan
        a0 = mat[:,Int((i-1)*8+8)]
        plot(t, a0)
    end
    xmin = t[1]
    xmax = t[end]
zmin = minimum(a0_all[:]) - 0.1*abs(minimum(a0_all[:]))
zmax = maximum(a0_all[:]) + 0.1*abs(maximum(a0_all[:]))
axis([xmin, xmax, zmin, zmax])
xlabel(L"$t^*$")
ylabel(L"$A_0$")
savefig("forcePlots/a0-spanwise.png")
close()

for i = 1:nspan
    a03d = mat[:,Int((i-1)*8+12)]
    plot(t, a03d)
end
xmin = t[1]
xmax = t[end]
zmin = minimum(a03d_all[:]) - 0.1*abs(minimum(a03d_all[:]))
zmax = maximum(a03d_all[:]) + 0.1*abs(maximum(a03d_all[:]))
axis([xmin, xmax, zmin, zmax])
xlabel(L"$t^*$")
ylabel(L"$A_0$")
savefig("forcePlots/a03d-spanwise.png")
close()

end

function makeVortPlots3D()

    dirvec = readdir()
    dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)
    #Determine axis limits
    dirmax = maximum(dirresults)
    cd("$(dirmax)")

    nspan = 0
    for i = 1:1000
        if isfile("boundv-$i") == true
            nspan += 1
        end
    end
    tip_ind = 1
    root_ind = nspan
    mid_ind= floor(Int, nspan/2)

    tev, _ = DelimitedFiles.readdlm("tev-1", '\t', Float64, header=true)
    if length(tev[:,1]) > 1
        tev = tev[2:end,:]
    end
    lev =   try
        DelimitedFiles.readdlm("lev-1", '\t', Float64, header=true)[1]
    catch
        zeros(0,3)
    end
    if length(lev[:,1]) > 1
        lev = lev[2:end,:]
    end
    bv, _ = DelimitedFiles.readdlm("boundv-1", '\t', Float64, header=true)
    for i in [mid_ind; root_ind]
        tevt, _ = DelimitedFiles.readdlm("tev-$i", '\t', Float64, header=true)
        if length(tevt[:,1]) > 1
            tevt = tevt[2:end,:]
        end
        tev = [tev; tevt]
        levt = try
            DelimitedFiles.readdlm("lev-$i", '\t', Float64, header=true)[1]
        catch
            zeros(0,3)
        end
        if length(levt[:,1]) > 1
            levt = levt[2:end,:]
        end
        lev = [lev; levt]
        bvt, _ = DelimitedFiles.readdlm("boundv-$i", '\t', Float64, header=true)
        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 "vortPlots" in dirvec
        rm("vortPlots", recursive=true)
    end
    mkdir("vortPlots")

    for i=1:length(dirresults)
        if dirresults[i] != 0
            dirstr="$(dirresults[i])"
            cd(dirstr)

            figure()

            tev, _ = DelimitedFiles.readdlm("tev-$root_ind", '\t', Float64, header=true)
            lev =   try
                DelimitedFiles.readdlm("lev-$root_ind", '\t', Float64, header=true)[1]
            catch
                zeros(0,3)
            end
            bv, _ = DelimitedFiles.readdlm("boundv-$root_ind", '\t', Float64, header=true)
            subplot(311)
            viewVort2D(tev, lev, bv)
            axis([xmin-1, xmax+1, zmin-1, zmax+1])
            ylabel("Root")

            tev, _ = DelimitedFiles.readdlm("tev-$mid_ind", '\t', Float64, header=true)
            lev =   try
                DelimitedFiles.readdlm("lev-$mid_ind", '\t', Float64, header=true)[1]
            catch
                zeros(0,3)
            end
            bv, _ = DelimitedFiles.readdlm("boundv-$mid_ind", '\t', Float64, header=true)
            subplot(312)
            viewVort2D(tev, lev, bv)
            axis([xmin-1, xmax+1, zmin-1, zmax+1])
            ylabel("Mid")

            tev, _ = DelimitedFiles.readdlm("tev-$tip_ind", '\t', Float64, header=true)
            lev =   try
                DelimitedFiles.readdlm("lev-$tip_ind", '\t', Float64, header=true)[1]
            catch
                zeros(0,3)
            end
            bv, _ = DelimitedFiles.readdlm("boundv-$tip_ind", '\t', Float64, header=true)
            subplot(313)
            viewVort2D(tev, lev, bv)
            axis([xmin-1, xmax+1, zmin-1, zmax+1])
            ylabel("Tip")

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

function makeInfoPlots3D()

    mat, _ = DelimitedFiles.readdlm("resultsSummary", '\t', Float64, header=true)
    nspan = Int((length(mat[1,:]) - 4)/8)

    t = mat[:,1]
    len = length(t)
    cl = mat[:,2]
    cdr = mat[:,3]
    cm = mat[:,4]

    dirvec = readdir()
    dirresults = map(x->(v = tryparse(Float64,x); typeof(v) == Nothing ? 0.0 : v),dirvec)

    #Determine axis limits
    dirmax = maximum(dirresults)

    cd("$(dirmax)")

    tip_ind = 1
    root_ind = nspan
    mid_ind= floor(Int, nspan/2)

    #Axis limits for vorticity plots based on last time step
    tev, _ = DelimitedFiles.readdlm("tev-1", '\t', Float64, header=true)
    if length(tev[:,1]) > 1
        tev = tev[2:end,:]
    end
    lev =   try
        DelimitedFiles.readdlm("lev-1", '\t', Float64, header=true)[1]
    catch
        zeros(0,3)
    end
    lev, _ = DelimitedFiles.readdlm("lev-1", '\t', Float64, header=true)
    if length(lev[:,1]) > 1
        lev = lev[2:end,:]
    end
    bv, _ = DelimitedFiles.readdlm("boundv-1", '\t', Float64, header=true)
    for i in [mid_ind; root_ind]
        tevt, _ = DelimitedFiles.readdlm("tev-$i", '\t', Float64, header=true)
        if length(tevt[:,1]) > 1
            tevt = tevt[2:end,:]
        end
        tev = [tev; tevt]
        levt = try
            DelimitedFiles.readdlm("lev-$i", '\t', Float64, header=true)[1]
        catch
            zeros(0,3)
        end
        if length(levt[:,1]) > 1
            levt = levt[2:end,:]
        end
        lev = [lev; levt]
        bvt, _ = DelimitedFiles.readdlm("boundv-$i", '\t', Float64, header=true)
        bv = [bv; bvt]
    end
    cd("..")

    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];])

    #Axis limits for spanwise variations based on all the time steps
    vorstr = Float64[]
    a0str = Float64[]
    gamstr = Float64[]
    for i=1:length(dirresults)
        if dirresults[i] != 0
            dirstr="$(dirresults[i])"
            cd(dirstr)

            matind = argmin(abs.(t .- dirresults[i]))

            spanvar, _  = DelimitedFiles.readdlm("spanwise-var", '\t', Float64, header=true)
            vorstr = [vorstr; spanvar[:,2]; spanvar[:,3]]
            a0str = [a0str; spanvar[:,5]; spanvar[:,6]; spanvar[:,7]]
            cl_load = mat[matind, ([1:nspan;].-1).*8 .+ 9]
            gamstr = [gamstr; spanvar[:,4]/spanvar[end,4]; cl_load]
            cd("..")
        end
    end
    vormin = minimum(vorstr) - 0.1*abs(minimum(vorstr))
    a0min = minimum(a0str) - 0.1*abs(minimum(a0str))
    gammin = minimum(gamstr) - 0.1*abs(minimum(gamstr))

    vormax = maximum(vorstr) + 0.1*abs(maximum(vorstr))
    a0max = maximum(a0str) + 0.1*abs(maximum(a0str))
    gammax = maximum(gamstr) + 0.1*abs(maximum(gamstr))

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

    for i=1:length(dirresults)
        if dirresults[i] != 0
            dirstr="$(dirresults[i])"
            cd(dirstr)

            figure(figsize=(12,8))

            #Kinematics
            alpha_all = Float64[]
            h_all = Float64[]
            for i = 1:nspan
                alpha = mat[:,Int((i-1)*8+5)]*180/pi
                alpha_all = [alpha_all; alpha]
                h = mat[:,Int((i-1)*8+6)]
                h_all = [h_all; h]
            end

            subplot2grid((4, 5), (0, 0))
            for i = 1:nspan
                alpha = mat[:,Int((i-1)*8+5)]*180/pi
                plot(t, alpha)
            end
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(alpha_all) - 0.1*abs(minimum(alpha_all))
            zmax = maximum(alpha_all) + 0.1*abs(maximum(alpha_all))
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$t^*$")
            ylabel(L"$\alpha$ (deg)")

            subplot2grid((4, 5), (0, 1))
            for i = 1:nspan
                h = mat[:,Int((i-1)*8+6)]
                plot(t, h)
            end
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(h_all) - 0.1*abs(minimum(h_all))
            zmax = maximum(h_all) + 0.1*abs(maximum(h_all))
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$t^*$")
            ylabel(L"$h/c$")

            subplot2grid((4, 5), (1, 0), colspan=4)
            tev, _ = DelimitedFiles.readdlm("tev-$root_ind", '\t', Float64, header=true)
            lev =   try
                DelimitedFiles.readdlm("lev-$root_ind", '\t', Float64, header=true)[1]
            catch
                zeros(0,3)
            end
            bv, _ = DelimitedFiles.readdlm("boundv-$root_ind", '\t', Float64, header=true)
            viewVort2D(tev, lev, bv)
            axis([xminv-1, xmaxv+1, zminv-1, zmaxv+1])
            ylabel("Root")

            subplot2grid((4, 5), (2, 0), colspan=4)
            tev, _ = DelimitedFiles.readdlm("tev-$mid_ind", '\t', Float64, header=true)
            lev =   try
                DelimitedFiles.readdlm("lev-$mid_ind", '\t', Float64, header=true)[1]
            catch
                zeros(0,3)
            end
            bv, _ = DelimitedFiles.readdlm("boundv-$mid_ind", '\t', Float64, header=true)
            viewVort2D(tev, lev, bv)
            axis([xminv-1, xmaxv+1, zminv-1, zmaxv+1])
            ylabel("Mid")

            subplot2grid((4, 5), (3, 0), colspan=4)
            tev, _ = DelimitedFiles.readdlm("tev-$tip_ind", '\t', Float64, header=true)
            lev =   try
                DelimitedFiles.readdlm("lev-$tip_ind", '\t', Float64, header=true)[1]
            catch
                zeros(0,3)
            end
            bv, _ = DelimitedFiles.readdlm("boundv-$tip_ind", '\t', Float64, header=true)
            viewVort2D(tev, lev, bv)
            axis([xminv-1, xmaxv+1, zminv-1, zmaxv+1])
            ylabel("Tip")

            #Force plots
            subplot2grid((4, 5), (0, 2))
            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(L"$t^*$")
            ylabel(L"$C_L$")

            subplot2grid((4, 5), (0, 3))
            plot(t, cdr)
            range = Int(round(0.05*len)):Int(round(0.95*len))
            xmin = t[1]
            xmax = t[end]
            zmin = minimum(cdr[range]) - 0.1*abs(minimum(cdr[range]))
            zmax = maximum(cdr[range]) + 0.1*abs(maximum(cdr[range]))
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$t^*$")
            ylabel(L"$C_D$")

            subplot2grid((4, 5), (0, 4))
            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(L"$t^*$")
            ylabel(L"$C_M$")

            #Spanwise variations
            spanvar, _ = DelimitedFiles.readdlm("spanwise-var", '\t', Float64, header=true)

            subplot2grid((4, 5), (1, 4))
            plot(spanvar[:,1], spanvar[:,2])
            plot(spanvar[:,1], spanvar[:,3])
            xmin = minimum(spanvar[:,1])
            xmax = maximum(spanvar[:,1])
            zmin = vormin
            zmax = vormax
            axis([xmin, xmax, zmin, zmax])
            xlabel(L"$y_{LE}$")
            ylabel(L"tev, lev")

            subplot2grid((4, 5), (2, 4))
            matind = argmin(abs.(t .- dirresults[i]))
            cl_load = [0.; mat[matind, ([1:nspan;].-1).*8 .+ 9]]
            plot(spanvar[:,1], spanvar[:,4]/spanvar[end,4])
            plot(spanvar[:,1], cl_load)
xmin = minimum(spanvar[:,1])
xmax = maximum(spanvar[:,1])
zmin = gammin
zmax = gammax
axis([xmin, xmax, zmin, zmax])
xlabel(L"$y_{LE}$")
ylabel(L"$\Gamma / \Gamma_r$, $C_l c/C_L c_{ref}$")

subplot2grid((4, 5), (3, 4))
plot(spanvar[:,1], spanvar[:,5])
plot(spanvar[:,1], spanvar[:,6])
plot(spanvar[:,1], spanvar[:,7])
xmin = minimum(spanvar[:,1])
xmax = maximum(spanvar[:,1])
zmin = a0min
zmax = a0max
axis([xmin, xmax, zmin, zmax])
xlabel(L"$y_{LE}$")
ylabel(L"$A0_{2D}$, $A0_{3D}$, $A0_{tot}$")

tight_layout()

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