text:
using DynamicAxisWarping
using DynamicAxisWarping.Datasets
using Plots
gr()

data, labels = ucr_traindata("50words")

plot(
    [fastdtw(data[:, i], data[:, i+10], radius)[1] for radius = 1:100, i = 1:10],
    show = true,
)
xlabel!("radius")


text:
module GetBlankets

using MLJ, DataFrames, ShapML
using StatsBase: sample
using CategoricalArrays: CategoricalArray
using Combinatorics: powerset
using Tables
xgc = @load XGBoostClassifier;
xgr = @load XGBoostRegressor;
rfc = MLJ.@load RandomForestClassifier  pkg = DecisionTree
rfr = MLJ.@load RandomForestRegressor   pkg = DecisionTree
include("ICPlibrary.jl")

export getBlanketAll, getBlanketBoosting, getBlanketRandomForest

function getBlanketAll(X::DataFrame)
    if ncol(X) >= 2
        s = collect(powerset(1:ncol(X)))
        return (s[2:end, 1])
    elseif ncol(X) == 1
        return ([[1]])
    else
        return ([[]])
    end
end # getBlanketAll

function getBlanketBoosting(
    X::DataFrame,
    target::DataFrame,
    selection::String="booster",
    maxNoVariables::Int64=10,
)
    if selection == "all" || ncol(X) <= maxNoVariables
        return ((blankets = getBlanketAll(X), featureImportances = missing))
    else
        indexFromPredictor = getIndexFromPredictor(X)
        (usingVariables, featureImportances) = selectXGBooster(X, target, maxNoVariables)

        if length(usingVariables) >= 1
            usingVariablesIndex = [indexFromPredictor[n] for n in usingVariables]
            s = sort(usingVariablesIndex)
            s = collect(powerset(s))
            return ((blankets = s[2:end, 1], featureImportances = featureImportances))
        else
            return ((blankets = [[]], featureImportances = missing))
        end
    end
end # getBlanketBoosting

function selectXGBooster(
    X::DataFrame,
    target::DataFrame,
    maxNoVariables::Int64,
)
    selectionVec = []
    featureImportances = trainXGBooster(X, target)
    println(featureImportances)

    gainSum = 0
    variableCount = 0
    for f in featureImportances
        gainSum += f.gain
        variableCount += 1
        push!(selectionVec, f.fname)
        if variableCount >= maxNoVariables || variableCount >= length(featureImportances) || gainSum >= 0.90
            break
        end
    end

    println("XGBooster selection:  ", selectionVec)
    GC.gc()
    return ((selectionVec = selectionVec, featureImportances = featureImportances))
end # selectXGBooster

function trainXGBooster(X::DataFrame, y::DataFrame)
    GC.gc()
    if length(unique(y[:, 1])) > 2
        pipeXGBoostRegressor = @pipeline XGBoostRegressor
        r1 = range(pipeXGBoostRegressor, :(xg_boost_regressor.max_depth), lower=3, upper=10)
        r2 = range(pipeXGBoostRegressor, :(xg_boost_regressor.num_round), lower=1, upper=25)

        tmXGBoostRegressor = TunedModel(
            model=pipeXGBoostRegressor,
            tuning=Grid(resolution=7),
            resampling=CV(rng=11),
            ranges=[r1, r2],
            measure=rms,
        )
        mtm = machine(
            tmXGBoostRegressor,
            setScientificTypes!(X),
            Float64.(y[:, 1]),
        )
        fit!(mtm)
        k = collect(keys(report(mtm).best_report.report_given_machine))[1]
        return (report(mtm).best_report.report_given_machine[k][1])
    else
        pipeXGBoostClassifier = @pipeline XGBoostClassifier
        r1 = range(pipeXGBoostClassifier, :(xg_boost_classifier.max_depth), lower=3, upper=10)
        r2 = range(pipeXGBoostClassifier, :(xg_boost_classifier.num_round), lower=1, upper=25)

        tmXGBoostClassifier = TunedModel(
            model=pipeXGBoostClassifier,
            tuning=Grid(resolution=7),
            resampling=CV(rng=11),
            ranges=[r1, r2],
            measure=cross_entropy, # don't use rms for probabilistic responses
        )
        mtm = machine(
            tmXGBoostClassifier,
            setScientificTypes!(X),
            categorical(y[:, 1]),
        )
        fit!(mtm)
        k = collect(keys(report(mtm).best_report.report_given_machine))[1]
        return (report(mtm).best_report.report_given_machine[k][1])
    end
end # trainXGBooster

#####################################################################################################################################################
function getBlanketRandomForest(
    X::DataFrame,
    target::DataFrame,
    selection::String="forest",
    maxNoVariables::Int64=10,
)
    if selection == "all"
        return ((blankets = getBlanketAll(X), featureImportances = missing))
    elseif selection == "booster"
        return (getBlanketBoosting(X, target, selection, maxNoVariables))
    else
        indexFromPredictor = getIndexFromPredictor(X)
        (usingVariables, featureImportances) = selectRandomForest(X, target, maxNoVariables)
        if length(usingVariables) >= 1
            usingVariablesIndex = [indexFromPredictor[String(n)] for n in usingVariables]
            s = sort(usingVariablesIndex)
            s = collect(powerset(s))
            return ((blankets = s[2:end, 1], featureImportances = featureImportances))
        else
            return ((blankets = [[]], featureImportances = missing))
        end
    end
end # getBlanketRandomForest

function selectRandomForest(
    X::DataFrame,
    target::DataFrame,
    maxNoVariables::Int64=10,
)
    selectionVec = []
    featureImportances = trainRandomForest(X, target)

    variableCount = 0
    meanEffectPercentSum = 0
    for f in featureImportances
        meanEffectPercentSum += f.meanEffectPercent
        variableCount += 1
        push!(selectionVec, f.feature_name)
        if variableCount >= maxNoVariables || variableCount >= length(featureImportances) || f.meanEffectPercent < 0.05 || meanEffectPercentSum >= 0.90
            break
        end
    end

    println("Random Forest with Shapley selection:  ", selectionVec)
    GC.gc()
    return ((selectionVec = selectionVec, featureImportances = featureImportances))
end # selectRandomForest

function predict_function(model, data)
    data_pred = DataFrame(y_pred=predict(model, data))
    return data_pred
end # predict_function

function predict_function_mode(model, data)
    ŷ = MLJ.predict(model, data)
    ŷMode = [convert(Int64, mode(ŷ[i])) for i in 1:length(ŷ)]
    data_pred = DataFrame(y_pred=ŷMode)
    return data_pred
end # predict_function_mode

function trainRandomForest(
    X::DataFrame,
    y::DataFrame
)
    GC.gc()
    if length(unique(y[:, 1])) > 2
        pipeRandomForestRegressor = @pipeline FeatureSelector RandomForestRegressor
        cases = [[Symbol(names(X)[j]) for j in 1:i] for i in 1:ncol(X)]
        r1 = range(pipeRandomForestRegressor, :(feature_selector.features), values=cases)
        tmRandomForestRegressor = TunedModel(
            model=pipeRandomForestRegressor,
            range=r1,
            measures=rms,
            resampling=CV(nfolds=5)
        )
        mtm = machine(tmRandomForestRegressor, setScientificTypes!(X), Float64.(y[:, 1]))
        Base.invokelatest(MLJ.fit!, mtm)

        predictor = predict_function
    else
        pipeRandomForestClassifier = @pipeline FeatureSelector RandomForestClassifier #prediction_type = :probabilistic
        cases = [[Symbol(names(X)[j]) for j in 1:i] for i in 1:ncol(X)]
        r1 = range(pipeRandomForestClassifier, :(feature_selector.features), values=cases)
        tmRandomForestClassifier = TunedModel(
            model=pipeRandomForestClassifier,
            range=r1,
            measures=[cross_entropy, BrierScore()],
            resampling=CV(nfolds=5)
        )
        mtm = machine(tmRandomForestClassifier, setScientificTypes!(X), categorical(y[:, 1]))
        Base.invokelatest(MLJ.fit!, mtm)

        predictor = predict_function_mode
    end

    dfShapMeanEffect = Shapley(X, mtm, predictor)
    println("Shapley Effect of Random Forest\n", dfShapMeanEffect, "\n")
    return (Tables.rowtable(dfShapMeanEffect))
end # trainRandomForest

function Shapley(
    X::DataFrame,
    mtm,
    predictor
)
    dataShap = Missing
    sample_size = 60            # Number of Monte Carlo samples for Shapley, must be > 30 to apply the Central Limit Theorem
    if ncol(X) > 12
    # ShapML is slow and runs out of memory on larger datasets
        if nrow(X) > 5_000
            lDivisors = [3,4,6,8,10,12]
        else
            lDivisors = [2,4,6,8,10,12]
        end
    else
        lDivisors = [1,2,3,4,6,8,10]
    end
    for d in lDivisors
        GC.gc()
        r = sample(1:nrow(X), Int(round(nrow(X) / d)), replace=false)
        explain = copy(X[r, :])     # Compute Shapley feature-level predictions
        reference = copy(X)         # An optional reference population to compute the baseline prediction.
        try
            println("Computing Shapley Effect of Random Forest using ", length(r), " random rows.  It might take a few minutes")
            dataShap = ShapML.shap( explain=explain,
                            reference=reference,
                            model=mtm,
                            predict_function=predictor,
                            sample_size=sample_size,
                            seed=20200629
            )
            GC.gc()
            break
        catch
            continue
        end
    end
    dfShapMean = DataFrames.by(dataShap, [:feature_name], mean_effect=[:shap_effect] => x -> mean(abs.(x.shap_effect)))
    dfShapMeanEffect = sort(dfShapMean, order(:mean_effect, rev=true))
    totalEffect = sum(dfShapMeanEffect.mean_effect)
    dfShapMeanEffect.meanEffectPercent = dfShapMeanEffect.mean_effect / totalEffect
    return dfShapMeanEffect
end # Shapley

end # model

ext: jl


text:
# Deprecated functionality

import Base: @deprecate_binding, @deprecate_moved false

@deprecate_binding SCHED_AUTO CTX_SCHED_AUTO false
@deprecate_binding SCHED_SPIN CTX_SCHED_SPIN false
@deprecate_binding SCHED_YIELD CTX_SCHED_YIELD false
@deprecate_binding SCHED_BLOCKING_SYNC CTX_SCHED_BLOCKING_SYNC false
@deprecate_binding MAP_HOST CTX_MAP_HOST false
@deprecate_binding LMEM_RESIZE_TO_MAX CTX_LMEM_RESIZE_TO_MAX false

@deprecate_binding MAX_THREADS_PER_BLOCK DEVICE_ATTRIBUTE_MAX_THREADS_PER_BLOCK false
@deprecate_binding MAX_BLOCK_DIM_X DEVICE_ATTRIBUTE_MAX_BLOCK_DIM_X false
@deprecate_binding MAX_BLOCK_DIM_Y DEVICE_ATTRIBUTE_MAX_BLOCK_DIM_Y false
@deprecate_binding MAX_BLOCK_DIM_Z DEVICE_ATTRIBUTE_MAX_BLOCK_DIM_Z false
@deprecate_binding MAX_GRID_DIM_X DEVICE_ATTRIBUTE_MAX_GRID_DIM_X false
@deprecate_binding MAX_GRID_DIM_Y DEVICE_ATTRIBUTE_MAX_GRID_DIM_Y false
@deprecate_binding MAX_GRID_DIM_Z DEVICE_ATTRIBUTE_MAX_GRID_DIM_Z false
@deprecate_binding MAX_SHARED_MEMORY_PER_BLOCK DEVICE_ATTRIBUTE_MAX_SHARED_MEMORY_PER_BLOCK false
@deprecate_binding TOTAL_CONSTANT_MEMORY DEVICE_ATTRIBUTE_TOTAL_CONSTANT_MEMORY false
@deprecate_binding WARP_SIZE DEVICE_ATTRIBUTE_WARP_SIZE false
@deprecate_binding MAX_PITCH DEVICE_ATTRIBUTE_MAX_PITCH false
@deprecate_binding MAX_REGISTERS_PER_BLOCK DEVICE_ATTRIBUTE_MAX_REGISTERS_PER_BLOCK false
@deprecate_binding CLOCK_RATE DEVICE_ATTRIBUTE_CLOCK_RATE false
@deprecate_binding TEXTURE_ALIGNMENT DEVICE_ATTRIBUTE_TEXTURE_ALIGNMENT false
@deprecate_binding GPU_OVERLAP DEVICE_ATTRIBUTE_GPU_OVERLAP false
@deprecate_binding MULTIPROCESSOR_COUNT DEVICE_ATTRIBUTE_MULTIPROCESSOR_COUNT false
@deprecate_binding KERNEL_EXEC_TIMEOUT DEVICE_ATTRIBUTE_KERNEL_EXEC_TIMEOUT false
@deprecate_binding INTEGRATED DEVICE_ATTRIBUTE_INTEGRATED false
@deprecate_binding CAN_MAP_HOST_MEMORY DEVICE_ATTRIBUTE_CAN_MAP_HOST_MEMORY false
@deprecate_binding COMPUTE_MODE DEVICE_ATTRIBUTE_COMPUTE_MODE false
@deprecate_binding MAXIMUM_TEXTURE1D_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE1D_WIDTH false
@deprecate_binding MAXIMUM_TEXTURE2D_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_WIDTH false
@deprecate_binding MAXIMUM_TEXTURE2D_HEIGHT DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_HEIGHT false
@deprecate_binding MAXIMUM_TEXTURE3D_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_WIDTH false
@deprecate_binding MAXIMUM_TEXTURE3D_HEIGHT DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_HEIGHT false
@deprecate_binding MAXIMUM_TEXTURE3D_DEPTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_DEPTH false
@deprecate_binding MAXIMUM_TEXTURE2D_LAYERED_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_LAYERED_WIDTH false
@deprecate_binding MAXIMUM_TEXTURE2D_LAYERED_HEIGHT DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_LAYERED_HEIGHT false
@deprecate_binding MAXIMUM_TEXTURE2D_LAYERED_LAYERS DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_LAYERED_LAYERS false
@deprecate_binding SURFACE_ALIGNMENT DEVICE_ATTRIBUTE_SURFACE_ALIGNMENT false
@deprecate_binding CONCURRENT_KERNELS DEVICE_ATTRIBUTE_CONCURRENT_KERNELS false
@deprecate_binding ECC_ENABLED DEVICE_ATTRIBUTE_ECC_ENABLED false
@deprecate_binding PCI_BUS_ID DEVICE_ATTRIBUTE_PCI_BUS_ID false
@deprecate_binding PCI_DEVICE_ID DEVICE_ATTRIBUTE_PCI_DEVICE_ID false
@deprecate_binding TCC_DRIVER DEVICE_ATTRIBUTE_TCC_DRIVER false
@deprecate_binding MEMORY_CLOCK_RATE DEVICE_ATTRIBUTE_MEMORY_CLOCK_RATE false
@deprecate_binding GLOBAL_MEMORY_BUS_WIDTH DEVICE_ATTRIBUTE_GLOBAL_MEMORY_BUS_WIDTH false
@deprecate_binding L2_CACHE_SIZE DEVICE_ATTRIBUTE_L2_CACHE_SIZE false
@deprecate_binding MAX_THREADS_PER_MULTIPROCESSOR DEVICE_ATTRIBUTE_MAX_THREADS_PER_MULTIPROCESSOR false
@deprecate_binding ASYNC_ENGINE_COUNT DEVICE_ATTRIBUTE_ASYNC_ENGINE_COUNT false
@deprecate_binding UNIFIED_ADDRESSING DEVICE_ATTRIBUTE_UNIFIED_ADDRESSING false
@deprecate_binding MAXIMUM_TEXTURE1D_LAYERED_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE1D_LAYERED_WIDTH false
@deprecate_binding MAXIMUM_TEXTURE1D_LAYERED_LAYERS DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE1D_LAYERED_LAYERS false
@deprecate_binding CAN_TEX2D_GATHER DEVICE_ATTRIBUTE_CAN_TEX2D_GATHER false
@deprecate_binding MAXIMUM_TEXTURE2D_GATHER_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_GATHER_WIDTH false
@deprecate_binding MAXIMUM_TEXTURE2D_GATHER_HEIGHT DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_GATHER_HEIGHT false
@deprecate_binding MAXIMUM_TEXTURE3D_WIDTH_ALTERNATE DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_WIDTH_ALTERNATE false
@deprecate_binding MAXIMUM_TEXTURE3D_HEIGHT_ALTERNATE DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_HEIGHT_ALTERNATE false
@deprecate_binding MAXIMUM_TEXTURE3D_DEPTH_ALTERNATE DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE3D_DEPTH_ALTERNATE false
@deprecate_binding PCI_DOMAIN_ID DEVICE_ATTRIBUTE_PCI_DOMAIN_ID false
@deprecate_binding TEXTURE_PITCH_ALIGNMENT DEVICE_ATTRIBUTE_TEXTURE_PITCH_ALIGNMENT false
@deprecate_binding MAXIMUM_TEXTURECUBEMAP_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURECUBEMAP_WIDTH false
@deprecate_binding MAXIMUM_TEXTURECUBEMAP_LAYERED_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURECUBEMAP_LAYERED_WIDTH false
@deprecate_binding MAXIMUM_TEXTURECUBEMAP_LAYERED_LAYERS DEVICE_ATTRIBUTE_MAXIMUM_TEXTURECUBEMAP_LAYERED_LAYERS false
@deprecate_binding MAXIMUM_SURFACE1D_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_SURFACE1D_WIDTH false
@deprecate_binding MAXIMUM_SURFACE2D_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_SURFACE2D_WIDTH false
@deprecate_binding MAXIMUM_SURFACE2D_HEIGHT DEVICE_ATTRIBUTE_MAXIMUM_SURFACE2D_HEIGHT false
@deprecate_binding MAXIMUM_SURFACE3D_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_SURFACE3D_WIDTH false
@deprecate_binding MAXIMUM_SURFACE3D_HEIGHT DEVICE_ATTRIBUTE_MAXIMUM_SURFACE3D_HEIGHT false
@deprecate_binding MAXIMUM_SURFACE3D_DEPTH DEVICE_ATTRIBUTE_MAXIMUM_SURFACE3D_DEPTH false
@deprecate_binding MAXIMUM_SURFACE1D_LAYERED_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_SURFACE1D_LAYERED_WIDTH false
@deprecate_binding MAXIMUM_SURFACE1D_LAYERED_LAYERS DEVICE_ATTRIBUTE_MAXIMUM_SURFACE1D_LAYERED_LAYERS false
@deprecate_binding MAXIMUM_SURFACE2D_LAYERED_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_SURFACE2D_LAYERED_WIDTH false
@deprecate_binding MAXIMUM_SURFACE2D_LAYERED_HEIGHT DEVICE_ATTRIBUTE_MAXIMUM_SURFACE2D_LAYERED_HEIGHT false
@deprecate_binding MAXIMUM_SURFACE2D_LAYERED_LAYERS DEVICE_ATTRIBUTE_MAXIMUM_SURFACE2D_LAYERED_LAYERS false
@deprecate_binding MAXIMUM_SURFACECUBEMAP_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_SURFACECUBEMAP_WIDTH false
@deprecate_binding MAXIMUM_SURFACECUBEMAP_LAYERED_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_SURFACECUBEMAP_LAYERED_WIDTH false
@deprecate_binding MAXIMUM_SURFACECUBEMAP_LAYERED_LAYERS DEVICE_ATTRIBUTE_MAXIMUM_SURFACECUBEMAP_LAYERED_LAYERS false
@deprecate_binding MAXIMUM_TEXTURE1D_LINEAR_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE1D_LINEAR_WIDTH false
@deprecate_binding MAXIMUM_TEXTURE2D_LINEAR_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_LINEAR_WIDTH false
@deprecate_binding MAXIMUM_TEXTURE2D_LINEAR_HEIGHT DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_LINEAR_HEIGHT false
@deprecate_binding MAXIMUM_TEXTURE2D_LINEAR_PITCH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_LINEAR_PITCH false
@deprecate_binding MAXIMUM_TEXTURE2D_MIPMAPPED_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_MIPMAPPED_WIDTH false
@deprecate_binding MAXIMUM_TEXTURE2D_MIPMAPPED_HEIGHT DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE2D_MIPMAPPED_HEIGHT false
@deprecate_binding COMPUTE_CAPABILITY_MAJOR DEVICE_ATTRIBUTE_COMPUTE_CAPABILITY_MAJOR false
@deprecate_binding COMPUTE_CAPABILITY_MINOR DEVICE_ATTRIBUTE_COMPUTE_CAPABILITY_MINOR false
@deprecate_binding MAXIMUM_TEXTURE1D_MIPMAPPED_WIDTH DEVICE_ATTRIBUTE_MAXIMUM_TEXTURE1D_MIPMAPPED_WIDTH false
@deprecate_binding STREAM_PRIORITIES_SUPPORTED DEVICE_ATTRIBUTE_STREAM_PRIORITIES_SUPPORTED false
@deprecate_binding GLOBAL_L1_CACHE_SUPPORTED DEVICE_ATTRIBUTE_GLOBAL_L1_CACHE_SUPPORTED false
@deprecate_binding LOCAL_L1_CACHE_SUPPORTED DEVICE_ATTRIBUTE_LOCAL_L1_CACHE_SUPPORTED false
@deprecate_binding MAX_SHARED_MEMORY_PER_MULTIPROCESSOR DEVICE_ATTRIBUTE_MAX_SHARED_MEMORY_PER_MULTIPROCESSOR false
@deprecate_binding MAX_REGISTERS_PER_MULTIPROCESSOR DEVICE_ATTRIBUTE_MAX_REGISTERS_PER_MULTIPROCESSOR false
@deprecate_binding MANAGED_MEMORY DEVICE_ATTRIBUTE_MANAGED_MEMORY false
@deprecate_binding MULTI_GPU_BOARD DEVICE_ATTRIBUTE_MULTI_GPU_BOARD false
@deprecate_binding MULTI_GPU_BOARD_GROUP_ID DEVICE_ATTRIBUTE_MULTI_GPU_BOARD_GROUP_ID false
@deprecate_binding HOST_NATIVE_ATOMIC_SUPPORTED DEVICE_ATTRIBUTE_HOST_NATIVE_ATOMIC_SUPPORTED false
@deprecate_binding SINGLE_TO_DOUBLE_PRECISION_PERF_RATIO DEVICE_ATTRIBUTE_SINGLE_TO_DOUBLE_PRECISION_PERF_RATIO false
@deprecate_binding PAGEABLE_MEMORY_ACCESS DEVICE_ATTRIBUTE_PAGEABLE_MEMORY_ACCESS false
@deprecate_binding CONCURRENT_MANAGED_ACCESS DEVICE_ATTRIBUTE_CONCURRENT_MANAGED_ACCESS false
@deprecate_binding COMPUTE_PREEMPTION_SUPPORTED DEVICE_ATTRIBUTE_COMPUTE_PREEMPTION_SUPPORTED false
@deprecate_binding CAN_USE_HOST_POINTER_FOR_REGISTERED_MEM DEVICE_ATTRIBUTE_CAN_USE_HOST_POINTER_FOR_REGISTERED_MEM false
@deprecate_binding CAN_USE_STREAM_MEM_OPS DEVICE_ATTRIBUTE_CAN_USE_STREAM_MEM_OPS false
@deprecate_binding CAN_USE_64_BIT_STREAM_MEM_OPS DEVICE_ATTRIBUTE_CAN_USE_64_BIT_STREAM_MEM_OPS false
@deprecate_binding CAN_USE_STREAM_WAIT_VALUE_NOR DEVICE_ATTRIBUTE_CAN_USE_STREAM_WAIT_VALUE_NOR false
@deprecate_binding COOPERATIVE_LAUNCH DEVICE_ATTRIBUTE_COOPERATIVE_LAUNCH false
@deprecate_binding COOPERATIVE_MULTI_DEVICE_LAUNCH DEVICE_ATTRIBUTE_COOPERATIVE_MULTI_DEVICE_LAUNCH false
@deprecate_binding MAX_SHARED_MEMORY_PER_BLOCK_OPTIN DEVICE_ATTRIBUTE_MAX_SHARED_MEMORY_PER_BLOCK_OPTIN false
@deprecate_binding CAN_FLUSH_REMOTE_WRITES DEVICE_ATTRIBUTE_CAN_FLUSH_REMOTE_WRITES false
@deprecate_binding HOST_REGISTER_SUPPORTED DEVICE_ATTRIBUTE_HOST_REGISTER_SUPPORTED false
@deprecate_binding PAGEABLE_MEMORY_ACCESS_USES_HOST_PAGE_TABLES DEVICE_ATTRIBUTE_PAGEABLE_MEMORY_ACCESS_USES_HOST_PAGE_TABLES false
@deprecate_binding DIRECT_MANAGED_MEM_ACCESS_FROM_HOST DEVICE_ATTRIBUTE_DIRECT_MANAGED_MEM_ACCESS_FROM_HOST false

@deprecate_binding MAX_REGISTERS JIT_MAX_REGISTERS false
@deprecate_binding THREADS_PER_BLOCK JIT_THREADS_PER_BLOCK false
@deprecate_binding WALL_TIME JIT_WALL_TIME false
@deprecate_binding INFO_LOG_BUFFER JIT_INFO_LOG_BUFFER false
@deprecate_binding INFO_LOG_BUFFER_SIZE_BYTES JIT_INFO_LOG_BUFFER_SIZE_BYTES false
@deprecate_binding ERROR_LOG_BUFFER JIT_ERROR_LOG_BUFFER false
@deprecate_binding ERROR_LOG_BUFFER_SIZE_BYTES JIT_ERROR_LOG_BUFFER_SIZE_BYTES false
@deprecate_binding OPTIMIZATION_LEVEL JIT_OPTIMIZATION_LEVEL false
@deprecate_binding TARGET_FROM_CUCONTEXT JIT_TARGET_FROM_CUCONTEXT false
@deprecate_binding TARGET JIT_TARGET false
@deprecate_binding FALLBACK_STRATEGY JIT_FALLBACK_STRATEGY false
@deprecate_binding GENERATE_DEBUG_INFO JIT_GENERATE_DEBUG_INFO false
@deprecate_binding LOG_VERBOSE JIT_LOG_VERBOSE false
@deprecate_binding GENERATE_LINE_INFO JIT_GENERATE_LINE_INFO false
@deprecate_binding CACHE_MODE JIT_CACHE_MODE false

@deprecate_binding CUBIN JIT_INPUT_CUBIN false
@deprecate_binding PTX JIT_INPUT_PTX false
@deprecate_binding FATBINARY JIT_INPUT_FATBINARY false
@deprecate_binding OBJECT JIT_INPUT_OBJECT false
@deprecate_binding LIBRARY JIT_INPUT_LIBRARY false

@deprecate_binding CuContext_t CUcontext false
@deprecate_binding CuStream_t CUstream false
@deprecate_binding CuEvent_t CUevent false
@deprecate_binding CuModule_t CUmodule false
@deprecate_binding CuDevice_t CUdevice false
@deprecate_binding CuFunction_t CUfunction false
@deprecate_binding CuLinkState_t CUlinkState false

# this function is broken, but let's at least keep downstream packages working
@deprecate synchronize(ctx::CuContext) synchronize()

import Base: get
Base.@deprecate get(var::CuGlobal) getindex(var)
Base.@deprecate set(var::CuGlobal, val) setindex!(var, val)

@deprecate CuStream(flags::CUstream_flags) CuStream(; flags=flags)


text:
function calibrate(x, poly)
    return sum(x.^i .* poly[i] for i in 1:length(poly))
end

function calibrate(x, shift, gain, gain²)
    return shift + x*gain + x^2*gain²
end

function calibrate(x::Real, shift, gain)
    return shift + (x+rand()*2-1)*gain
end

function ucalibrate(x, shift, gain)
    return shift + x*gain
end

function ucalibrate(x, shift, gain, gain²)
    return shift + gain*x + gain²*x^2
end

function ucalibrate(x, shift, gain, gain², gain³)
    return shift + x*gain + x^2 * gain² + x^3*gain³
end

function gainshift!(m::Array{T, 1}, shift, gain, gain²) where T <: Integer
    tmp = zeros(T, size(m))
    first = firstindex(m)
    last = lastindex(m)
    @inbounds @fastmath for i in eachindex(m)
        j = unsafe_trunc.(T, calibrate.(i .+ 2 .*rand(m[i]) .- 1, shift, gain, gain²))
        for _j in j[first .≤ j .≤ last]
            tmp[_j] += 1
        end
    end
    m .= tmp
end

function gainshift!(m::Array{T, 1}, shift, gain) where T <: Integer
    tmp = zeros(T, size(m))
    first = firstindex(m)
    last = lastindex(m)
    @inbounds @fastmath for i in eachindex(m)
        j = unsafe_trunc.(T, ucalibrate.(i .+ 2 .*rand(m[i]) .- 1, shift, gain))
        for _j in j[first .≤ j .≤ last]
            tmp[_j] += 1
        end
    end
    m .= tmp
end

function gainshift!(m::OMatrix{Array{T, 1}}, shift, gain, gain²) where T <: Integer
    gainshift!(m.matrix, shift, gain, gain²)
end

function gainshift(m::Array{T, 1}, shift, gain, gain²) where T <: Integer
    tmp = zeros(T, size(m))
    first = firstindex(m)
    last = lastindex(m)
    @inbounds @fastmath for i in eachindex(m)
        j = unsafe_trunc.(T, ucalibrate.(i .+ 2 .*rand(m[i]) .- 1, shift, gain, gain²))
        for _j in j[first .≤ j .≤ last]
            tmp[_j] += 1
        end
    end
    tmp
end

function gainshift(m::Array{T, 1}, shift, gain) where T <: Integer
    tmp = zeros(T, size(m))
    first = firstindex(m)
    last = lastindex(m)
    @inbounds @fastmath for i in eachindex(m)
        j = unsafe_trunc.(T, ucalibrate.(i .+ 2 .*rand(m[i]) .- 1, shift, gain))
        for _j in j[first .≤ j .≤ last]
            tmp[_j] += 1
        end
    end
    tmp
end

function gainshift(m::Array{T, 1}, shift, gain, gain², gain³) where T <: Integer
    tmp = zeros(T, size(m))
    first = firstindex(m)
    last = lastindex(m)
    @inbounds @fastmath for i in eachindex(m)
        j = unsafe_trunc.(T, ucalibrate.(i .+ 2 .*rand(m[i]) .- 1, shift, gain, gain², gain³))
        for _j in j[first .≤ j .≤ last]
            tmp[_j] += 1
        end
    end
    tmp
end

function gainshift(m::OMatrix{Array{T, 2}}, coeffsx, coeffxy) where T <: Integer
    throw("Gainshift for OMatrix not yet implemented")
end

function calibrate(x, event::MiniEvent, value::Symbol, parameters::Parameters)::Float64
    if value == :e || value == :Δe
        coeffs = parameters.gainshift[value][event.back, event.front, :]::Array{Float64, 1}
        return calibrate(x, coeffs...)
    elseif value == :ex
        coeffs = parameters.exfromeΔe[event.front, :]::Array{Float64, 1}
        theoretical = calibrate(x, coeffs...)
        return parameters.excorr[event.front, 1] + theoretical*parameters.excorr[event.front, 2]
    elseif value == :range

    end
    error("Unsupported symbol for calibration")
end

function calibrate(event::MiniEvent, value::Symbol, parameters::Parameters)::Float64
    coeffs = parameters.gainshift[value][event.back, event.front, :]::Array{Float64, 1}
    return calibrate(getfield(event, value), coeffs...)
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
module DarcyFEOperatorsTests

##
using Test
using Gridap
import Gridap: ∇

import LinearAlgebra: det
using Base: transpose
##

# model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,2))
# model.cgrid
# ps = points(model.cgrid)
# typeof(ps)
# collect(ps)
#
#
#
# f(x) = Point(x[1],5*x[2])
#
# broadcast(f,ps)
#
#
# struct MappedCartesianGrid <: CartesianGrid
#   map
#
# end
# with mapping
#
#
# trian = Triangulation(model)
#
#
# phi = CellGeomap(trian)
# jac = gradient(phi)
# jact = transpose(jac)
#
#
# piola_map = inv(jact)
#
#
# graph = GridGraph(model)
# pols = CellPolytopes(Grid(model))
# pt = pols.value
#
# order = 1
# _reffe = NedelecRefFE(pt,Float64,order)
#
# _labels = FaceLabels(model)
#
# # V = ConformingFESpace(_reffe,trian,graph,_labels,[5,6])
# V = CurlConformingFESpace(_reffe,trian,graph,_labels,[5,6])
##

u(x) = VectorValue(x[1]*x[2], -0.5*x[2]^2)
∇u(x) = TensorValue(x[2],0.0,x[1],-x[2])
divu(x) = 0.0
∇(::typeof(u)) = ∇u

p(x) = x[1] + x[2]
∇p(x) = VectorValue(1.0,1.0)
∇(::typeof(p)) = ∇p

f(x) = divu(x)
g(x) = p(x)

const kinv_1 = TensorValue(1.0,0.0,0.0,1.0)
r(x) = kinv_1*u(x) + ∇p(x)

gmap(x) = 2*x
model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(50,50), map = gmap)
order = 2
V = FESpace( reffe=:RaviartThomas, conformity=:HDiv, order=2, model=model, diritags = [5,6])
_Q = FESpace( reffe=:QLagrangian, conformity=:L2, valuetype = Float64, order = 1,
                    model = model)
# If all faces of the domain are zero values, we must fix zero mean value
# _Q = FESpace( reffe=:QLagrangian, conformity=:L2, valuetype = Float64, order = 1,
                    # model = model, constraint = :zeromean)

V_0 = TestFESpace(V)
Q = TestFESpace(_Q)
Y = [V_0, Q]

V_g = TrialFESpace(V,u)
X = [V_g, Q]

trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)

@law σ(x,u) = kinv_1*u

function a(y,x)
  v, q = y
  u, p = x
  inner(v,σ(u)) - inner(div(v),p) + inner(q,div(u))
end

function b_Ω(y)
  v, q = y
  inner(v,r) + inner(q,f)
end

t_Ω = AffineFETerm(a,b_Ω,trian,quad)

neumanntags = [7,8]
btrian = BoundaryTriangulation(model,neumanntags)
bquad = CellQuadrature(btrian,degree=order*2)
nb = NormalVector(btrian)

function b_Γ(y)
  v, q = y
  - inner(v*nb,g)
end

t_Γ = FESource(b_Γ,btrian,bquad)

op = LinearFEOperator(Y,X,t_Ω,t_Γ)

uh, ph = solve(op)

e1 = u - uh
e2 = p - ph

# Define norms to measure the error
l2(u) = inner(u,u)
hdiv(u) = inner(div(u),div(u)) + l2(u)

# Compute errors
e1l2 = sqrt(sum( integrate(l2(e1),trian,quad) ))
e1hdiv = sqrt(sum( integrate(hdiv(e1),trian,quad) ))
e2l2 = sqrt(sum( integrate(l2(e2),trian,quad) ))

# writevtk(trian,"/home/santiago/github-repos/Gridap/tmp/darcyresults",
         # cellfields=["uh"=>uh,"ph"=>ph, "euh"=> e1, "eph"=> e2])

@test e1l2 < 1.e-8
@test e1hdiv < 1.e-8
@test e2l2 < 1.e-8
##
end # module


text:
import Base.showerror

struct NMRToolsException <: Exception
    message::String
end

Base.showerror(io::IO, e::NMRToolsException) = printstyled(io, "(NMRTools) " * e.message, color=:red)


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using Test

Test.@testset "plotting" begin
end


text:
###
  Name     : c10_31_timeUsed.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 4/24/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
###

using QuantEcon
nodes, weights = qnwlege(65, -2pi, 2pi);
integral = do_quad(x -> cos(x), nodes, weights)
@time quadgk(x -> cos.(x), -2pi, 2pi)
@time do_quad(x -> cos.(x), nodes, weights)


text:
module NonlinearSolverInterface

abstract type AbstNonlinearSolver{Tvec} end

"""
A concrete type for nonlinear solver should at least include the following members:
  residual! -> function residual!(res::Tvec, u::Tvec)
  jacobian! -> function jacobian!(jac::SparseMatrix, u::Tvec)
  linsolver -> type of AbstLinearSolver
  precond   -> a preconditioner object of type AbstPreconditioner
  debug     -> a function used only for debugging, which should be left undefined
               in a productive version

  These memebers are set using the following functions
"""
function setResidual!(solver::AbstNonlinearSolver{Tvec},
                      residual!::AbstFunction{Tvec}) where Tvec
  solver.residual! = residual!
end

function setLinearSolver!(solver::AbstNonlinearSolver{Tvec},
                          linsolver::AbstLinearSolver{Tvec})
  solver.linsolver = linsolver
end

function setJacobian!(solver::AbstNonlinearSolver{Tvec},
                      jacobian!::AbstLinearOperator{Tvec})
  solver.jacobian = jacobian!
end

function setPreconditioner!(precond::AbstLinearOperator{Tvec},
                            solver::AbstNonlinearSolver{Tvec})
  solver.precond = precond
end

struct InexactNewton{Tvec} <: AbstNonlinearSolver{Tvec}
  res::Tvec
  dx::Tvec
  outprint::Function  # for information output
  histprint::Function # for convergence history plotting
  debug::Function     # for debugging
  vecnorm::Function   # to define non-standard norms
  state::Dict         # data at each iteration   Anthony wants this changed to a different name
  residual!::AbstFunction{Tvec} # to evaluate the residual
  linsolver::AbstLinearSolver{Tvec} # to solve linear system
  jacobian::AbstLinearOperator{Tvec}
  precond::AbstLinearOperator{Tvec}

  function InexactNewton{Tvec}(;outprint::Function=(AbstractString)->nothing,
                               histprint::Function=(AbstractString)->nothing,
                               vecnorm::Function=norm) where Tvec
    solver = new()
    solver.res = Tvec()
    solver.dx  = Tvec()
    solver.outprint = outprint
    solver.histprint = histprint
    solver.debug = debug
    solver.vecnorm = vecnorm
    solver.state = Dict()

    return solver
  end
end

"""
Update the jacobian matrix and preconditioner.

Input/Output:
jac :: jacobian matrix
prec:: preconditioner

Input:
x :: at which value \p jac and \p prec are evaluated
freeze_jac :: don't update jac
freeze_prec:: don't update preconditioner
"""
function updateJacAndPrec!(jac::AbstMatrix{Tvec},
                           precond::AbstPreconditioner{Tve},
                           x::Tvec;
                           freeze_jac::Bool=false,
                           freeze_prec::Bool=false)
end

function solve!(x::Tvec, solver::InexactNewton{Tvec}, options::Dict)
  # clear state dictionary
  empty!(solver.state)

  # write some data to file
  solver.outprint(repeat('=',80))
  solver.outprint("Nonlinear solve method: Inexact Newton")
  solver.histprint("# " * rpad("iter",5) * rpad("norm",10) * rpad("rel. norm",10))

  # initial residual and update, then begin Newton iterations
  fill!(solver.res, 0.0) # need to define fill!() for our parallel vectors
  fill!(solver.dx, 0.0)  # neither of these fill! statements may be needed

  for k = 1:get!(options, "max_iter", 10)
    solver.state["iter"] = k

    # might want some generic timing calls around the iteration and parts of it;
    # this could also be accomplished with a member function like vecnorm or
    # histprint.

    debug(sovler, res, x)

    # get the residual and its norm
    # Probably will have to use a try catch statement here
    if solver.residual!(solver.res, x) == false
      solver.outprint("failure: residual! invalid at solution")
      return false
    end
    resnorm = norm(solver, solver.res)
    # if we store the whole history in state, we can avoid storing resnorm0
    solver.state["resnorm"] = resnorm # do we forego the variable resnorm for state?
    if k == 1
      resnorm0 = resnorm
      state["resnorm0"] = resnorm0
    end

    # check for convergence
    # This should be an optional function
    if resnorm < options["abs_tol"]
      solver.outprint("success: absolute tolerance satisfied")
      return true
    elseif resnorm < resnorm0*options["rel_tol"]
      solver.outprint("success: relative tolerance satisfied")
      return true
    end

    # Update the matrices, if necessary
    # This is where "policy" decision can be made.  Where should updateMatrices! be defined?
    solver.updateMatrices!(solver.jacobian, solver.precond, x, solver.state)

    # solve for the Newton update; this will call a linear solver that is
    # similar in spirit to this nonlinear solver.
    solve!(solver.dx, solver.linsolver, solver.res, solver.jacobian,
           solver.precond, options["linsolver_options"])

    # update the solution; would need a line search here in general
    axpy!(1.0, solver.dx, x) # this function needs to be generalized for parallel vectors
  end
  # if we get here, then we failed to converge in max_iter iterations
  solver.outprint("Inexact Newton failed to converge in "*
                  string(options["max_iter"])*" iterations")
  return false
end

function debug(solver, x, res, state)
  # if mpi_rank == 3
  # open your file, using state["iter"] to help name it
  # convert x to q format

end

end # of module NonlinearSolverInterface

max_stars_repo_name: OptimalDesignLab/CSEInterfaces

max_issues_repo_name: OptimalDesignLab/CSEInterfaces

max_forks_repo_name: OptimalDesignLab/CSEInterfaces


text:
module ESS

function help_categories()
    Base.Help.init_help()
    for cat = Base.Help.CATEGORY_LIST
        if !isempty(Base.Help.CATEGORY_DICT[cat])
            println("\"$cat\" ")
        end
    end
end

function all_help_topics()
    Base.Help.init_help()
    ## show all categories
    for cat = Base.Help.CATEGORY_LIST
        if !isempty(Base.Help.CATEGORY_DICT[cat])
            println()
            show(cat); println();
            for func = Base.Help.CATEGORY_DICT[cat]
                print("  ")
                show(func)
            end
        end
    end
end

function help(topic::String)
    Base.Help.init_help()
    if !haskey(Base.Help.CATEGORY_DICT, topic)
        # if it's not a category, try another named thing
        try
            obj = eval(current_module(), parse(topic))
            Base.Help.help(obj)
        catch
            print("No help information found")
        end
    else
        ln = length(topic)
        bar = " ===" * "="^ln * "==="
        println(bar)
        println(" =  ", topic, "  =")
        println(bar, "\n")
        for func = Base.Help.CATEGORY_DICT[topic]
            Base.Help.help(func)
            println()
            println("-"^72)
        end
    end
end

## modified version of function show(io::IO, m::Method)
function fun_args(m::Method)
        tv = m.tvars
        io = STDOUT::IO
        if !isa(tv,Tuple)
            tv = (tv,)
        end
        if !isempty(tv)
            Base.show_delim_array(io, tv, '{', ',', '}', false)
        end
        li = m.func.code
        e = Base.uncompressed_ast(li)
        argnames = e.args[1]
        decls = map(Base.argtype_decl_string, argnames, {m.sig...})
        print(io, "(")
        print_joined(io, decls, ",", ",")
        print(io, ")")
end

## modified versionof show(io::IO, mt::MethodTable)
function fun_args(f::Function)
    mt = f.env
    mod = Base.function_module(f)
    if mod == Main
        mod = "nil"
    end
    print("(list \"$mod\" nil '(")
    d = mt.defs
    while !is(d,())
        print("\"")
        fun_args(d)
        print("\" ")
        d = d.next
        if is(d,())
            print("))")
        end
    end
end

function fun_args(s::ASCIIString)
    try
        m = eval(current_module(), parse(s))
        if typeof(m) != ASCIIString
            fun_args(m)
        end
    catch
        print("(list nil nil nil)")
    end
end

function fun_args(t::DataType)
    print("(list nil nil '(")
    for d = names(t)
        print("\"$d\" ")
    end
    print("))")
end

### OBJECT COMPLETION
function components(m::Module)
    for v in sort(names(m))
        s = string(v)
        if isdefined(m,v)
            println(rpad(s, 30), summary(eval(m,v)))
        end
    end
end

function components(t::DataType)
    for v in sort(names(t))
        println(rpad(string(v), 30), "field")
    end
end

function components(v)
    t = typeof(v)
    if isa(t, DataType)
        return components(t)
    end
end

### MISC
function main_modules()
    mainmod = current_module()
    for nm in names(mainmod)
        if isdefined(mainmod, nm)
            mod = eval(mainmod, nm)
            if isa(mod, Module)
                print("\"$nm\" ")
            end
        end
    end
end

end

max_stars_repo_name: samrain/emacs-config

max_issues_repo_name: samrain/emacs-config

max_forks_repo_name: samrain/emacs-config


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

#################################################################################
# Exports

# Base Models
export Stage
export Simulation
export OperationsProblem
export OperationsProblemTemplate
export InitialCondition
export SimulationSequence

#Network Relevant Exports
export StandardPTDFModel
export CopperPlatePowerModel
export AreaBalancePowerModel

######## Device Models ########
export DeviceModel
export FixedOutput
######## Service Models ########
export ServiceModel
export RangeReserve
export StepwiseCostReserve
export PIDSmoothACE
export GroupReserve
######## Branch Models ########
export StaticLine
export StaticLineBounds
export StaticLineUnbounded
export StaticTransformer
export StaticTransformerBounds
export StaticTransformerUnbounded
export FlowMonitoredLine
export HVDCLossless
export HVDCDispatch
export HVDCUnbounded
#export VoltageSourceDC
######## Load Models ########
export StaticPowerLoad
export InterruptiblePowerLoad
export DispatchablePowerLoad
######## Renewable Formulations ########
export RenewableFullDispatch
export RenewableConstantPowerFactor
######## Hydro Formulations ########
export HydroDispatchRunOfRiver
export HydroDispatchReservoirBudget
export HydroDispatchReservoirStorage
export HydroCommitmentRunOfRiver
export HydroCommitmentReservoirBudget
export HydroCommitmentReservoirStorage
export HydroDispatchPumpedStorage
export HydroDispatchPumpedStoragewReservation
######## Renewable Formulations ########
export BookKeeping
export BookKeepingwReservation
######## Thermal Formulations ########
export ThermalStandardUnitCommitment
export ThermalBasicUnitCommitment
export ThermalDispatch
export ThermalRampLimited
export ThermalDispatchNoMin
export ThermalMultiStartUnitCommitment

###### Regulation Device Formulation #######
export DeviceLimitedRegulation
export ReserveLimitedRegulation

# feedforward chrons
export RecedingHorizon
export Synchronize
export Consecutive
export FullHorizon
export Range

# feedforward models
export UpperBoundFF
export SemiContinuousFF
export RangeFF
export IntegralLimitFF
export ParameterFF

# InitialConditions chrons
export InterStageChronology
export IntraStageChronology

# Initial Conditions Quantities
export DevicePower
export DeviceStatus
export TimeDurationON
export TimeDurationOFF
export EnergyLevel

# cache_models
export TimeStatusChange
export StoredEnergy

#operation_models
export GenericOpProblem
export UnitCommitmentProblem
export EconomicDispatchProblem
#export OptimalPowerFlow

# Functions
## Construction Exports
export construct_device!
export construct_network!
## Op Model Exports
export solve!
export get_initial_conditions
export set_transmission_model!
export set_devices_template!
export set_branches_template!
export set_services_template!
export set_device_model!
export set_branch_model!
export set_device_model!
export set_model!
## Sim Model Exports
export build!
export execute!
export make_references
## Template Exports
export template_economic_dispatch
export template_unit_commitment
export template_agc_reserve_deployment
export EconomicDispatchProblem
export UnitCommitmentProblem
export AGCReserveDeployment
export run_economic_dispatch
export run_unit_commitment
## Results interfaces
export get_duals

## Utils Exports
export SimulationResultsReference
export write_results
export check_file_integrity
export load_results
export load_operation_results
export load_simulation_results
export write_to_CSV
export get_all_constraint_index
export get_all_var_index
export get_con_index
export get_var_index
export get_result_variable
export get_variable_names
export configure_logging
export show_recorder_events
export list_simulation_events
export show_simulation_events

# Variables / Parameters
export ACTIVE_POWER
export ENERGY
export ENERGY_BUDGET
export FLOW_ACTIVE_POWER
export ON
export REACTIVE_POWER
export ACTIVE_POWER_IN
export ACTIVE_POWER_OUT
export RESERVE
export SERVICE_REQUIREMENT
export START
export STOP
export THETA
export VM
export INFLOW
export SPILLAGE
export SLACK_UP
export SLACK_DN

# Constraints
export ACTIVE
export ACTIVE_RANGE
export ACTIVE_RANGE_LB
export ACTIVE_RANGE_UB
export COMMITMENT
export DURATION
export DURATION_DOWN
export DURATION_UP
export ENERGY_CAPACITY
export ENERGY_LIMIT
export FEEDFORWARD
export FEEDFORWARD_UB
export FEEDFORWARD_BIN
export FEEDFORWARD_INTEGRAL_LIMIT
export FLOW_LIMIT
export FLOW_LIMIT_FROM_TO
export FLOW_LIMIT_TO_FROM
export FLOW_REACTIVE_POWER_FROM_TO
export FLOW_REACTIVE_POWER_TO_FROM
export FLOW_ACTIVE_POWER_FROM_TO
export FLOW_ACTIVE_POWER_TO_FROM
export FLOW_ACTIVE_POWER
export FLOW_REACTIVE_POWER
export INPUT_POWER_RANGE
export OUTPUT_POWER_RANGE
export RAMP
export RAMP_DOWN
export RAMP_UP
export RATE_LIMIT
export RATE_LIMIT_FT
export RATE_LIMIT_TF
export REACTIVE
export REACTIVE_RANGE
export REQUIREMENT
export INFLOW_RANGE

#################################################################################
# Imports
import Logging
import Serialization
#Modeling Imports
import JuMP
# so that users do not need to import JuMP to use a solver with PowerModels
import JuMP: optimizer_with_attributes
export optimizer_with_attributes
import MathOptInterface
import ParameterJuMP
import LinearAlgebra
import PowerSystems
import InfrastructureSystems
# so that users have access to IS.Results interfaces
import InfrastructureSystems:
    get_variables, get_total_cost, get_optimizer_log, write_results, get_timestamp, get_name
export get_name
export get_model_base_power
export get_variables
export get_dual_values
export get_total_cost
export get_optimizer_log
export get_timestamp
export write_results
import PowerModels
import TimerOutputs

#TimeStamp Management Imports
import Dates
import TimeSeries

#I/O Imports
import DataFrames
import Feather
import JSON
import CSV
import SHA

include("core/definitions.jl")

################################################################################
# Includes

include("logging.jl")
include("utils.jl")

#Models and constructors
include("core/results.jl")
include("core/abstract_types.jl")
include("core/aux_structs.jl")

include("services_models/services_model.jl")
include("devices_models/device_model.jl")
include("network_models/networks.jl")

include("core/parameters.jl")
include("core/variables.jl")
include("core/constraints.jl")
include("core/cache.jl")
include("core/initial_condition_types.jl")
include("core/initial_condition.jl")
include("core/initial_conditions.jl")
include("core/operations_problem_template.jl")
include("core/psi_settings.jl")
include("core/psi_container.jl")
include("core/update_initial_conditions.jl")
include("core/operations_problem_results.jl")
include("core/operations_problem.jl")
include("core/simulation_stages.jl")
include("core/simulation_sequence.jl")
include("core/simulation.jl")
include("devices_models/devices/common.jl")
include("core/feedforward.jl")
include("core/simulation_results.jl")
include("core/recorder_events.jl")

#Device Modeling components
include("devices_models/devices/common/device_range_constraints.jl")
include("devices_models/devices/common/nodal_expression.jl")
include("devices_models/devices/renewable_generation.jl")
include("devices_models/devices/thermal_generation.jl")
include("devices_models/devices/electric_loads.jl")
include("devices_models/devices/AC_branches.jl")
include("devices_models/devices/DC_branches.jl")
include("devices_models/devices/storage.jl")
include("devices_models/devices/hydro_generation.jl")
include("devices_models/devices/regulation_device.jl")

#Services Models
include("services_models/agc.jl")
include("services_models/reserves.jl")
include("services_models/GroupReserve.jl")
include("services_models/service_slacks.jl")
include("services_models/services_constructor.jl")

#Network models
include("network_models/copperplate_model.jl")
include("network_models/powermodels_interface.jl")
include("network_models/ptdf_model.jl")
include("network_models/network_slack_variables.jl")
include("network_models/area_balance_model.jl")

#Device constructors
include("devices_models/device_constructors/common/constructor_validations.jl")
include("devices_models/device_constructors/thermalgeneration_constructor.jl")
include("devices_models/device_constructors/hydrogeneration_constructor.jl")
include("devices_models/device_constructors/branch_constructor.jl")
include("devices_models/device_constructors/renewablegeneration_constructor.jl")
include("devices_models/device_constructors/load_constructor.jl")
include("devices_models/device_constructors/storage_constructor.jl")
include("devices_models/device_constructors/regulationdevice_constructor.jl")

#Network constructors
include("network_models/network_constructor.jl")

#Templates
include("operations_problems_templates.jl")
# Printing
include("printing.jl")

end


text:
# Data structure and method for cross-links
# Alex Tam, 12/10/2020

"Data structure for a cross-linker"
struct Cross_Link
    f1::Actin_Filament # 1st cross-linked actin filament
    f2::Actin_Filament # 2nd cross-linked actin filament
    seg1_index::Int # Index of relevant segment on 1st filament
    seg2_index::Int # Index of relevant segment on 2nd filament
    relative_pos1::Float64 # Relative position along the 1st segment
    relative_pos2::Float64 # Relative position along the 2nd segment
    t1::Vector{Int} # Translation of the 1st segment
    t2::Vector{Int} # Translation of the 2nd segment
end

"Obtain dimensionless cross-linker position"
function get_xl_pos(l::Cross_Link, s::State)
    # Find actin segment nodes
    m1x, m1y, p1x, p1y = get_segment_nodes(l.f1, l.f1.segments[l.seg1_index], s);
    m2x, m2y, p2x, p2y = get_segment_nodes(l.f2, l.f2.segments[l.seg2_index], s);
    # Compute the cross-link site
    x1 = m1x + l.relative_pos1*(p1x - m1x);
    y1 = m1y + l.relative_pos1*(p1y - m1y);
    x2 = m2x + l.relative_pos2*(p2x - m2x);
    y2 = m2y + l.relative_pos2*(p2y - m2y);
    return x1, y1, x2, y2
end


text:
using OpenQuantumTools, Test

u0 = PauliVec[1][2]
ρ0 = u0*u0'

@test u0 == OpenQuantumTools.build_u0(u0, :v, vectorize=false)
@test ρ0 == OpenQuantumTools.build_u0(u0, :m, vectorize=false)
@test ρ0 == OpenQuantumTools.build_u0(ρ0, :m, vectorize=false)
@test ρ0[:] == OpenQuantumTools.build_u0(ρ0, :m, vectorize=true)


text:
"""
Constructs LTI system with quadratic reward and uniformly sampled initial conditions.
"""
function LTIQuadraticEnv(
    A::AbstractVecOrMat,
    B::AbstractVecOrMat,
    C::AbstractVecOrMat,
    D::AbstractVecOrMat,
    Q::AbstractVecOrMat,
    R::AbstractVecOrMat,
    s0_lb::Union{Real,Vector{<:Real}},
    s0_ub::Union{Real,Vector{<:Real}},
    tspan::Tuple{<:Real,<:Real},
    dt::Real;
    #= Keyword arguments =#
    o_ub::Union{Nothing,Real,Vector{<:Real}}=nothing, # upper bound for observation space
    o_lb::Union{Nothing,Real,Vector{<:Real}}=nothing, # lower bound for observation space
    a_ub::Union{Nothing,Real,Vector{<:Real}}=nothing, # upper bound for action space
    a_lb::Union{Nothing,Real,Vector{<:Real}}=nothing, # lower bound for action space
    solver::DiffEqBase.AbstractODEAlgorithm=Tsit5(),
    reltol::Real=1e-8,
    abstol::Real=1e-8, # TODO: add more kwargs for integrator
    T=Float32
    )

    n_states, n_actions, _ = state_space_validation(A, B, C, D, Continuous())
    ic_sampler = UniformSampler(s0_lb, s0_ub)

    # Sample random IC for dimension check and to define ODEProblem
    s0 = ic_sampler()
    length(s0) == n_states || throw(ArgumentError("Length $(length(s0)) of s0 doesn't match state dimension $(n_states)"))

    function ode(s, a, t)
        ṡ = A * s + B * a
    end

    problem = ODEProblem(ode, s0, tspan)
    reward_fn = QuadraticReward(Q, R)
    observation_fn = LinearObservation(C, D)

    return DiffEqEnv(problem, reward_fn, n_actions, dt, s0_lb, s0_ub;
        observation_fn=observation_fn,
        o_ub=o_ub, o_lb=o_lb, a_ub=a_ub, a_lb=a_lb,
        solver=solver, reltol=reltol, abstol=abstol,
        T=T
    )
end

"""
Constructs LTI system with quadratic reward and constant initial conditions.
"""
function LTIQuadraticEnv(
    A::AbstractVecOrMat,
    B::AbstractVecOrMat,
    C::AbstractVecOrMat,
    D::AbstractVecOrMat,
    Q::AbstractVecOrMat,
    R::AbstractVecOrMat,
    s0::Union{Real,Vector{<:Real}},
    tspan::Tuple{<:Real,<:Real},
    dt::Real;
    #= Keyword arguments =#
    o_ub::Union{Nothing,Real,Vector{<:Real}}=nothing, # upper bound for observation space
    o_lb::Union{Nothing,Real,Vector{<:Real}}=nothing, # lower bound for observation space
    a_ub::Union{Nothing,Real,Vector{<:Real}}=nothing, # upper bound for action space
    a_lb::Union{Nothing,Real,Vector{<:Real}}=nothing, # lower bound for action space
    solver::DiffEqBase.AbstractODEAlgorithm=Tsit5(),
    reltol::Real=1e-8,
    abstol::Real=1e-8, # TODO: add more kwargs for integrator
    T=Float32
    )

    s0_lb = s0
    s0_ub = s0

    return LTIQuadraticEnv(A, B, C, D, Q, R,
        s0_lb, s0_ub, tspan, dt;
        o_ub=o_ub, o_lb=o_lb, a_ub=a_ub, a_lb=a_lb,
        solver=solver, reltol=reltol, abstol=abstol, T=T)
end


text:
module TestMeans
using GaussianProcesses, Calculus
using Test, Statistics, Random

Random.seed!(1)

@testset "Means" begin
    n = 5 # number of observations
    d = 4 # dimension
    D = 3 # degree of polynomial mean function

    means = [MeanZero(),
             MeanConst(3.0),
             MeanLin(rand(d)),
             MeanPoly(rand(d, D)),
             MeanPeriodic(randn(d), randn(d), randn(d)),
             MeanConst(3.0) * MeanLin(rand(d)),
             MeanLin(rand(d)) + MeanPoly(rand(d, D))]

    X = rand(d, n)
    x = view(X, :, 1)

    @testset "Mean $(typeof(m))" for m in means
        println("\tTesting ", nameof(typeof(m)), "...")

        means = mean(m, X)
        params = GaussianProcesses.get_params(m)

        @testset "Run" begin
            mean(m, x)
            GaussianProcesses.grad_mean(m, x)
            GaussianProcesses.grad_stack(m, X)
            show(devnull, m)
            GaussianProcesses.get_param_names(m)
            set_params!(m, params)
        end

        @testset "Consistency" begin
            @test means ≈ invoke(mean, Tuple{GaussianProcesses.Mean, Matrix{Float64}}, m, X)
        end

        @testset "Gradients" begin
            @testset "Observation #$i" for i in 1:n
                Xi = view(X, :, i)
                theor_grad = GaussianProcesses.grad_mean(m, Xi)
                num_grad = Calculus.gradient(params) do params
                    set_params!(m, params)
                    mean(m, Xi)
                end
                @test theor_grad ≈ num_grad rtol=1e-5 atol=1e-5
            end
        end
    end
end
end


text:
#=

 Test cases for interval exceptions from IEEE Std 1788-2015

 Copyright 2016 Oliver Heimlich (oheim@posteo.de)

 Copying and distribution of this file, with or without modification,
 are permitted in any medium without royalty provided the copyright
 notice and this notice are preserved.  This file is offered as-is,
 without any warranty.

=#
#Language imports

#Test library imports
using Test

#Arithmetic library imports
using IntervalArithmetic

#Preamble
setprecision(53)
setprecision(Interval, Float64)
setrounding(Interval, :tight)
# Set full format, and show decorations
@format full

@testset "exceptions" begin
    @test @interval("[+infinity]") == ∅
    @test Interval(Inf, -Inf) == ∅
    @test @interval("[1.0000000000000001, 1.0000000000000002]") == Interval(1.0, 0x1.0000000000001p+0)
end
# FactCheck.exitstatus()


text:
"""
    drawdowns(pnl; compound)

Calculates the drawdown based on a returns time series.

# Arguments
- `returns`:    Vector of asset returns (usually log-returns).
- `geometric`:  Use geometric compounding (cumprod) if set to true, otherwise simple arithmetic sum (cumsum), e.g. for log-returns (default).
"""
function drawdowns(returns; geometric::Bool=false)
    if geometric
        returns = cumprod(1.0 .+ returns)
    else
        returns = 1.0 .+ cumsum(returns)
    end

    # cumulative returns needs to start at 1.0 for cumulative max to be correct
    returns_max = accumulate(max, [1.0; returns])[2:end]
    returns ./ returns_max .- 1.0
end

"""
    drawdowns_pnl(pnl)

Calculates the drawdown based on a Profit-and-Loss (PnL) time series, e.g. daily equity changes in USD.

# Arguments
- `pnl`:    Vector of Profit-and-Loss (PnL) values, e.g. daily equity changes in USD.
"""
function drawdowns_pnl(pnl)
    pnl = cumsum(pnl)

    # cumulative PnL needs to start at 0.0 for cumulative max to be correct
    pnl_max = accumulate(max, [0.0; pnl])[2:end]
    pnl .- pnl_max
end


text: include("universal_functions_parameters.jl")

max_stars_repo_name: climate-machine/CLIMAParameters

max_issues_repo_name: climate-machine/CLIMAParameters

max_forks_repo_name: climate-machine/CLIMAParameters


text:
module LDBPlots

export plot_emap,
       plot_emaps

using Plots, Wavelets, WaveletsExt

function plot_emap(emap::AbstractArray{T,2}, tree::BitVector;
                   clim::Tuple{T,T} = extrema(emap)) where T<:Real
    start = 0                           # start index
    ny, nx = size(emap)                 # emap dimensions
    emap = emap'

    plt = heatmap(start:(ny+start-1),             # x-axis labels
                  0:(nx-1),                       # y-axis labels
                  emap,                           # heatmap values
                  c = :matter,                    # color gradient
                  background_color = :white,
                  legend = false,
                  clim = clim)                    # color bar limit
    # Add tfbdry lines on top
    plot_tfbdry!(plt, tree, node_color = :transparent, line_color = :black,
                 background_color = :transparent)
    return plt
end

function plot_emap(emap::AbstractArray{T,3}, tree::BitVector;
                   clim::Tuple{T,T} = extrema(emap)) where T<:Real
    @assert isvalidtree(emap[:,:,1], tree)
    nx, ny, _ = size(emap)
    data = getbasiscoef(emap, tree)             # get basis coefficients

    plt = heatmap(0:(ny-1),                     # x-axis labels
                  0:(nx-1),                     # y-axis labels
                  data,                         # heatmap values
                  c = :matter,                  # color gradient
                  background_color = :white,
                  legend = false,
                  clim = clim)                  # color bar limit
    plot_tfbdry2!(plt, tree, line_color = :black, background_color = :transparent)
    return plt
end

function plot_emaps(emap::AbstractArray{T}, tree::BitVector;
                    clim::Tuple{T,T} = extrema(emap)) where T<:Real
    @assert 3 ≤ ndims(emap) ≤ 4
    test_signal = ndims(emap) == 3 ? emap[:,1,1] : emap[:,:,1,1]
    @assert isvalidtree(test_signal, tree)            # Check if tree is valid for 1D signals
    plts = []                                         # Collect all plots
    # Plot energy map for each class
    for (emapᵢ, i) in zip(eachslice(emap, dims=ndims(emap)), axes(emap,ndims(emap)))
        plt = plot_emap(emapᵢ, tree, clim = clim)
        plot!(plt, ylabel = "Class $i")
        push!(plts, plt)
    end
    plts = plot(plts..., layout = (size(emap,ndims(emap)), 1))  # Display all plots
    return plts
end

end # module


text:
Base.@kwdef mutable struct Axis
    isenabled::Bool = true
    thickness::Cfloat = Cfloat(1)
end

function isenabled(axis::Axis)
    axis.isenabled
end

function enable!(axis::Axis)
    axis.isenabled = true
end

function disable!(axis::Axis)
    axis.isenabled = false
end

function get_thickness(axis::Axis)
    axis.thickness
end

function set_thickness!(axis::Axis, thickness::Cfloat)
    axis.thickness = thickness
end


text:
module SpringCollab2020CustomizableGlassBlock

using ..Ahorn, Maple

@mapdef Entity "SpringCollab2020/CustomizableGlassBlock" CustomizableGlassBlock(x::Integer, y::Integer, width::Integer=Maple.defaultBlockWidth, height::Integer=Maple.defaultBlockHeight, behindFgTiles::Bool=false)

const placements = Ahorn.PlacementDict(
    "Customizable Glass Block (Spring Collab 2020)" => Ahorn.EntityPlacement(
        CustomizableGlassBlock,
        "rectangle"
    )
)

Ahorn.minimumSize(entity::CustomizableGlassBlock) = 8, 8
Ahorn.resizable(entity::CustomizableGlassBlock) = true, true

Ahorn.selection(entity::CustomizableGlassBlock) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::CustomizableGlassBlock, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (1.0, 1.0, 1.0, 0.5), (1.0, 1.0, 1.0, 0.5))
end

end


text:
module BellComm

# include external modules
include("./ConvexPolytope.jl")

# include internal modules
include("./LocalPolytope.jl")
include("./Behavior.jl")
include("./Degeneracy.jl")
include("./PrepareAndMeasure.jl")
include("./DichotomicLocalPolytope.jl")
include("./QuantumBehavior.jl")
include("./QuantumOpt.jl")

export QBase, ConvexPolytope, Degeneracy, LocalPolytope, Behavior, QuantumBehavior, QuantumOpt, PrepareAndMeasure

# function local_polytope()
#
# end
#
# function local_vertices()
# end
#
# function local_facets()
#
# end
#
# function optimize_quantum_nonlocality()
#
# end

end