text:
module WishartTest

using Test
using ForneyLab
using ForneyLab: prod!, unsafeMean, unsafeVar, unsafeDetLogMean, outboundType, isApplicable, dims, isProper, logPdf
using ForneyLab: SPWishartOutNPP, VBWishartOut
using SpecialFunctions: digamma

@testset "dims" begin
    @test dims(ProbabilityDistribution(MatrixVariate, Wishart, v=diageye(3), nu=4.0)) == (3, 3)
end

@testset "vague" begin
    @test vague(Wishart, 3) == ProbabilityDistribution(MatrixVariate, Wishart, v=huge*diageye(3), nu=3.0)
    @test vague(Wishart, (3,3)) == ProbabilityDistribution(MatrixVariate, Wishart, v=huge*diageye(3), nu=3.0)
    @test vague(Union{Gamma, Wishart}, (3,3)) == ProbabilityDistribution(MatrixVariate, Wishart, v=huge*diageye(3), nu=3.0)
    @test vague(Union{Gamma, Wishart}, ()) == ProbabilityDistribution(Univariate, Gamma, a=1.0, b=tiny)
end

@testset "isProper" begin
    @test isProper(ProbabilityDistribution(MatrixVariate, Wishart, v=mat(1.0), nu=1.0)) == true
    @test isProper(ProbabilityDistribution(MatrixVariate, Wishart, v=transpose([-1.0]), nu=2.0)) == false
    @test isProper(ProbabilityDistribution(MatrixVariate, Wishart, v=mat(1.0), nu=0.0)) == false
end

@testset "prod!" begin
    @test ProbabilityDistribution(MatrixVariate, Wishart, v=mat(1.0), nu=2.0) * ProbabilityDistribution(MatrixVariate, Wishart, v=mat(1.0),nu=2.0) == ProbabilityDistribution(MatrixVariate, Wishart, v=transpose([0.4999999999999999]),nu=2.0)
    @test ProbabilityDistribution(MatrixVariate, Wishart, v=mat(1.0), nu=2.0) * ProbabilityDistribution(MatrixVariate, PointMass, m=mat(1.0)) == ProbabilityDistribution(MatrixVariate, PointMass, m=mat(1.0))
    @test ProbabilityDistribution(MatrixVariate, PointMass, m=mat(1.0)) * ProbabilityDistribution(MatrixVariate, Wishart, v=mat(1.0), nu=2.0) == ProbabilityDistribution(MatrixVariate, PointMass, m=mat(1.0))
    @test_throws Exception ProbabilityDistribution(MatrixVariate, PointMass, m=transpose([-1.0])) * ProbabilityDistribution(MatrixVariate, Wishart, v=mat(1.0), nu=2.0)
end

@testset "unsafe mean and variance" begin
    @test unsafeMean(ProbabilityDistribution(MatrixVariate, Wishart, v=[2.0 1.0; 1.0 2.0], nu=3.0)) == 3.0*[2.0 1.0; 1.0 2.0]
    @test unsafeVar(ProbabilityDistribution(MatrixVariate, Wishart, v=[2.0 1.0; 1.0 2.0], nu=3.0)) == [24.0 15.0; 15.0 24.0]
    @test unsafeVar(ProbabilityDistribution(MatrixVariate, Wishart, v=diageye(3), nu=3.0)) == [6.0 3.0 3.0; 3.0 6.0 3.0; 3.0 3.0 6.0]
    @test unsafeDetLogMean(ProbabilityDistribution(MatrixVariate, Wishart, v=mat(1.0), nu=1.0)) == digamma(0.5) + log(2)
    @test unsafeDetLogMean(ProbabilityDistribution(MatrixVariate, Wishart, v=eye(2), nu=2.0)) == digamma(0.5) + digamma(1) + 2*log(2)
end

@testset "log pdf" begin
    @test isapprox(logPdf(ProbabilityDistribution(MatrixVariate, Wishart, v=[3.0 1.0; 1.0 1.2], nu=6.0),[2.0 1.0; 1.0 2.0]), -8.15846321016661)
end

#-------------
# Update rules
#-------------

@testset "SPWishartOutNPP" begin
    @test SPWishartOutNPP <: SumProductRule{Wishart}
    @test outboundType(SPWishartOutNPP) == Message{Wishart}
    @test isApplicable(SPWishartOutNPP, [Nothing, Message{PointMass}, Message{PointMass}])

    @test ruleSPWishartOutNPP(nothing, Message(MatrixVariate, PointMass, m=transpose([1.0])), Message(Univariate, PointMass, m=2.0)) == Message(MatrixVariate, Wishart, v=transpose([1.0]), nu=2.0)
end

@testset "VBWishartOut" begin
    @test VBWishartOut <: NaiveVariationalRule{Wishart}
    @test outboundType(VBWishartOut) == Message{Wishart}
    @test isApplicable(VBWishartOut, [Nothing, ProbabilityDistribution, ProbabilityDistribution])
    @test !isApplicable(VBWishartOut, [ProbabilityDistribution, ProbabilityDistribution, Nothing])

    @test ruleVBWishartOut(nothing, ProbabilityDistribution(MatrixVariate, PointMass, m=transpose([1.5])), ProbabilityDistribution(Univariate, PointMass, m=3.0)) == Message(MatrixVariate, Wishart, v=transpose([1.5]), nu=3.0)
end

@testset "averageEnergy and differentialEntropy" begin
    @test differentialEntropy(ProbabilityDistribution(MatrixVariate, Wishart, v=mat(1.0), nu=2.0)) == averageEnergy(Wishart, ProbabilityDistribution(MatrixVariate, Wishart, v=mat(1.0), nu=2.0), ProbabilityDistribution(MatrixVariate, PointMass, m=mat(1.0)), ProbabilityDistribution(Univariate, PointMass, m=2.0))
    @test differentialEntropy(ProbabilityDistribution(MatrixVariate, Wishart, v=mat(1.0), nu=2.0)) == differentialEntropy(ProbabilityDistribution(Univariate, Gamma, a=1.0, b=0.5))
end

end #module


text:
# Autogenerated wrapper script for CUDA_full_jll for x86_64-w64-mingw32
JLLWrappers.@generate_wrapper_header("CUDA_full")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
using Images

img = load("sample-images/cats-3061372_640.jpg");
img_channel_view = channelview(img);
img_channel_view = permuteddimsview(img_channel_view, (2, 3, 1));
# img_channel_view[:, :, 2] *= 0.9;
img_channel_view[img_channel_view .> 0.7] *= 0.9
imshow(img_channel_view)

img_channel_view_1d = view(img_channel_view, 1, :, :)
img_channel_view_1d = img_channel_view_1d * 2

img_float = Float16.(img_channel_view)
img_float[2, :, :] *= 0.9
colorview(RGB, img_float)
imshow(img)

img_adjusted = copy(img)
img_adjusted_channel_view = channelview(img_adjusted)
img_adjusted *= 2

img_channel_view[img_channel_view[:, :, :] .> 0.75N0f8] = img_channel_view[img_channel_view[:, :, :] .> 0.75N0f8] .+ 1N0f8
imshow(img)

img = load("sample-images/cats-3061372_640.jpg");
img_channel_view = channelview(img);
img_channel_view = permuteddimsview(img_channel_view, (2, 3, 1));

x_coords = 320:640

img_channel_view[:, x_coords, 1] = min.(img_channel_view[:, x_coords, 1] .* 1.1, 1);
img_channel_view[:, x_coords, 2] = min.(img_channel_view[:, x_coords, 2] .* 1.2, 1);
img_channel_view[:, x_coords, 3] = min.(img_channel_view[:, x_coords, 3] .* 1.4, 1);
# img_channel_view[img_channel_view .> 0.7] *= 0.9
imshow(img)

img = load("sample-images/cats-3061372_640.jpg");
img_area_range = 320:640

img_area = img[:, img_area_range]
img_area = imfilter(img_area, Kernel.gaussian(5))
img[:, img_area_range] = img_area
imshow(img)

img = load("sample-images/cats-3061372_640.jpg");

height_range = 1:size(img, 1)
width_range = 1:size(img, 2)
noise_color = RGB4{N0f8}(0.,0.,0.)

for i = 1:2000
    img[rand(height_range, 1)[], rand(width_range, 1)[]] = noise_color
end

img = imfilter(img, Kernel.gaussian(3))
img[:, img_area_range] = img_area
imshow(img)

border_size = 50
blur_distance = 10

img = load("sample-images/cats-3061372_640.jpg");
img = padarray(img, Pad(:reflect, border_size, 0))
img = parent(img) # reset the indices

img_area_top = 1:border_size
img[img_area_top, :] = imfilter(img[img_area_top, :], Kernel.gaussian(blur_distance))

img_area_bottom = size(img, 1)-border_size:size(img, 1)
img[img_area_bottom, :] = imfilter(img[img_area_bottom, :], Kernel.gaussian(blur_distance))
imshow(img)

img_large = fill(RGB4{N0f8}(0.,0.,0.), 640, 640)
img_large[141:500, :] = img
img_large[1:140, :] = imfilter(img[1:140, :], Kernel.gaussian(10))
imshow(img_large)

colorview(RGB, Float16.(channelview(img)))
img = padarray(img, Pad(:replicate, 0, 50))
img_area_range = 1:50
img_area = img[:, img_area_range]
img_area = imfilter(img_area, Kernel.gaussian(5))
img[:, img_area_range] = img_area
imshow(img)

# // sharpen image using "unsharp mask" algorithm
# Mat blurred; double sigma = 1, threshold = 5, amount = 1;
# GaussianBlur(img, blurred, Size(), sigma, sigma);
# Mat lowContrastMask = abs(img - blurred) < threshold;
# Mat sharpened = img*(1+amount) + blurred*(-amount);
# img.copyTo(sharpened, lowContrastMask);

# parameters
gaussian_smoothing = 2
intensity = 1

# load an image and apply Gaussian filter
img = load("sample-images/cats-3061372_640.jpg");
imgb = imfilter(img, Kernel.gaussian(gaussian_smoothing));

# convert images to Float to perform mathematical operations
img_array = Float16.(channelview(img));
imgb_array = Float16.(channelview(imgb));

sharpened = img_array .* (1 + intensity) .+ imgb_array .* (-intensity);
sharpened = max.(sharpened, 0);
sharpened = min.(sharpened, 1);

img[:, 1:321] = img[:, 320:640];
img[:, 320:640] = colorview(RGB, sharpened)[:, 320:640];
imshow(img)

using Images, ImageView

pad_top = 25
pad_bottom = 25
pad_color = RGB4{N0f8}(.5,0.5,0.) # color

img = load("sample-images/cats-3061372_640.jpg");
img = padarray(img, Fill(pad_color, (pad_top, pad_bottom)))
img = parent(img)

# load an image and create a grayscale copy
img = load("sample-images/busan-night-scene.jpg");
img_gray = RGB.(Gray.(img));

# get channel representation
img_channel_view = channelview(img);
img_gray_channel_view = channelview(img_gray);

# make channel dimension last
img_arr = permuteddimsview(img_channel_view, (2, 3, 1));
img_grayscale_arr = permuteddimsview(img_gray_channel_view, (2, 3, 1));

# create a mask with all values beiung true
img_mask = fill(true, size(img));

# spot are with colours to retain
img_spot_height = 430:460
img_spot_width = 430:460

for channel_id = 1:3

    # select current channel and crop are of interest
    current_channel = view(img_arr, :, :, channel_id)
    current_channel_area = current_channel[img_spot_height, img_spot_width, :]

    # identify min and max values in a cropped area
    channel_min = minimum(current_channel_area)
    channel_max = maximum(current_channel_area)

    # merge existing mask with a channel specific mask

    channel_mask = channel_min .< current_channel .< channel_max
    img_mask = img_mask .& channel_mask

    #img_mask = img_mask .& (current_channel .> channel_min) .& (current_channel .< channel_max)
end

# apply mask to original and grayscale images
img_masked = img_arr .* img_mask .+ img_grayscale_arr .* (~img_mask);
imshow(colorview(RGB, permutedims(img_masked, (3, 1, 2))))


text:
using Typeclass
import ..Types.Molecule
import ..Abstract

@instance MolecularContainer Molecule begin
    clone(x::Molecule, new_props::Dict=Dict()) = Abstract.clone(x, new_props)
    obtain(x::Molecule, k::Symbol) = Abstract.obtain(x, k)
    features(x::Molecule) = Abstract.features(x)
    ==(x::Molecule, y::Molecule) = Abstract.==(x, y)
end


text:
"""
    provider_standalone_step1(config::ProviderStandaloneConfig; kwargs...)

## Optional Keyword Arguments:
- `scope::AbstractString`. Default value: `"launch"`.
"""
function provider_standalone_step1(config::ProviderStandaloneConfig;
                                   scope::AbstractString = "launch")
    uri = URIs.URI(
        URIs.URI(config.authorize_endpoint);
        query = Dict(
            "client_id" => config.client_id,
            "scope" => scope,
            "response_type" => "code",
            "redirect_uri" => config.redirect_uri,
        )
    )
    return uri
end

"""
    provider_standalone_step2(config::ProviderStandaloneConfig, uri_string::AbstractString)
"""
function provider_standalone_step2(config::ProviderStandaloneConfig,
                                   uri_string::AbstractString)
    uri = URIs.URI(uri_string)
    return provider_standalone_step2(config, uri)
end

"""
    provider_standalone_step2(config::ProviderStandaloneConfig, uri::URIs.URI)
"""
function provider_standalone_step2(config::ProviderStandaloneConfig,
                                   uri::URIs.URI)
    queryparams = URIs.queryparams(uri)
    authorization_code = queryparams["code"]

    authz_code_info = _AuthorizationCodeInformation(;
        authorization_code = authorization_code,
        client_id          = config.client_id,
        redirect_uri       = config.redirect_uri,
        token_endpoint     = config.token_endpoint,
    )

    access_token_info = authorization_code_to_access_token(authz_code_info)

    standalone_result = ProviderStandaloneResult(;
        authorization_code             = access_token_info.authorization_code,
        authorization_code_is_jwt      = access_token_info.authorization_code_is_jwt,
        authorization_code_jwt_decoded = access_token_info.authorization_code_jwt_decoded,
        access_token                   = access_token_info.access_token,
        access_token_is_jwt            = access_token_info.access_token_is_jwt,
        access_token_jwt_decoded       = access_token_info.access_token_jwt_decoded,
        access_token_response          = access_token_info.access_token_response,
    )

    return standalone_result
end


text:
# Automatically-generated headers for cupti
#
# DO NOT EDIT THIS FILE DIRECTLY.

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

@checked function cuptiGetResultString(result, str)
    @runtime_ccall((:cuptiGetResultString, libcupti()), CUptiResult,
                   (CUptiResult, Ptr{Cstring}),
                   result, str)
end
# Julia wrapper for header: cupti_version.h
# Automatically generated using Clang.jl

@checked function cuptiGetVersion(version)
    @runtime_ccall((:cuptiGetVersion, libcupti()), CUptiResult,
                   (Ptr{UInt32},),
                   version)
end
# Julia wrapper for header: cupti_activity.h
# Automatically generated using Clang.jl

@checked function cuptiGetTimestamp(timestamp)
    initialize_api()
    @runtime_ccall((:cuptiGetTimestamp, libcupti()), CUptiResult,
                   (Ptr{UInt64},),
                   timestamp)
end

@checked function cuptiGetContextId(context, contextId)
    initialize_api()
    @runtime_ccall((:cuptiGetContextId, libcupti()), CUptiResult,
                   (CUcontext, Ptr{UInt32}),
                   context, contextId)
end

@checked function cuptiGetStreamId(context, stream, streamId)
    initialize_api()
    @runtime_ccall((:cuptiGetStreamId, libcupti()), CUptiResult,
                   (CUcontext, CUstream, Ptr{UInt32}),
                   context, stream, streamId)
end

@checked function cuptiGetStreamIdEx(context, stream, perThreadStream, streamId)
    initialize_api()
    @runtime_ccall((:cuptiGetStreamIdEx, libcupti()), CUptiResult,
                   (CUcontext, CUstream, UInt8, Ptr{UInt32}),
                   context, stream, perThreadStream, streamId)
end

@checked function cuptiGetDeviceId(context, deviceId)
    initialize_api()
    @runtime_ccall((:cuptiGetDeviceId, libcupti()), CUptiResult,
                   (CUcontext, Ptr{UInt32}),
                   context, deviceId)
end

@checked function cuptiActivityEnable(kind)
    initialize_api()
    @runtime_ccall((:cuptiActivityEnable, libcupti()), CUptiResult,
                   (CUpti_ActivityKind,),
                   kind)
end

@checked function cuptiActivityDisable(kind)
    initialize_api()
    @runtime_ccall((:cuptiActivityDisable, libcupti()), CUptiResult,
                   (CUpti_ActivityKind,),
                   kind)
end

@checked function cuptiActivityEnableContext(context, kind)
    initialize_api()
    @runtime_ccall((:cuptiActivityEnableContext, libcupti()), CUptiResult,
                   (CUcontext, CUpti_ActivityKind),
                   context, kind)
end

@checked function cuptiActivityDisableContext(context, kind)
    initialize_api()
    @runtime_ccall((:cuptiActivityDisableContext, libcupti()), CUptiResult,
                   (CUcontext, CUpti_ActivityKind),
                   context, kind)
end

@checked function cuptiActivityGetNumDroppedRecords(context, streamId, dropped)
    initialize_api()
    @runtime_ccall((:cuptiActivityGetNumDroppedRecords, libcupti()), CUptiResult,
                   (CUcontext, UInt32, Ptr{Csize_t}),
                   context, streamId, dropped)
end

@checked function cuptiActivityGetNextRecord(buffer, validBufferSizeBytes, record)
    initialize_api()
    @runtime_ccall((:cuptiActivityGetNextRecord, libcupti()), CUptiResult,
                   (Ptr{UInt8}, Csize_t, Ptr{Ptr{CUpti_Activity}}),
                   buffer, validBufferSizeBytes, record)
end

@checked function cuptiActivityRegisterCallbacks(funcBufferRequested, funcBufferCompleted)
    initialize_api()
    @runtime_ccall((:cuptiActivityRegisterCallbacks, libcupti()), CUptiResult,
                   (CUpti_BuffersCallbackRequestFunc, CUpti_BuffersCallbackCompleteFunc),
                   funcBufferRequested, funcBufferCompleted)
end

@checked function cuptiActivityFlush(context, streamId, flag)
    initialize_api()
    @runtime_ccall((:cuptiActivityFlush, libcupti()), CUptiResult,
                   (CUcontext, UInt32, UInt32),
                   context, streamId, flag)
end

@checked function cuptiActivityFlushAll(flag)
    initialize_api()
    @runtime_ccall((:cuptiActivityFlushAll, libcupti()), CUptiResult,
                   (UInt32,),
                   flag)
end

@checked function cuptiActivityGetAttribute(attr, valueSize, value)
    initialize_api()
    @runtime_ccall((:cuptiActivityGetAttribute, libcupti()), CUptiResult,
                   (CUpti_ActivityAttribute, Ptr{Csize_t}, Ptr{Cvoid}),
                   attr, valueSize, value)
end

@checked function cuptiActivitySetAttribute(attr, valueSize, value)
    initialize_api()
    @runtime_ccall((:cuptiActivitySetAttribute, libcupti()), CUptiResult,
                   (CUpti_ActivityAttribute, Ptr{Csize_t}, Ptr{Cvoid}),
                   attr, valueSize, value)
end

@checked function cuptiActivityConfigureUnifiedMemoryCounter(config, count)
    initialize_api()
    @runtime_ccall((:cuptiActivityConfigureUnifiedMemoryCounter, libcupti()), CUptiResult,
                   (Ptr{CUpti_ActivityUnifiedMemoryCounterConfig}, UInt32),
                   config, count)
end

@checked function cuptiGetAutoBoostState(context, state)
    initialize_api()
    @runtime_ccall((:cuptiGetAutoBoostState, libcupti()), CUptiResult,
                   (CUcontext, Ptr{CUpti_ActivityAutoBoostState}),
                   context, state)
end

@checked function cuptiActivityConfigurePCSampling(ctx, config)
    initialize_api()
    @runtime_ccall((:cuptiActivityConfigurePCSampling, libcupti()), CUptiResult,
                   (CUcontext, Ptr{CUpti_ActivityPCSamplingConfig}),
                   ctx, config)
end

@checked function cuptiGetLastError()
    initialize_api()
    @runtime_ccall((:cuptiGetLastError, libcupti()), CUptiResult, ())
end

@checked function cuptiSetThreadIdType(type)
    initialize_api()
    @runtime_ccall((:cuptiSetThreadIdType, libcupti()), CUptiResult,
                   (CUpti_ActivityThreadIdType,),
                   type)
end

@checked function cuptiGetThreadIdType(type)
    initialize_api()
    @runtime_ccall((:cuptiGetThreadIdType, libcupti()), CUptiResult,
                   (Ptr{CUpti_ActivityThreadIdType},),
                   type)
end

@checked function cuptiComputeCapabilitySupported(major, minor, support)
    initialize_api()
    @runtime_ccall((:cuptiComputeCapabilitySupported, libcupti()), CUptiResult,
                   (Cint, Cint, Ptr{Cint}),
                   major, minor, support)
end

@checked function cuptiDeviceSupported(dev, support)
    initialize_api()
    @runtime_ccall((:cuptiDeviceSupported, libcupti()), CUptiResult,
                   (CUdevice, Ptr{Cint}),
                   dev, support)
end

@checked function cuptiDeviceVirtualizationMode(dev, mode)
    initialize_api()
    @runtime_ccall((:cuptiDeviceVirtualizationMode, libcupti()), CUptiResult,
                   (CUdevice, Ptr{CUpti_DeviceVirtualizationMode}),
                   dev, mode)
end

@checked function cuptiFinalize()
    initialize_api()
    @runtime_ccall((:cuptiFinalize, libcupti()), CUptiResult, ())
end

@checked function cuptiActivityPushExternalCorrelationId(kind, id)
    initialize_api()
    @runtime_ccall((:cuptiActivityPushExternalCorrelationId, libcupti()), CUptiResult,
                   (CUpti_ExternalCorrelationKind, UInt64),
                   kind, id)
end

@checked function cuptiActivityPopExternalCorrelationId(kind, lastId)
    initialize_api()
    @runtime_ccall((:cuptiActivityPopExternalCorrelationId, libcupti()), CUptiResult,
                   (CUpti_ExternalCorrelationKind, Ptr{UInt64}),
                   kind, lastId)
end

@checked function cuptiActivityEnableLatencyTimestamps(enable)
    initialize_api()
    @runtime_ccall((:cuptiActivityEnableLatencyTimestamps, libcupti()), CUptiResult,
                   (UInt8,),
                   enable)
end
# Julia wrapper for header: cupti_callbacks.h
# Automatically generated using Clang.jl

@checked function cuptiSupportedDomains(domainCount, domainTable)
    initialize_api()
    @runtime_ccall((:cuptiSupportedDomains, libcupti()), CUptiResult,
                   (Ptr{Csize_t}, Ptr{CUpti_DomainTable}),
                   domainCount, domainTable)
end

@checked function cuptiSubscribe(subscriber, callback, userdata)
    initialize_api()
    @runtime_ccall((:cuptiSubscribe, libcupti()), CUptiResult,
                   (Ptr{CUpti_SubscriberHandle}, CUpti_CallbackFunc, Ptr{Cvoid}),
                   subscriber, callback, userdata)
end

@checked function cuptiUnsubscribe(subscriber)
    initialize_api()
    @runtime_ccall((:cuptiUnsubscribe, libcupti()), CUptiResult,
                   (CUpti_SubscriberHandle,),
                   subscriber)
end

@checked function cuptiGetCallbackState(enable, subscriber, domain, cbid)
    initialize_api()
    @runtime_ccall((:cuptiGetCallbackState, libcupti()), CUptiResult,
                   (Ptr{UInt32}, CUpti_SubscriberHandle, CUpti_CallbackDomain,
                    CUpti_CallbackId),
                   enable, subscriber, domain, cbid)
end

@checked function cuptiEnableCallback(enable, subscriber, domain, cbid)
    initialize_api()
    @runtime_ccall((:cuptiEnableCallback, libcupti()), CUptiResult,
                   (UInt32, CUpti_SubscriberHandle, CUpti_CallbackDomain, CUpti_CallbackId),
                   enable, subscriber, domain, cbid)
end

@checked function cuptiEnableDomain(enable, subscriber, domain)
    initialize_api()
    @runtime_ccall((:cuptiEnableDomain, libcupti()), CUptiResult,
                   (UInt32, CUpti_SubscriberHandle, CUpti_CallbackDomain),
                   enable, subscriber, domain)
end

@checked function cuptiEnableAllDomains(enable, subscriber)
    initialize_api()
    @runtime_ccall((:cuptiEnableAllDomains, libcupti()), CUptiResult,
                   (UInt32, CUpti_SubscriberHandle),
                   enable, subscriber)
end

@checked function cuptiGetCallbackName(domain, cbid, name)
    initialize_api()
    @runtime_ccall((:cuptiGetCallbackName, libcupti()), CUptiResult,
                   (CUpti_CallbackDomain, UInt32, Ptr{Cstring}),
                   domain, cbid, name)
end
# Julia wrapper for header: cupti_events.h
# Automatically generated using Clang.jl

@checked function cuptiSetEventCollectionMode(context, mode)
    initialize_api()
    @runtime_ccall((:cuptiSetEventCollectionMode, libcupti()), CUptiResult,
                   (CUcontext, CUpti_EventCollectionMode),
                   context, mode)
end

@checked function cuptiDeviceGetAttribute(device, attrib, valueSize, value)
    initialize_api()
    @runtime_ccall((:cuptiDeviceGetAttribute, libcupti()), CUptiResult,
                   (CUdevice, CUpti_DeviceAttribute, Ptr{Csize_t}, Ptr{Cvoid}),
                   device, attrib, valueSize, value)
end

@checked function cuptiDeviceGetTimestamp(context, timestamp)
    initialize_api()
    @runtime_ccall((:cuptiDeviceGetTimestamp, libcupti()), CUptiResult,
                   (CUcontext, Ptr{UInt64}),
                   context, timestamp)
end

@checked function cuptiDeviceGetNumEventDomains(device, numDomains)
    initialize_api()
    @runtime_ccall((:cuptiDeviceGetNumEventDomains, libcupti()), CUptiResult,
                   (CUdevice, Ptr{UInt32}),
                   device, numDomains)
end

@checked function cuptiDeviceEnumEventDomains(device, arraySizeBytes, domainArray)
    initialize_api()
    @runtime_ccall((:cuptiDeviceEnumEventDomains, libcupti()), CUptiResult,
                   (CUdevice, Ptr{Csize_t}, Ptr{CUpti_EventDomainID}),
                   device, arraySizeBytes, domainArray)
end

@checked function cuptiDeviceGetEventDomainAttribute(device, eventDomain, attrib,
                                                     valueSize, value)
    initialize_api()
    @runtime_ccall((:cuptiDeviceGetEventDomainAttribute, libcupti()), CUptiResult,
                   (CUdevice, CUpti_EventDomainID, CUpti_EventDomainAttribute,
                    Ptr{Csize_t}, Ptr{Cvoid}),
                   device, eventDomain, attrib, valueSize, value)
end

@checked function cuptiGetNumEventDomains(numDomains)
    initialize_api()
    @runtime_ccall((:cuptiGetNumEventDomains, libcupti()), CUptiResult,
                   (Ptr{UInt32},),
                   numDomains)
end

@checked function cuptiEnumEventDomains(arraySizeBytes, domainArray)
    initialize_api()
    @runtime_ccall((:cuptiEnumEventDomains, libcupti()), CUptiResult,
                   (Ptr{Csize_t}, Ptr{CUpti_EventDomainID}),
                   arraySizeBytes, domainArray)
end

@checked function cuptiEventDomainGetAttribute(eventDomain, attrib, valueSize, value)
    initialize_api()
    @runtime_ccall((:cuptiEventDomainGetAttribute, libcupti()), CUptiResult,
                   (CUpti_EventDomainID, CUpti_EventDomainAttribute, Ptr{Csize_t},
                    Ptr{Cvoid}),
                   eventDomain, attrib, valueSize, value)
end

@checked function cuptiEventDomainGetNumEvents(eventDomain, numEvents)
    initialize_api()
    @runtime_ccall((:cuptiEventDomainGetNumEvents, libcupti()), CUptiResult,
                   (CUpti_EventDomainID, Ptr{UInt32}),
                   eventDomain, numEvents)
end

@checked function cuptiEventDomainEnumEvents(eventDomain, arraySizeBytes, eventArray)
    initialize_api()
    @runtime_ccall((:cuptiEventDomainEnumEvents, libcupti()), CUptiResult,
                   (CUpti_EventDomainID, Ptr{Csize_t}, Ptr{CUpti_EventID}),
                   eventDomain, arraySizeBytes, eventArray)
end

@checked function cuptiEventGetAttribute(event, attrib, valueSize, value)
    initialize_api()
    @runtime_ccall((:cuptiEventGetAttribute, libcupti()), CUptiResult,
                   (CUpti_EventID, CUpti_EventAttribute, Ptr{Csize_t}, Ptr{Cvoid}),
                   event, attrib, valueSize, value)
end

@checked function cuptiEventGetIdFromName(device, eventName, event)
    initialize_api()
    @runtime_ccall((:cuptiEventGetIdFromName, libcupti()), CUptiResult,
                   (CUdevice, Cstring, Ptr{CUpti_EventID}),
                   device, eventName, event)
end

@checked function cuptiEventGroupCreate(context, eventGroup, flags)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupCreate, libcupti()), CUptiResult,
                   (CUcontext, Ptr{CUpti_EventGroup}, UInt32),
                   context, eventGroup, flags)
end

@checked function cuptiEventGroupDestroy(eventGroup)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupDestroy, libcupti()), CUptiResult,
                   (CUpti_EventGroup,),
                   eventGroup)
end

@checked function cuptiEventGroupGetAttribute(eventGroup, attrib, valueSize, value)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupGetAttribute, libcupti()), CUptiResult,
                   (CUpti_EventGroup, CUpti_EventGroupAttribute, Ptr{Csize_t}, Ptr{Cvoid}),
                   eventGroup, attrib, valueSize, value)
end

@checked function cuptiEventGroupSetAttribute(eventGroup, attrib, valueSize, value)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupSetAttribute, libcupti()), CUptiResult,
                   (CUpti_EventGroup, CUpti_EventGroupAttribute, Csize_t, Ptr{Cvoid}),
                   eventGroup, attrib, valueSize, value)
end

@checked function cuptiEventGroupAddEvent(eventGroup, event)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupAddEvent, libcupti()), CUptiResult,
                   (CUpti_EventGroup, CUpti_EventID),
                   eventGroup, event)
end

@checked function cuptiEventGroupRemoveEvent(eventGroup, event)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupRemoveEvent, libcupti()), CUptiResult,
                   (CUpti_EventGroup, CUpti_EventID),
                   eventGroup, event)
end

@checked function cuptiEventGroupRemoveAllEvents(eventGroup)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupRemoveAllEvents, libcupti()), CUptiResult,
                   (CUpti_EventGroup,),
                   eventGroup)
end

@checked function cuptiEventGroupResetAllEvents(eventGroup)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupResetAllEvents, libcupti()), CUptiResult,
                   (CUpti_EventGroup,),
                   eventGroup)
end

@checked function cuptiEventGroupEnable(eventGroup)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupEnable, libcupti()), CUptiResult,
                   (CUpti_EventGroup,),
                   eventGroup)
end

@checked function cuptiEventGroupDisable(eventGroup)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupDisable, libcupti()), CUptiResult,
                   (CUpti_EventGroup,),
                   eventGroup)
end

@checked function cuptiEventGroupReadEvent(eventGroup, flags, event,
                                           eventValueBufferSizeBytes, eventValueBuffer)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupReadEvent, libcupti()), CUptiResult,
                   (CUpti_EventGroup, CUpti_ReadEventFlags, CUpti_EventID, Ptr{Csize_t},
                    Ptr{UInt64}),
                   eventGroup, flags, event, eventValueBufferSizeBytes, eventValueBuffer)
end

@checked function cuptiEventGroupReadAllEvents(eventGroup, flags,
                                               eventValueBufferSizeBytes, eventValueBuffer,
                                               eventIdArraySizeBytes, eventIdArray,
                                               numEventIdsRead)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupReadAllEvents, libcupti()), CUptiResult,
                   (CUpti_EventGroup, CUpti_ReadEventFlags, Ptr{Csize_t}, Ptr{UInt64},
                    Ptr{Csize_t}, Ptr{CUpti_EventID}, Ptr{Csize_t}),
                   eventGroup, flags, eventValueBufferSizeBytes, eventValueBuffer,
                   eventIdArraySizeBytes, eventIdArray, numEventIdsRead)
end

@checked function cuptiEventGroupSetsCreate(context, eventIdArraySizeBytes, eventIdArray,
                                            eventGroupPasses)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupSetsCreate, libcupti()), CUptiResult,
                   (CUcontext, Csize_t, Ptr{CUpti_EventID}, Ptr{Ptr{CUpti_EventGroupSets}}),
                   context, eventIdArraySizeBytes, eventIdArray, eventGroupPasses)
end

@checked function cuptiEventGroupSetsDestroy(eventGroupSets)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupSetsDestroy, libcupti()), CUptiResult,
                   (Ptr{CUpti_EventGroupSets},),
                   eventGroupSets)
end

@checked function cuptiEventGroupSetEnable(eventGroupSet)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupSetEnable, libcupti()), CUptiResult,
                   (Ptr{CUpti_EventGroupSet},),
                   eventGroupSet)
end

@checked function cuptiEventGroupSetDisable(eventGroupSet)
    initialize_api()
    @runtime_ccall((:cuptiEventGroupSetDisable, libcupti()), CUptiResult,
                   (Ptr{CUpti_EventGroupSet},),
                   eventGroupSet)
end

@checked function cuptiEnableKernelReplayMode(context)
    initialize_api()
    @runtime_ccall((:cuptiEnableKernelReplayMode, libcupti()), CUptiResult,
                   (CUcontext,),
                   context)
end

@checked function cuptiDisableKernelReplayMode(context)
    initialize_api()
    @runtime_ccall((:cuptiDisableKernelReplayMode, libcupti()), CUptiResult,
                   (CUcontext,),
                   context)
end

@checked function cuptiKernelReplaySubscribeUpdate(updateFunc, customData)
    initialize_api()
    @runtime_ccall((:cuptiKernelReplaySubscribeUpdate, libcupti()), CUptiResult,
                   (CUpti_KernelReplayUpdateFunc, Ptr{Cvoid}),
                   updateFunc, customData)
end
# Julia wrapper for header: cupti_metrics.h
# Automatically generated using Clang.jl

@checked function cuptiGetNumMetrics(numMetrics)
    initialize_api()
    @runtime_ccall((:cuptiGetNumMetrics, libcupti()), CUptiResult,
                   (Ptr{UInt32},),
                   numMetrics)
end

@checked function cuptiEnumMetrics(arraySizeBytes, metricArray)
    initialize_api()
    @runtime_ccall((:cuptiEnumMetrics, libcupti()), CUptiResult,
                   (Ptr{Csize_t}, Ptr{CUpti_MetricID}),
                   arraySizeBytes, metricArray)
end

@checked function cuptiDeviceGetNumMetrics(device, numMetrics)
    initialize_api()
    @runtime_ccall((:cuptiDeviceGetNumMetrics, libcupti()), CUptiResult,
                   (CUdevice, Ptr{UInt32}),
                   device, numMetrics)
end

@checked function cuptiDeviceEnumMetrics(device, arraySizeBytes, metricArray)
    initialize_api()
    @runtime_ccall((:cuptiDeviceEnumMetrics, libcupti()), CUptiResult,
                   (CUdevice, Ptr{Csize_t}, Ptr{CUpti_MetricID}),
                   device, arraySizeBytes, metricArray)
end

@checked function cuptiMetricGetAttribute(metric, attrib, valueSize, value)
    initialize_api()
    @runtime_ccall((:cuptiMetricGetAttribute, libcupti()), CUptiResult,
                   (CUpti_MetricID, CUpti_MetricAttribute, Ptr{Csize_t}, Ptr{Cvoid}),
                   metric, attrib, valueSize, value)
end

@checked function cuptiMetricGetIdFromName(device, metricName, metric)
    initialize_api()
    @runtime_ccall((:cuptiMetricGetIdFromName, libcupti()), CUptiResult,
                   (CUdevice, Cstring, Ptr{CUpti_MetricID}),
                   device, metricName, metric)
end

@checked function cuptiMetricGetNumEvents(metric, numEvents)
    initialize_api()
    @runtime_ccall((:cuptiMetricGetNumEvents, libcupti()), CUptiResult,
                   (CUpti_MetricID, Ptr{UInt32}),
                   metric, numEvents)
end

@checked function cuptiMetricEnumEvents(metric, eventIdArraySizeBytes, eventIdArray)
    initialize_api()
    @runtime_ccall((:cuptiMetricEnumEvents, libcupti()), CUptiResult,
                   (CUpti_MetricID, Ptr{Csize_t}, Ptr{Cint}),
                   metric, eventIdArraySizeBytes, eventIdArray)
end

@checked function cuptiMetricGetNumProperties(metric, numProp)
    initialize_api()
    @runtime_ccall((:cuptiMetricGetNumProperties, libcupti()), CUptiResult,
                   (CUpti_MetricID, Ptr{UInt32}),
                   metric, numProp)
end

@checked function cuptiMetricEnumProperties(metric, propIdArraySizeBytes, propIdArray)
    initialize_api()
    @runtime_ccall((:cuptiMetricEnumProperties, libcupti()), CUptiResult,
                   (CUpti_MetricID, Ptr{Csize_t}, Ptr{CUpti_MetricPropertyID}),
                   metric, propIdArraySizeBytes, propIdArray)
end

@checked function cuptiMetricGetRequiredEventGroupSets(context, metric, eventGroupSets)
    initialize_api()
    @runtime_ccall((:cuptiMetricGetRequiredEventGroupSets, libcupti()), CUptiResult,
                   (CUcontext, CUpti_MetricID, Ptr{Ptr{Cint}}),
                   context, metric, eventGroupSets)
end

@checked function cuptiMetricCreateEventGroupSets(context, metricIdArraySizeBytes,
                                                  metricIdArray, eventGroupPasses)
    initialize_api()
    @runtime_ccall((:cuptiMetricCreateEventGroupSets, libcupti()), CUptiResult,
                   (CUcontext, Csize_t, Ptr{CUpti_MetricID}, Ptr{Ptr{Cint}}),
                   context, metricIdArraySizeBytes, metricIdArray, eventGroupPasses)
end

@checked function cuptiMetricGetValue(device, metric, eventIdArraySizeBytes, eventIdArray,
                                      eventValueArraySizeBytes, eventValueArray,
                                      timeDuration, metricValue)
    initialize_api()
    @runtime_ccall((:cuptiMetricGetValue, libcupti()), CUptiResult,
                   (CUdevice, CUpti_MetricID, Csize_t, Ptr{Cint}, Csize_t, Ptr{UInt64},
                    UInt64, Ptr{CUpti_MetricValue}),
                   device, metric, eventIdArraySizeBytes, eventIdArray,
                   eventValueArraySizeBytes, eventValueArray, timeDuration, metricValue)
end

@checked function cuptiMetricGetValue2(metric, eventIdArraySizeBytes, eventIdArray,
                                       eventValueArraySizeBytes, eventValueArray,
                                       propIdArraySizeBytes, propIdArray,
                                       propValueArraySizeBytes, propValueArray, metricValue)
    initialize_api()
    @runtime_ccall((:cuptiMetricGetValue2, libcupti()), CUptiResult,
                   (CUpti_MetricID, Csize_t, Ptr{Cint}, Csize_t, Ptr{UInt64}, Csize_t,
                    Ptr{CUpti_MetricPropertyID}, Csize_t, Ptr{UInt64},
                    Ptr{CUpti_MetricValue}),
                   metric, eventIdArraySizeBytes, eventIdArray, eventValueArraySizeBytes,
                   eventValueArray, propIdArraySizeBytes, propIdArray,
                   propValueArraySizeBytes, propValueArray, metricValue)
end
# Julia wrapper for header: cupti_profiler_target.h
# Automatically generated using Clang.jl

@checked function cuptiProfilerInitialize(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerInitialize, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_Initialize_Params},),
                   pParams)
end

@checked function cuptiProfilerDeInitialize(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerDeInitialize, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_DeInitialize_Params},),
                   pParams)
end

@checked function cuptiProfilerCounterDataImageCalculateSize(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerCounterDataImageCalculateSize, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_CounterDataImage_CalculateSize_Params},),
                   pParams)
end

@checked function cuptiProfilerCounterDataImageInitialize(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerCounterDataImageInitialize, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_CounterDataImage_Initialize_Params},),
                   pParams)
end

@checked function cuptiProfilerCounterDataImageCalculateScratchBufferSize(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerCounterDataImageCalculateScratchBufferSize, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_CounterDataImage_CalculateScratchBufferSize_Params},
                    ),
                   pParams)
end

@checked function cuptiProfilerCounterDataImageInitializeScratchBuffer(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerCounterDataImageInitializeScratchBuffer, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_CounterDataImage_InitializeScratchBuffer_Params},),
                   pParams)
end

@checked function cuptiProfilerBeginSession(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerBeginSession, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_BeginSession_Params},),
                   pParams)
end

@checked function cuptiProfilerEndSession(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerEndSession, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_EndSession_Params},),
                   pParams)
end

@checked function cuptiProfilerSetConfig(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerSetConfig, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_SetConfig_Params},),
                   pParams)
end

@checked function cuptiProfilerUnsetConfig(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerUnsetConfig, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_UnsetConfig_Params},),
                   pParams)
end

@checked function cuptiProfilerBeginPass(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerBeginPass, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_BeginPass_Params},),
                   pParams)
end

@checked function cuptiProfilerEndPass(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerEndPass, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_EndPass_Params},),
                   pParams)
end

@checked function cuptiProfilerEnableProfiling(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerEnableProfiling, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_EnableProfiling_Params},),
                   pParams)
end

@checked function cuptiProfilerDisableProfiling(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerDisableProfiling, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_DisableProfiling_Params},),
                   pParams)
end

@checked function cuptiProfilerIsPassCollected(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerIsPassCollected, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_IsPassCollected_Params},),
                   pParams)
end

@checked function cuptiProfilerFlushCounterData(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerFlushCounterData, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_FlushCounterData_Params},),
                   pParams)
end

@checked function cuptiProfilerPushRange(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerPushRange, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_PushRange_Params},),
                   pParams)
end

@checked function cuptiProfilerPopRange(pParams)
    initialize_api()
    @runtime_ccall((:cuptiProfilerPopRange, libcupti()), CUptiResult,
                   (Ptr{CUpti_Profiler_PopRange_Params},),
                   pParams)
end


text:
"""
    alias_setup(probs)

Computes alias probabilities.
"""
alias_setup(probs::AbstractVector{<:Real}) = alias_setup(sparse(probs))

function alias_setup(probs::SparseVector{<:Real})
    K = length(probs)
    J = spzeros(Int, K)
    q = probs * K

    smaller = Int[] # prob idxs < 1/K
    larger = Int[]  # prob idxs >= 1/k

    for i in 1:length(probs)
        if q[i] < 1.0  # equivalent to prob < 1/K but saves the division
            push!(smaller, i)
        else
            push!(larger, i)
        end
    end

    while length(smaller) > 0 && length(larger) > 0
        small = pop!(smaller)
        large = pop!(larger)
        J[small] = large
        q[large] = q[large] + q[small] - 1.0
        if q[large] < 1.0
            push!(smaller, large)
        else
            push!(larger, large)
        end
    end

    return J, q
end

"""
    alias_sample(J, q)

Alias Sampling first described in [1]. [2] might be a helpful resource to understand alias sampling.

[1] A. Kronmal and A. V. Peterson. On the alias method for generating random variables from a
    discrete distribution. The American Statistician, 33(4):214-218, 1979.
[2] https://lips.cs.princeton.edu/the-alias-method-efficient-sampling-with-many-discrete-outcomes/
"""
function alias_sample(J::AbstractVector{<:Integer}, q::AbstractVector{<:Real})
    small_index = ceil(Int, rand() * length(J))
    if rand() < q[small_index]
        return small_index
    else
        return J[small_index]
    end
end


text:
# Tests related to ModelConstructors and making sure
# its code related to parameters work w/DSGE.jl

# subspecs
function sstest(m::AnSchorfheide)

    # Change all the fields of an unfixed parameter
    m <= parameter(:ι_w, 0.000, (0.0, .9999), (0.0,0.9999), ModelConstructors.Untransformed(), Normal(0.0,1.0), fixed=false,
                   description="ι_w: A new parameter.",
                   tex_label="\\iota_w")

    # Change an unfixed parameter to be fixed
    m <= parameter(:ι_p, 0.000, fixed=true,
                   description= "ι_p: Another new parameter",
                   tex_label="\\iota_p")

    # Change a fixed parameter
    m <= parameter(:δ, 0.02,  fixed=true,
                   description="δ: The capital depreciation rate.", tex_label="\\delta" )

    # Overwrite a fixed parameter with an unfixed parameter
    m <= parameter(:ϵ_p, 0.750, (1e-5, 10.),   (1e-5, 10.),     ModelConstructors.Exponential(),    GammaAlt(0.75, 0.4),
                   fixed=false,  scaling = x -> 1 + x/100,
                   description="ϵ_p: No description available.",
                   tex_label="\\varepsilon_{p}")

    steadystate!(m)
end

# Test update! with regime switching
m = Model1002("ss10")
m <= Setting(:regime_switching, true)
m <= Setting(:n_regimes, 2)
m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m),
                                            2 => Date(2012, 3, 31)))

ModelConstructors.set_regime_val!(m[:σ_g], 1, 2.523)
ModelConstructors.set_regime_val!(m[:σ_g], 2, 0.; override_bounds = true)
ModelConstructors.set_regime_val!(m[:σ_b], 1, 0.0292)
ModelConstructors.set_regime_val!(m[:σ_b], 2, 0.0292)
ModelConstructors.set_regime_val!(m[:α], 1, 0.1596)
ModelConstructors.set_regime_val!(m[:α], 2, 0.0292)

draw = vec(rand(m.parameters, 1, regime_switching = true))
ModelConstructors.update(m.parameters, draw)
σ_g_ind = findfirst(x->x.key==:σ_g, m.parameters)
σ_b_ind = findfirst(x->x.key==:σ_b, m.parameters)
α_ind = findfirst(x->x.key==:α, m.parameters)
@testset "Test regime switching with two regimes" begin
    @test draw[α_ind] == m.parameters[α_ind].regimes[:value][1]
    @test draw[σ_g_ind] ==  m.parameters[σ_g_ind].regimes[:value][1]
    @test draw[σ_b_ind] ==  m.parameters[σ_b_ind].regimes[:value][1]
    @test draw[end-2] ==  m.parameters[α_ind].regimes[:value][2]
    @test draw[end-1] ==  m.parameters[σ_g_ind].regimes[:value][2]
    @test draw[end] ==  m.parameters[σ_b_ind].regimes[:value][2]
end

m = Model1002("ss10")
m <= Setting(:regime_switching, true)
m <= Setting(:n_regimes, 2)
m <= Setting(:regime_dates, Dict{Int, Date}(1 => date_presample_start(m),
                                            2 => Date(2012, 3, 31),
                                            3 => Date(2015, 3, 31)))

ModelConstructors.set_regime_val!(m[:σ_g], 1, 2.523)
ModelConstructors.set_regime_val!(m[:σ_g], 2, 0.; override_bounds = true)
ModelConstructors.set_regime_val!(m[:σ_g], 3, 2.523)

draw = vec(rand(m.parameters, 1, regime_switching = true))
ModelConstructors.update(m.parameters, draw)

@testset "Test regime switching with two regimes" begin
    @test draw[σ_g_ind] ==  m.parameters[σ_g_ind].regimes[:value][1]
    @test draw[end-1] ==  m.parameters[σ_g_ind].regimes[:value][2]
    @test draw[end] ==  m.parameters[σ_g_ind].regimes[:value][3]
end

m = AnSchorfheide()
sstest(m)

@testset "Test steady-state parameters" begin
    @test m[:ι_w].value == 0.0
    @test m[:ι_w].valuebounds == (0.0, .9999)
    @test m[:ι_w].transform == ModelConstructors.Untransformed()
    @test m[:ι_w].transform_parameterization == (0.0,0.9999)
    @test isa(m[:ι_w].prior.value, Normal)

    @test m[:ι_p].value == 0.0
    @test m[:ι_p].valuebounds == (0.0, 0.0)
    @test isnull(m[:ι_p].prior)
    @test m[:ι_p].fixed == true

    @test m[:δ].value == 0.02
    @test m[:δ].valuebounds == (0.02, 0.02)
    @test isnull(m[:δ].prior)
    @test m[:δ].fixed == true

    @test m[:ϵ_p].value == 0.750
    @test m[:ϵ_p].transform == ModelConstructors.Exponential()
    @test isa(m[:ϵ_p].prior.value, Gamma)
    @test m[:ϵ_p].fixed==false
end

m = AnSchorfheide()
let lastparam = parameter(:p, 0.0)
    for θ in m.parameters
        isa(θ, Parameter) && (lastparam = θ)
    end
    @testset "Check AnSchorfheide last parameter" begin
        @test isa(lastparam, Parameter)
        @test lastparam.value == 0.20*2.237937
    end
end

# transform_to_real_line and transform_to_model_space, acting on the entire parameter vector. they should be inverses!
pvec = m.parameters
vals = transform_to_real_line(pvec)
transform_to_model_space!(m, vals)
@testset "Check parameter transformations for optimization part 2" begin
    @test pvec == m.parameters
end

# all fixed parameters should be unchanged by both transform_to_real_line and transform_to_model_space
@testset "Check fixed parameters are unchanged by optimization transformations" begin
    for θ in m.parameters
        if θ.fixed
            @test θ.value == transform_to_real_line(θ)
            @test θ.value == transform_to_model_space(θ, θ.value)
        end
    end
end


text:
function padarray(img::AbstractArray, prepad::Union{Vector{Int},Dims}, postpad::Union{Vector{Int},Dims}, border::AbstractString, value)
    if border == "value"
        depwarn("string-valued borders are deprecated, use `padarray(img, Fill(value, prepad, postpad))` instead, where the padding entries are Dims-tuples", :padarray)
        return padarray(img, Fill(value, (prepad...,), (postpad...,)))
    end
    padarray(img, prepad, postpad, border)
end

function padarray(img::AbstractArray, prepad::Union{Vector{Int},Dims}, postpad::Union{Vector{Int},Dims}, border::AbstractString)
    if border ∈ ["replicate", "circular", "reflect", "symmetric"]
        depwarn("string-valued borders are deprecated, use `padarray(img, Pad(:$border, prepad, postpad)` instead, where the padding entries are Dims-tuples", :padarray)
        return padarray(img, Pad(Symbol(border), (prepad...,), (postpad...,)))
    elseif border == "inner"
        depwarn("string-valued borders are deprecated, use `padarray(img, Inner)` instead", :padarray)
        return padarray(img, Inner())
    else
        throw(ArgumentError("$border not a recognized border"))
    end
end

padarray(img::AbstractArray, padding::Union{Vector{Int},Dims}, border::AbstractString = "replicate") = padarray(img, padding, padding, border)
padarray{T<:Number}(img::AbstractArray{T}, padding::Union{Vector{Int},Dims}, value::T) = padarray(img, padding, padding, "value", value)

function padarray{T,n}(img::AbstractArray{T,n}, padding::Union{Vector{Int},Dims}, border::AbstractString, direction::AbstractString)
    if direction == "both"
        return padarray(img, padding, padding, border)
    elseif direction == "pre"
        return padarray(img, padding, zeros(Int, n), border)
    elseif direction == "post"
        return padarray(img, zeros(Int, n), padding, border)
    end
end

function padarray{T<:Number,n}(img::AbstractArray{T,n}, padding::Vector{Int}, value::T, direction::AbstractString)
    if direction == "both"
        return padarray(img, padding, padding, "value", value)
    elseif direction == "pre"
        return padarray(img, padding, zeros(Int, n), "value", value)
    elseif direction == "post"
        return padarray(img, zeros(Int, n), padding, "value", value)
    end
end

function imfilter(img::AbstractArray, kern, border::AbstractString, value)
    if border == "value"
        depwarn("string-valued borders are deprecated, use `imfilter(img, kern, Fill(value))` instead", :imfilter)
        return imfilter(img, kern, Fill(value))
    end
    imfilter(img, kern, border)
end

export imfilter_fft
function imfilter_fft(img, kern, border::AbstractString, value)
    if border == "value"
        depwarn("string-valued 'fill' borders are deprecated, use `imfilter(img, kern, Fill(value), Algorithm.FFT())` instead", :imfilter_fft)
        return imfilter(img, kern, Fill(value), Algorithm.FFT())
    elseif border ∈ ["replicate", "circular", "reflect", "symmetric"]
        return imfilter(img, kern, border, Algorithm.FFT())
    elseif border == "inner"
        depwarn("specifying \"inner\" as a string is deprecated, use `imfilter(img, kern, Inner(), Algorithm.FFT())` instead", :imfilter_fft)
        return imfilter(img, kern, Inner(), Algorithm.FFT())
    else
        throw(ArgumentError("$border not a recognized border"))
    end
end

imfilter_fft(img, filter) = imfilter_fft(img, filter, "replicate", 0)
imfilter_fft(img, filter, border) = imfilter_fft(img, filter, border, 0)

export imfilter_gaussian
function imfilter_gaussian(img, sigma; emit_warning=true, astype=nothing)
    if astype != nothing
        depwarn("imfilter_gaussian(img, sigma; astype=$astype, kwargs...) is deprecated; use `imfilter($astype, img, IIRGaussian(sigma; kwargs...))` instead, possibly with `NA()`", :imfilter_gaussian)
        factkernel = KernelFactors.IIRGaussian(astype, sigma; emit_warning=emit_warning)
        return imfilter(astype, img, factkernel, NA())
    end
    depwarn("imfilter_gaussian(img, sigma; kwargs...) is deprecated; use `imfilter(img, IIRGaussian(sigma; kwargs...))` instead, possibly with `NA()`", :imfilter_gaussian)
    factkernel = KernelFactors.IIRGaussian(sigma; emit_warning=emit_warning)
    imfilter(_eltype(Float64, eltype(img)), img, factkernel, NA())
end

_eltype{T,C<:Colorant}(::Type{T}, ::Type{C}) = base_colorant_type(C){T}
_eltype{T,R<:Real}(::Type{T}, ::Type{R}) = T


text:
#=
makemask.jl
methods relating to forming mask/support array
=#

export makemask, circle, ellipse
#using ImageGeoms: ImageGeom, Mask*

makemask(ig::ImageGeom, ::MaskAll) = trues(ig)
makemask(ig::ImageGeom, ::MaskCircle) = circle(ig)
makemask(ig::ImageGeom, ::MaskEllipse) = ellipse(ig)
makemask(ig::ImageGeom, ::MaskAllButEdge) = all_but_edge(ig)

"""
    ellipse(ig::ImageGeom{2} ; ...)
Ellipse that just inscribes the rectangle,
but keeping a 1 pixel border due to some regularizer limitations.
"""
function ellipse(ig::ImageGeom{2} ;
    dx::RealU = ig.deltas[1],
    dy::RealU = ig.deltas[2],
    rx::RealU = abs((ig.dims[1]/2-1)*dx),
    ry::RealU = abs((ig.dims[2]/2-1)*dy),
    cx::RealU = zero(dx),
    cy::RealU = zero(dy),
#   over::Int=2,
)

    x,y = axes(ig)
    fun = (x,y) -> ((x - cx) / rx)^2 + ((y - cy) / ry)^2 < 1
    return fun.(x, y')
end

function ellipse(ig::ImageGeom{3} ; kwargs...)
    i2 = ImageGeom(ig.dims[1:2], ig.deltas[1:2], ig.offsets[1:2])
    return repeat(ellipse(i2 ; kwargs...), 1, 1, ig.dims[3])
end

# default is a circle that just inscribes the square
# but keeping a 1 pixel border due to ASPIRE regularization restriction
function circle(ig::ImageGeom{2} ;
    rx::RealU = abs((ig.dims[1]/2-1)*ig.deltas[1]),
    ry::RealU = abs((ig.dims[2]/2-1)*ig.deltas[2]),
    r::RealU = min(rx,ry),
    kwargs...,
#   over::Int=2,
)
    return ellipse(ig ; rx=r, ry=r, kwargs...)
end

function circle(ig::ImageGeom{3} ; kwargs...)
    i2 = ImageGeom(ig.dims[1:2], ig.deltas[1:2], ig.offsets[1:2])
    return repeat(circle(i2 ; kwargs...), 1, 1, ig.dims[3])
end

function all_but_edge(ig::ImageGeom{2})
    mask = trues(ig.dims)
    mask[  1,   :] .= false
    mask[end,   :] .= false
    mask[  :,   1] .= false
    mask[  :, end] .= false
    return mask
end

function all_but_edge(ig::ImageGeom{3})
    i2 = ImageGeom( dims=ig.dims[1:2] )
    return repeat(all_but_edge(i2), 1, 1, ig.dims[3])
end


text:
export Line, Wire, StickNetwork, WireProp, NWN_circuit, len

using StaticArrays, LinearAlgebra, LightGraphs, SimpleWeightedGraphs

"""
```julia
Line(p₁,p₂)
```
Represents a line segment with end points `p₁` and `p₂` given
by `SVector{N,T}`. Parameterized by the datatype of the
coordinates, `T`, and the length (i.e. dimension) `N`.
"""
struct Line{T<:Real,N}
    p₁::SVector{N,T}
    p₂::SVector{N,T}
end

"""
```julia
WireProp(ρ,D)
```
Contains the resistivity, `ρ`, and diameter, `D`, data for a nanowire.
"""
struct WireProp{T<:Number}
    ρ::T # Resistivity
    D::T # Wire diameter
end

"""
```julia
StickNetwork(lines,dims)
```
Stores a network of sticks, `lines`, bounded in the region `[0,0]` to `dims`.
"""
struct StickNetwork{T<:Real,N}
    lines::Vector{Line{T,N}}
    dims::SVector{N,T}
end

"""
```julia
NWN_circuit{M}(lines::Vector{Line{T,N}},props::Vector{WireProp{T}},elecs::Vector{Line{T,N}},volts::Vector{T},grnds::Vector{Line{T,N}},dims::SVector{N,T},Rⱼ::T,Rₑ::T)
```
Contains the data to perform circuit analysis of a nanowire
network including all geometric, electrical, and structural
data.
"""
mutable struct NWN_circuit{M,T<:Number,N} # T is numeric type, N is number of spatial dimensions, M is Model parameter
    lines::Vector{Line{T,N}}
    props::Vector{WireProp{T}}
    elecs::Vector{Line{T,N}}
    volts::Vector{T}
    grnds::Vector{Line{T,N}}
    dims::SVector{N,T}
    Rⱼ::T
    Rₑ::T
    graph::SimpleWeightedGraph{Int, T}
    ls::Vector{Int}
    ns::Vector{Int}
    function NWN_circuit{M}(lines::Vector{Line{T,N}},props::Vector{WireProp{T}},elecs::Vector{Line{T,N}},volts::Vector{T},grnds::Vector{Line{T,N}},dims::SVector{N,T},Rⱼ::T,Rₑ::T) where {M,T,N}
        update_graph(new{M,T,N}(lines,props,elecs,volts,grnds,dims,Rⱼ,Rₑ))
    end
end

NWN_circuit{M2}(o::NWN_circuit{M1}) where {M1,M2} = NWN_circuit{M2}(o.lines, o.props, o.elecs, o.volts, o.grnds, o.dims, o.Rⱼ, o.Rₑ)

# Conversion methods
import Base.convert
Line(p₁::Vector{T},p₂::Vector{S}) where {T,S} = Line(promote(SVector{length(p₁)}(p₁), SVector{length(p₂)}(p₂))...)
convert(::Type{Line{T,N}},x::Line{S,N}) where {T,S,N} = Line(convert(SVector{N,T},x.p₁), convert(SVector{N,T},x.p₂))

"""
```julia
len(ℓ::Line{T,N})
```
Determines the Euclidean length of the line segment ℓ.
"""
len(ℓ::Line{T,N}) where {T,N} = sqrt(sum((ℓ.p₁[i]-ℓ.p₂[i])^2 for i ∈ 1:N))


text:
"""
    gmtset(cmd0::String="", kwargs...)

Adjust individual GMT defaults settings in the current directory’s gmt.conf file.

Full option list at [`gmtset`]($(GMTdoc)gmtset.html)

Parameters
----------

- **D** | **units** :: [Type => Str | []]

    Modify the GMT defaults based on the system settings. Append u for US defaults or s for SI defaults.
    ($(GMTdoc)gmtinfo.html#d)
- **G** | **defaultsfile** :: [Type => Str]

    Name of specific gmt.conf file to read and modify.
    ($(GMTdoc)gmtinfo.html#g)
- $(GMT.opt_V)
- $(GMT.opt_write)
"""
function gmtset(cmd0::String=""; kwargs...)

	length(kwargs) == 0 && return monolitic("gmtset", cmd0)

	d = KW(kwargs)
	help_show_options(d)			# Check if user wants ONLY the HELP mode
	cmd = parse_V(d, "")

	cmd = add_opt(d, cmd, 'D', [:D :units], nothing, true)
	cmd = add_opt(d, cmd, 'G', [:G :defaultsfile], nothing, true)

	key = collect(keys(d))
	for k = 1:length(d)
		(key[k] == :Vd)	&& continue
		cmd *= " " * string(key[k]) * " " * string(d[key[k]])
		delete!(d, key[k])
	end

	cmd = "gmtset " * cmd
	if (dbg_print_cmd(d, cmd) !== nothing)  return cmd  end
	gmt(cmd)
	gmt("destroy")
end


text:
using HSA
using FactCheck

type Fable
end

facts("The Macros") do
    context("can provide deterministic destruction") do
        finalized = false

        function finalize_fable(f)
            finalized = true
        end

        @use f = Fable() begin
            finalizer(f, finalize_fable)

            @fact finalized --> false
        end

        @fact finalized --> true
    end
end


text:
using Rasters, Test, Dates, DimensionalData
using Rasters.LookupArrays, Rasters.Dimensions

# RasterSeries from Raster/RasterStack components

data1 = [1 2 3 4
         5 6 7 8]
data2 = 2 * data1
data3 = 3 * data1
data4 = 4 * data1
dimz = X([30, 40]), Y(-10.0:10.0:20.0)
ga1 = Raster(data1, dimz)
ga2 = Raster(data2, dimz)
ga1a = Raster(data3, dimz)
ga2a = Raster(data4, dimz)
stack1 = RasterStack(ga1, ga2; keys=(:ga1, :ga2))
stack2 = RasterStack(ga1a, ga2a; keys=(:ga1, :ga2))
dates =[DateTime(2017), DateTime(2018)]
ser = RasterSeries([stack1, stack2], Ti(dates))
@test issorted(dates)

@testset "getindex returns the currect types" begin
    @test ser[Ti(1)] isa RasterStack{<:NamedTuple}
    @test ser[Ti(1)][:ga2] isa Raster{Int,2}
    @test ser[Ti(1)][:ga2, 1, 1] isa Int
    @test ser[Ti(1)][:ga2][1, 1] isa Int
end

@testset "map" begin
    @test map(x -> x[1], ser) == Raster([(ga1=1, ga2=2), (ga1=3, ga2=4)], dims(ser))
    @test map(x -> x, ser) == ser;
end

@testset "properties" begin
    @test refdims(ser) === ()
    # Should these be real fields? what is the use-case?
    @test metadata(ser) === NoMetadata()
    @test name(ser) === DimensionalData.NoName()
    @test label(ser) === ""
end

@testset "getindex returns the currect results" begin
    @test ser[Ti(Near(DateTime(2017)))][:ga1][X(1), Y(3)] === 3
    @test ser[Ti(At(DateTime(2017)))][:ga1, X(1), Y(3)] === 3
    @test ser[Ti(At(DateTime(2018)))][:ga2][X(2), Y(4)] === 32
    @test ser[Ti(At(DateTime(2018)))][:ga2, X(2), Y(4)] === 32
    @test ser[Ti(1)][:ga1, X(1), Y(2)] == 2
    @test ser[Ti(1)][:ga2, X(2), Y(3:4)] == [14, 16]
end

@testset "getindex is type stable all the way down" begin
    # @inferred ser[Ti(At(DateTime(2017)))][:ga1, X(1), Y(2)]
    @inferred ser[Ti(1)][:ga1][X(1), Y(2)]
    # @inferred ser[Ti(1)][:ga1, X(1), Y(2:4)]
    # @inferred ser[Ti(1)][:ga1][X(1), Y(2:4)]
    # @inferred ser[1][:ga1, X(1:2), Y(:)]
    # @inferred ser[1][:ga1][X(1:2), Y(:)]
end

@testset "setindex!" begin
    ser[1] = ser[2]
    @test typeof(ser[1]) == typeof(ser[2]) == eltype(parent(ser))
    typeof(parent(ser)[1]) == eltype(ser)
    ser[Ti(1)] = ser[Ti(2)]
    @test ser[Ti(1)] == ser[Ti(2)]
end

@testset "rebuild" begin
    @test rebuild(ser, parent(ser)) === ser
    @test rebuild(ser; dims=dims(ser)) === ser
    @test rebuild(ser; dims=(X(),)) !== ser
    @test rebuild(ser; name=nothing) === ser
    @test rebuild(ser; metadata=nothing) === ser
end

@testset "slice, combine" begin
    ga1 = Raster(ones(4, 5, 10), (X(), Y(), Ti(10:10:100))) .* reshape(1:10, (1, 1, 10))
    ga2 = ga1 .* 2
    ser = slice(ga1, Ti)
    @test size(ser) == (10,)
    combined = Rasters.combine(ser, Ti())
    @test combined == ga1
    @test dims(combined) === dims(ga1)
    ser = slice(ga1, (X, Ti))
    @test size(ser) == (4, 10)
    ser = slice(ga1, (X, Y, Ti))
    combined2 = Rasters.combine(ser, (X, Y, Ti))
    @test combined == ga1 == permutedims(combined2, (X, Y, Ti))
    @test dims(combined) === dims(ga1) == dims(permutedims(combined2, (X, Y, Ti)))
    stack = RasterStack((ga1=ga1, ga2=ga2))
    ser = slice(stack, Ti)
    @test size(ser) == (10,)
    combined = Rasters.combine(ser, Ti)
    ser = slice(stack, (Y, Ti))
    @test size(ser) == (5, 10,)
    combined = Rasters.combine(ser, (Y, Ti))
end

@testset "show" begin
    # 2d
    ser2 = slice(ga1, (X, Y))
    sh = sprint(show, MIME("text/plain"), ser2)
    # Test but don't lock this down too much
    @test occursin("RasterSeries", sh)
    @test occursin("Raster", sh)
    @test occursin("X", sh)
    @test occursin("Y", sh)
    # 1d
    ser1 = slice(ga1, X)
    sh = sprint(show, MIME("text/plain"), ser1)
    # Test but don't lock this down too much
    @test occursin("RasterSeries", sh)
    @test occursin("Raster", sh)
    @test occursin("X", sh)
end


text:
function collatz_steps(t::Int64)
    t < 1 && throw(DomainError(t))
    s::Array = []
    while true
        if t % 2 == 0
            t = t ÷ 2
        else
            t = 3t + 1
        end
        push!(s,t)
        t == 1 && break
    end
    return s
end

function collatz_steps_02(t::Int64)
    t < 1 && throw(DomainError(t))
    s = []
    collatz = (x) -> iseven(x) ? x ÷ 2 : 3x + 1
    while true
        t = collatz(t)
        push!(s,t)
        t == 1 && break
    end
    return s
end

function max_stop_time(t::Int64)
    t < 1 && throw(DomainError(t))
    l = []
    w = []
    for i in 1:t
        push!(l, collatz_steps(i))
    end
    for j in l
        push!(w, length(j))
    end
    return findmax(w)
 end

function max_stop_time_02(t::Int64)
    t < 1 && throw(DomainError(t))
    l::Tuple{Int, Array} = (0,[])
    Threads.@threads for i in 1:t
        m = collatz_steps(i)
        if length(m) > length(l[2])
             l = (i,m)
        end
    end
    return l[1],length(l[2])
end

function max_stop_time_03(t::Int64)::Tuple{Int64,Int64}
    collatz = (x::Int64) -> iseven(x) ? x ÷ 2 : 3x + 1
    eval = (s) ->
    begin
        x = 0
        while s > 1
            s = collatz(s)
            x += 1
        end
        x
    end
    m::Tuple{Int64,Int64} = (0,0)
    Threads.@threads for i in 1:t
        if (s = eval(i)) > m[2] m = (i,s) end
    end
    return m
end


text:
@testset "structures.jl" begin
    include(joinpath("structures", "gauging_test.jl"))
    include(joinpath("structures", "ratingcurve_test.jl"))
    include(joinpath("structures", "compoundratingcurve_test.jl"))
end


text:
"""
Struct for SubordinatedBrownianMotion

		subordinatedBrownianMotion=SubordinatedBrownianMotion(sigma::num,drift::num1)

Where:\n
		sigma       =	Volatility.
		drift       = 	drift.
		underlying  = 	underlying.
"""
mutable struct SubordinatedBrownianMotionVec{num <: Number, num1 <: Number, num4 <: Number, Distr <: Distribution{Univariate, Continuous}, numtype <: Number} <: AbstractMonteCarloProcess{numtype}
    sigma::num
    drift::Curve{num1, num4}
    subordinator_::Distr
    function SubordinatedBrownianMotionVec(sigma::num, drift::Curve{num1, num4}, dist::Distr) where {num <: Number, num1 <: Number, num4 <: Number, Distr <: Distribution{Univariate, Continuous}}
        if sigma <= 0
            error("volatility must be positive")
        else
            zero_typed = zero(num) + zero(num1) + zero(num4)
            return new{num, num1, num4, Distr, typeof(zero_typed)}(sigma, drift, dist)
        end
    end
end

function SubordinatedBrownianMotion(σ::num, drift::Curve{num1, num4}, subordinator_::Distr) where {num <: Number, num1 <: Number, num4 <: Number, Distr <: Distribution{Univariate, Continuous}}
    return SubordinatedBrownianMotionVec(σ, drift, subordinator_)
end

export SubordinatedBrownianMotionVec;

function simulate!(X, mcProcess::SubordinatedBrownianMotionVec, mcBaseData::MonteCarloConfiguration{type1, type2, type3, SerialMode, type5}, T::numb) where {numb <: Number, type1 <: Number, type2 <: Number, type3 <: StandardMC, type5 <: Random.AbstractRNG}
    Nsim = mcBaseData.Nsim
    Nstep = mcBaseData.Nstep

    drift = mcProcess.drift
    sigma = mcProcess.sigma
    @assert T > 0
    type_sub = typeof(quantile(mcProcess.subordinator_, 0.5))
    dt_s = Array{type_sub}(undef, Nsim)
    t_s = zeros(type_sub, Nsim)
    dt = T / Nstep
    zero_drift = drift(zero(type_sub), zero(type_sub) + dt) * 0
    isDualZero = sigma * zero(type_sub) * 0 * zero_drift
    @views X[:, 1] .= isDualZero
    Z = Array{Float64}(undef, Nsim)
    for i = 1:Nstep
        randn!(mcBaseData.rng, Z)
        rand!(mcBaseData.rng, mcProcess.subordinator_, dt_s)
        tmp_drift = [drift(t_, dt_) for (t_, dt_) in zip(t_s, dt_s)]
        t_s += dt_s
        # SUBORDINATED BROWNIAN MOTION (dt_s=time change)
        @views X[:, i+1] = X[:, i] + tmp_drift + sigma * sqrt.(dt_s) .* Z
    end

    return
end

function simulate!(X, mcProcess::SubordinatedBrownianMotionVec, mcBaseData::MonteCarloConfiguration{type1, type2, type3, SerialMode, type5}, T::numb) where {numb <: Number, type1 <: Number, type2 <: Number, type3 <: AntitheticMC, type5 <: Random.AbstractRNG}
    Nsim = mcBaseData.Nsim
    Nstep = mcBaseData.Nstep
    type_sub = typeof(quantile(mcProcess.subordinator_, 0.5))
    dt_s = Array{type_sub}(undef, Nsim)
    drift = mcProcess.drift
    sigma = mcProcess.sigma
    @assert T > 0
    t_s = zeros(type_sub, Nsim)
    dt = T / Nstep
    zero_drift = drift(zero(type_sub), zero(type_sub) + dt)
    isDualZero = sigma * zero(type_sub) * 0 * zero_drift
    #X=Matrix{typeof(isDualZero)}(undef,Nsim,Nstep+1);
    @views X[:, 1] .= isDualZero
    for i = 1:Nstep
        NsimAnti = div(Nsim, 2)
        Z = randn(mcBaseData.rng, NsimAnti)
        #this is bad
        Z = [Z; -Z]
        rand!(mcBaseData.rng, mcProcess.subordinator_, dt_s)
        tmp_drift = [drift(t_, dt_) for (t_, dt_) in zip(t_s, dt_s)]
        t_s += dt_s
        # SUBORDINATED BROWNIAN MOTION (dt_s=time change)
        @views X[:, i+1] = X[:, i] .+ tmp_drift .+ sigma .* sqrt.(dt_s) .* Z
    end

    return
end


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Unit tests for !=(x, y)
"""

# --- Imports

# Standard library
import InteractiveUtils.subtypes
using Test

# GeometricAlgebra.jl
using GeometricAlgebra

# --- Tests

# ------ M::Multivector

@testset "!=(M::Multivector, N::Blade)" begin
    @test_skip 1
end

@testset "!=(M::Multivector, N::Pseudoscalar)" begin
    @test_skip 1
end

@testset "!=(M::Multivector, N::Scalar)" begin
    @test_skip 1
end

@testset "!=(M::Multivector, N::One)" begin
    @test_skip 1
end

@testset "!=(M::Multivector, N::Zero)" begin
    @test_skip 1
end

@testset "!=(M::Multivector, N::Real)" begin
    @test_skip 1
end

@testset "!=(M::Multivector, N::Vector)" begin
    @test_skip 1
end

# ------ B::Blade

@testset "!=(B::Blade, M::Multivector)" begin
    @test_skip 1
end

# ------ B::Pseudoscalar

@testset "!=(B::Pseudoscalar, M::Multivector)" begin
    @test_skip 1
end

# ------ B::Scalar

@testset "!=(B::Scalar, M::Multivector)" begin
    @test_skip 1
end

# ------ B::One

@testset "!=(B::One, M::Multivector)" begin
    @test_skip 1
end

# ------ B::Zero

@testset "!=(B::Zero, M::Multivector)" begin
    @test_skip 1
end

# ------ B::Real

@testset "!=(B::Real, M::Multivector)" begin
    @test_skip 1
end

# ------ B::Vector

@testset "!=(B::Vector, M::Multivector)" begin
    @test_skip 1
end


text:
using Test
using MPI

if get(ENV,"JULIA_MPI_TEST_ARRAYTYPE","") == "CuArray"
    import CUDA
    ArrayType = CUDA.CuArray
else
    ArrayType = Array
end

provided = MPI.Init_thread(MPI.THREAD_MULTIPLE)

@test MPI.THREAD_SINGLE <= provided <= MPI.THREAD_MULTIPLE
@test MPI.Query_thread() == provided
@test MPI.Is_thread_main()

comm = MPI.COMM_WORLD
size = MPI.Comm_size(comm)
rank = MPI.Comm_rank(comm)

const N = 10

dst = mod(rank+1, size)
src = mod(rank-1, size)

if provided == MPI.THREAD_MULTIPLE
    send_arr = collect(1.0:N)
    recv_arr = zeros(N)

    reqs = Array{MPI.Request}(undef, 2N)

    Threads.@threads for i = 1:N
        reqs[N+i] = MPI.Irecv!(@view(recv_arr[i:i]), src, i, comm)
        reqs[i] = MPI.Isend(@view(send_arr[i:i]), dst, i, comm)
    end

    MPI.Waitall!(reqs)

    @test recv_arr == send_arr
end

MPI.Finalize()


text:
# v0.6.0

using Requests

str = readstring(get("https://sourceforge.net/"))

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Metatheory

@testset "Reduction Basics" begin
	t = @theory begin
	    a + a => 2a
	    x / x => 1
	    x * 1 => x
	end

    # basic theory to check that everything works
    @test rewrite(:(a + a), t) == :(2a)
    @test rewrite(:(a + a), t) == :(2a)
    @test rewrite(:(a + (x * 1)), t) == :(a + x)
	@test rewrite(:(a + (a * 1)), t; order=:inner) == :(2a)
end

import Base.(+)
@testset "Extending Algebra Operators" begin
    t = @theory begin
        a + a => 2a
    end

    # Let's extend an operator from base, for sake of example
    function +(x::Symbol, y)
        rewrite(:($x + $y), t)
    end

    @test (:x + :x) == :(2x)
end

## Free Monoid

@testset "Free Monoid - Overriding identity" begin
    # support symbol literals
	symbol_monoid = @theory begin
		a ⋅ :ε => a
		:ε ⋅ a => a
		a::Symbol => a
		a::Symbol ⋅ b::Symbol |> Symbol(String(a) * String(b))
		# i |> error("unsupported ", i)
	end;

    @test rewrite(:(ε ⋅ a ⋅ ε ⋅ b ⋅ c ⋅ (ε ⋅ ε ⋅ d) ⋅ e), symbol_monoid; order=:inner) == :abcde
end

## Interpolation should be possible at runtime

@testset "Calculator" begin
	calculator = @theory begin
		x::Number ⊕ y::Number |> x + y
		x::Number ⊗ y::Number |> x * y
		x::Number ⊖ y::Number |> x ÷ y
		x::Symbol => x
		x::Number => x
	end;
	a = 10

	@test rewrite(:(3 ⊕ 1 ⊕ $a), calculator; order=:inner) == 14
end

# FIXME add support for splatting
# @testset "Type assertions and destructuring" begin
#     # let's try type assertions and destructuring
#     t = @theory begin
#         f(a::Number) => a
#         f(a...) => a
#     end

#     @test rewrite(:(f(3)), t) == 3
#     @test rewrite(:(f(2, 3)), t) == [2, 3]

#     # destructuring in right hand
#     n = @theory begin
#         f(a...) => +(a...)
#     end

#     @test rewrite(:(f(2, 3)), n) == :(2 + 3)
#     @test rewrite(:(f(a, b, c, d)), n) == :(((a + b) + c) + d)
# end

## Direct rules
@testset "Direct Rules" begin
    t = @theory begin
        # maps
        a * b |> ((a isa Number && b isa Number) ? a * b : _lhs_expr)
    end
    @test rewrite(:(3 * 1), t) == 3

    t = @theory begin
        # maps
        a::Number * b::Number |> a * b
    end
    @test rewrite(:(3 * 1), t) == 3
end

## Take advantage of subtyping.
# Subtyping in Julia has been formalized in this paper
# [Julia Subtyping: A Rational Reconstruction](https://benchung.github.io/papers/jlsub.pdf)

abstract type Vehicle end
abstract type GroundVehicle <: Vehicle end
abstract type AirVehicle <: Vehicle end
struct Airplane <: AirVehicle end
struct Car <: GroundVehicle end

airpl = Airplane()
car = Car()

@testset "Subtyping" begin
	t = @theory begin
		a::AirVehicle * b => "flies"
		a::GroundVehicle * b => "doesnt_fly"
	end

	sf = rewrite(:($airpl * c), t)
	df = rewrite(:($car * c), t)

    @test sf == "flies"
    @test df == "doesnt_fly"
end

airpl = Airplane()
car = Car()

@testset "Interpolation" begin
	t = @theory begin
		$airpl * b => "flies"
		$car * b => "doesnt_fly"
	end

	sf = rewrite(:($airpl * c), t)
	df = rewrite(:($car * c), t)

    @test sf == "flies"
    @test df == "doesnt_fly"
end

## Multiplicative Monoid

# TODO consider autocompiling from src/Library. This is very ugly.

# FIXME TYPE VARIABLES IN CLASSICAL MATCHCORE ALGORITHM
macro red_commutative_monoid(t, op, id, func)
	zop = Meta.quot(op)
   quote
	   @theory begin
			($op)($id, a) => a
			($op)(a, $id) => a
			# closure
			($op)(a::$t, b::$t) |> ($func)(a,b)
			# associativity reductions implying commutativity
			($op)(a::$t, ($op)(b::$t, c)) |>
				(z = ($func)(a, b); Expr(:call, ($zop), z, c))

			($op)(($op)(a::$t, c), b::$t) |>
				(z = ($func)(a, b); Expr(:call, ($zop), z, c))

			($op)(a::$t, ($op)(c, b::$t)) |>
				(z = ($func)(a, b); Expr(:call, ($zop), z, c))

			($op)(($op)(c, a::$t), b::$t) |>
				(z = ($func)(a, b); Expr(:call, ($zop), z, c))
	   end
   end
end

mult_monoid = @red_commutative_monoid Int (*) 1 (*)

@testset "Multiplication Monoid" begin
	res_macro = rewrite(:(3 * (1 * 2)), mult_monoid)
	res_sym = rewrite(:(3 * (1 * 2)), mult_monoid; order=:inner)
	res_macro_2 = rewrite(:(3 * (1 * 2)), mult_monoid; order=:inner)
	res_macro_3 = rewrite(:(2a * (3 * 2)), mult_monoid; order=:inner)

	@test res_macro == 6
	@test res_macro_2 == 6
	@test res_sym == 6
	@test res_macro_3 == :(12a)
end

# constructs a semantic theory about a an abelian group
# The definition of a group does not require that a ⋅ b = b ⋅ a
# for all elements a and b in G. If this additional condition holds,
# then the operation is said to be commutative, and the group is called an abelian group.
macro red_abelian_group(t, op, id, invop, func)
	zop = Meta.quot(op)
   quote
	   @theory begin
			# identity element
			($op)($id, a) => a
			($op)(a, $id) => a
			# inversibility
			($op)(a, ($invop)(a)) => $id
			($op)(($invop)(a), a) => $id
			# closure
			($op)(a::$t, b::$t) |> ($func)(a,b)
			# inverse
			# associativity reductions
			($op)(a::$t, ($op)(b::$t, c)) |> (z = ($func)(a, b); Expr(:call, ($zop), z, c))
			($op)(($op)(a::$t, c), b::$t) |> (z = ($func)(a, b); Expr(:call, ($zop), z, c))
			($op)(a::$t, ($op)(c, b::$t)) |> (z = ($func)(a, b); Expr(:call, ($zop), z, c))
			($op)(($op)(c, a::$t), b::$t) |> (z = ($func)(a, b); Expr(:call, ($zop), z, c))
	   end
   end
end

addition_group = @red_abelian_group Int (+) 0 (-) (+) ;
@testset "Addition Abelian Group" begin
	zero = rewrite(:((x+y) +  -(x+y)), addition_group)
	@test zero == 0
end

# distributivity of two operations
# example: `@distrib (⋅) (⊕)`
macro red_distrib(outop, inop)
	@assert Base.isbinaryoperator(outop)
	@assert Base.isbinaryoperator(inop)
	quote
		@theory begin
			($outop)(a, ($inop)(b,c)) => ($inop)(($outop)(a,b),($outop)(a,c))
			($outop)(($inop)(a,b), c) => ($inop)(($outop)(a,c),($outop)(b,c))
		end
	end
end
distr = @red_distrib (*) (+)
Z = mult_monoid ∪ addition_group ∪ distr;
@testset "Composing Theories, distributivity" begin
	res_1 = rewrite(:(((2 + (b + -b)) * 3) * (a + b)), Z; order=:inner)
	e = rewrite(:((2 * (3 + b + (4 * 2)))), Z; order=:inner)

	@test_skip res_1 == :(6a+6b)
	@test_skip e == :(22+2b)
end

logids = @theory begin
	x * x => x^2
	x^n * x |> :($x^($(n+1)))
	x * x^n |> :($x^($(n+1)))
	a * (a * b) => a^2 * b
	a * (b * a) => a^2 * b
	(a * b) * a => a^2 * b
	(b * a) * a => a^2 * b
	# catch e as ℯ
	:e => ℯ
	log(a^n) => n * log(a)
	log(x * y) => log(x) * log(y)
	log(1) => 0
	log(:ℯ) => 1
	:ℯ^(log(x)) => x
end;

t = logids ∪ Z;
@testset "Symbolic Logarithmic Identities, Composing Theories" begin
    @test rewrite(:(log(ℯ)), t) == 1
    @test rewrite(:(log(x)), t) == :(log(x))
    @test rewrite(:(log(ℯ^3)), t) == 3
    @test rewrite(:(log(a^3)), t) == :(3 * log(a))
    # Reduce.jl wtf u doing? log(:x^-2 * ℯ^3) = :(log(5021384230796917 / (250000000000000 * x ^ 2)))

    @test rewrite(:(log(x^2 * ℯ^3)), t) == :((6 * log(x)))
    @test_skip rewrite(:(log(x^2 * ℯ^(log(x)))), t; order=:inner) == :(3 * log(x))
end

using Calculus: differentiate
diff = @theory begin
	∂(y)/∂(x) |> differentiate(y, x)
	a * (1/x) => a/x
	a * 0 => 0
	0 * a => 0
end
finalt = t ∪ diff
@testset "Symbolic Differentiation, Accessing Outer Variables" begin
	@test rewrite(:(∂( log(x^2 * ℯ^(log(x))) )/∂(x)), finalt) == :(3/x)
end;

## let's loop this
# @testset "Reduction loop should return." begin
#     t = @theory begin
#         a + b => b + a
#     end

#     @test :(a+b) == rewrite(:(a + b), t)

# 	t = @theory begin
# 	    a + b => b + a
# 	    b + a => a + b
# 	end

# 	@test :(a+b) == rewrite(:(a + b), t)
# end

# safe_var = 0
# ref_var = Ref{Real}(0)

# reft = @theory begin
# 	:safe |> (safe_var = π)
# 	:ref |> (ref_var[] = π)
# end

# @testset "Evil Assignment" begin
# 	rewrite(:(safe), reft; order=:inner)
# 	rewrite(:(ref), reft; order=:inner)

# 	@test safe_var == 0
# 	@test ref_var[] == π

# 	rewrite(ref, reft)
# 	rewrite(safe, reft)

# 	@test safe_var == 0
# 	@test ref_var[] == π
# end;

# dollart = @theory begin
# 	# dollars are interpolated at rule creation time?
# 	$(Expr(:call, :foo, :bar)) => :match
# 	:x + $(3+2) => :weeee
# end

# @testset "Dollar Escapes" begin
# 	r1 = rewrite(:(foo(bar)), dollart; m=@__MODULE__)
# 	r2 = rewrite(:(x + 5), dollart; m=@__MODULE__)

# 	@test r1 == :match
# 	@test r2 == :weeee
# end;


text:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ DYNAMICS FUNCTIONS ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~#

"Propagate the dynamics forward, storing the result in the next knot point"
function propagate_dynamics(::Type{Q}, model::AbstractModel, z_::KnotPoint, z::KnotPoint) where Q<:Implicit
    x_next = discrete_dynamics(Q, model, z)
    set_state!(z_, x_next)
end

"Evaluate the discrete dynamics for all knot points"
function discrete_dynamics!(f, model, Z::Traj)
    for k in eachindex(Z)
        f[k] = discrete_dynamics(model, Z[k])
    end
end

"Evaluate the discrete dynamics Jacobian for all knot points"
function discrete_jacobian!(∇f, model, Z::Traj)
    for k in eachindex(∇f)
        ∇f[k] .= discrete_jacobian(model, Z[k])
    end
end

function discrete_jacobian!(∇f::Vector{<:SMatrix}, model, Z::Traj)
    for k in eachindex(∇f)
        ∇f[k] = discrete_jacobian(model, Z[k])
    end
end


text:
using Revise
using POMDPs
using GridInterpolations
using POMDPModelTools
using POMDPPolicies
using Parameters
using StaticArrays
using DiscreteValueIteration
using Distributed
using AutomotiveDrivingModels
using AutoViz
using AutomotiveSensors
using AutomotivePOMDPs
using PedestrianAvoidancePOMDP

using FileIO
using JLD2

using ArgParse

s = ArgParseSettings()

@add_arg_table s begin
    "--policy_name"
        arg_type = String
        default = "default_name"
        help = "Name of the policy"
    "--policy_type"
        arg_type = String
        default = "longitudinal"
        help = "choose among 'longitudinal', 'lateral'"

end

parsed_args = parse_args(ARGS, s)

policy_name = string(parsed_args["policy_name"], ".jld2")
policy_type = parsed_args["policy_type"]
println("Policy name: ", policy_name)
println("Policy type: ", policy_type)

pomdp = SingleOCFPOMDP()
solver = SparseValueIterationSolver(max_iterations=200, belres=1e-4, include_Q=true, verbose=true)
mdp = UnderlyingMDP(pomdp);

vi_policy = solve(solver, mdp)

policy = AlphaVectorPolicy(pomdp, vi_policy.qmat, vi_policy.action_map)

# save policy!
FileIO.save(policy_name, "policy", policy)


text:
using TailRec
using Test

@testset "TailRec.jl" begin
    # Write your tests here.
end


text:
#Example lowpass/bandpass real/quadrature modulator design
include("dsexample1_fn.jl")

opt=2
#dsm = RealDSM(order=5, OSR=32, form=:CRFB, Hinf=1.5, opt=2)
dsm = RealDSM(order=5, OSR=32, form=:CRFB, Hinf=1.5, opt=opt)
results = dsexample1(dsm, LiveDemo=true)

:END_OF_EXAMPLE


text:
# run tests in parallel (np=4)

"""
  Test the parallel communication primatives in the Utils module
"""
function test_parallel2_comm()

  fname = "input_vals_parallel2.jl"
  mesh, sbp, eqn, opts = solvePDE(fname)

  myrank = mesh.myrank
  commsize = mesh.commsize

  AdvectionEquationMod.ICDict["ICp1"](mesh, sbp, eqn, opts, eqn.q_vec)
  array1DTo3D(mesh, sbp, eqn, opts, eqn.q_vec, eqn.q)
  # add offset so we can tell where data came from
  for i=1:length(eqn.q)
    eqn.q[i] += myrank
  end

#=
  function wrap(i, commsize)
    if i > (commsize-1)
      i = 0
    elseif i < 0
      i = commsize - 1
    end

    return i
  end

  function populate_buffer(mesh, sbp, eqn, opts, data_i)

    return nothing
  end
=#
  @testset "----- Testing Parallel Communication -----" begin

    # parallel_data == PARALLEL_DATA_FACE
    populate_buffer = Utils.getSendDataFace
    shared_data = getSharedFaceData(Float64, mesh, sbp, opts, opts["parallel_data"])
    shared_data_bar = getSharedFaceData(Float64, mesh, sbp, opts, opts["parallel_data"])

    exchangeData(mesh, sbp, eqn, opts, shared_data, populate_buffer)
    calc_func = (mesh, sbp, eqn, opts, data) -> return nothing
    finishExchangeData(mesh, sbp, eqn, opts, shared_data, calc_func)

    testSharedDataFace(mesh, sbp, eqn, opts, shared_data)

    # test exchanging element data
    setParallelData(shared_data, PARALLEL_DATA_ELEMENT)
    populate_buffer = Utils.getSendDataElement

    exchangeData(mesh, sbp, eqn, opts, shared_data, populate_buffer)
    calc_func = (mesh, sbp, eqn, opts, data) -> return nothing
    finishExchangeData(mesh, sbp, eqn, opts, shared_data, calc_func)

    testSharedDataElement(mesh, sbp, eqn, opts, shared_data)

    # test finishExchangeData_rev2

    fill!(eqn.q_bar, 0)
    for i=1:mesh.npeers
      rand!(eqn.shared_data_bar[i].q_recv)
    end
    #TODO: fill recv buffers with random values to test zeroing out
    startSolutionExchange(mesh, sbp, eqn, opts)
    exchangeData_rev(mesh, sbp, eqn, opts, eqn.shared_data, eqn.shared_data_bar, calc_func_rev)
    finishSolutionBarExchange(mesh, sbp, eqn, opts)
    testSharedDataElement_rev(mesh, sbp, eqn, opts)

    calc_func2 = (mesh, sbp, eqn, opts, data, data2) -> return nothing
    setParallelData(shared_data, PARALLEL_DATA_FACE)
    setParallelData(shared_data_bar, PARALLEL_DATA_FACE)
    copy!(eqn.res_bar, eqn.q)
    for i=1:length(eqn.res_bar)
      eqn.res_bar[i] += 1   # add offset to distinguish q and res_bar
    end

    populate_buffer = Utils.getSendDataFace
    populate_buffer_rev = Utils.getSendDataFace_resbar

    #TODO: repalce with startSolutionExchange_rev2
    exchangeData(mesh, sbp, eqn, opts, shared_data, populate_buffer)
    exchangeData(mesh, sbp, eqn, opts, shared_data_bar, populate_buffer_rev)

    finishExchangeData_rev2(mesh, sbp, eqn, opts, shared_data, shared_data_bar,
                           calc_func2)
    testSharedDataFace(mesh, sbp, eqn, opts, shared_data)
    testSharedDataFace(mesh, sbp, eqn, opts, shared_data_bar, 1)

    # make parallel_data different
    setParallelData(shared_data, PARALLEL_DATA_ELEMENT)
    setParallelData(shared_data_bar, PARALLEL_DATA_FACE)

    populate_buffer = Utils.getSendDataElement
    populate_buffer_rev = Utils.getSendDataFace_resbar

    exchangeData(mesh, sbp, eqn, opts, shared_data, populate_buffer)
    exchangeData(mesh, sbp, eqn, opts, shared_data_bar, populate_buffer_rev)

    finishExchangeData_rev2(mesh, sbp, eqn, opts, shared_data, shared_data_bar,
                           calc_func2)

    testSharedDataElement(mesh, sbp, eqn, opts, shared_data)
    testSharedDataFace(mesh, sbp, eqn, opts, shared_data_bar, 1)

  end  # end facts block

  return nothing
end

function calc_func_rev(mesh, sbp, eqn, opts, data::SharedFaceData, data_bar::SharedFaceData)

  # use Remote metrics to add the polynomial + comm_rank to q_recv
  idx = data_bar.peeridx
  coords = mesh.remote_metrics[idx].coords
  myrank = mesh.myrank

  for i=1:size(coords, 3)
    for j=1:size(coords, 2)
      coords_j = sview(coords, :, j, i)
      q_exact = AdvectionEquationMod.calc_p1(eqn.params, coords_j, 0.0) + myrank
      data_bar.q_recv[1, j, i] += q_exact
    end
  end

  return nothing
end

function testSharedDataElement_rev(mesh, sbp, eqn, opts)

  # test that each element of eqn.q_bar has the sum of the contributions of
  # each remote process

  q_bar2 = zeros(eqn.q_bar)
  for i=1:mesh.npeers
    peer = mesh.peer_parts[i]
    for el in mesh.local_element_lists[i]
      for j=1:mesh.numNodesPerElement
        coords_j = sview(mesh.coords, :, j, el)
        q_exact = AdvectionEquationMod.calc_p1(eqn.params, coords_j, 0.0) + peer
        q_bar2[1, j, el] += q_exact
      end
    end
  end

  @test maximum(abs.(q_bar2 - eqn.q_bar)) < 1e-13

  return nothing
end

"""
  Test that face parallel data was sent correctly.  This requires the solution
  on each process is linear and has a constant added to it equal to the
  MPI rank of the process.
"""
function testSharedDataFace(mesh, sbp, eqn, opts, shared_data, offset=0)

  # test the interpolation was exact + the data came from the right place
  # using the peer rank offset
  for i=1:mesh.npeers
    peer_rank = mesh.peer_parts[i]
    coords_peer = mesh.coords_sharedface[i]
    data_i = shared_data[i]
    for j=1:length(data_i.bndries_local)
      for k=1:mesh.numNodesPerFace
        coords_k = sview(coords_peer, :, k, j)
        q_exact = AdvectionEquationMod.calc_p1(eqn.params, coords_k, 0.0) + peer_rank + offset
        @test abs(data_i.q_recv[1, k, j] - q_exact) < 1e-13
      end
    end
  end

  return nothing
end

"""
  Tests element parallel data was sent correctly.  Same requirements as the
  above function
"""
function testSharedDataElement(mesh, sbp, eqn, opts, shared_data)

  q_faceL = zeros(eltype(eqn.q), mesh.numDofPerNode, mesh.numNodesPerFace)
  q_faceR = zeros(q_faceL)
  for i=1:mesh.npeers
    peer_rank = mesh.peer_parts[i]
    data_i = shared_data[i]
    el_offset = mesh.shared_element_offsets[i]
    for j=1:length(data_i.interfaces)
      iface_i = data_i.interfaces[j]

      qL = sview(eqn.q, :, :, iface_i.elementL)
      elR = iface_i.elementR - el_offset + 1
      qR = sview(data_i.q_recv, :, :, elR)

      interiorFaceInterpolate!(mesh.sbpface, iface_i, qL, qR, q_faceL, q_faceR)

      for k=1:mesh.numNodesPerFace
        @test abs(q_faceR[1, k] - peer_rank - (q_faceL[1, k] - mesh.myrank)) < 1e-13
      end
    end
  end

  return nothing
end

add_func1!(AdvectionTests, test_parallel2_comm, [TAG_SHORTTEST, TAG_TMP])


text:
"""
	sample1d(cmd0::String="", arg1=nothing, kwargs...)

Resample 1-D table data using splines

Full option list at [`sample1d`]($(GMTdoc)sample1d.html)

Parameters
----------

- **A** | **resamp** :: [Type => Str]        ``Arg = f|p|m|r|R``

    For track resampling (if -T…unit is set) we can select how this is to be performed.
    ($(GMTdoc)sample1d.html#a)
- **F** | **interp_type** :: [Type => Str]   ``Arg = l|a|c|n[+1|+2]``

    Choose from l (Linear), a (Akima spline), c (natural cubic spline), and n (no interpolation:
    nearest point) [Default is Akima].
    ($(GMTdoc)sample1d.html#f)
- **N** | **time_col** :: [Type => Int]      ``Arg = t_col``

    Indicates which column contains the independent variable (time). The left-most column
    is # 0, the right-most is # (n_cols - 1). [Default is 0].
    ($(GMTdoc)sample1d.html#n)
- **T** | **equi_space** :: [Type => List | Str]     ``Arg = [min/max/]inc[+a|n]] or file|list``

    Evaluate the best-fit regression model at the equidistant points implied by the arguments.
    ($(GMTdoc)sample1d.html#t)
- $(GMT.opt_V)
- $(GMT.opt_write)
- $(GMT.opt_append)
- $(GMT.opt_b)
- $(GMT.opt_d)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_o)
- $(GMT.opt_swap_xy)
"""
function sample1d(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("sample1d", cmd0, arg1)

	d = KW(kwargs)
	cmd = parse_common_opts(d, "", [:V_params :b :d :e :f :g :h :i :o :yx])
	cmd = parse_these_opts(cmd, d, [[:A :resamp], [:F :interp_type], [:N :time_col], [:T :equi_space]])

	common_grd(d, cmd0, cmd, "sample1d ", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
sample1d(arg1, cmd0::String=""; kw...) = sample1d(cmd0, arg1; kw...)


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Compute the beta angle of a satellite.
#
#   The algorithm was based on:
#
#       Mortari, D., Wilkins, M. P., and Bruccoleri, C.  On Sun-Synchronous
#       Orbits and Associated Constellations
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

export satellite_beta_angle

"""
    function satellite_beta_angle(JD0::Number, a::Number, e::Number, i::Number, RAAN::Number, numDays::Integer)

Compute the beta angle of a satellite.

# Args

* `JD0`: Initial instant for the analysis [Julian day].
* `a`: Semi-major axis of the orbit [m].
* `e`: Orbit eccentricity.
* `i`: Orbit inclination [rad].
* `RAAN`: Right ascension of the ascending node at `JD0` [rad].
* `numDays`: Number of days of the analysis.

# Returns

The beta angle [deg] computed for each day.

"""
function satellite_beta_angle(JD0::Number,
                              a::Number,
                              e::Number,
                              i::Number,
                              RAAN::Number,
                              numDays::Integer)
    # Constants
    rad2deg = 180.0/pi

    # Initialization of variables.
    theta = 0.0                   # Sun angle relative to the inertial
                                  # coordinate frame.

    days = collect(0:1:numDays-1) # Vector of the days in which the beta angle
                                  # will be computed.

    N = [0.0; 0.0; 0.0]           # Versor normal to the orbital plane,
                                  # represented in the inertial coordinate
                                  # frame.

    S = [0.0; 0.0; 0.0]           # Versor that points to the Sun, represented
                                  # in the inertial coordinate frame.

    # Output vector.
    beta = Vector{Float64}(undef,numDays)

    # RAAN rotation rate [rad/day].
    dOmega = dRAAN(a, e, i, :J2)*24.0*3600.0

    # Loop
    for t in days
        # Compute the RAAN at the day d.
        RAAN_d = RAAN + dOmega*t

        # Compute the versor N represented in the Inertial ref. frame.
        Dio = angle_to_dcm(-i, -RAAN_d, 0.0, :XZX)
        N_i = Dio*SVector{3}(0,0,1)

        # Compute the Sun position at noon (UT) represented in the Inertial ref.
        # frame.
        S_i = sun_position_i(JD0+t)
        S_i = S_i/norm(S_i)

        # Get the angle between N_i and S_i [deg].
        beta[t+1] = 90.0-acos(dot(N_i,S_i))*rad2deg
    end

    [days beta]
end


text:
mutable struct BlackCalculator{S <: StrikedTypePayoff}
  payoff::S
  strike::Float64
  forward::Float64
  stdDev::Float64
  discount::Float64
  variance::Float64
  d1::Float64
  d2::Float64
  alpha::Float64
  beta::Float64
  DalphaDd1::Float64
  DbetaDd2::Float64
  n_d1::Float64
  cum_d1::Float64
  n_d2::Float64
  cum_d2::Float64
  x::Float64
  DxDs::Float64
  DxDstrike::Float64
end

function BlackCalculator(p::S, fwd::Float64, stdDev::Float64, disc::Float64) where {S <: StrikedTypePayoff}
  calc = BlackCalculator{S}(p, p.strike, fwd, stdDev, disc, stdDev * stdDev, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0)
  initialize!(calc, p)

  return calc
end

function gen_addl_black_vars!(calc::BlackCalculator, optType::Call)
  calc.alpha = calc.cum_d1
  calc.DalphaDd1 = calc.n_d1
  calc.beta = -calc.cum_d2
  calc.DbetaDd2 = -calc.n_d2

  return calc
end

function gen_addl_black_vars!(calc::BlackCalculator, optType::Put)
  calc.alpha = -1.0 + calc.cum_d1
  calc.DalphaDd1 = calc.n_d1
  calc.beta = 1.0 - calc.cum_d2
  calc.DbetaDd2 = -calc.n_d2

  return calc
end

function initialize!(calc::BlackCalculator, p::StrikedTypePayoff)
  calc.strike >= 0.0 || error("strike $strike must be non-negative")
  calc.forward > 0.0 || error("forward $forward must be positive")
  calc.stdDev >= 0.0 || error("stdDev $stdDev must be non-negative")
  calc.discount > 0.0 || error("discount $discount must be positive")

  if calc.stdDev >= eps()
    if is_close(calc.strike, 0.0)
      calc.d1 = JQ_MAX
      calc.d2 = JQ_MAX
      calc.cum_d1 = 1.0
      calc.cum_d2 = 1.0
      calc.n_d1 = 0.0
      calc.n_d2 = 0.0
    else
      calc.d1 = log(calc.forward / calc.strike) / calc.stdDev + 0.5 * calc.stdDev
      calc.d2 = calc.d1 - calc.stdDev
      w = Normal() # Distributions
      calc.cum_d1 = cdf(w, calc.d1)
      calc.cum_d2 = cdf(w, calc.d2)
      calc.n_d1 = distribution_derivative(w, calc.d1)
      calc.n_d2 = distribution_derivative(w, calc.d2)
    end
  else
    if is_close(calc.forward, calc.strike)
      calc.d1 = 0.0
      calc.d2 = 0.0
      calc.cum_d1 = 0.5
      calc.cum_d2 = 0.5
      calc.n_d1 = M_SQRT_2 * M_1_SQRTPI
      calc.n_d2 = M_SQRT_2 * M_1_SQRTPI
    elseif calc.forward > calc.strike
      calc.d1 = JQ_MAX
      calc.d2 = JQ_MAX
      calc.cum_d1 = 1.0
      calc.cum_d2 = 1.0
      calc.n_d1 = 0.0
      calc.n_d2 = 0.0
    else
      calc.d1 = JQ_MIN
      calc.d2 = JQ_MIN
      calc.cum_d1 = 0.0
      calc.cum_d2 = 0.0
      calc.n_d1 = 0.0
      calc.n_d2 = 0.0
    end
  end

  calc.x = calc.strike
  calc.DxDstrike = 1.0

  calc.DxDs = 0.0

  gen_addl_black_vars!(calc, p.optionType)

  return calc
end

value(calc::BlackCalculator) = calc.discount * (calc.forward * calc.alpha + calc.x * calc.beta)

itm_cash_probability(calc::BlackCalculator) = calc.cum_d2
itm_asset_probability(calc::BlackCalculator) = calc.cum_d1

function delta(calc::BlackCalculator, spot::Float64)
  spot > 0.0 || error("positive spot value required.  $spot not allowed")

  DfowardDs = calc.forward / spot

  temp = calc.stdDev * spot
  DalphaDs = calc.DalphaDd1 / temp
  DbetaDs = calc.DbetaDd2 / temp
  temp2 = DalphaDs * calc.forward + calc.alpha * DfowardDs + DbetaDs * calc.x + calc.beta * calc.DxDs

  return calc.discount * temp2
end

function delta_forward(calc::BlackCalculator)
  temp = calc.stdDev * calc.forward
  DalphaDforward = calc.DalphaDd1 / temp
  DbetaDforward = calc.DbetaDd2 / temp
  temp2 = DalphaDforward * calc.forward + calc.alpha + DbetaDforward * calc.x

  return calc.discount * temp2
end

function elasticity(calc::BlackCalculator, spot::Float64)
  val = value(calc)
  del = delta(calc, spot)

  if val > eps()
    return del / val * spot
  elseif abs(del) < eps()
    return 0.0
  elseif del > 0.0
    return JQ_MAX
  else
    return JQ_MIN
  end
end

function gamma(calc::BlackCalculator, spot::Float64)
  spot > 0.0 || error("positive spot value required.  $spot not allowed")

  DforwardDs = calc.forward / spot

  temp = calc.stdDev * spot
  DalphaDs = calc.DalphaDd1 / temp
  DbetaDs = calc.DbetaDd2 / temp

  D2alphaDs2 = -DalphaDs / spot * (1 + calc.d1 / calc.stdDev)
  D2betaDs2 = -DbetaDs / spot * (1 + calc.d2 / calc.stdDev)

  temp2 = D2alphaDs2 * calc.forward + 2.0 * DalphaDs * DforwardDs + D2betaDs2 * calc.x + 2.0 * DbetaDs * calc.DxDs

  return calc.discount * temp2
end

function rho(calc::BlackCalculator, mat::Float64)
  mat >= 0.0 || error("negative maturity not allowed")

  DalphaDr = calc.DalphaDd1 / calc.stdDev
  DbetaDr = calc.DbetaDd2 / calc.stdDev

  temp = DalphaDr * calc.forward + calc.alpha * calc.forward + DbetaDr * calc.x

  return mat * (calc.discount * temp - value(calc))
end

function dividend_rho(calc::BlackCalculator, mat::Float64)
  mat >= 0.0 || error("negative maturity not allowed")

  DalphaDq = -calc.DalphaDd1 / calc.stdDev
  DbetaDq = -calc.DbetaDd2 / calc.stdDev

  temp = DalphaDq * calc.forward - calc.alpha * calc.forward + DbetaDq * calc.x

  return mat * calc.discount * temp
end

function vega(calc::BlackCalculator, mat::Float64)
  temp = log(calc.strike / calc.forward) / calc.variance
  DalphaDsigma = calc.DalphaDd1 * (temp + 0.5)
  DbetaDsigma = calc.DbetaDd2 * (temp - 0.5)

  temp2 = DalphaDsigma * calc.forward + DbetaDsigma * calc.x

  return calc.discount * sqrt(mat) * temp2
end

theta_per_day(calc::BlackCalculator, spot::Float64, mat::Float64) = theta(calc, spot, mat) / 365.0

function theta(calc::BlackCalculator, spot::Float64, mat::Float64)
  if mat < 0.0
    println("WARN: Maturity is negative")
    return -1.0
  end

  is_close(mat, 0.0) && return 0.0

  return -( log(calc.discount) * value(calc) + log(calc.forward / spot) * spot * delta(calc, spot)
          + 0.5 * calc.variance * spot * spot * gamma(calc, spot)) / mat
end

function strike_sensitivity(calc::BlackCalculator)
  temp = calc.stdDev * calc.strike
  DalphaDstrike = -calc.DalphaDd1 / temp
  DbetaDstrike = -calc.DbetaDd2 / temp

  temp2 = DalphaDstrike * calc.forward + DbetaDstrike * calc.x + calc.beta * calc.DxDstrike

  return calc.discount * temp2
end


text:
module TestDoctest

using PerformanceTestTools
using Test

function test(; skip = true)
    if skip
        if lowercase(get(ENV, "JULIA_PKGEVAL", "false")) == "true"
            @info "Skipping doctests on PkgEval."
            @test_skip nothing
            return
        elseif VERSION < v"1.6"
            @info "Skipping doctests on Julia $VERSION."
            @test_skip nothing
            return
        end
    end
    PerformanceTestTools.@include_foreach("__test_doctest.jl", [[]])
end

end  # module