text:
# Julia wrapper for header: stb_image.h
# Automatically generated using Clang.jl

function stbi_load_from_memory(buffer, len, x, y, channels_in_file, desired_channels)
    ccall((:stbi_load_from_memory, libstbimage), Ptr{stbi_uc}, (Ptr{stbi_uc}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), buffer, len, x, y, channels_in_file, desired_channels)
end

function stbi_load_from_callbacks(clbk, user, x, y, channels_in_file, desired_channels)
    ccall((:stbi_load_from_callbacks, libstbimage), Ptr{stbi_uc}, (Ptr{stbi_io_callbacks}, Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), clbk, user, x, y, channels_in_file, desired_channels)
end

function stbi_load(filename, x, y, channels_in_file, desired_channels)
    ccall((:stbi_load, libstbimage), Ptr{stbi_uc}, (Cstring, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), filename, x, y, channels_in_file, desired_channels)
end

function stbi_load_from_file(f, x, y, channels_in_file, desired_channels)
    ccall((:stbi_load_from_file, libstbimage), Ptr{stbi_uc}, (Ptr{FILE}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), f, x, y, channels_in_file, desired_channels)
end

function stbi_load_gif_from_memory(buffer, len, delays, x, y, z, comp, req_comp)
    ccall((:stbi_load_gif_from_memory, libstbimage), Ptr{stbi_uc}, (Ptr{stbi_uc}, Cint, Ptr{Ptr{Cint}}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), buffer, len, delays, x, y, z, comp, req_comp)
end

function stbi_load_16_from_memory(buffer, len, x, y, channels_in_file, desired_channels)
    ccall((:stbi_load_16_from_memory, libstbimage), Ptr{stbi_us}, (Ptr{stbi_uc}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), buffer, len, x, y, channels_in_file, desired_channels)
end

function stbi_load_16_from_callbacks(clbk, user, x, y, channels_in_file, desired_channels)
    ccall((:stbi_load_16_from_callbacks, libstbimage), Ptr{stbi_us}, (Ptr{stbi_io_callbacks}, Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), clbk, user, x, y, channels_in_file, desired_channels)
end

function stbi_load_16(filename, x, y, channels_in_file, desired_channels)
    ccall((:stbi_load_16, libstbimage), Ptr{stbi_us}, (Cstring, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), filename, x, y, channels_in_file, desired_channels)
end

function stbi_load_from_file_16(f, x, y, channels_in_file, desired_channels)
    ccall((:stbi_load_from_file_16, libstbimage), Ptr{stbi_us}, (Ptr{FILE}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), f, x, y, channels_in_file, desired_channels)
end

function stbi_loadf_from_memory(buffer, len, x, y, channels_in_file, desired_channels)
    ccall((:stbi_loadf_from_memory, libstbimage), Ptr{Cfloat}, (Ptr{stbi_uc}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), buffer, len, x, y, channels_in_file, desired_channels)
end

function stbi_loadf_from_callbacks(clbk, user, x, y, channels_in_file, desired_channels)
    ccall((:stbi_loadf_from_callbacks, libstbimage), Ptr{Cfloat}, (Ptr{stbi_io_callbacks}, Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), clbk, user, x, y, channels_in_file, desired_channels)
end

function stbi_loadf(filename, x, y, channels_in_file, desired_channels)
    ccall((:stbi_loadf, libstbimage), Ptr{Cfloat}, (Cstring, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), filename, x, y, channels_in_file, desired_channels)
end

function stbi_loadf_from_file(f, x, y, channels_in_file, desired_channels)
    ccall((:stbi_loadf_from_file, libstbimage), Ptr{Cfloat}, (Ptr{FILE}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}, Cint), f, x, y, channels_in_file, desired_channels)
end

function stbi_hdr_to_ldr_gamma(gamma)
    ccall((:stbi_hdr_to_ldr_gamma, libstbimage), Cvoid, (Cfloat,), gamma)
end

function stbi_hdr_to_ldr_scale(scale)
    ccall((:stbi_hdr_to_ldr_scale, libstbimage), Cvoid, (Cfloat,), scale)
end

function stbi_ldr_to_hdr_gamma(gamma)
    ccall((:stbi_ldr_to_hdr_gamma, libstbimage), Cvoid, (Cfloat,), gamma)
end

function stbi_ldr_to_hdr_scale(scale)
    ccall((:stbi_ldr_to_hdr_scale, libstbimage), Cvoid, (Cfloat,), scale)
end

function stbi_is_hdr_from_callbacks(clbk, user)
    ccall((:stbi_is_hdr_from_callbacks, libstbimage), Cint, (Ptr{stbi_io_callbacks}, Ptr{Cvoid}), clbk, user)
end

function stbi_is_hdr_from_memory(buffer, len)
    ccall((:stbi_is_hdr_from_memory, libstbimage), Cint, (Ptr{stbi_uc}, Cint), buffer, len)
end

function stbi_is_hdr(filename)
    ccall((:stbi_is_hdr, libstbimage), Cint, (Cstring,), filename)
end

function stbi_is_hdr_from_file(f)
    ccall((:stbi_is_hdr_from_file, libstbimage), Cint, (Ptr{FILE},), f)
end

function stbi_failure_reason()
    ccall((:stbi_failure_reason, libstbimage), Cstring, ())
end

function stbi_image_free(retval_from_stbi_load)
    ccall((:stbi_image_free, libstbimage), Cvoid, (Ptr{Cvoid},), retval_from_stbi_load)
end

function stbi_info_from_memory(buffer, len, x, y, comp)
    ccall((:stbi_info_from_memory, libstbimage), Cint, (Ptr{stbi_uc}, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), buffer, len, x, y, comp)
end

function stbi_info_from_callbacks(clbk, user, x, y, comp)
    ccall((:stbi_info_from_callbacks, libstbimage), Cint, (Ptr{stbi_io_callbacks}, Ptr{Cvoid}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), clbk, user, x, y, comp)
end

function stbi_is_16_bit_from_memory(buffer, len)
    ccall((:stbi_is_16_bit_from_memory, libstbimage), Cint, (Ptr{stbi_uc}, Cint), buffer, len)
end

function stbi_is_16_bit_from_callbacks(clbk, user)
    ccall((:stbi_is_16_bit_from_callbacks, libstbimage), Cint, (Ptr{stbi_io_callbacks}, Ptr{Cvoid}), clbk, user)
end

function stbi_info(filename, x, y, comp)
    ccall((:stbi_info, libstbimage), Cint, (Cstring, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), filename, x, y, comp)
end

function stbi_info_from_file(f, x, y, comp)
    ccall((:stbi_info_from_file, libstbimage), Cint, (Ptr{FILE}, Ptr{Cint}, Ptr{Cint}, Ptr{Cint}), f, x, y, comp)
end

function stbi_is_16_bit(filename)
    ccall((:stbi_is_16_bit, libstbimage), Cint, (Cstring,), filename)
end

function stbi_is_16_bit_from_file(f)
    ccall((:stbi_is_16_bit_from_file, libstbimage), Cint, (Ptr{FILE},), f)
end

function stbi_set_unpremultiply_on_load(flag_true_if_should_unpremultiply)
    ccall((:stbi_set_unpremultiply_on_load, libstbimage), Cvoid, (Cint,), flag_true_if_should_unpremultiply)
end

function stbi_convert_iphone_png_to_rgb(flag_true_if_should_convert)
    ccall((:stbi_convert_iphone_png_to_rgb, libstbimage), Cvoid, (Cint,), flag_true_if_should_convert)
end

function stbi_set_flip_vertically_on_load(flag_true_if_should_flip)
    ccall((:stbi_set_flip_vertically_on_load, libstbimage), Cvoid, (Cint,), flag_true_if_should_flip)
end

function stbi_zlib_decode_malloc_guesssize(buffer, len, initial_size, outlen)
    ccall((:stbi_zlib_decode_malloc_guesssize, libstbimage), Cstring, (Cstring, Cint, Cint, Ptr{Cint}), buffer, len, initial_size, outlen)
end

function stbi_zlib_decode_malloc_guesssize_headerflag(buffer, len, initial_size, outlen, parse_header)
    ccall((:stbi_zlib_decode_malloc_guesssize_headerflag, libstbimage), Cstring, (Cstring, Cint, Cint, Ptr{Cint}, Cint), buffer, len, initial_size, outlen, parse_header)
end

function stbi_zlib_decode_malloc(buffer, len, outlen)
    ccall((:stbi_zlib_decode_malloc, libstbimage), Cstring, (Cstring, Cint, Ptr{Cint}), buffer, len, outlen)
end

function stbi_zlib_decode_buffer(obuffer, olen, ibuffer, ilen)
    ccall((:stbi_zlib_decode_buffer, libstbimage), Cint, (Cstring, Cint, Cstring, Cint), obuffer, olen, ibuffer, ilen)
end

function stbi_zlib_decode_noheader_malloc(buffer, len, outlen)
    ccall((:stbi_zlib_decode_noheader_malloc, libstbimage), Cstring, (Cstring, Cint, Ptr{Cint}), buffer, len, outlen)
end

function stbi_zlib_decode_noheader_buffer(obuffer, olen, ibuffer, ilen)
    ccall((:stbi_zlib_decode_noheader_buffer, libstbimage), Cint, (Cstring, Cint, Cstring, Cint), obuffer, olen, ibuffer, ilen)
end
# Julia wrapper for header: stb_image_write.h
# Automatically generated using Clang.jl

function stbi_write_png(filename, w, h, comp, data, stride_in_bytes)
    ccall((:stbi_write_png, libstbimage), Cint, (Cstring, Cint, Cint, Cint, Ptr{Cvoid}, Cint), filename, w, h, comp, data, stride_in_bytes)
end

function stbi_write_bmp(filename, w, h, comp, data)
    ccall((:stbi_write_bmp, libstbimage), Cint, (Cstring, Cint, Cint, Cint, Ptr{Cvoid}), filename, w, h, comp, data)
end

function stbi_write_tga(filename, w, h, comp, data)
    ccall((:stbi_write_tga, libstbimage), Cint, (Cstring, Cint, Cint, Cint, Ptr{Cvoid}), filename, w, h, comp, data)
end

function stbi_write_hdr(filename, w, h, comp, data)
    ccall((:stbi_write_hdr, libstbimage), Cint, (Cstring, Cint, Cint, Cint, Ptr{Cfloat}), filename, w, h, comp, data)
end

function stbi_write_jpg(filename, x, y, comp, data, quality)
    ccall((:stbi_write_jpg, libstbimage), Cint, (Cstring, Cint, Cint, Cint, Ptr{Cvoid}, Cint), filename, x, y, comp, data, quality)
end

function stbi_write_png_to_func(func, context, w, h, comp, data, stride_in_bytes)
    ccall((:stbi_write_png_to_func, libstbimage), Cint, (Ptr{stbi_write_func}, Ptr{Cvoid}, Cint, Cint, Cint, Ptr{Cvoid}, Cint), func, context, w, h, comp, data, stride_in_bytes)
end

function stbi_write_bmp_to_func(func, context, w, h, comp, data)
    ccall((:stbi_write_bmp_to_func, libstbimage), Cint, (Ptr{stbi_write_func}, Ptr{Cvoid}, Cint, Cint, Cint, Ptr{Cvoid}), func, context, w, h, comp, data)
end

function stbi_write_tga_to_func(func, context, w, h, comp, data)
    ccall((:stbi_write_tga_to_func, libstbimage), Cint, (Ptr{stbi_write_func}, Ptr{Cvoid}, Cint, Cint, Cint, Ptr{Cvoid}), func, context, w, h, comp, data)
end

function stbi_write_hdr_to_func(func, context, w, h, comp, data)
    ccall((:stbi_write_hdr_to_func, libstbimage), Cint, (Ptr{stbi_write_func}, Ptr{Cvoid}, Cint, Cint, Cint, Ptr{Cfloat}), func, context, w, h, comp, data)
end

function stbi_write_jpg_to_func(func, context, x, y, comp, data, quality)
    ccall((:stbi_write_jpg_to_func, libstbimage), Cint, (Ptr{stbi_write_func}, Ptr{Cvoid}, Cint, Cint, Cint, Ptr{Cvoid}, Cint), func, context, x, y, comp, data, quality)
end

function stbi_flip_vertically_on_write(flip_boolean)
    ccall((:stbi_flip_vertically_on_write, libstbimage), Cvoid, (Cint,), flip_boolean)
end
# Julia wrapper for header: stb_image_resize.h
# Automatically generated using Clang.jl

function stbir_resize_uint8(input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels)
    ccall((:stbir_resize_uint8, libstbimage), Cint, (Ptr{Cuchar}, Cint, Cint, Cint, Ptr{Cuchar}, Cint, Cint, Cint, Cint), input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels)
end

function stbir_resize_float(input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels)
    ccall((:stbir_resize_float, libstbimage), Cint, (Ptr{Cfloat}, Cint, Cint, Cint, Ptr{Cfloat}, Cint, Cint, Cint, Cint), input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels)
end

function stbir_resize_uint8_srgb(input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels, alpha_channel, flags)
    ccall((:stbir_resize_uint8_srgb, libstbimage), Cint, (Ptr{Cuchar}, Cint, Cint, Cint, Ptr{Cuchar}, Cint, Cint, Cint, Cint, Cint, Cint), input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels, alpha_channel, flags)
end

function stbir_resize_uint8_srgb_edgemode(input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels, alpha_channel, flags, edge_wrap_mode)
    ccall((:stbir_resize_uint8_srgb_edgemode, libstbimage), Cint, (Ptr{Cuchar}, Cint, Cint, Cint, Ptr{Cuchar}, Cint, Cint, Cint, Cint, Cint, Cint, stbir_edge), input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels, alpha_channel, flags, edge_wrap_mode)
end

function stbir_resize_uint8_generic(input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels, alpha_channel, flags, edge_wrap_mode, filter, space, alloc_context)
    ccall((:stbir_resize_uint8_generic, libstbimage), Cint, (Ptr{Cuchar}, Cint, Cint, Cint, Ptr{Cuchar}, Cint, Cint, Cint, Cint, Cint, Cint, stbir_edge, stbir_filter, stbir_colorspace, Ptr{Cvoid}), input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels, alpha_channel, flags, edge_wrap_mode, filter, space, alloc_context)
end

function stbir_resize_uint16_generic(input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels, alpha_channel, flags, edge_wrap_mode, filter, space, alloc_context)
    ccall((:stbir_resize_uint16_generic, libstbimage), Cint, (Ptr{stbir_uint16}, Cint, Cint, Cint, Ptr{stbir_uint16}, Cint, Cint, Cint, Cint, Cint, Cint, stbir_edge, stbir_filter, stbir_colorspace, Ptr{Cvoid}), input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels, alpha_channel, flags, edge_wrap_mode, filter, space, alloc_context)
end

function stbir_resize_float_generic(input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels, alpha_channel, flags, edge_wrap_mode, filter, space, alloc_context)
    ccall((:stbir_resize_float_generic, libstbimage), Cint, (Ptr{Cfloat}, Cint, Cint, Cint, Ptr{Cfloat}, Cint, Cint, Cint, Cint, Cint, Cint, stbir_edge, stbir_filter, stbir_colorspace, Ptr{Cvoid}), input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, num_channels, alpha_channel, flags, edge_wrap_mode, filter, space, alloc_context)
end

function stbir_resize(input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, datatype, num_channels, alpha_channel, flags, edge_mode_horizontal, edge_mode_vertical, filter_horizontal, filter_vertical, space, alloc_context)
    ccall((:stbir_resize, libstbimage), Cint, (Ptr{Cvoid}, Cint, Cint, Cint, Ptr{Cvoid}, Cint, Cint, Cint, stbir_datatype, Cint, Cint, Cint, stbir_edge, stbir_edge, stbir_filter, stbir_filter, stbir_colorspace, Ptr{Cvoid}), input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, datatype, num_channels, alpha_channel, flags, edge_mode_horizontal, edge_mode_vertical, filter_horizontal, filter_vertical, space, alloc_context)
end

function stbir_resize_subpixel(input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, datatype, num_channels, alpha_channel, flags, edge_mode_horizontal, edge_mode_vertical, filter_horizontal, filter_vertical, space, alloc_context, x_scale, y_scale, x_offset, y_offset)
    ccall((:stbir_resize_subpixel, libstbimage), Cint, (Ptr{Cvoid}, Cint, Cint, Cint, Ptr{Cvoid}, Cint, Cint, Cint, stbir_datatype, Cint, Cint, Cint, stbir_edge, stbir_edge, stbir_filter, stbir_filter, stbir_colorspace, Ptr{Cvoid}, Cfloat, Cfloat, Cfloat, Cfloat), input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, datatype, num_channels, alpha_channel, flags, edge_mode_horizontal, edge_mode_vertical, filter_horizontal, filter_vertical, space, alloc_context, x_scale, y_scale, x_offset, y_offset)
end

function stbir_resize_region(input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, datatype, num_channels, alpha_channel, flags, edge_mode_horizontal, edge_mode_vertical, filter_horizontal, filter_vertical, space, alloc_context, s0, t0, s1, t1)
    ccall((:stbir_resize_region, libstbimage), Cint, (Ptr{Cvoid}, Cint, Cint, Cint, Ptr{Cvoid}, Cint, Cint, Cint, stbir_datatype, Cint, Cint, Cint, stbir_edge, stbir_edge, stbir_filter, stbir_filter, stbir_colorspace, Ptr{Cvoid}, Cfloat, Cfloat, Cfloat, Cfloat), input_pixels, input_w, input_h, input_stride_in_bytes, output_pixels, output_w, output_h, output_stride_in_bytes, datatype, num_channels, alpha_channel, flags, edge_mode_horizontal, edge_mode_vertical, filter_horizontal, filter_vertical, space, alloc_context, s0, t0, s1, t1)
end


text:
if v"1.0" <= VERSION < v"1.1"
    isnothing(x) = x===nothing
end

ext: jl


text:
using PowerFlowData
using Documenter
using Tables

DocMeta.setdocmeta!(PowerFlowData, :DocTestSetup, :(using PowerFlowData); recursive=true)

makedocs(;
    modules=[PowerFlowData],
    authors="Nick Robinson <npr251@gmail.com> and contributors",
    repo="https://github.com/nickrobinson251/PowerFlowData.jl/blob/{commit}{path}#{line}",
    sitename="PowerFlowData.jl",
    format=Documenter.HTML(
        canonical="https://nickrobinson251.github.io/PowerFlowData.jl",
        prettyurls=false,
    ),
    pages=[
        "Home" => "index.md",
        "API" => "api.md",
        "Alternatives" => "alternatives.md",
        "Implementation" => "implementation.md",
    ],
    strict=true,
    checkdocs=:exports,
)

deploydocs(;
    repo="github.com/nickrobinson251/PowerFlowData.jl",
    devbranch="main",
    push_preview=true,
)


text:
module CirculantMatrices
using FFTW
using FFTW: Plan
using LinearAlgebra

import Base: convert, eltype, *, \, /, getindex, print_matrix, size, Matrix, +, copy, -,
             real, imag, conj

#import Base: copyto!
import LinearAlgebra: mul!, ldiv!, inv, norm, eigvals, eigvecs, eigen, adjoint, cond,
                     issymmetric, ishermitian, rank, axpy!, adjoint!, transpose!

import LinearAlgebra: rmul!, lmul!

# import Base: setindex!,
#               Nothing, copyto!,
#              , ^, conj!, promote_rule, similar, fill!,

export Circulant

mutable struct Circulant{T <: Number, S <: Number} <: AbstractMatrix{T}
    vc::Vector{T}
    vc_dft::Vector{S}
    fftplan::Plan
    ifftplan::Plan
end

function Circulant{T}(vc::AbstractVector{T}) where T
    tmp = plan_fft(vc)
    return Circulant(copy(convert(Array{T,1},vc)), tmp*vc, tmp, plan_ifft(vc))
end

Circulant(vc::AbstractVector) = Circulant{ eltype(vc) <: Integer ?
                                            promote_type(eltype(vc), Float64) :
                                            promote_type(eltype(vc), Float32)}(vc)
Circulant{T}(A::AbstractMatrix) where T = Circulant{T}(A[:,1])
Circulant(A::AbstractMatrix) = Circulant(A[:,1])

## From Base
size(A::Circulant) = (size(A, 1), size(A, 2))

size(A::Circulant, _) = length(A.vc)

eltype(A::Circulant) = eltype(A.vc)

function getindex(C::Circulant, i::Integer, j::Integer)
    n = size(C, 1)
    if i > n || j > n
        error("BoundsError()")
    end
    return C.vc[mod(i - j, length(C.vc)) + 1]
end

function copy(A::Circulant)
    #println("We are copying it")
    return Circulant(copy(A.vc))
end

#function (=)(A::Circulant)
# Convert an Circulant matrix to a full matrix
function Matrix(A::Circulant{T}) where T
    m, n = size(A)
    #println("Generating a matrix of dimension: " * string(m) * "x" * string(n))
    Af = Matrix{T}(undef, m, n)
    for j = 1:n
        for i = 1:m
            Af[i,j] = A[i,j]
        end
    end
    return Af
end

function real(A::Circulant)
    return Circulant(real(A.vc))
end

function imag(A::Circulant)
    return Circulant(imag(A.vc))
end

function isreal(A::Circulant)
    if all(imag(A.vc) .!= 0)
        return false
    else
        return true
    end
end

function conj(A::Circulant)
    return Circulant(conj(A.vc))
end

function mul!(Y::Circulant, A::Circulant, B::Circulant)
    if (size(A, 2) != size(B, 1)) || (size(A, 1) != size(Y, 1))
        error("DimensionsMismatch()")
    end
    #println("Here we are 4")
    Y.vc_dft[:] = A.vc_dft.*B.vc_dft
    Y.vc[:] = convert(Vector{eltype(Y)},
              eltype(Y) <: Real ? real(Y.ifftplan*Y.vc_dft) : Y.ifftplan*Y.vc_dft )
end

function (*)(A::Circulant, B::Circulant)
    Y = Circulant(Vector{promote_type(eltype(A), eltype(B))}(undef, size(A,1)))
    mul!(Y, A, B)
    return Y
end

function (+)(A::Circulant, B::Circulant)
    if (size(A, 1) != size(B, 1))
        error("DimensionsMismatch()")
    end
    return Circulant(
            Vector{promote_type(eltype(A), eltype(B))}(A.vc .+ B.vc),
            A.vc_dft .+ B.vc_dft,
            A.fftplan,
            A.ifftplan)
end

function (-)(A::Circulant, B::Circulant)
    if (size(A, 1) != size(B, 1))
        error("DimensionsMismatch()")
    end
    return Circulant(
            Vector{promote_type(eltype(A), eltype(B))}(A.vc .- B.vc),
            A.vc_dft .- B.vc_dft,
            A.fftplan,
            A.ifftplan)
end

function (-)(A::Circulant)
    return Circulant(
            -A.vc,
            -A.vc_dft,
            A.fftplan,
            A.ifftplan)
end

function mul!(Y::AbstractMatrix, A::Circulant, B::AbstractMatrix)
    if (size(A, 2) != size(B, 1)) || (size(A, 1) != size(Y, 1)) ||
        (size(B, 2) != size(Y, 2))
        error("DimensionsMismatch()")
    end

    rettype = eltype(Y)
    if rettype <: Real && ( eltype(A) <: Complex || eltype(B) <:  Complex )
        println("Output matrix is real, while input(s) are complex")
    end
    #println("Here we are 5")

    for i=1:size(B, 2)
        Y[:, i] = convert(Vector{eltype(Y)},
                eltype(Y) <: Real ?
                real(A.ifftplan*(A.vc_dft.*(A.fftplan*B[:, i]))) :
                A.ifftplan*(A.vc_dft.*(A.fftplan*B[:, i])))
    end
end

function mul!(Y::AbstractMatrix, A::Circulant, B::Adjoint{<:Any,<:AbstractMatrix})
    #println("Adjoint mul R")
    mul!(Y, A, copy(B))
end

function mul!(Y::AbstractMatrix, A::Circulant, B::Transpose{<:Any,<:AbstractMatrix})
    #println("Transpose mul R")
    mul!(Y, A, copy(B))
end

function (*)(A::Circulant, B::AbstractMatrix)
    Y = Matrix{promote_type(eltype(A), eltype(B))}(undef, size(A, 1), size(B, 2))
    mul!(Y, A, B)
    return Y
end

function mul!(y::AbstractVector, A::Circulant, x::AbstractVector)
    if length(y) != length(x) || size(A, 2) != length(x)
        throw(DimensionMismatch(""))
    end

    rettype = eltype(y)
    if rettype <: Real && ( eltype(A) <: Complex || eltype(x) <:  Complex )
        println("Output matrix is real, while input(s) are complex")
    end

    #println("Here we are 222")
    y[:] = convert(Vector{eltype(y)},
            eltype(y) <: Real ?
            real(A.ifftplan*(A.vc_dft.*(A.fftplan*x))) :
            A.ifftplan*(A.vc_dft.*(A.fftplan*x)))
end

function (*)(A::Circulant, x::AbstractVector)
    y = Vector{promote_type(eltype(A), eltype(x))}(undef, size(A, 1))
    #println("Here we are, what about 222")
    mul!(y, A, x)
    return y
end

function mul!(Y::AbstractMatrix, A::AbstractMatrix, B::Circulant)
    if (size(A, 2) != size(B, 1)) || (size(A, 1) != size(Y, 1)) ||
        (size(B, 2) != size(Y, 2))
        error("DimensionsMismatch()")
    end

    rettype = eltype(Y)
    if rettype <: Real && ( eltype(A) <: Complex || eltype(B) <:  Complex )
        println("Output matrix is real, while input(s) are complex")
    end
    #println("Here we are 6")

    for i=1:size(Y, 1)
        Y[i, :] = convert(Vector{eltype(Y)},
                eltype(Y) <: Real ?
                real(B.fftplan*(B.vc_dft.*(B.ifftplan*A[i, :]))) :
                B.fftplan*(B.vc_dft.*(B.ifftplan*A[i, :])))
    end
end

function mul!(Y::AbstractMatrix, A::Adjoint{<:Any,<:AbstractMatrix}, B::Circulant)
    #println("Adjoint mul L")
    mul!(Y, copy(A), B)
end

function mul!(Y::AbstractMatrix, A::Transpose{<:Any,<:AbstractMatrix}, B::Circulant)
    #println("Transpose mul 1")
    mul!(Y, copy(A), B)
end

function (*)(A::AbstractMatrix, B::Circulant)
    Y = Matrix{promote_type(eltype(A), eltype(B))}(undef, size(A, 1), size(B, 2))
    mul!(Y, A, B)
    return Y
end

function mul!(Y::Circulant, A::Circulant, b::Number)
    if (size(Y, 1) != size(A, 1)) || (size(Y, 2) != size(A, 2))
        error("DimensionsMismatch()")
    end

    rettype = eltype(Y)
    if rettype <: Real && ( eltype(A) <: Complex)
        error("Output matrix is real, while input(s) are complex")
    end
    #println("Here we are 7")

    Y.vc = A.vc .* b
    Y.vc_dft = A.vc_dft .* b
    Y.fftplan = A.fftplan
    Y.ifftplan = A.ifftplan
end

function mul!(Y::Circulant, b::Number, A::Circulant)
    mul!(Y, A, b)
end

function (*)(A::Circulant, b::Number)
    Y = Circulant(Vector{promote_type(eltype(A), eltype(b))}(undef, size(A,1)))
    mul!(Y, A, b)
    return Y
end

function (/)(A::Circulant, b::Number)
    Y = Circulant(Vector{promote_type(eltype(A), eltype(b))}(undef, size(A,1)))
    mul!(Y, A, 1/b)
    return Y
end

# function .(*)(A::Circulant,b::Number)
#     return A*b
# end

function (*)(b::Number, A::Circulant)
    Y = Circulant(Vector{promote_type(eltype(A), eltype(b))}(undef, size(A,1)))
    mul!(Y, A, b)
    return Y
end

function (/)(b::Number, A::Circulant)
    error("ERROR #84616: Not implemented yet")
end

function lmul!(a::Number, B::Circulant)
    error("ERROR #78913: Not implemented yet")
end

function rmul!(A::Circulant, b::Number)
    error("ERROR #123564: Not implemented yet")
end

function axpy!(a, X::Circulant, Y::Circulant)
    Y.vc[:] = Y.vc[:] + X.vc[:]*a
    Y.vc_dft[:] = Y.vc_dft[:] + X.vc_dft[:]*a
end

#
# function .(*)(b::Number, A::Circulant)
#     return A*b
# end

function ldiv!(y::AbstractVector, A::Circulant, x::AbstractVector)
    @assert length(x) == length(y)
    if T <: Real
        y[:] = real(A.ifftplan*((1 ./ A.vc_dft).*(A.fftplan*x)))
    else
        y[:] = A.ifftplan*((1 ./ A.vc_dft).*(A.fftplan*x))
    end
end

function ldiv!(Y::AbstractMatrix, A::Circulant, X::AbstractMatrix)
    @assert size(Y) == size(X)
    for i=1:size(X, 2)
        if T <: Real
            Y[:,i] = real(A.ifftplan*((1 ./ A.vc_dft).*(A.fftplan*X[:,i])))
        else
            Y[:,i] = A.ifftplan*((1 ./ A.vc_dft).*(A.fftplan*X[:,i]))
        end
    end
end

# function copy(A::Adjoint{Circulant})
#     return(\))
# end

function adjoint!(A::Circulant)
    error("Error #78617")
end

function transpose!(A::Circulant)
    error("Error #123489")
end

function adjoint(A::Circulant)
    n = size(A,1)
    Adj = zeros(eltype(A),n)
    Adj[1] = conj(A.vc[1])

    for i = 2:n
        Adj[i] = conj(A.vc[n-i+2])
    end

    return Circulant(Adj)
end

function transpose(A::Circulant)
    n = size(A,1)
    Trsp = zeros(eltype(A), n)
    Trsp[1] = A.vc[1]

    for i = 2:n
        Trsp[i] = A.vc[n-i+2]
    end

    return Circulant(Trsp)
end

function inv(A::Circulant)
    if eltype(A) <: Real
        return Circulant(convert(Vector{eltype(A)}, real(A.ifftplan*(1 ./ A.vc_dft))), 1 ./ A.vc_dft, A.fftplan, A.ifftplan)
    else
        return Circulant(convert(Vector{eltype(A)}, A.ifftplan*(1 ./ A.vc_dft)), 1 ./ A.vc_dft, A.fftplan, A.ifftplan)
    end
end

function norm(A::Circulant)
    return maximum(abs.(A.vc_dft))
end

function eigvals(A::Circulant, vl=1, vr=size(A,1))
    return A.vc_dft[vl:vr]
end

function eigvecs(A::Circulant, vl=1, vr=size(A,1))
    n = size(A, 1)
    E = zeros(Complex{Float64}, n, vr-vl+1)
    for i = 1:size(A,1)
        for j = vl:vr
            E[i, j-vl+1] = (exp(im*2*π*(j-1)/n))^(i-1)/sqrt(n)
        end
    end
    return E
end

function eigen(A::Circulant, vl=1, vr=size(A,1))
    return (eigvals(A, vl, vr), eigvecs(A, vl, vr))
end

function cond(A::Circulant)
    σ = abs.(eigvals(A))
    return maximum(σ)/minimum(σ)
end

function (\)(A::Circulant, b::AbstractVector)
    return inv(A)*b
end
#
# function ishermitian(A::Circulant)
#     if all((imag(A.vc_dft) ./ real(A.vc_dft)) .< eps()*100)
#         println("Matrix is Hermitian")
#         return true
#     else
#         println("Matrix is not Hermitian")
#         return false
#     end
# end

function ishermitian(A::Circulant)
    diffvec = abs.(A.vc[2:end] - conj(reverse(A.vc[2:end])))
    if all(diffvec .< eps()*10)
        #println("Matrix is Hermitian")
        return true
    else
        #println("Matrix is not Hermitian")
        return false
    end
end

function issymmetric(A::Circulant)
    diffvec = abs.(A.vc[2:end] - reverse(A.vc[2:end]))
    if all(diffvec .< eps()*10)
        #println("Matrix is symmetric")
        return true
    else
        #println("Matrix is not symmetric")
        return false
    end
end

function rank(A::Circulant; rtol=eps())
    myrank = 0
    atol = norm(A)*rtol
    σ  = abs.(A.vc_dft)
    for i=1:size(A,1)
        if σ[i] > atol
            myrank = myrank + 1
        end
    end
    return myrank
end
#function
end


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

export threadid, nthreads, @threads, @spawn,
       threadpool, nthreadpools

"""
    Threads.threadid() -> Int

Get the ID number of the current thread of execution. The master thread has
ID `1`.
"""
threadid() = Int(ccall(:jl_threadid, Int16, ())+1)

"""
    Threads.nthreads([:default|:interactive]) -> Int

Get the number of threads (across all thread pools or within the specified
thread pool) available to Julia. The number of threads across all thread
pools is the inclusive upper bound on [`threadid()`](@ref).

See also: `BLAS.get_num_threads` and `BLAS.set_num_threads` in the
[`LinearAlgebra`](@ref man-linalg) standard library, and `nprocs()` in the
[`Distributed`](@ref man-distributed) standard library.
"""
function nthreads end

nthreads() = Int(unsafe_load(cglobal(:jl_n_threads, Cint)))
function nthreads(pool::Symbol)
    if pool == :default
        tpid = Int8(0)
    elseif pool == :interactive
        tpid = Int8(1)
    else
        error("invalid threadpool specified")
    end
    return _nthreads_in_pool(tpid)
end
function _nthreads_in_pool(tpid::Int8)
    p = unsafe_load(cglobal(:jl_n_threads_per_pool, Ptr{Cint}))
    return Int(unsafe_load(p, tpid + 1))
end

"""
    Threads.threadpool(tid = threadid()) -> Symbol

Returns the specified thread's threadpool; either `:default` or `:interactive`.
"""
function threadpool(tid = threadid())
    tpid = ccall(:jl_threadpoolid, Int8, (Int16,), tid-1)
    return tpid == 0 ? :default : :interactive
end

"""
    Threads.nthreadpools() -> Int

Returns the number of threadpools currently configured.
"""
nthreadpools() = Int(unsafe_load(cglobal(:jl_n_threadpools, Cint)))

function threading_run(fun, static)
    ccall(:jl_enter_threaded_region, Cvoid, ())
    n = nthreads()
    tasks = Vector{Task}(undef, n)
    for i = 1:n
        t = Task(() -> fun(i)) # pass in tid
        t.sticky = static
        static && ccall(:jl_set_task_tid, Cint, (Any, Cint), t, i-1)
        tasks[i] = t
        schedule(t)
    end
    for i = 1:n
        Base._wait(tasks[i])
    end
    ccall(:jl_exit_threaded_region, Cvoid, ())
    failed_tasks = filter(istaskfailed, tasks)
    if !isempty(failed_tasks)
        throw(CompositeException(map(TaskFailedException, failed_tasks)))
    end
end

function _threadsfor(iter, lbody, schedule)
    lidx = iter.args[1]         # index
    range = iter.args[2]
    quote
        local threadsfor_fun
        let range = $(esc(range))
        function threadsfor_fun(tid = 1; onethread = false)
            r = range # Load into local variable
            lenr = length(r)
            # divide loop iterations among threads
            if onethread
                tid = 1
                len, rem = lenr, 0
            else
                len, rem = divrem(lenr, nthreads())
            end
            # not enough iterations for all the threads?
            if len == 0
                if tid > rem
                    return
                end
                len, rem = 1, 0
            end
            # compute this thread's iterations
            f = firstindex(r) + ((tid-1) * len)
            l = f + len - 1
            # distribute remaining iterations evenly
            if rem > 0
                if tid <= rem
                    f = f + (tid-1)
                    l = l + tid
                else
                    f = f + rem
                    l = l + rem
                end
            end
            # run this thread's iterations
            for i = f:l
                local $(esc(lidx)) = @inbounds r[i]
                $(esc(lbody))
            end
        end
        end
        if $(schedule === :dynamic || schedule === :default)
            threading_run(threadsfor_fun, false)
        elseif ccall(:jl_in_threaded_region, Cint, ()) != 0 # :static
            error("`@threads :static` cannot be used concurrently or nested")
        else # :static
            threading_run(threadsfor_fun, true)
        end
        nothing
    end
end

"""
    Threads.@threads [schedule] for ... end

A macro to execute a `for` loop in parallel. The iteration space is distributed to
coarse-grained tasks. This policy can be specified by the `schedule` argument. The
execution of the loop waits for the evaluation of all iterations.

See also: [`@spawn`](@ref Threads.@spawn) and
`pmap` in [`Distributed`](@ref man-distributed).

# Extended help

## Semantics

Unless stronger guarantees are specified by the scheduling option, the loop executed by
`@threads` macro have the following semantics.

The `@threads` macro executes the loop body in an unspecified order and potentially
concurrently. It does not specify the exact assignments of the tasks and the worker threads.
The assignments can be different for each execution. The loop body code (including any code
transitively called from it) must not make any assumptions about the distribution of
iterations to tasks or the worker thread in which they are executed. The loop body for each
iteration must be able to make forward progress independent of other iterations and be free
from data races. As such, invalid synchronizations across iterations may deadlock while
unsynchronized memory accesses may result in undefined behavior.

For example, the above conditions imply that:

- The lock taken in an iteration *must* be released within the same iteration.
- Communicating between iterations using blocking primitives like `Channel`s is incorrect.
- Write only to locations not shared across iterations (unless a lock or atomic operation is
  used).
- The value of [`threadid()`](@ref Threads.threadid) may change even within a single
  iteration.

## Schedulers

Without the scheduler argument, the exact scheduling is unspecified and varies across Julia
releases. Currently, `:dynamic` is used when the scheduler is not specified.

!!! compat "Julia 1.5"
    The `schedule` argument is available as of Julia 1.5.

### `:dynamic` (default)

`:dynamic` scheduler executes iterations dynamically to available worker threads. Current
implementation assumes that the workload for each iteration is uniform. However, this
assumption may be removed in the future.

This scheduling option is merely a hint to the underlying execution mechanism. However, a
few properties can be expected. The number of `Task`s used by `:dynamic` scheduler is
bounded by a small constant multiple of the number of available worker threads
([`nthreads()`](@ref Threads.nthreads)). Each task processes contiguous regions of the
iteration space. Thus, `@threads :dynamic for x in xs; f(x); end` is typically more
efficient than `@sync for x in xs; @spawn f(x); end` if `length(xs)` is significantly
larger than the number of the worker threads and the run-time of `f(x)` is relatively
smaller than the cost of spawning and synchronizing a task (typically less than 10
microseconds).

!!! compat "Julia 1.8"
    The `:dynamic` option for the `schedule` argument is available and the default as of Julia 1.8.

### `:static`

`:static` scheduler creates one task per thread and divides the iterations equally among
them, assigning each task specifically to each thread. In particular, the value of
[`threadid()`](@ref Threads.threadid) is guranteed to be constant within one iteration.
Specifying `:static` is an error if used from inside another `@threads` loop or from a
thread other than 1.

!!! note
    `:static` scheduling exists for supporting transition of code written before Julia 1.3.
    In newly written library functions, `:static` scheduling is discouraged because the
    functions using this option cannot be called from arbitrary worker threads.

## Example

To illustrate of the different scheduling strategies, consider the following function
`busywait` containing a non-yielding timed loop that runs for a given number of seconds.

```julia-repl
julia> function busywait(seconds)
            tstart = time_ns()
            while (time_ns() - tstart) / 1e9 < seconds
            end
        end

julia> @time begin
            Threads.@spawn busywait(5)
            Threads.@threads :static for i in 1:Threads.nthreads()
                busywait(1)
            end
        end
6.003001 seconds (16.33 k allocations: 899.255 KiB, 0.25% compilation time)

julia> @time begin
            Threads.@spawn busywait(5)
            Threads.@threads :dynamic for i in 1:Threads.nthreads()
                busywait(1)
            end
        end
2.012056 seconds (16.05 k allocations: 883.919 KiB, 0.66% compilation time)
```

The `:dynamic` example takes 2 seconds since one of the non-occupied threads is able
to run two of the 1-second iterations to complete the for loop.
"""
macro threads(args...)
    na = length(args)
    if na == 2
        sched, ex = args
        if sched isa QuoteNode
            sched = sched.value
        elseif sched isa Symbol
            # for now only allow quoted symbols
            sched = nothing
        end
        if sched !== :static && sched !== :dynamic
            throw(ArgumentError("unsupported schedule argument in @threads"))
        end
    elseif na == 1
        sched = :default
        ex = args[1]
    else
        throw(ArgumentError("wrong number of arguments in @threads"))
    end
    if !(isa(ex, Expr) && ex.head === :for)
        throw(ArgumentError("@threads requires a `for` loop expression"))
    end
    if !(ex.args[1] isa Expr && ex.args[1].head === :(=))
        throw(ArgumentError("nested outer loops are not currently supported by @threads"))
    end
    return _threadsfor(ex.args[1], ex.args[2], sched)
end

"""
    Threads.@spawn [:default|:interactive] expr

Create a [`Task`](@ref) and [`schedule`](@ref) it to run on any available
thread in the specified threadpool (`:default` if unspecified). The task is
allocated to a thread once one becomes available. To wait for the task to
finish, call [`wait`](@ref) on the result of this macro, or call
[`fetch`](@ref) to wait and then obtain its return value.

Values can be interpolated into `@spawn` via `\$`, which copies the value
directly into the constructed underlying closure. This allows you to insert
the _value_ of a variable, isolating the asynchronous code from changes to
the variable's value in the current task.

!!! note
    See the manual chapter on [multi-threading](@ref man-multithreading)
    for important caveats. See also the chapter on [threadpools](@ref man-threadpools).

!!! compat "Julia 1.3"
    This macro is available as of Julia 1.3.

!!! compat "Julia 1.4"
    Interpolating values via `\$` is available as of Julia 1.4.

!!! compat "Julia 1.9"
    A threadpool may be specified as of Julia 1.9.
"""
macro spawn(args...)
    tpid = Int8(0)
    na = length(args)
    if na == 2
        ttype, ex = args
        if ttype isa QuoteNode
            ttype = ttype.value
        elseif ttype isa Symbol
            # TODO: allow unquoted symbols
            ttype = nothing
        end
        if ttype === :interactive
            tpid = Int8(1)
        elseif ttype !== :default
            throw(ArgumentError("unsupported threadpool in @spawn: $ttype"))
        end
    elseif na == 1
        ex = args[1]
    else
        throw(ArgumentError("wrong number of arguments in @spawn"))
    end

    letargs = Base._lift_one_interp!(ex)

    thunk = esc(:(()->($ex)))
    var = esc(Base.sync_varname)
    quote
        let $(letargs...)
            local task = Task($thunk)
            task.sticky = false
            ccall(:jl_set_task_threadpoolid, Cint, (Any, Int8), task, $tpid)
            if $(Expr(:islocal, var))
                put!($var, task)
            end
            schedule(task)
            task
        end
    end
end

# This is a stub that can be overloaded for downstream structures like `Channel`
function foreach end

# Scheduling traits that can be employed for downstream overloads
abstract type AbstractSchedule end
struct StaticSchedule <: AbstractSchedule end
struct FairSchedule <: AbstractSchedule end

max_stars_repo_name: TechPenguineer/julia

max_issues_repo_name: TechPenguineer/julia

max_forks_repo_name: TechPenguineer/julia


text:
export BlankRobot

mutable struct BlankRobot <: Robot
  btCollisionObject
end
BlankRobot() = BlankRobot(BT.convex_hull([zeros(3)]))


text:
#  Copyright 2016, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at http://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modelling langauge for Julia
# See http://github.com/JuliaOpt/JuMP.jl
#############################################################################
# basic.jl
#
# Solves a simple LP:
# max 5x + 3y
#  st 1x + 5y <= 3
#      0 <= x <= 2
#      0 <= y <= 20
#############################################################################

using JuMP

m = Model()

@variable(m, 0 <= x <= 2)
@variable(m, 0 <= y <= 30)

@objective(m, Max, 5x + 3y)
@constraint(m, 1x + 5y <= 3.0)

print(m)

status = solve(m)

println("Objective value: ", getobjectivevalue(m))
println("x = ", getvalue(x))
println("y = ", getvalue(y))


text:
# The purpose of this script is to validate first order hold derivation (including minimum time)
using ForwardDiff
using LinearAlgebra
using Test
model, obj = TrajectoryOptimization.Dynamics.dubinscar
n = model.n
m = model.m
dt = 1.0
h = sqrt(dt)
solver = Solver(model,obj,integration=:rk3_foh,dt=dt)
results = UnconstrainedVectorResults(model.n,model.m,solver.N)
U = ones(solver.model.m,solver.N)
copyto!(results.U,U)
rollout!(results,solver)
update_jacobians!(results,solver)

# Results from backward pass
k = 2
fcx, fcu = results.fcx[k], results.fcu[k]
fcy, fcv = results.fcx[k+1], results.fcu[k+1]
Q = obj.Q
R = getR(solver)
c = .0075
xf = obj.xf
X = results.X
U = results.U
xm = results.xm[k]
um = (results.U[k] + results.U[k+1])/2.

Lx = dt/6*Q*(X[k] - xf) + 4*dt/6*(I/2 + dt/8*fcx)'*Q*(xm - xf)
Lu = dt/6*R*U[k] + 4*dt/6*((dt/8*fcu)'*Q*(xm - xf) + 0.5*R*um)
Ly = dt/6*Q*(X[k+1] - xf) + 4*dt/6*(I/2 - dt/8*fcy)'*Q*(xm - xf)
Lv = dt/6*R*U[k+1] + 4*dt/6*((-dt/8*fcv)'*Q*(xm - xf) + 0.5*R*um)

Lxx = dt/6.0*Q + 4.0*dt/6.0*(I/2.0 + dt/8.0*fcx)'*Q*(I/2.0 + dt/8.0*fcx)
Luu = dt/6*R + 4*dt/6*((dt/8*fcu)'*Q*(dt/8*fcu) + 0.5*R*0.5)
Lyy = dt/6*Q + 4*dt/6*(I/2 - dt/8*fcy)'*Q*(I/2 - dt/8*fcy)
Lvv = dt/6*R + 4*dt/6*((-dt/8*fcv)'*Q*(-dt/8*fcv) + 0.5*R*0.5)

Lxu = 4*dt/6*(I/2 + dt/8*fcx)'*Q*(dt/8*fcu)
Lxy = 4*dt/6*(I/2 + dt/8*fcx)'*Q*(I/2 - dt/8*fcy)
Lxv = 4*dt/6*(I/2 + dt/8*fcx)'*Q*(-dt/8*fcv)
Luy = 4*dt/6*(dt/8*fcu)'*Q*(I/2 - dt/8*fcy)
Luv_ = 4*dt/6*((dt/8*fcu)'*Q*(-dt/8*fcv) + 0.5*R*0.5)  # note the name change; workspace conflict
Lyv = 4*dt/6*(I/2 - dt/8*fcy)'*Q*(-dt/8*fcv)

# test functions
function el(x,u)
    0.5*(x - xf)'*Q*(x - xf) + 0.5*u'*R*u
end

function el(s)
    x = s[1:n]
    u = s[n+1:n+m]
    el(x,u)
end

function x_midpoint(x,y,xdot,ydot,h)
    0.5*x + 0.5*y + (h^2)/8*xdot - (h^2)/8*ydot
end

function u_midpoint(u,v)
    0.5*u + 0.5*v
end

function fc_(x,u)
    xdot = zero(x)
    model.f(xdot,x,u)
    xdot
end

function fc_(s)
    x = s[1:n]
    u = s[n+1:n+m]
    fc_(x,u)
end

function stage_cost(x,y,u,v,h)
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    (h^2)/6*(el(x,u) + 4*el(xm,um) + el(y,v)) + c*h^2
end

function stage_cost(s)
    x = s[1:n]
    u = s[n+1:n+m]
    y = s[n+m+1:n+m+n]
    v = s[n+m+n+1:n+m+n+m]
    h = s[n+m+n+m+1]
    stage_cost(x,y,u,v,h)
end

x = results.X[k]
y = results.X[k+1]
u = results.U[k]
v = results.U[k+1]

s = [x;u;y;v;h]
stage_cost(x,y,u,v,h)
stage_cost(s)

# Get ForwardDiff gradients and Jacobians
L_gradient = ForwardDiff.gradient(stage_cost,s)

@test isapprox(L_gradient[1:n],Lx)
@test isapprox(L_gradient[n+1:n+m],Lu)
@test isapprox(L_gradient[n+m+1:n+m+n],Ly)
@test isapprox(L_gradient[n+m+n+1:n+m+n+m],Lv)

# NOTE: Because we neglect second-order derivatives of the dynamics (continous or discrete) we can't simply call ForwardDiff's Hessian and get results that match the equations we use in backwardpass_foh, we need to make functions that match the gradients first, then take the jacobians of the gradients without getting second order dynamics derivatives
function Lx_func(x,u,y,v,h)
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    (h^2)/6*Q*(x - xf) + 4*(h^2)/6*(I/2 + (h^2)/8*fcx)'*Q*(xm - xf)
end
function Lx_func(s)
    x = s[1:n]
    u = s[n+1:n+m]
    y = s[n+m+1:n+m+n]
    v = s[n+m+n+1:n+m+n+m]
    h = s[n+m+n+m+1]
    Lx_func(x,u,y,v,h)
end

function Lu_func(x,u,y,v,h)
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    (h^2)/6*R*u + 4*(h^2)/6*(((h^2)/8*fcu)'*Q*(xm - xf) + 0.5*R*um)
end
function Lu_func(s)
    x = s[1:n]
    u = s[n+1:n+m]
    y = s[n+m+1:n+m+n]
    v = s[n+m+n+1:n+m+n+m]
    h = s[n+m+n+m+1]
    Lu_func(x,u,y,v,h)
end

function Ly_func(x,u,y,v,h)
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    Ly = (h^2)/6*Q*(y - xf) + 4*(h^2)/6*(I/2 - (h^2)/8*fcy)'*Q*(xm - xf)
end
function Ly_func(s)
    x = s[1:n]
    u = s[n+1:n+m]
    y = s[n+m+1:n+m+n]
    v = s[n+m+n+1:n+m+n+m]
    h = s[n+m+n+m+1]
    Ly_func(x,u,y,v,h)
end

function Lv_func(x,u,y,v,h)
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    Lv = (h^2)/6*R*v + 4*(h^2)/6*((-(h^2)/8*fcv)'*Q*(xm - xf) + 0.5*R*um)
end
function Lv_func(s)
    x = s[1:n]
    u = s[n+1:n+m]
    y = s[n+m+1:n+m+n]
    v = s[n+m+n+1:n+m+n+m]
    h = s[n+m+n+m+1]
    Lv_func(x,u,y,v,h)
end

# Check that foh Hessians match ForwardDiff
@test isapprox(ForwardDiff.jacobian(Lx_func,s)[1:n,1:n],Lxx)
@test isapprox(ForwardDiff.jacobian(Lu_func,s)[1:m,n+1:n+m],Luu)
@test isapprox(ForwardDiff.jacobian(Ly_func,s)[1:n,n+m+1:n+m+n],Lyy)
@test isapprox(ForwardDiff.jacobian(Lv_func,s)[1:m,n+m+n+1:n+m+n+m],Lvv)

@test isapprox(ForwardDiff.jacobian(Lx_func,s)[1:n,n+1:n+m],Lxu)
@test isapprox(ForwardDiff.jacobian(Lx_func,s)[1:n,n+m+1:n+m+n],Lxy)
@test isapprox(ForwardDiff.jacobian(Lx_func,s)[1:n,n+m+n+1:n+m+n+m],Lxv)
@test isapprox(ForwardDiff.jacobian(Lu_func,s)[1:m,n+m+1:n+m+n],Luy)
@test isapprox(ForwardDiff.jacobian(Lu_func,s)[1:m,n+m+n+1:n+m+n+m],Luv_)
@test isapprox(ForwardDiff.jacobian(Ly_func,s)[1:n,n+m+n+1:n+m+n+m],Lyv)

function L2_func(s)
    x = s[1:n]
    u = s[n+1:n+m]
    y = s[n+m+1:n+m+n]
    v = s[n+m+n+1:n+m+n+m]
    h = s[n+m+n+m+1]
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    4/6*(h^2)*el(xm,um)
end

function L2x_func(s)
    x = s[1:n]
    u = s[n+1:n+m]
    y = s[n+m+1:n+m+n]
    v = s[n+m+n+1:n+m+n+m]
    h = s[n+m+n+m+1]
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    4*(h^2)/6*(I/2 + (h^2)/8*fcx)'*Q*(xm - xf)
end

function L2u_func(s)
    x = s[1:n]
    u = s[n+1:n+m]
    y = s[n+m+1:n+m+n]
    v = s[n+m+n+1:n+m+n+m]
    h = s[n+m+n+m+1]
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    4*(h^2)/6*(((h^2)/8*fcu)'*Q*(xm - xf) + 0.5*R*um)
end

function L2y_func(s)
    x = s[1:n]
    u = s[n+1:n+m]
    y = s[n+m+1:n+m+n]
    v = s[n+m+n+1:n+m+n+m]
    h = s[n+m+n+m+1]
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    4*(h^2)/6*(I/2 - (h^2)/8*fcy)'*Q*(xm - xf)
end

function L2v_func(s)
    x = s[1:n]
    u = s[n+1:n+m]
    y = s[n+m+1:n+m+n]
    v = s[n+m+n+1:n+m+n+m]
    h = s[n+m+n+m+1]
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    4*(h^2)/6*(-((h^2)/8*fcv)'*Q*(xm - xf) + 0.5*R*um)
end

function L2h_func(s)
    x = s[1:n]
    u = s[n+1:n+m]
    y = s[n+m+1:n+m+n]
    v = s[n+m+n+1:n+m+n+m]
    h = s[n+m+n+m+1]
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    xmh = 2/8*h*(fc_(x,u) - fc_(y,v))
    xmu = (h^2)/8*fcu
    xmy = -(h^2)/8*fcy
    xmv = -(h^2)/8*fcv
    ℓ2h = xmh'*Q*(xm-xf)
    return 4/6*(2*h*el(xm,um) + (h^2)*ℓ2h)
end

## L(x,u,y,v) = L(x,u) + L(xm,um) + L(y,v) = L1 + L2 + L3
# Assembling ℓ(x,u) expansion
ℓ1 = el(x,u)
ℓ1x = Q*(x - xf)
ℓ1u = R*u
ℓ1y = zeros(n)
ℓ1v = zeros(m)

ℓ1xx = Q
ℓ1uu = R
ℓ1yy = zeros(n,n)
ℓ1vv = zeros(m,m)

ℓ1xu = zeros(n,m)
ℓ1xy = zeros(n,n)
ℓ1xv = zeros(n,m)
ℓ1uy = zeros(m,n)
ℓ1uv = zeros(m,m)
ℓ1yv = zeros(n,m)

# Assembling ℓ(xm,um) expansion
ℓ2 = el(xm,um)
ℓ2x = (I/2 + dt/8*fcx)'*Q*(xm - xf)
ℓ2u = ((dt/8*fcu)'*Q*(xm - xf) + 0.5*R*um)
ℓ2y = (I/2 - dt/8*fcy)'*Q*(xm - xf)
ℓ2v = ((-dt/8*fcv)'*Q*(xm - xf) + 0.5*R*um)

ℓ2xx = (I/2.0 + dt/8.0*fcx)'*Q*(I/2.0 + dt/8.0*fcx)
ℓ2uu = ((dt/8*fcu)'*Q*(dt/8*fcu) + 0.5*R*0.5)
ℓ2yy = (I/2 - dt/8*fcy)'*Q*(I/2 - dt/8*fcy)
ℓ2vv = ((-dt/8*fcv)'*Q*(-dt/8*fcv) + 0.5*R*0.5)

ℓ2xu = (I/2 + dt/8*fcx)'*Q*(dt/8*fcu)
ℓ2xy = (I/2 + dt/8*fcx)'*Q*(I/2 - dt/8*fcy)
ℓ2xv = (I/2 + dt/8*fcx)'*Q*(-dt/8*fcv)
ℓ2uy = (dt/8*fcu)'*Q*(I/2 - dt/8*fcy)
ℓ2uv = ((dt/8*fcu)'*Q*(-dt/8*fcv) + 0.5*R*0.5)  # note the name change; workspace conflict
ℓ2yv = (I/2 - dt/8*fcy)'*Q*(-dt/8*fcv)

# Confirm that partials of L2 still match ForwardDiff
@test isapprox(ForwardDiff.gradient(L2_func,s)[1:n],4*(h^2)/6*ℓ2x)
@test isapprox(ForwardDiff.gradient(L2_func,s)[n+1:n+m],4*(h^2)/6*ℓ2u)
@test isapprox(ForwardDiff.gradient(L2_func,s)[n+m+1:n+m+n],4*(h^2)/6*ℓ2y)
@test isapprox(ForwardDiff.gradient(L2_func,s)[n+m+n+1:n+m+n+m],4*(h^2)/6*ℓ2v)

@test isapprox(ForwardDiff.jacobian(L2x_func,s)[1:n,1:n],4*(h^2)/6*ℓ2xx)
@test isapprox(ForwardDiff.jacobian(L2u_func,s)[1:m,n+1:n+m],4*(h^2)/6*ℓ2uu)
@test isapprox(ForwardDiff.jacobian(L2y_func,s)[1:n,n+m+1:n+m+n],4*(h^2)/6*ℓ2yy)
@test isapprox(ForwardDiff.jacobian(L2v_func,s)[1:m,n+m+n+1:n+m+n+m],4*(h^2)/6*ℓ2vv)

@test isapprox(ForwardDiff.jacobian(L2x_func,s)[1:n,n+1:n+m],4*(h^2)/6*ℓ2xu)
@test isapprox(ForwardDiff.jacobian(L2x_func,s)[1:n,n+m+1:n+m+n],4*(h^2)/6*ℓ2xy)
@test isapprox(ForwardDiff.jacobian(L2x_func,s)[1:n,n+m+n+1:n+m+n+m],4*(h^2)/6*ℓ2xv)
@test isapprox(ForwardDiff.jacobian(L2u_func,s)[1:m,n+m+1:n+m+n],4*(h^2)/6*ℓ2uy)
@test isapprox(ForwardDiff.jacobian(L2u_func,s)[1:m,n+m+n+1:n+m+n+m],4*(h^2)/6*ℓ2uv)
@test isapprox(ForwardDiff.jacobian(L2y_func,s)[1:n,n+m+n+1:n+m+n+m],4*(h^2)/6*ℓ2yv)

# Assembling ℓ(y,v) expansion
ℓ3 = el(y,v)
ℓ3x = zeros(n)
ℓ3u = zeros(m)
ℓ3y = Q*(y - xf)
ℓ3v = R*v

ℓ3xx = zeros(n,n)
ℓ3uu = zeros(m,m)
ℓ3yy = Q
ℓ3vv = R

ℓ3xu = zeros(n,m)
ℓ3xy = zeros(n,n)
ℓ3xv = zeros(n,m)
ℓ3uy = zeros(m,n)
ℓ3uv = zeros(m,m)
ℓ3yv = zeros(n,m)

## Assemble action-value function Q(x,̂u,y,̂v) expansion where ̂u = [u;h], ̂v = [v;w]
Lx = dt/6*ℓ1x + 4*dt/6*ℓ2x + dt/6*ℓ3x
Lu = dt/6*ℓ1u + 4*dt/6*ℓ2u + dt/6*ℓ3u
Ly = dt/6*ℓ1y + 4*dt/6*ℓ2y + dt/6*ℓ3y
Lv = dt/6*ℓ1v + 4*dt/6*ℓ2v + dt/6*ℓ3v

Lxx = dt/6*ℓ1xx + 4*dt/6*ℓ2xx + dt/6*ℓ3xx
Luu = dt/6*ℓ1uu + 4*dt/6*ℓ2uu + dt/6*ℓ3uu
Lyy = dt/6*ℓ1yy + 4*dt/6*ℓ2yy + dt/6*ℓ3yy
Lvv = dt/6*ℓ1vv + 4*dt/6*ℓ2vv + dt/6*ℓ3vv

Lxu = dt/6*ℓ1xu + 4*dt/6*ℓ2xu + dt/6*ℓ3xu
Lxy = dt/6*ℓ1xy + 4*dt/6*ℓ2xy + dt/6*ℓ3xy
Lxv = dt/6*ℓ1xv + 4*dt/6*ℓ2xv + dt/6*ℓ3xv
Luy = dt/6*ℓ1uy + 4*dt/6*ℓ2uy + dt/6*ℓ3uy
Luv_ = dt/6*ℓ1uv + 4*dt/6*ℓ2uv + dt/6*ℓ3uv
Lyv = dt/6*ℓ1yv + 4*dt/6*ℓ2yv + dt/6*ℓ3yv

dx = results.dx[k]
dy = results.dx[k+1]
dxm = 2*h/8*dx - 2*h/8*dy

# Simplified
xmh = 2/8*h*(dx - dy)
xmu = (h^2)/8*fcu
xmy = 0.5*Matrix(I,n,n) - (h^2)/8*fcy
xmv = -(h^2)/8*fcv
ℓ2h = xmh'*Q*(xm-xf)
L2h = 4/6*(2*h*ℓ2 + (h^2)*ℓ2h)
ℓ2hh = 2/8*((dx - dy)'*Q*(xm - xf) + h*(dx - dy)'*Q*xmh)
L2hh = 4/6*(2*h*ℓ2h + 2*ℓ2 + (h^2)*ℓ2hh + 2*h*ℓ2h)
L2xh = 4/6*(2*h*ℓ2x + (h^2)*(0.5*Matrix(I,n,n) + (h^2)/8*fcx)'*Q*xmh + 2/8*h*fcx'*Q*(xm - xf))
L2uh = 4/6*(2*h*ℓ2u + (h^2)*((h^2)/8*fcu)'*Q*xmh + 2/8*h*fcu'*Q*(xm - xf))
L2hu = 4/6*(2*h*ℓ2u + 2/8*(h^3)*(fcu'*Q*(xm - xf) + xmu'*Q*(dx - dy)))
L2hy = 4/6*(2*h*ℓ2y + 2/8*(h^3)*(-fcy'*Q*(xm - xf) + xmy'*Q*(dx - dy)))
L2hv = 4/6*(2*h*ℓ2v + 2/8*(h^3)*(-fcv'*Q*(xm - xf) + xmv'*Q*(dx - dy)))
# Gradients/Jacobians/Hessians for minimum time

# L2h = 4/6*((h^2)*dxm'*Q*(xm - xf) + 2*h*el(xm,um))
# L2hh = 4/6*(2/8*((h^3)*dxm'*Q*dx + 3*(h^2)*dx'*Q*xm) - 2/8*((h^3)*dxm'*Q*dy + 3*(h^2)*dy'*Q*xm) - 6*(h^2)/8*dx'*Q*xf + 6*(h^2)/8*dy'*Q*xf + 2*(h*dxm'*Q*(xm - xf) + el(xm,um)))
# L2hu = 4*(h^2)/6*(2*(h^3)/8*(fcu'*Q*xm + (h^2)/8*fcu'*Q*dx) - 2*(h^3)/8*fcu'*Q*xf - 2*(h^5)/64*fcu'*Q*dy + 2*h*ℓ2u)

# L2xh = 2/6*Q*(h*x + h*y + (h^3)/2*dx - (h^3)/2*dy) - 4/6*h*Q*xf + 1/12*fcx'*Q*(2*(h^3)*x + 2*(h^3)*y + 6/8*(h^5)*dx -6/8*(h^5)*dy) - (h^3)/3*fcx'*Q*xf
# L2uh = 1/12*fcu'*Q*(2*(h^3)*x + 2*(h^3)*y + 6/8*(h^5)*dx - 6/8*(h^5)*dy) - 1/3*(h^3)*fcu'*Q*xf + 4/6*h*R*um
# L2yh = 2/6*Q*(h*x + h*y + (h^3)/2*dx - (h^3)/2*dy) - 4/6*h*Q*xf - 1/12*fcy'*Q*(2*(h^3)*x + 2*(h^3)*y + 6/8*(h^5)*dx -6/8*(h^5)*dy) + (h^3)/3*fcy'*Q*xf
# L2vh = -1/12*fcv'*Q*(2*(h^3)*x + 2*(h^3)*y + 6/8*(h^5)*dx - 6/8*(h^5)*dy) + 1/3*(h^3)*fcv'*Q*xf + 4/6*h*R*um

# Confirm that new partials match ForwardDiff
@test isapprox(ForwardDiff.gradient(L2_func,s)[end],L2h)
@test isapprox(ForwardDiff.hessian(L2_func,s)[end,end],L2hh)
@test isapprox(ForwardDiff.jacobian(L2x_func,s)[:,end],L2xh)
@test isapprox(ForwardDiff.jacobian(L2u_func,s)[:,end],L2uh)
@test isapprox(ForwardDiff.gradient(L2h_func,s)[n+m+1:n+m+n],L2hy)
@test isapprox(ForwardDiff.gradient(L2h_func,s)[n+m+n+1:n+m+n+m],L2hv)
@test isapprox(L2h_func(s),L2h)

# Augmented expansions for minimum time
_Lx = Lx
_Lu = [Lu;2*h/6*ℓ1 + L2h + 2*h/6*ℓ3 + 2*c*h]
_Ly = Ly
_Lv = [Lv;0]

_Lxx = Lxx
_Luu = [Luu (2/6*h*ℓ1u + L2uh);(2/6*h*ℓ1u + L2hu)' (2/6*ℓ1 + L2hh + 2/6*ℓ3 + 2*c)]
_Lyy = Lyy
_Lvv = [Lvv zeros(m); zeros(m)' 0]

_Lxu = [Lxu (2/6*h*ℓ1x + L2xh)]
_Lxy = Lxy
_Lxv = [Lxv zeros(n)]
_Luy = [Luy; (L2hy + 2*h/6*ℓ3y)']
_Luv = [Luv_ zeros(m); (L2hv + 2*h/6*ℓ3v)' 0]
_Lyv = [Lyv zeros(n)]

###
function stage_cost_alt(s)
    x = s[1:n]
    u = s[n+1:n+m]
    h = s[n+m+1]
    y = s[n+m+1+1:n+m+1+n]
    v = s[n+m+1+n+1:n+m+1+n+m]
    w = s[n+m+1+n+m+1]
    stage_cost(x,y,u,v,h)
end

function Lx_func_alt(s)
    x = s[1:n]
    u = s[n+1:n+m]
    h = s[n+m+1]
    y = s[n+m+1+1:n+m+1+n]
    v = s[n+m+1+n+1:n+m+1+n+m]
    w = s[n+m+1+n+m+1]
    Lx_func(x,u,y,v,h)
end

function Lu_func_alt(s)
    x = s[1:n]
    u = s[n+1:n+m]
    h = s[n+m+1]
    y = s[n+m+1+1:n+m+1+n]
    v = s[n+m+1+n+1:n+m+1+n+m]
    w = s[n+m+1+n+m+1]
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    dxm = 2*h/8*fc_(x,u) - 2*h/8*fc_(y,v)
    L2h = 4/6*((h^2)*dxm'*Q*(xm - xf) + 2*h*el(xm,um))
    [Lu_func(x,u,y,v,h); 2*h/6*el(x,u) + L2h + 2*h/6*el(y,v) + 2*c*h]
end

function Ly_func_alt(s)
    x = s[1:n]
    u = s[n+1:n+m]
    h = s[n+m+1]
    y = s[n+m+1+1:n+m+1+n]
    v = s[n+m+1+n+1:n+m+1+n+m]
    w = s[n+m+1+n+m+1]
    Ly_func(x,u,y,v,h)
end

function Lv_func_alt(s)
    x = s[1:n]
    u = s[n+1:n+m]
    h = s[n+m+1]
    y = s[n+m+1+1:n+m+1+n]
    v = s[n+m+1+n+1:n+m+1+n+m]
    w = s[n+m+1+n+m+1]
    [Lv_func(x,u,y,v,h);0]
end

function Lh_func_alt(s)
    x = s[1:n]
    u = s[n+1:n+m]
    h = s[n+m+1]
    y = s[n+m+1+1:n+m+1+n]
    v = s[n+m+1+n+1:n+m+1+n+m]
    w = s[n+m+1+n+m+1]
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    dxm = 2*h/8*fc_(x,u) - 2*h/8*fc_(y,v)
    L2h = 4/6*((h^2)*dxm'*Q*(xm - xf) + 2*h*el(xm,um))
    2/6*h*el(x,u) + L2h + 2/6*h*el(y,v) + 2*c*h
end

s_alt = [x;u;h;y;v;0]
L_gradient_alt = ForwardDiff.gradient(stage_cost_alt,s_alt)

# Confirm that analytical minimum time expansions match ForwardDiff
@test isapprox(L_gradient_alt[1:n],_Lx)
@test isapprox(L_gradient_alt[n+1:n+m+1],_Lu)
@test isapprox(L_gradient_alt[n+m+1+1:n+m+1+n],_Ly)
@test isapprox(L_gradient_alt[n+m+1+n+1:n+m+1+n+m+1],_Lv)

@test isapprox(ForwardDiff.jacobian(Lx_func_alt,s_alt)[1:n,1:n],_Lxx)
@test isapprox(ForwardDiff.jacobian(Lu_func_alt,s_alt)[:,n+1:n+m+1],_Luu)
@test isapprox(ForwardDiff.jacobian(Ly_func_alt,s_alt)[:,n+m+1+1:n+m+1+n],_Lyy)
@test isapprox(ForwardDiff.jacobian(Lv_func_alt,s_alt)[:,n+m+1+n+1:n+m+1+n+m+1],_Lvv)

@test isapprox(ForwardDiff.jacobian(Lx_func_alt,s_alt)[1:n,n+1:n+m+1],_Lxu)
@test isapprox(ForwardDiff.jacobian(Lx_func_alt,s_alt)[1:n,n+m+1+1:n+m+1+n],_Lxy)
@test isapprox(ForwardDiff.jacobian(Lx_func_alt,s_alt)[1:n,n+m+1+n+1:n+m+1+n+m+1],_Lxv)
@test isapprox(ForwardDiff.jacobian(Lu_func_alt,s_alt)[1:m+1,n+m+1+1:n+m+1+n],_Luy)
@test isapprox(ForwardDiff.jacobian(Lu_func_alt,s_alt)[1:m+1,n+m+1+n+1:n+m+1+n+m+1],_Luv)
@test isapprox(ForwardDiff.jacobian(Ly_func_alt,s_alt)[1:n,n+m+1+n+1:n+m+1+n+m+1],_Lyv)

# Final component check
@test _Lx == (h^2)/6*ℓ1x + 4/6*(h^2)*ℓ2x
@test _Lu == [(h^2)/6*ℓ1u + 4/6*(h^2)*ℓ2u; (2/6*h*ℓ1 + L2h + 2/6*ℓ3 + 2*c*h)]
@test _Ly == 4/6*(h^2)*ℓ2y + (h^2)/6*ℓ3y
@test _Lv == [4/6*(h^2)*ℓ2v + (h^2)/6*ℓ3v; 0]

@test _Lxx == (h^2)/6*ℓ1xx + 4/6*(h^2)*ℓ2xx
@test _Luu == [((h^2)/6*ℓ1uu + 4/6*(h^2)*ℓ2uu) (2/6*h*ℓ1u + L2uh); (2/6*h*ℓ1u + L2hu)' (2/6*ℓ1 + L2hh + 2/6*ℓ3 + 2*c)]
@test _Lyy ==  4/6*(h^2)*ℓ2yy + (h^2)/6*ℓ3yy
@test _Lvv == [(4/6*(h^2)*ℓ2vv + (h^2)/6*ℓ3vv) zeros(m); zeros(m)' 0]

@test _Lxu == [4/6*(h^2)*ℓ2xu (2/6*h*ℓ1x + L2xh)]
@test _Lxy == 4/6*(h^2)*ℓ2xy
@test _Lxv == [4/6*(h^2)*ℓ2xv zeros(n)]
@test _Luy == [4/6*(h^2)*ℓ2uy; (L2hy + 2/6*h*ℓ3y)']
@test _Luv == [4/6*(h^2)*ℓ2uv zeros(m); (L2hv + 2/6*h*ℓ3v)' 0]
@test _Lyv == [4/6*(h^2)*ℓ2yv zeros(n)]

# Confirm that changing dt in discrete dynamics works
x01 = zeros(model.n)
x02 = zeros(model.n)
x1 = zeros(model.n)
x2 = zeros(model.n)
u1 = rand(model.m)
u2 = rand(model.m)

solver.fd(x1,x01,u1,u2,dt) == solver.fd(x2,x02,u1,u2,dt)
solver.fd(x1,x01,u1,u2,dt) != solver.fd(x2,x02,u1,u2,2*dt)

# Infeasible
Ri = .77*Matrix(I,n,n)
ui = rand(n)
vi = rand(n)
w = 0.4
u_bar = [u; h; ui]
v_bar = [v; w; vi]
m_bar = m+1
mm = m_bar + n

si = [x;u_bar;y;v_bar]

function stage_cost_infeasible(x,y,u,v,h,w,ui,vi)
    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    (h^2)/6*(el(x,u) + 4*el(xm,um) + el(y,v)) + c*h^2 + 0.5*ui'*Ri*ui
end

function stage_cost_infeasible(s)
    x_ = s[1:n]
    u_bar_ = s[n+1:n+mm]
    y_ = s[n+mm+1:n+mm+n]
    v_bar_ = s[n+mm+n+1:n+mm+n+mm]

    u_ = u_bar_[1:m]
    h_ = u_bar_[m_bar]
    ui_ = u_bar_[m_bar+1:mm]

    v_ = v_bar_[1:m]
    w_ = v_bar_[m_bar]
    vi_ = v_bar_[m_bar+1:mm]

    @test x_ == x
    @test u_bar_ == u_bar
    @test y_ == y
    @test v_bar_ == v_bar
    @test u_ == u
    @test h_ == h
    @test ui_ == ui
    @test v_ == v
    @test w_ == w
    @test vi_ == vi

    stage_cost_infeasible(x_,y_,u_,v_,h_,w_,ui_,vi_)
end

# Gradients
Gz = ForwardDiff.gradient(stage_cost_infeasible,si)

@test isapprox(Gz[1:n], _Lx)
@test isapprox(Gz[n+1:n+m_bar+n], [_Lu; Ri*ui])
@test isapprox(Gz[n+m_bar+n+1:n+m_bar+n+n], _Ly)
@test isapprox(Gz[n+m_bar+n+n+1:n+m_bar+n+n+m_bar+n],[_Lv; zeros(n)])

# Hessians
function Lz_func_inf(s)
    x = s[1:n]
    u_bar_ = s[n+1:n+mm]
    y = s[n+mm+1:n+mm+n]
    v_bar_ = s[n+mm+n+1:n+mm+n+mm]

    u = u_bar_[1:m]
    h = u_bar_[m_bar]
    ui = u_bar_[m_bar+1:mm]

    v = v_bar_[1:m]
    w = v_bar_[m_bar]
    vi = v_bar_[m_bar+1:mm]

    xm = x_midpoint(x,y,fc_(x,u),fc_(y,v),h)
    um = u_midpoint(u,v)
    dxm = 2*h/8*fc_(x,u) - 2*h/8*fc_(y,v)
    L2h = 4/6*((h^2)*dxm'*Q*(xm - xf) + 2*h*el(xm,um))

    [Lx_func(x,u,y,v,h);
     Lu_func(x,u,y,v,h);
     2*h/6*el(x,u) + L2h + 2*h/6*el(y,v) + 2*c*h;
     Ri*ui;
     Ly_func(x,u,y,v,h);
     Lv_func(x,u,y,v,h);
     0;
     zeros(n)]
end

Lz_func_inf(si)

@test isapprox(Gz,Lz_func_inf(si))

Hz = ForwardDiff.jacobian(Lz_func_inf,si)
@test isapprox(Hz[1:n,1:n],_Lxx)
@test isapprox(Hz[n+1:n+m_bar+n,n+1:n+m_bar+n],[_Luu zeros(m_bar,n); zeros(n,m_bar) Ri])
@test isapprox(Hz[n+m_bar+n+1:n+m_bar+n+n,n+m_bar+n+1:n+m_bar+n+n],_Lyy)
@test isapprox(Hz[n+m_bar+n+n+1:n+m_bar+n+n+m_bar+n,n+m_bar+n+n+1:n+m_bar+n+n+m_bar+n],[_Lvv zeros(m_bar,n); zeros(n,m_bar) zeros(n,n)])

@test isapprox(Hz[1:n,n+1:n+m_bar+n],[_Lxu zeros(n,n)])
@test isapprox(Hz[1:n,n+m_bar+n+1:n+m_bar+n+n],_Lxy)
@test isapprox(Hz[1:n,n+m_bar+n+n+1:n+m_bar+n+n+m_bar+n],[_Lxv zeros(n,n)])
@test isapprox(Hz[n+1:n+m_bar+n,n+m_bar+n+1:n+m_bar+n+n],[_Luy; zeros(n,n)'])
@test isapprox(Hz[n+1:n+m_bar+n,n+m_bar+n+n+1:n+m_bar+n+n+m_bar+n],[_Luv zeros(m_bar,n); zeros(n,m_bar) zeros(n,n)])
@test isapprox(Hz[n+m_bar+n+1:n+m_bar+n+n,n+m_bar+n+n+1:n+m_bar+n+n+m_bar+n],[_Lyv zeros(n,n)])

# Simplified notation
xmh = 2/8*h*(dx - dy)
xmu = (h^2)/8*fcu
xmy = 0.5*Matrix(I,n,n) - (h^2)/8*fcy
xmv = -(h^2)/8*fcv
_ℓ2h = xmh'*Q*(xm-xf)
_ℓ2hh = 2/8*((dx - dy)'*Q*(xm - xf) + h*(dx - dy)'*Q*xmh)
_L2hh = 4/6*(2*h*_ℓ2h + 2*ℓ2 + (h^2)*_ℓ2hh + 2*h*_ℓ2h)
_L2xh = 4/6*(2*h*ℓ2x + (h^2)*(0.5*Matrix(I,n,n) + (h^2)/8*fcx)'*Q*xmh + 2/8*h*fcx'*Q*(xm - xf))
_L2uh = 4/6*(2*h*ℓ2u + (h^2)*((h^2)/8*fcu)'*Q*xmh + 2/8*h*fcu'*Q*(xm - xf))
_L2hu = 4/6*(2*h*ℓ2u + 2/8*(h^3)*(fcu'*Q*(xm - xf) + xmu'*Q*(dx - dy)))
_L2hy = 4/6*(2*h*ℓ2y + 2/8*(h^3)*(-fcy'*Q*(xm - xf) + xmy'*Q*(dx - dy)))
_L2hv = 4/6*(2*h*ℓ2v + 2/8*(h^3)*(-fcv'*Q*(xm - xf) + xmv'*Q*(dx - dy)))

_ℓ2h == ℓ2h
_ℓ2hh == ℓ2hh
_L2hh == L2hh
_L2xh == L2xh
_L2uh == L2uh
_L2hu == L2hu
_L2hy == L2hy
_L2hv == L2hv


text:
# plot gaussian probabilistic distribution
# using probability density function

module GaussProbDistPdf
    using Plots, Distributions
    pyplot()

    function main(is_test=false)
        norm_dist = Normal(209.7, 4.84) # mu, std_dev

        zs = Array(193:229)
        ys = [pdf(norm_dist, z) for z in zs]
        plot(zs, ys, label="PDF")

        if is_test == false
            save_path = joinpath(split(@__FILE__, "src")[1], "src/prob_stats/prob_model/gauss_prob_dist_pdf/gauss_prob_dist_pdf.png")
            savefig(save_path)
        end
    end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
#=
run examples tests from the examples folder
=#

using Test
import DataFrames
include(joinpath(@__DIR__, "../examples/Examples.jl"))
using Main.Examples
import Hypatia
import Hypatia.Solvers

# script verbosity
script_verbose = false

# default options to solvers
default_options = (
    verbose = false,
    # verbose = true,
    default_tol_relax = 10,
    # stepper = Solvers.CombinedStepper{Float64}(),
    # stepper = Solvers.PredOrCentStepper{Float64}(),
    iter_limit = 250,
    )

# instance sets and real types to run and corresponding time limits (seconds)
inst_sets = [
    ("minimal", Float64, 60),
    # ("minimal", Float32, 60),
    # ("minimal", BigFloat, 60),
    # ("fast", Float64, 60),
    # ("various", Float64, 120),
    ]

perf = Examples.setup_benchmark_dataframe()

@testset "examples tests" begin
test_insts = Examples.get_test_instances()

@testset "$mod, $ex" for (mod, mod_insts) in test_insts,
    (ex, (ex_type, ex_insts)) in mod_insts
@testset "$inst_set, $T, $time_limit" for (inst_set, T, time_limit) in inst_sets
    haskey(ex_insts, inst_set) || continue
    inst_subset = ex_insts[inst_set]
    isempty(inst_subset) && continue

    info_perf = (; inst_set, :real_T => T, :example => ex, :model_type => mod)
    new_default_options = (; default_options..., time_limit = time_limit)
    ex_type_T = ex_type{T}

    str = "$mod $ex $T $inst_set"
    println("\nstarting $str tests")
    @testset "$str" begin
        Examples.run_instance_set(inst_subset, ex_type_T, info_perf,
            new_default_options, script_verbose, perf)
    end
end
end

# println("\n")
# DataFrames.show(perf, allrows = true, allcols = true)
# println("\n")
end
;


text:
using Sundials, Base.Test

h1 = Sundials.Handle(Sundials.CVodeCreate(Sundials.CV_BDF, Sundials.CV_NEWTON))
h2 = h1

@test !isempty(h1)
@test !isempty(h2)

empty!(h1)

@test isempty(h1)
@test isempty(h2) # Check that the pointers are correct, and h1 empties h2

empty!(h1) # Make sure this does not throw

h = Sundials.Handle(h1.ptr_ref[]) # Check construction with null pointers
@test isempty(h)

neq = 3
h3 = Sundials.MatrixHandle(Sundials.SUNDenseMatrix(neq,neq),Sundials.DenseMatrix())
h3 = Sundials.MatrixHandle(Sundials.SUNDenseMatrix(neq,neq),Sundials.DenseMatrix())
empty!(h3)
@test isempty(h3)
empty!(h3)
@test isempty(h3)

h3 = Sundials.MatrixHandle(Sundials.SUNBandMatrix(100,3,3,6),Sundials.BandMatrix())
h3 = Sundials.MatrixHandle(Sundials.SUNBandMatrix(100,3,3,6),Sundials.BandMatrix())
empty!(h3)
@test isempty(h3)
empty!(h3)
@test isempty(h3)

h3 = Sundials.MatrixHandle(Sundials.SUNSparseMatrix(neq,neq,neq,Sundials.CSC_MAT),Sundials.SparseMatrix())
h3 = Sundials.MatrixHandle(Sundials.SUNSparseMatrix(neq,neq,neq,Sundials.CSC_MAT),Sundials.SparseMatrix())
empty!(h3)
@test isempty(h3)
empty!(h3)
@test isempty(h3)

A = Sundials.SUNDenseMatrix(neq,neq)
u0 = rand(neq)
Sundials.SUNDenseLinearSolver(u0,A)
h3 = Sundials.LinSolHandle(Sundials.SUNDenseLinearSolver(u0,A),Sundials.Dense())
h3 = Sundials.LinSolHandle(Sundials.SUNDenseLinearSolver(u0,A),Sundials.Dense())
empty!(h3)
@test isempty(h3)
empty!(h3)
@test isempty(h3)


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -12.548880784274985 <= q <= 15.06837792907793)

                     add_NL_constraint(m, :(1/(1 + exp(-(0.4386422028907555 + 0.03049853723546203*1/(1 + exp(-(-0.33806544253006265 + -0.5419050015026476*1/(1 + exp(-(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])))) + -0.12420368558845318*1/(1 + exp(-(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])))) + 0.008939690927983968*1/(1 + exp(-(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))))))) + -0.7963936693468932*1/(1 + exp(-(0.35455126051631725 + -0.3141650319667155*1/(1 + exp(-(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])))) + 0.26311214619974344*1/(1 + exp(-(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])))) + -0.4178451532515153*1/(1 + exp(-(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))))))) + -0.19871150900740808*1/(1 + exp(-(-0.9712553762218641 + -0.9788234903921222*1/(1 + exp(-(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])))) + -0.6598676192535482*1/(1 + exp(-(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])))) + 0.7900118209668339*1/(1 + exp(-(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4])))))))))) + 1/(1 + exp(-(0.039843038969352484 + 0.2885956651060697*1/(1 + exp(-(-0.33806544253006265 + -0.5419050015026476*1/(1 + exp(-(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])))) + -0.12420368558845318*1/(1 + exp(-(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])))) + 0.008939690927983968*1/(1 + exp(-(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))))))) + 0.8477392827604406*1/(1 + exp(-(0.35455126051631725 + -0.3141650319667155*1/(1 + exp(-(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])))) + 0.26311214619974344*1/(1 + exp(-(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])))) + -0.4178451532515153*1/(1 + exp(-(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))))))) + -0.38122333260621444*1/(1 + exp(-(-0.9712553762218641 + -0.9788234903921222*1/(1 + exp(-(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])))) + -0.6598676192535482*1/(1 + exp(-(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])))) + 0.7900118209668339*1/(1 + exp(-(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4])))))))))) + 1/(1 + exp(-(0.5234531332153054 + -0.26529753551452195*1/(1 + exp(-(-0.33806544253006265 + -0.5419050015026476*1/(1 + exp(-(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])))) + -0.12420368558845318*1/(1 + exp(-(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])))) + 0.008939690927983968*1/(1 + exp(-(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))))))) + 0.29155449146294243*1/(1 + exp(-(0.35455126051631725 + -0.3141650319667155*1/(1 + exp(-(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])))) + 0.26311214619974344*1/(1 + exp(-(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])))) + -0.4178451532515153*1/(1 + exp(-(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4]))))))) + -0.357141816330512*1/(1 + exp(-(-0.9712553762218641 + -0.9788234903921222*1/(1 + exp(-(0.3281182697751017 + -0.8099128015858463*$(x[1]) + 0.13317347784439315*$(x[2]) + 0.4405761556905756*$(x[3]) + 0.9454727306320025*$(x[4])))) + -0.6598676192535482*1/(1 + exp(-(-0.11273025364322331 + 0.7707180324977863*$(x[1]) + -0.4030548048256195*$(x[2]) + 0.6515536701116291*$(x[3]) + -0.7265975512301988*$(x[4])))) + 0.7900118209668339*1/(1 + exp(-(-0.1658555453462811 + -0.8928983125267513*$(x[1]) + 0.8779481001349092*$(x[2]) + -0.7642847980322278*$(x[3]) + -0.7283572839362114*$(x[4])))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# implements vanilla value iteration solver

type SerialSolver <: Solver

    maxIterations::Int64

    tolerance::Float64

    gaussSiedel::Bool

    includeV::Bool

    includeQ::Bool

    includeA::Bool

end

function SerialSolver(;maxIterations::Int64=1000, tolerance::Float64=1e-3, gaussSiedel::Bool=true,
                       includeV::Bool=true, includeQ::Bool=true, includeA::Bool=true)
    return SerialSolver(maxIterations, tolerance, gaussSiedel, includeV, includeQ, includeA)
end

function solve(solver::SerialSolver, mdp::DiscreteMDP; verbose::Bool=false)

    gs = solver.gaussSiedel
    u = []
    q = []

    if gs
        u, q = solveGS(solver, mdp, verbose=verbose)
    else
        u, q = solveRegular(solver, mdp, verbose=verbose)
    end

    p = []
    # check policy flag
    if solver.includeA
        p = computePolicy(mdp, u)
    end

    policy = DiscretePolicy(V=u, Q=q, P=p)

    return policy
end

function solveGS(solver::SerialSolver, mdp::DiscreteMDP; verbose::Bool=false)

    nStates  = numStates(mdp)
    nActions = numActions(mdp)

    maxIter = solver.maxIterations
    tol     = solver.tolerance

    valU = zeros(nStates)
    valQ = zeros(nActions, nStates)

    totalTime = 0.0
    iterTime  = 0.0

    for i = 1:maxIter
        tic()
        for si = 1:nStates
            qHi = 0.0

            for ai = 1:nActions
                states, probs = nextStates(mdp, si, ai)
                qNow = reward(mdp, si, ai)

                for sp = 1:length(states)
                    spi = states[sp]
                    qNow += probs[sp] * valU[spi]
                end # sp loop
                valQ[ai, si] = qNow
                if ai == 1 || qNow > qHi
                    qHi = qNow
                    valU[si] = qHi
                end
            end # action loop
        end # state loop

        iterTime = toq();
        totalTime += iterTime
        verbose ? (println("Iteration : $i, iteration run-time: $iterTime, total run-time: $totalTime")) : nothing
    end # main loop

    return valU, valQ

end

function solveRegular(solver::SerialSolver, mdp::DiscreteMDP; verbose::Bool=false)

    nStates  = numStates(mdp)
    nActions = numActions(mdp)

    maxIter = alg.maxIterations
    tol     = alg.tolerance

    valOld = zeros(nStates)
    valNew = zeros(nStates)
    valQ = zeros(nActions, nStates)

    totalTime = 0.0
    iterTime  = 0.0

    residual = 0.0

    for i = 1:maxIter
        tic()
        for si = 1:nStates
            qHi = 0.0

            for ai = 1:nActions
                states, probs = nextStates(mdp, si, ai)
                qNow = reward(mdp, si, ai)

                for sp = 1:length(states)
                    spi = states[sp]
                    qNow += probs[sp] * valOld[spi]
                end # sp loop
                valQ[ai, si] = qNow
                if ai == 1 || qNow > qHi
                    qHi = qNow
                    valNew[si] = qHi
                end
            end # action loop
            newResidual = (valOld[si] - valNew[si])^2
            newResidual > residual? (residual = newResidual) : (nothing)
            valOld[si] = valNew[si]
        end # state loop

        iterTime = toq();
        totalTime += iterTime
        verbose ? (println("Iteration : $i, iteration run-time: $iterTime, total run-time: $totalTime")) : nothing

        # terminate if tolerance value is reached
        if residual < tol; lastIdx = uIdx; break; end
    end # main loop

    return valNew, valQ
end


text:
using SRKGenerator
#Standard Driver Script

if isempty(ARGS)

  ## Called from test. Use a default setup
  imin = -12
  jmax = 1

  ## Parallelize correctly
  host = gethostname()
  if host == "ChrisRackBed"
    #numCores = 4
    #addprocs(numCores-length(workers()))
    ptx_str = joinpath(Pkg.dir("SRKGenerator"),"deps","integrationWin.ptx")
    cudaCores=[1664,640]
  end

  if host == "crackauc2"
    #numCores = 16
    #addprocs(numCores-length(workers()))
    ptx_str = joinpath(Pkg.dir("SRKGenerator"),"deps","integrationLinux.ptx")
    cudaCores = [2816]
  end

  if host == "ChrisRackTV"
    #numCores = 4
    #addprocs(numCores-length(workers()))
    ptx_str = joinpath(Pkg.dir("SRKGenerator"),"deps","integrationWin.ptx")
    cudaCores=[2560]
  end

else

  ## Called from commandline. Assume Linux.
  ptx_str = joinpath(Pkg.dir("SRKGenerator"),"deps","integrationLinux.ptx")
  ## Arg Form:
  cudaCores = [parse(Int,ARGS[1])]
  imin = parse(Int,ARGS[2])
  jmax = parse(Int,ARGS[3]); jmin = -jmax

end

initCon =  [0.13803991182809405,0.5764629771061757,0.4235370228938241,0.45988142444172125,0.8084859191118128,-0.2683673489095204,0.45482884244918836,0.7519704601202507,0.24802953987974916,0.6961685917464052,0.3460790565637937,-0.04224764831019894,0.08513289389320751,1.0272930820405226,0.4565894293733811,1.7266402546369712,-0.46040416640520754,-0.9623742513556953,0.6744100174205603,-0.07377591084107533,-0.4985455217486592,-0.5623853411970345,0.023606441329882558,-0.8953332413687689,-0.14925077940925455,0.7532260416025327,0.6911225838401098,-0.2950978480741201,-0.4550650633914512,0.8347196075716,0.3790963153683881,0.24124914164997122,-0.5014671170389827,0.9198342527636846,-0.2556662881012024,-0.16270085567402398,1.4550650304633204,-0.8347196116875014,-0.37909625532888075,-0.24124916721822443,-0.4998591615011572,0.9168848205568053,-1.4115003346202473,0.994474673400943]
x_L = initCon .- 0.01
x_U = initCon .+ 0.01
resString = srk_optimize(:LN_AUGLAG_EQ,1/100,0,imin,1,-jmax,jmax,
                         rand_minmax = rand(1:3),tol=1e-15,
                         initCon = initCon,x_L=x_L,x_U=x_U,
                         NLoptRandSeed=rand(1:Int(1e8)),gpuEnabled = true,
                         cudaCores=cudaCores,ptx_str=ptx_str)


text:
import Plots, Plots.plot, Plots.plot!, Plots.current

function plot!(tree::Tree{1}; filter::Function = active)
    X = Vector()
    Y = Vector()
    for cell ∈ cells(tree, filter=filter)
        push!(X, cell.position[1])
        push!(Y, cell.state)
    end

    max_level = levels(root(tree)) - 1
    plot(X, Y, legend = false)
    plot!(X, zeros(length(X)), seriestype=:scatter, markersize = 20. /  (1 << max_level))

    return current()
end

function plot(tree::Tree{2}; markers::Bool = false, max_marker_level::Int = 5, path::Bool = false, filter::Function = active)
    plot()
    plot!(tree, markers = markers, max_marker_level = max_marker_level, path = path, filter = filter)
end

function plot!(tree::Tree{2}; markers::Bool = false, max_marker_level::Int = 5, path::Bool = false, filter::Function = active)
    X = Vector()
    Y = Vector()

    for cell ∈ cells(tree, filter=filter)
        append!(X, NaN)
        append!(Y, NaN)

        append!(X, cell.position[1] .+ [-1., 1., 1., -1., -1.] / (2 << (level(cell))))
        append!(Y, cell.position[2] .+ [-1., -1., 1., 1., -1.] / (2 << (level(cell))))
    end
    max_level = levels(root(tree)) - 1
    linewidth0 = 4.
    linewidth = linewidth0 / (1 << max_level)
    plot!(X, Y, legend = false, color=:black, linewidth = linewidth)

    if path
        X = Vector()
        Y = Vector()
        for cell ∈ cells(tree, filter=filter)
            pos = centroid(cell)
            push!(X, pos[1])
            push!(Y, pos[2])
        end
        plot!(X, Y, linewidth = linewidth, color = :lightgray)
    end

    if markers
        markerstrokewidth = linewidth0 /  (1 << min(max_marker_level, max_level))
        markersize = 20. /  (1 << min(max_marker_level, max_level))

        X = Vector()
        Y = Vector()
        for cell ∈ cells(tree, filter=filter, max_level = max_marker_level)
            pos = centroid(cell)
            push!(X, pos[1])
            push!(Y, pos[2])
        end
        plot!(X, Y, seriestype = :scatter, markersize = markersize, markerstrokewidth = markerstrokewidth, marker = :circle)

        X = Vector()
        Y = Vector()
        for face ∈ all_faces(tree, filter = filter)
            if level(face) < max_marker_level
                pos = centroid(face)
                push!(X, pos[1])
                push!(Y, pos[2])
            end
        end
        plot!(X, Y, seriestype = :scatter, markersize = markersize, markerstrokewidth = markerstrokewidth, marker = :square)
    end

    return Plots.current()
end

function plot!(tree::Tree{3}; markers::Bool = false, max_marker_level::Int = 5, path::Bool = false, wireframe=true)
    X = Vector()
    Y = Vector()
    Z = Vector()

    max_level = levels(root(tree)) - 1
    linewidth0 = 4.
    linewidth = linewidth0 / (1 << max_level)
    if wireframe
        xloop = [-1., 1., 1., -1., -1.]
        yloop = [-1., -1., 1., 1., -1.]
        for cell ∈ cells(tree, filter=filter)
            append!(X, NaN)
            append!(Y, NaN)
            append!(Z, NaN)

            append!(X, cell.position[1] .+ xloop / (2 << (level(cell))))
            append!(Y, cell.position[2] .+ yloop / (2 << (level(cell))))
            append!(Z, cell.position[3] .+ -ones(5) / (2 << (level(cell))))
            append!(X, cell.position[1] .+ xloop / (2 << (level(cell))))
            append!(Y, cell.position[2] .+ yloop / (2 << (level(cell))))
            append!(Z, cell.position[3] .+ ones(5) / (2 << (level(cell))))
            for i=2:4
                append!(X, NaN)
                append!(Y, NaN)
                append!(Z, NaN)
                append!(X, cell.position[1] .+ [xloop[i], xloop[i]] / (2 << (level(cell))))
                append!(Y, cell.position[2] .+ [yloop[i], yloop[i]] / (2 << (level(cell))))
                append!(Z, cell.position[3] .+ [-1., 1.] / (2 << (level(cell))))
            end
        end
        plot(X, Y, Z, legend = false, color = :black, linewidth = linewidth)
    end

    if path
        X = Vector()
        Y = Vector()
        Z = Vector()
        for cell ∈ cells(tree, filter=filter)
            pos = centroid(cell)
            push!(X, pos[1])
            push!(Y, pos[2])
            push!(Z, pos[3])
        end
        plot!(X, Y, Z, linewidth = linewidth, color = :lightgray)
    end

    if markers
        markerstrokewidth = linewidth0 /  (1 << min(max_marker_level, max_level))
        markersize = 20. /  (1 << min(max_marker_level, max_level))

        X = Vector()
        Y = Vector()
        Z = Vector()
        for cell ∈ cells(tree, filter=filter)
            if level(cell) < max_marker_level
                pos = centroid(cell)
                push!(X, pos[1])
                push!(Y, pos[2])
                push!(Z, pos[3])
            end
        end
        plot!(X, Y, Z, seriestype = :scatter, markersize = markersize, markerstrokewidth = markerstrokewidth, marker = :circle)

        X = Vector()
        Y = Vector()
        Z = Vector()
        for face ∈ all_faces(tree, filter = filter)
            if level(face) < max_marker_level
                pos = centroid(face)
                push!(X, pos[1])
                push!(Y, pos[2])
                push!(Z, pos[3])
            end
        end
        plot!(X, Y, Z, seriestype = :scatter, markersize = markersize, markerstrokewidth = markerstrokewidth, marker = :square)
    end

    return current()
end


text:
isdefined(Base, :__precompile__) && __precompile__()

module Graphs
using DataStructures
using SparseArrays

import Base: show, ==, <
import Base: length, isempty, size, getindex, isless

export
    AbstractGraph,

    # concept checking
    implements_vertex_list,
    implements_edge_list,
    implements_vertex_map,
    implements_edge_map,
    implements_adjacency_list,
    implements_incidence_list,
    implements_bidirectional_adjacency_list,
    implements_bidirectional_incidence_list,
    implements_adjacency_matrix,
    @graph_implements, @graph_requires,

    # common
    make_vertex, make_edge,
    vertex_type, edge_type, source, target, revedge,
    is_directed, is_mutable, vertex_index, edge_index,
    num_vertices, vertices, num_edges, edges,
    out_degree, out_neighbors, out_edges,
    in_degree, in_neighbors, in_edges,
    attributes,

    KeyVertex, Edge, WeightedEdge, ExVertex, ExEdge, AttributeDict,
    collect_edges, collect_weighted_edges,

    add_edge!, add_vertex!,
    delete_vertex!,

    AbstractEdgePropertyInspector, VectorEdgePropertyInspector,
    ConstantEdgePropertyInspector, AttributeEdgePropertyInspector,
    edge_property, edge_property_requirement,

    # edge_list
    GenericEdgeList, EdgeList, simple_edgelist, edgelist,

    # adjacency_list
    GenericAdjacencyList, SimpleAdjacencyList, AdjacencyList,
    simple_adjlist, adjlist,

    # incidence_list
    GenericIncidenceList, SimpleIncidenceList, VectorIncidenceList, IncidenceList,
    simple_inclist, inclist,

    # dict based graphs
    IncidenceDict, incdict,

    # graph
    GenericGraph, SimpleGraph, Graph, simple_graph, graph,

    # gmatrix
    adjacency_matrix, weight_matrix, distance_matrix, laplacian_matrix,
    adjacency_matrix_sparse, weight_matrix_sparse, laplacian_matrix_sparse,
    sparse2adjacencylist,

    # graph_visit
    AbstractGraphVisitor, TrivialGraphVisitor, LogGraphVisitor,
    discover_vertex!, open_vertex!, close_vertex!,
    examine_neighbor!, examine_edge!,
    visited_vertices, traverse_graph, traverse_graph_withlog,

    # breadth_first_visit
    BreadthFirst, gdistances, gdistances!,

    # depth_first_visit
    DepthFirst, test_cyclic_by_dfs, topological_sort_by_dfs,

    # maximum_adjacency_visit
    MaximumAdjacency, AbstractMASVisitor, min_cut, maximum_adjacency_visit,

    # connected_components
    connected_components, strongly_connected_components,

    # cliques
    maximal_cliques,

    # dijkstra_spath
    DijkstraStates, create_dijkstra_states, AbstractDijkstraVisitor,
    dijkstra_shortest_paths!, dijkstra_shortest_paths,
    dijkstra_shortest_paths_withlog,
    enumerate_paths, enumerate_indices,

    # bellmanford
    BellmanFordStates, create_bellman_ford_states, NegativeCycleError,
    bellman_ford_shortest_paths!, bellman_ford_shortest_paths,
    has_negative_edge_cycle,

    # a_star_spath
    shortest_path,

    # prim_mst
    PrimStates, create_prim_states, AbstractPrimVisitor,
    prim_minimum_spantree!, prim_minimum_spantree, prim_minimum_spantree_withlog,

    # kruskal_mst
    kruskal_select, kruskal_minimum_spantree,

    # floyd_warshall
    floyd_warshall!, floyd_warshall,

    # Graphviz
    to_dot, plot,

    # Random Graph Generation
    erdos_renyi_graph, watts_strogatz_graph,

    # Static Graph Generation
    simple_complete_graph, simple_star_graph, simple_path_graph,
    simple_wheel_graph, simple_diamond_graph, simple_bull_graph,
    simple_chvatal_graph, simple_cubical_graph, simple_desargues_graph,
    simple_dodecahedral_graph, simple_frucht_graph, simple_heawood_graph,
    simple_house_graph, simple_house_x_graph, simple_icosahedral_graph,
    simple_krackhardt_kite_graph, moebius_kantor_graph, simple_octahedral_graph,
    simple_pappus_graph, simple_petersen_graph, simple_sedgewick_maze_graph,
    simple_tetrahedral_graph, simple_truncated_cube_graph,
    simple_truncated_tetrahedron_graph, simple_tutte_graph

## source files

include("concepts.jl")
include("common.jl")

include("edge_list.jl")
include("adjacency_list.jl")
include("incidence_list.jl")
include("graph.jl")
include("show.jl")
include("gmatrix.jl")

include("graph_visit.jl")
include("breadth_first_visit.jl")
include("depth_first_visit.jl")
include("maximum_adjacency_visit.jl")

include("connected_components.jl")
include("dijkstra_spath.jl")
include("bellmanford.jl")
include("a_star_spath.jl")
include("prim_mst.jl")
include("kruskal_mst.jl")
include("floyd_warshall.jl")
include("cliques.jl")

include("dot.jl")

include("random.jl")
include("generators.jl")
end


text:
# Abstract agent definitions
using Agents

abstract type PolicyAgent <: AbstractAgent end

mutable struct SimplePolicyAgent <: PolicyAgent
    id::Int # The identifier number of the agent
    pos::Int # node id in the network
    spos::Tuple{Int,Int} # The x, y location of the agent on a 2D grid
end

# Agent Perception
# SEIIHDR in their local 1kmGrid
# SEIIHDR in the global environment
# SEIIHDR from Network
# Opinion of other Agents
# Media Response

# Agent Opinions
# Vector of opinions
# Media (TV, RADIO, NEWSPAPER, Social Media)

# Description of Agent Behavior
# Opinion shapes behavior
# Mask Wearing
# Social Distancing
# keep 1.5 Distance
# Work/School
#

# Dampening on Virus output
# Function of Behavior.

## MynotesAndré

## EndnotsAndré


text:
## -- Load external packages
    if VERSION>=v"0.7"
        using Statistics
        using DelimitedFiles
        using SpecialFunctions
    end

    using Plots; gr(); default(fmt=:svg);
    using StatGeochem

## --- Read Ronov volume data
    ronov = readdlm("RonovVolumes.csv",'|')
    ronov = elementify(ronov)

    x = vcat(ronov["t_age"]',ronov["b_age"]')[:]
    dt = ronov["b_age"]-ronov["t_age"]
    dt[1] = NaN # Exclude quaternary alluvium

    h = plot(x, repeat(ronov["Ronov_globe"]./dt,inner=2),label="Global")
    plot!(h, x, repeat(ronov["Ronov_NAm"]./dt,inner=2),label="N.Am.")
    plot!(h, x, repeat((ronov["Ronov_globe"] .- ronov["Ronov_NAm"])./dt,inner=2),label="Globe-N.Am.")
    plot!(h, xlims=(0,1600), xlabel="Age (Ma)", ylabel="Volume (km3/yr)", xflip=true, legend=:topleft, fg_color_legend=:white, framestyle=:box)
    savefig(h,"RonovVolume[Globe-NAm].pdf")
    display(h);

## --- Plot each continent by itself
    cont = ["Ronov_Afr", "Ronov_Eur", "Ronov_NAm", "Ronov_SAm", "Ronov_Aus"]
    h = plot(xlims=(0,1600),xlabel="Age (Ma)", ylims=(0,2),ylabel="Volume (km3/yr)")
    for i=1:5
        plot!(h, x, repeat(ronov[cont[i]]./dt,inner=2),label=continents[i], color=continentcolors[i])
    end
    plot!(h, xflip=true, legend=:topleft, fg_color_legend=:white, framestyle=:box)
    savefig(h,"RonovVolumebyContinent.pdf")
    display(h);

## --- Stacked histogram
    cont = ["Ronov_Afr", "Ronov_Eur", "Ronov_NAm", "Ronov_SAm", "Ronov_Aus"]
    vol_cumulative = Array{Float64}(undef, length(x),length(cont))
    total=fill(0.0, size(x))
    for i=5:-1:1
        total .+= repeat(ronov[cont[i]]./dt, inner=2)
        vol_cumulative[:,i] = total
    end

    h = plot()
    for i=1:5
        plot!(h,x,vol_cumulative[:,i],fill=0,label=continents[i],color=continentcolors[i])
    end
    plot!(h, xflip=true, legend=:topleft, fg_color_legend=:white, framestyle=:box)
    plot!(h, xlims=(0,1600), xlabel="Age (Ma)", ylims=(0,3.1), ylabel="Volume (km3/yr)")
    savefig(h,"RonovVolumebyContinentCumulative.pdf")
    display(h);

## --- Plot each continent by itself, normalized by area
    cont = ["Ronov_Afr", "Ronov_Eur", "Ronov_NAm", "Ronov_SAm", "Ronov_Aus"]
    contarea = [30.3E6, 53.4E6, 24.228E6, 18.28E6, 8.8015E6, 1.4E7]
    h = plot(xlims=(0,1600),xlabel="Age (Ma)", ylims=(0,4),ylabel="Volume flux (E-8 km3/yr/km2)")
    for i=1:5
        plot!(h, x, repeat(ronov[cont[i]]./dt/contarea[i]*1E8,inner=2), label=continents[i], color=continentcolors[i])
    end
    plot!(h, xflip=true, legend=:topleft, fg_color_legend=:white, framestyle=:box)
    savefig(h,"RonovVolumeFluxbyContinent.pdf")
    display(h);

## --- Compare Ronov and Macrostrat, for North America and globally

    # Load Macrostrat
    data = readdlm("../macrostrat/macrostrat.csv",',')
    macrostrat = elementify(data)

    # Compare Ronov and Macrostrat for North America
    h = plot(xlims=(0,1600),xlabel="Age (Ma)", ylims=(0,0.5),ylabel="Volume (km3/yr)")
    plot!(h, x, repeat(ronov["Ronov_NAm"]./dt,inner=2),label="Ronov (North America)", color=continentcolors[3])
    plot!(h, macrostrat["Time"], macrostrat["Volume"]/1E6,label="Macrostrat (North America)",color=continentcolors[3])
    plot!(h, xflip=true, legend=:topleft, fg_color_legend=:white, framestyle=:box)

    plot!(h,ylims=(0.006,5),yscale=:log)

    # Compare Ronov and Macrostrat scaled
    plot!(h, x, repeat(ronov["Ronov_globe"]./dt,inner=2),label="Ronov (global)", color=:black)
    plot!(h, macrostrat["Time"], macrostrat["Volume"]/1E6*sum(contarea)/contarea[3],label="Macrostrat scaled", color=:black)
    plot!(h, xflip=true, legend=:topleft, fg_color_legend=:white, framestyle=:box)
    display(h);

    savefig(h,"RonovVsMacrostrat.pdf")
    display(h);

## --- Plot glaciations
    h = plot(xlims=(0,1600),xflip=true,ylims=(0,1))
    plot!(h, [717,717],[0,1], label="")
    plot!(h, [635,635],[0,1], label="")
    plot!(h, [580,580],[0,1], label="")
    savefig(h, "glaciations.pdf")

## ---

max_stars_repo_name: brenhinkeller/GreatUnconformity

max_issues_repo_name: brenhinkeller/GreatUnconformity

max_forks_repo_name: brenhinkeller/GreatUnconformity


text:
# This script steps through the calculation of the gradient of the log-likelihood
# for the multivariate skew t distribution. For debugging purposes.

using SkewDist

n, p = size(X)
k = size(Y,2)
(β, ρ, A, η, ν) = SkewDist.read_params(params, p, k)
# print_params(β, ρ, A, η, ν)
# println("————————————–")
U = SkewDist._U(X,Y,β)
Q = SkewDist._Q(U, A, ρ)
L = SkewDist._L(U,η)
t = SkewDist._t(L, Q, ν, k)

D = diagm(exp(-2*ρ))
Dinv = diagm(exp(2*ρ))
Ωinv = A'D*A
sf = SkewDist._sf(Q,ν,k)
∂logT₁∂t = SkewDist._∂logT₁∂t(t, ν, k)
g_Q = SkewDist._gQ(sf)
t_Q = SkewDist._tQ(L,Q,sf,ν)

# Calculate derivatives
∂ℓ∂β = -2X'diagm(g_Q + ∂logT₁∂t.*t_Q)*U*Ωinv - X'diagm(∂logT₁∂t.*SkewDist._tL(sf))*ones(n)*η'
∂ℓ∂D = eye(k).*(A*U'diagm(g_Q + ∂logT₁∂t.*t_Q)*U*A') + 0.5 * n * Dinv
∂ℓ∂ρ = diag(∂ℓ∂D).*(-2*diag(D))

∂ℓ∂A = 2 * triu(D*A*U'diagm(g_Q + ∂logT₁∂t.*t_Q)*U)

∂ℓ∂η = U'diagm(∂logT₁∂t.*SkewDist._tL(sf))*ones(n)
∂ℓ∂ν = sum(SkewDist._∂log_g∂ν(Q,ν,k) + SkewDist._∂logT₁∂ν(L,Q,ν,k))
∂ℓ∂logν = ∂ℓ∂ν *

#print_params(∂ℓ∂β, ∂ℓ∂ρ, ∂ℓ∂A , ∂ℓ∂η, ∂ℓ∂logν)
ℓ = n * ( log(2) + 0.5 * logdet(D) ) + sum(SkewDist._log_g(Q, ν, k)) + sum(SkewDist._logT₁(t, ν + k))

grad = -2*SkewDist.write_params(∂ℓ∂β, ∂ℓ∂ρ, ∂ℓ∂A, ∂ℓ∂η, ∂ℓ∂logν)
grad


text:
# This file generates types and function overloads for specified functions acting on VectorizedExpressions
# The list of functions is not complete, and is restricted to those that naturally operate
# index-by-index on an an expression.
using Compat

import Base.exp, Base.expm1, Base.log, Base.log10, Base.log1p, Base.sqrt, Base.cbrt, Base.exponent
import Base.significand, Base.sin, Base.sinpi, Base.cos, Base.cospi, Base.tan, Base.sec, Base.cot, Base.csc
import Base.sinh, Base.cosh, Base.tanh, Base.coth, Base.sech, Base.csch
import Base.asin, Base.acos, Base.atan, Base.acot, Base.asec, Base.acsc
import Base.asinh, Base.acosh, Base.atanh, Base.acoth, Base.asech, Base.acsch, Base.sinc, Base.cosc

# single-argument functions
univar_funcs = [:exp, :expm1, :log, :log10, :log1p, :sqrt, :cbrt, :exponent, :significand,
         :sin, :sinpi, :cos, :cospi, :tan, :sec, :cot, :csc,
         :sinh, :cosh, :tanh, :coth, :sech, :csch,
         :asin, :acos, :atan, :acot, :asec, :acsc,
         :asinh, :acosh, :atanh, :acoth, :asech, :acsch, :sinc, :cosc]

# Generate types and methods for all functions
for i = 1:length(univar_funcs)

    op = univar_funcs[i]

    # generate type name from function name
    opname = Symbol("ET_", op)

    eval(quote

         immutable $opname{T,R,ET<:VectorizedExpression}<:VectorizedExpression{T,R}
         VE::ET
         end

         @inline function getindex(A::$opname, i...)
         return $op(A.VE[i...])
         end

         @inline function length(A::$opname)
         return length(A.VE)
         end

         @inline function size(A::$opname)
         return size(A.VE)
         end

        @inline function ($op){T,R}(v::VectorizedExpression{T,R})
            return $opname{T,R, typeof(v)}(v)
        end

        export $op
    end)

end


text:
module MyExample

using ForwardDiff

println("entering MyExample")

include("extra_file.jl")

export my_f, derivative_of_my_f

end