text:
# Device type and auxiliary functions

export
    CuDevice, name, totalmem, attribute

const CuDevice_t = Cint

"""
    CuDevice(i::Integer)

Get a handle to a compute device.
"""
struct CuDevice
    handle::CuDevice_t

    # CuDevice is just an integer, but we need (?) to call cuDeviceGet to make sure this
    # integer is valid. to avoid ambiguity, add a bogus argument (cfr. `checkbounds`)
    CuDevice(::Type{Bool}, handle::CuDevice_t) = new(handle)
end

function CuDevice(ordinal::Integer)
    device_ref = Ref{CuDevice_t}()
    @apicall(:cuDeviceGet, (Ptr{CuDevice_t}, Cint), device_ref, ordinal)
    CuDevice(Bool, device_ref[])
end

Base.convert(::Type{CuDevice_t}, dev::CuDevice) = dev.handle

Base.:(==)(a::CuDevice, b::CuDevice) = a.handle == b.handle
Base.hash(dev::CuDevice, h::UInt) = hash(dev.handle, h)

Base.show(io::IO, dev::CuDevice) = print(io, "CuDevice($(dev.handle)): $(name(dev))")

"""
    name(dev::CuDevice)

Returns an identifier string for the device.
"""
function name(dev::CuDevice)
    buflen = 256
    buf = Vector{Cchar}(undef, buflen)
    @apicall(:cuDeviceGetName, (Ptr{Cchar}, Cint, CuDevice_t),
                               buf, buflen, dev)
    buf[end] = 0
    return unsafe_string(pointer(buf))
end

"""
    totalmem(dev::CuDevice)

Returns the total amount of memory (in bytes) on the device.
"""
function totalmem(dev::CuDevice)
    mem_ref = Ref{Csize_t}()
    @apicall(:cuDeviceTotalMem, (Ptr{Csize_t}, CuDevice_t), mem_ref, dev)
    return mem_ref[]
end

@enum(CUdevice_attribute, MAX_THREADS_PER_BLOCK = Cint(1),
                          MAX_BLOCK_DIM_X,
                          MAX_BLOCK_DIM_Y,
                          MAX_BLOCK_DIM_Z,
                          MAX_GRID_DIM_X,
                          MAX_GRID_DIM_Y,
                          MAX_GRID_DIM_Z,
                          MAX_SHARED_MEMORY_PER_BLOCK,
                          TOTAL_CONSTANT_MEMORY,
                          WARP_SIZE,
                          MAX_PITCH,
                          MAX_REGISTERS_PER_BLOCK,
                          CLOCK_RATE,
                          TEXTURE_ALIGNMENT,
                          GPU_OVERLAP,
                          MULTIPROCESSOR_COUNT,
                          KERNEL_EXEC_TIMEOUT,
                          INTEGRATED,
                          CAN_MAP_HOST_MEMORY,
                          COMPUTE_MODE,
                          MAXIMUM_TEXTURE1D_WIDTH,
                          MAXIMUM_TEXTURE2D_WIDTH,
                          MAXIMUM_TEXTURE2D_HEIGHT,
                          MAXIMUM_TEXTURE3D_WIDTH,
                          MAXIMUM_TEXTURE3D_HEIGHT,
                          MAXIMUM_TEXTURE3D_DEPTH,
                          MAXIMUM_TEXTURE2D_LAYERED_WIDTH,
                          MAXIMUM_TEXTURE2D_LAYERED_HEIGHT,
                          MAXIMUM_TEXTURE2D_LAYERED_LAYERS,
                          SURFACE_ALIGNMENT,
                          CONCURRENT_KERNELS,
                          ECC_ENABLED,
                          PCI_BUS_ID,
                          PCI_DEVICE_ID,
                          TCC_DRIVER,
                          MEMORY_CLOCK_RATE,
                          GLOBAL_MEMORY_BUS_WIDTH,
                          L2_CACHE_SIZE,
                          MAX_THREADS_PER_MULTIPROCESSOR,
                          ASYNC_ENGINE_COUNT,
                          UNIFIED_ADDRESSING,
                          MAXIMUM_TEXTURE1D_LAYERED_WIDTH,
                          MAXIMUM_TEXTURE1D_LAYERED_LAYERS,
                          CAN_TEX2D_GATHER,
                          MAXIMUM_TEXTURE2D_GATHER_WIDTH,
                          MAXIMUM_TEXTURE2D_GATHER_HEIGHT,
                          MAXIMUM_TEXTURE3D_WIDTH_ALTERNATE,
                          MAXIMUM_TEXTURE3D_HEIGHT_ALTERNATE,
                          MAXIMUM_TEXTURE3D_DEPTH_ALTERNATE,
                          PCI_DOMAIN_ID,
                          TEXTURE_PITCH_ALIGNMENT,
                          MAXIMUM_TEXTURECUBEMAP_WIDTH,
                          MAXIMUM_TEXTURECUBEMAP_LAYERED_WIDTH,
                          MAXIMUM_TEXTURECUBEMAP_LAYERED_LAYERS,
                          MAXIMUM_SURFACE1D_WIDTH,
                          MAXIMUM_SURFACE2D_WIDTH,
                          MAXIMUM_SURFACE2D_HEIGHT,
                          MAXIMUM_SURFACE3D_WIDTH,
                          MAXIMUM_SURFACE3D_HEIGHT,
                          MAXIMUM_SURFACE3D_DEPTH,
                          MAXIMUM_SURFACE1D_LAYERED_WIDTH,
                          MAXIMUM_SURFACE1D_LAYERED_LAYERS,
                          MAXIMUM_SURFACE2D_LAYERED_WIDTH,
                          MAXIMUM_SURFACE2D_LAYERED_HEIGHT,
                          MAXIMUM_SURFACE2D_LAYERED_LAYERS,
                          MAXIMUM_SURFACECUBEMAP_WIDTH,
                          MAXIMUM_SURFACECUBEMAP_LAYERED_WIDTH,
                          MAXIMUM_SURFACECUBEMAP_LAYERED_LAYERS,
                          MAXIMUM_TEXTURE1D_LINEAR_WIDTH,
                          MAXIMUM_TEXTURE2D_LINEAR_WIDTH,
                          MAXIMUM_TEXTURE2D_LINEAR_HEIGHT,
                          MAXIMUM_TEXTURE2D_LINEAR_PITCH,
                          MAXIMUM_TEXTURE2D_MIPMAPPED_WIDTH,
                          MAXIMUM_TEXTURE2D_MIPMAPPED_HEIGHT,
                          COMPUTE_CAPABILITY_MAJOR,
                          COMPUTE_CAPABILITY_MINOR,
                          MAXIMUM_TEXTURE1D_MIPMAPPED_WIDTH,
                          STREAM_PRIORITIES_SUPPORTED,
                          GLOBAL_L1_CACHE_SUPPORTED,
                          LOCAL_L1_CACHE_SUPPORTED,
                          MAX_SHARED_MEMORY_PER_MULTIPROCESSOR,
                          MAX_REGISTERS_PER_MULTIPROCESSOR,
                          MANAGED_MEMORY,
                          MULTI_GPU_BOARD,
                          MULTI_GPU_BOARD_GROUP_ID)
@assert Cint(MULTI_GPU_BOARD_GROUP_ID) == 85

"""
    attribute(dev::CuDevice, code)

Returns information about the device.
"""
function attribute(dev::CuDevice, code::CUdevice_attribute)
    value_ref = Ref{Cint}()
    @apicall(:cuDeviceGetAttribute, (Ptr{Cint}, Cint, CuDevice_t),
                                    value_ref, code, dev)
    return value_ref[]
end

## device iteration

export devices

struct DeviceSet end

"""
    devices()

Get an iterator for the compute devices.
"""
devices() = DeviceSet()

Base.eltype(::DeviceSet) = CuDevice

function Base.iterate(iter::DeviceSet, i=1)
    i >= length(iter) + 1 ? nothing : (CuDevice(i-1), i+1)
end

function Base.length(::DeviceSet)
    count_ref = Ref{Cint}()
    @apicall(:cuDeviceGetCount, (Ptr{Cint},), count_ref)
    return count_ref[]
end

Base.IteratorSize(::DeviceSet) = Base.HasLength()

## convenience attribute getters

export warpsize, capability

"""
    warpsize(dev::CuDevice)

Returns the warp size (in threads) of the device.
"""
warpsize(dev::CuDevice) = attribute(dev, WARP_SIZE)

"""
    capability(dev::CuDevice)

Returns the compute capability of the device.
"""
function capability(dev::CuDevice)
    return VersionNumber(attribute(dev, COMPUTE_CAPABILITY_MAJOR),
                         attribute(dev, COMPUTE_CAPABILITY_MINOR))
end


text:
using NeuralVerification
using Test

macro no_error(ex)
    quote
        try $(esc(ex))
            true
        catch
            false
        end
    end
end

at = @__DIR__

small_nnet = read_nnet("$at/../examples/networks/small_nnet.txt")
A = Matrix{Float64}(undef, 2,1)
A[1:2] = [1, -1]

# TODO can we unify all of the reachability test conditions?
### reverify
inputSet  = HPolytope(A, [1.0, 1.0])
outputSet = HPolytope(A, [1.0, 1.0])
problem_reverify = Problem(small_nnet, inputSet, outputSet)
solver_reverify = NSVerify(GLPKSolverMIP(), 1000.0)

### maxSens
inputSet  = HPolytope(A, [1.0, 1.0])
outputSet = HPolytope(A, [100.0, 1.0])
problem_maxSens = Problem(small_nnet, inputSet, outputSet)
solver_maxSens = MaxSens(0.3)

### exactReach
inputSet  = HPolytope(A,           [0.0, 2.0])
outputSet = HPolytope(zeros(2, 1), [100.0, 1.0])
problem_exactReach = Problem(small_nnet, inputSet, outputSet)
solver_exactReach = ExactReach()

#reluVal
inputSet  = Hyperrectangle(low = [-1.0], high = [1.0])
outputSet = Hyperrectangle(low = [-1.0], high = [50.0])
problem_reluVal = Problem(small_nnet, inputSet, outputSet)
solver_reluVal = ReluVal(2)

### Ai2
input  = HPolytope(A, [0.0, 0.0])
output = HPolytope(A, [100.0, 0.0])
problem_ai2 = Problem(small_nnet, input, output)
solver_ai2 = Ai2()

### BaB
#inputSet  = Hyperrectangle([-1.0], [0.5])
#outputSet = HPolytope(ones(1,1), [18.5])
#problem_bab = Problem(small_nnet, inputSet, outputSet)
#solver_bab = BaB(0.1, GLPKSolverMIP())

### Certify
inputSet = Hyperrectangle([1.0], [1.0])
outputSet = HPolytope(ones(1,1), [2.1])
problem_certify = Problem(small_nnet, inputSet, outputSet)
solver_certify = Certify(SCSSolver())

### ConvDual
inputSet  = Hyperrectangle([0.0], [1.0])
outputSet = HPolytope(-ones(1,1), [1.0])
problem_convdual = Problem(small_nnet, inputSet, outputSet)
solver_convdual = ConvDual()

### Duality
inputSet  = Hyperrectangle([1.0], [1.0])
outputSet = HPolytope(ones(1,1), [120.0])
problem_duality = Problem(small_nnet, inputSet, outputSet)
solver_duality = Duality(GLPKSolverMIP())

### iLP
inputSet  = Hyperrectangle([-2.0], [.5])
outputSet = HPolytope(ones(1,1),[72.5])
problem_ilp = Problem(small_nnet, inputSet, outputSet)
solver_ilp = ILP(GLPKSolverMIP(), 1)

### MIPVerify
inputSet  = Hyperrectangle([0.0], [.5])
outputSet = HPolytope(ones(1,1), [102.5])
problem_mipverify = Problem(small_nnet, inputSet, outputSet)
solver_mipverify = MIPVerify(GLPKSolverMIP())

### Planet
inputSet  = Hyperrectangle([-1.0], [0.5])
outputSet = HPolytope(ones(1,1), [2.5])
problem_planet = Problem(small_nnet, inputSet, outputSet)
solver_planet = Planet(GLPKSolverMIP())

### Reluplex
inputSet  = Hyperrectangle([1.0],[.2])
outputSet = Hyperrectangle([2.0,], [100.0])
problem_reluplex = Problem(small_nnet, inputSet, outputSet)
solver_reluplex = Reluplex()

### Sherlock
inputSet  = Hyperrectangle([2.0], [.5])
outputSet = Hyperrectangle([10.0], [10.0])
problem_sherlock = Problem(small_nnet, inputSet, outputSet)
solver_sherlock = Sherlock(GLPKSolverMIP(), 0.1)

#@test @no_error solve(solver_reverify,   problem_reverify)
#@test @no_error solve(solver_maxSens,    problem_maxSens)
@test @no_error solve(solver_exactReach, problem_exactReach)
@test @no_error solve(solver_reluVal,    problem_reluVal)
@test @no_error solve(solver_ai2,        problem_ai2)
@test @no_error solve(solver_bab,        problem_bab)
@test @no_error solve(solver_certify,    problem_certify)
@test @no_error solve(solver_convdual,   problem_convdual)
@test @no_error solve(solver_duality,    problem_duality)
@test @no_error solve(solver_ilp,        problem_ilp)
@test @no_error solve(solver_mipverify,  problem_mipverify)
@test @no_error solve(solver_planet,     problem_planet)
@test @no_error solve(solver_reluplex,   problem_reluplex)
@test @no_error solve(solver_sherlock,   problem_sherlock)

ext: jl


text:
# function get_surface_grid_fixed_brick_size(segmentation::AbstractArray, block_size::ArrayOrTuple; return_all::Bool=false)
#
#     Fchar = __grid_get_surface_Fchar_per_block(segmentation, block_size)
#     #
#     # data_size = size_as_array(size(segmentation))
#     # bigV, FVreduced = grid_Fchar_to_V_FVreduced(Fchar, data_size)
#     #
#     # # return filtered_bigFV, Fchar, (bigV, model)
#     #
#     # if return_all
#     #     return (bigV,[FVreduced]), Fchar, (bigV, [FVreduced])
#     # else
#     #     return (bigV,[FVreduced])
#     # end
# end


text:
module Pendulum2

using  Modia

Pendulum = Model(
    world = Object3D(feature=Scene(animationFile="Pendulum2.json")),
    obj1  = Object3D(feature=Solid(shape=Box(lengthX=1.0, lengthY=0.2, lengthZ=0.2),
                solidMaterial="Steel", visualMaterial=VisualMaterial(color="Blue"))),
    obj2  = Object3D(parent=:obj1, feature=Visual(shape=Cylinder(diameter=0.1, length=0.21),
                visualMaterial=VisualMaterial(color="Red")), translation=[-0.5, 0.0, 0.0]),
    rev   = Revolute(obj1=:world, obj2=:obj2)
)

pendulum = @instantiateModel(buildModia3D(Pendulum), unitless=true)
simulate!(pendulum, stopTime=3.0)

@usingModiaPlot
plot(pendulum, "rev.phi")
end


text:
# given an 1D interval, return quadrature nodes and weights

using LinearAlgebra
using FastGaussQuadrature

function quadgrid(xmin, xmax, n)
    # get nodes and weights in reference coord [-1, 1]
    x, w = gausslegendre(n)
    # map nodes and weights to current interval
    dist = (xmax - xmin)/2
    center = (xmax + xmin)/2
    weights = dist .* w
    nodes = dist .* x .+ center
    return nodes, weights
end

function quadgrid_old(xmin, xmax, n)
    dist = (xmax - xmin)/2
    center = (xmax + xmin)/2
    if n == 1
        w = 2
        x = 0
    elseif n == 2
        w = [1, 1]
        x =  1/sqrt(3) * [1, -1]
    elseif n == 3
        w = [5/9, 8/9, 5/9]
        x = -sqrt(3/5) * [-1, 0, 1]
    elseif n == 4
        w = [(18+sqrt(30))/36, (18+sqrt(30))/36,
             (18-sqrt(30))/36, (18-sqrt(30))/36]
        x = [sqrt(3/7 - 2/7*sqrt(6/5)), -sqrt(3/7 - 2/7*sqrt(6/5)),
             sqrt(3/7 + 2/7*sqrt(6/5)), -sqrt(3/7 + 2/7*sqrt(6/5))]
    elseif n == 5
        w = [128/225, (322+13*sqrt(70))/900, (322+13*sqrt(70))/900,
                      (322-13*sqrt(70))/900, (322-13*sqrt(70))/900]
        x = [0, sqrt(5-2*sqrt(10/7))/3, -sqrt(5-2*sqrt(10/7))/3,
                sqrt(5+2*sqrt(10/7))/3, -sqrt(5+2*sqrt(10/7))/3]
    end

    w_new = dist .* w
    x_new = dist .* x .+ center
    # fx = f.(x_new)
    # int = dot(w_new, fx)

    return x_new, w_new

end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module CLIMA

include("Utilities/ParametersType/src/ParametersType.jl")
include("Utilities/PlanetParameters/src/PlanetParameters.jl")
include("Utilities/RootSolvers/src/RootSolvers.jl")
include("Utilities/MoistThermodynamics/src/MoistThermodynamics.jl")
include("Atmos/Parameterizations/SurfaceFluxes/src/SurfaceFluxes.jl")
include("Atmos/Parameterizations/TurbulenceConvection/src/TurbulenceConvection.jl")
include("Mesh/Topologies.jl")
include("Mesh/Grids.jl")
include("Arrays/MPIStateArrays.jl")
include("Atmos/Dycore/src/AtmosDycore.jl")
include("DGmethods/SpaceMethods.jl")
include("DGmethods/DGBalanceLawDiscretizations.jl")
include("ODESolvers/ODESolvers.jl")
include("ODESolvers/LowStorageRungeKuttaMethod.jl")
include("ODESolvers/GenericCallbacks.jl")

end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
###############################################################
# Sample that caches properties returned by an AbstractSample #
###############################################################

mutable struct Sample
    # The sample itself
    value::Any
    # Properties of the sample
    properties::Dict{Symbol,Any}

    function Sample(
        value::Any = nothing;
        properties::Dict{Symbol,Any} = Dict{Symbol,Any}(),
        sample_rate=get_job().sample_rate,
        total_memory_usage=nothing
    )
        newsample = new(value, properties)

        # We compute the `memory_usage` lazily

        # Fill in properties if possible
        if !isnothing(total_memory_usage)
            setsample!(newsample, :memory_usage, round(total_memory_usage / sample_rate))
        end
        setsample!(newsample, :rate, sample_rate)

        newsample
    end
    # TODO: Un-comment if needed
    # Sample(sample::Sample, properties::Vector{String}) =
    #     new(sample.value, sample.stale, Dict(
    #         sample.sample_properties[prop]
    #         for prop in properties
    #     ))
end

ExactSample(value::Any = nothing; kwargs...) = Sample(value; sample_rate=1, kwargs...)

# sample* functions always return a concrete value or a dict with properties.
# To get a `Sample`, access the property.

# TODO: Lazily compute samples by storing sample computation in a DAG if its
# getting too expensive
sample(fut::AbstractFuture) = sample(get_location(fut))
sample(sample::Sample) = sample.value

sample(fut::AbstractFuture, propertykeys...) = sample(get_location(fut).sample, propertykeys...)
function sample(s::Sample, propertykeys...)
    properties = s.properties
    for (i, propertykey) in enumerate(propertykeys)
        properties = get!(
            properties,
            propertykey,
            if i < length(propertykeys)
                Dict()
            else
                sample(s.value, propertykeys...)
            end
        )
    end
    properties
end

setsample!(fut::AbstractFuture, value) = setsample!(get_location(fut).sample, value)
function setsample!(sample::Sample, value)
    sample.value = value
end

setsample!(fut::AbstractFuture, propertykeys...) = setsample!(get_location(fut).sample, propertykeys...)
function setsample!(sample::Sample, propertykeys...)
    if length(propertykeys) == 1
        setsample!(sample, first(propertykeys))
    else
        properties = sample.properties
        propertyvalue = last(propertykeys)
        propertykeys = propertykeys[1:end-1]
        for (i, propertykey) in enumerate(propertykeys)
            if i < length(propertykeys)
                properties = get!(properties, propertykey, Dict())
            end
        end
        properties[last(propertykeys)] = propertyvalue
    end
end

# TODO: For futures with locations like Size, "scale up" the computed sample
# for a useful approximation of things like length of an array

####################################################################
# AbstractSample to be implemented by anything that can be sampled #
####################################################################

# NOTE: We use strings for things that will be serialized to JSON and symbols
# for everything else

# NOTE: We use upper-camel-case for user-facing names (like names of PTs) and
# all-caps-snake-case for anything internal (like names of constraints)

# The purpose of the `sample` function is to allow for computing various
# properties of the sample by property key name instead of an explicit
# function call. This makes it easier for the `sample` and `setsample!`
# functions for `Future`s to compute and cache samples.

sample(as::Any, properties...) =
    if length(properties) <= 2 && first(properties) == :statistics
        # If the statistic is not cached in `Sample.properties`, then we just
        # return an empty dictionary
        Dict()
    elseif length(properties) == 1
        if first(properties) == :memory_usage
            sample_memory_usage(as)
        elseif first(properties) == :rate
            # This is the default but the `Sample` constructor overrides this
            # before-hand to allow some samples to be "exact" with a sample
            # rate of 1
            get_job().sample_rate
        elseif first(properties) == :keys
            sample_keys(as)
        elseif first(properties) == :axes
            sample_axes(as)
        elseif first(properties) == :groupingkeys
            # This is just the initial value for grouping keys. Calls to
            # `keep_*` functions will expand it.
            []
        else
            # println(typeof(as))
            # println(typeof(as) <: Any)
            throw(ArgumentError("Invalid sample properties: $properties"))
        end
    elseif length(properties) == 3 && first(properties) == :statistics
        key = properties[2]
        query = properties[3]
        if query == :max_ngroups
            sample_max_ngroups(as, key)
        elseif query == :divisions
            sample_divisions(as, key)
        elseif query == :min
            sample_min(as, key)
        elseif query == :max
            sample_max(as, key)
        else
            throw(ArgumentError("Invalid sample properties: $properties"))
        end
    elseif length(properties) == 5 && first(properties) == :statistics
        key, query, minvalue, maxvalue = properties[2:end]
        if query == :percentile
            sample_percentile(as, key, minvalue, maxvalue)
        else
            throw(ArgumentError("Invalid sample properties: $properties"))
        end
    else
        throw(ArgumentError("Invalid sample properties: $properties"))
    end

sample_memory_usage(as::Any) = total_memory_usage(as)

# Implementation error
impl_error(fn_name, as) = error("$fn_name not implemented for $(typeof(as))")

# Functions to implement for Any (e.g., for DataFrame or
# Array)
sample_axes(as::Any) = impl_error("sample_axes", as)
sample_keys(as::Any) = impl_error("sample_keys", as)
sample_divisions(as::Any, key) = impl_error("sample_divisions", as)
sample_percentile(as::Any, key, minvalue, maxvalue) = impl_error("sample_percentile", as)
sample_max_ngroups(as::Any, key) = impl_error("sample_max_ngroups", as)
sample_min(as::Any, key) = impl_error("sample_min", as)
sample_max(as::Any, key) = impl_error("sample_max", as)

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
#Dummy module used to install C-Data suite
module CData
using Pkg

parsestring(::Type{String}, s::AbstractString) = s
parsestring(::Type{T}, s::AbstractString) where T = parse(T, s)

function input(::Type{T}, prompt::String, default) where T
	println(prompt)
	print("[$default]: ")
	result = parsestring(T, strip(readline(stdin)))
	if "" == result
		result = default
	end
	return result
end

function isinstalled(modname::String)
	return in(modname, keys(Pkg.installed()))
end

#Install module using Pkg.clone(), if not installed:
function cond_install_clone(path::String, modname::String; jl_suffix::Bool = true)
	suffix = jl_suffix ? ".jl.git" : ""
	println("\nInstalling module $modname...")

	if isinstalled(modname)
		println("        already present.")
	else
		pkgpath = joinpath(path, "$modname$suffix")
		Pkg.add(PackageSpec(url=pkgpath))
		println("done.")
	end
end

#Vectorize cond_install_clone:
function cond_install_clone(path::String, v::Vector; jl_suffix::Bool = true)
	for modname in v
		cond_install_clone(path, modname, jl_suffix=jl_suffix)
	end
end

function install()
	path_laforge = "git://github.com/ma-laforge" #WANTCONST

	modlist = [
		"FileIO2", "MDDatasets", "NumericIO",
		"InspectDR", "GracePlot",
		"EasyPlot", "EasyData", "EasyPlotInspect", "EasyPlotGrace", #"EasyPlotPlots", #Not working
		"CircuitAnalysis", "SignalProcessing", "NetwAnalysis",
      "SpiceData", "LibPSF", "PSFWrite", "EDAData",
		"CppSimData", "LibPSFC", #NOTE: install CppSimData/LibPSFC for sample data files only.
	]
	cond_install_clone(path_laforge, modlist)

	if "Y" == uppercase(input(String, "Install PyPlot/PyCall dependencies (y/n)?", "y"))
		modlist = [
			"EasyPlotMPL", #"EasyPlotQwt" #Not working
		]
		cond_install_clone(path_laforge, modlist)
	end

	println("\n\nPackage installation complete.")
	println("Refresh out-of-date packages using \"Pkg.update()\".")
end

end #module CData

max_stars_repo_name: ma-laforge/CData

max_issues_repo_name: ma-laforge/CData

max_forks_repo_name: ma-laforge/CData


text:
using Base.Test
cmp{T}(x::T,y::T) = all(Bool[getfield(x,fname) == getfield(y,fname) for fname in fieldnames(x)])

#@perpulse a,b,c,d = do_filter(filter, pulse_record)
#@threshold calibration = calibrate_nofit(filt_value_dc_hist, known_energies) when counted(filt_value_dc_hist)>3000
#macroexpand( quote @histogram update_histogram!(filt_value_hist, selection_good, filt_value) end)
@test cmp(PerPulseStep(sum,[:filter,:pulse_record],[:filt_value]), @perpulse filt_value = sum(filter, pulse_record))
@test cmp(ThresholdStep(calibrate_nofit, [:filt_value_dc_hist,:known_energies],[:calibration],:filt_value_dc_hist, counted, 3000, true),
@threshold calibration = calibrate_nofit(filt_value_dc_hist, known_energies) when counted(filt_value_dc_hist)>3000 )
@test cmp(HistogramStep(update_histogram!, [:selection_good, :filt_value], [:filt_value_hist]),
@histogram update_histogram!(filt_value_hist, selection_good, filt_value))


text: foo() = @fact 10

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
module StatsLH

# using Base.Test
using StatsBase

export std_w, std_std, density_weighted, discretize

include("stats/std_w.jl")
include("stats/std_std.jl")
include("stats/density_weighted.jl")
include("stats/discretize.jl")

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule SymEngine_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("SymEngine")
JLLWrappers.@generate_main_file("SymEngine", UUID("3428059b-622b-5399-b16f-d347a77089a4"))
end  # module SymEngine_jll


text:
using GridapFSI
using DifferentialEquations

function get_cylinder_motion()
  # Parameters
  A₀ = 1.5
  A₁ = 8.5
  A₂ = -11.1
  A₃ = 2.6
  B₀ = 4.2
  B₁ = 11.3
  B₂ = -68.7
  B₃ = 50.5

  # Physics Properties
  k = 1.0
  m = 1.0
  D = 0.1
  L = 1.0
  ρ = 1.0e3
  b = 1.0
  V = 1.0

  # Parameters
  q₁ = 2
  St = 0.2
  Cx0 = 2.0
  Cy0 = 0.3
  A = 12
  ε = 0.3

  # Auxiliar variables
  mₐ = π*ρ*D^2*L/4
  mʳ = m/mₐ
  ωₙ = sqrt(k/(m+mₐ))
  ωₛ = St*V/(2*π*D)
  Ωₙ = ωₙ/ωₛ
  ζ = b/(2*sqrt((m+mₐ)*k))

  function oscillator(x,p,t)
    y, q, dy, dq = x
    y_dot = dy
    q_dot = dq
    Cvy= (-2^π*St*dy*Cx0 + Cy0/q₁*q) * sqrt(1+4*π^2*St^2*dy^2)
    dy_dot = ρ*D^2*L/(m+mₐ)*1.0/(8*π^2*St^2)*Cvy - 2*ζ*Ωₙ*dy - Ωₙ^2*y
    dq_dot = A*dy_dot - ε*(q^2-1)*dq - q
    dx = [y_dot; q_dot; dy_dot; dq_dot]
  end

  prob = ODEProblem(oscillator,[0.0;2.0;0.0;0.0],(0.0,100.0))
  sol=DifferentialEquations.solve(prob)
  y(t) = sol(t)[1]
  dy(t) = sol(t)[3]

  return y,dy
end

# Define cylinder motion
u, du = get_cylinder_motion()

output = main(problemName="oscillator",dt=0.025,tf=10.0,u_y=u, du_y=du)


text:
module TestConvert
using Compat
using Compat.Test
using CategoricalArrays
using CategoricalArrays: DefaultRefType, level, reftype, leveltype, catvalue, iscatvalue, CatValue

@testset "convert() for CategoricalPool{Int, DefaultRefType} and values" begin
    pool = CategoricalPool([1, 2, 3])
    @test convert(CategoricalPool{Int, DefaultRefType}, pool) === pool
    @test convert(CategoricalPool{Int}, pool) === pool
    @test convert(CategoricalPool, pool) === pool
    convert(CategoricalPool{Float64, UInt8}, pool)
    convert(CategoricalPool{Float64}, pool)
    convert(CategoricalPool, pool)

    v1 = catvalue(1, pool)
    v2 = catvalue(2, pool)
    v3 = catvalue(3, pool)
    @test iscatvalue(v1)
    @test iscatvalue(typeof(v1))
    @test eltype(v1) === Any
    @test eltype(typeof(v1)) === Any
    @test leveltype(v1) === Int
    @test leveltype(typeof(v1)) === Int
    @test reftype(v1) === DefaultRefType
    @test reftype(typeof(v1)) === DefaultRefType
    @test v1 isa CategoricalArrays.CategoricalValue{Int, DefaultRefType}

    @test convert(Int, v1) === 1
    @test convert(Int, v2) === 2
    @test convert(Int, v3) === 3

    @test convert(Int32, v1) === Int32(1)
    @test convert(Int32, v2) === Int32(2)
    @test convert(Int32, v3) === Int32(3)

    @test convert(UInt8, v1) === 0x01
    @test convert(UInt8, v2) === 0x02
    @test convert(UInt8, v3) === 0x03

    @test convert(CategoricalValue, v1) === v1
    @test convert(CategoricalValue{Int}, v1) === v1
    @test convert(CategoricalValue{Int, DefaultRefType}, v1) === v1

    @test convert(Any, v1) === v1
    @test convert(Any, v2) === v2
    @test convert(Any, v3) === v3

    for T in (typeof(v1), CatValue, CategoricalValue), U in (Missing, Nothing)
        @test convert(Union{T, U}, v1) === v1
        @test convert(Union{T, U}, v2) === v2
        @test convert(Union{T, U}, v3) === v3
    end

    @test get(v1) === 1
    @test get(v2) === 2
    @test get(v3) === 3

    @test promote(1, v1) === (1, 1)
    @test promote(1.0, v1) === (1.0, 1.0)
    @test promote(0x1, v1) === (1, 1)
end

@testset "convert() for CategoricalPool{String, DefaultRefType} and values" begin
    pool = CategoricalPool(["a", "b", "c"])
    @test convert(CategoricalPool{String, DefaultRefType}, pool) === pool
    @test convert(CategoricalPool{String}, pool) === pool
    @test convert(CategoricalPool, pool) === pool
    convert(CategoricalPool{String, UInt8}, pool)

    v1 = catvalue(1, pool)
    v2 = catvalue(2, pool)
    v3 = catvalue(3, pool)
    @test iscatvalue(v1)
    @test iscatvalue(typeof(v1))
    @test eltype(v1) === Char
    @test eltype(typeof(v1)) === Char
    @test leveltype(v1) === String
    @test leveltype(typeof(v1)) === String
    @test reftype(v1) === DefaultRefType
    @test reftype(typeof(v1)) === DefaultRefType
    @test v1 isa CategoricalArrays.CategoricalString{DefaultRefType}

    @test convert(String, v1) == "a"
    @test convert(String, v2) == "b"
    @test convert(String, v3) == "c"

    @test convert(AbstractString, v1) == "a"
    @test convert(AbstractString, v2) == "b"
    @test convert(AbstractString, v3) == "c"

    @test convert(CategoricalString, v1) === v1
    @test convert(CategoricalString{DefaultRefType}, v1) === v1

    @test convert(Any, v1) === v1
    @test convert(Any, v2) === v2
    @test convert(Any, v3) === v3

    for T in (typeof(v1), CatValue, CategoricalString, CategoricalString{DefaultRefType}),
        U in (Missing, Nothing)
        @test convert(Union{T, U}, v1) === v1
        @test convert(Union{T, U}, v2) === v2
        @test convert(Union{T, U}, v3) === v3
    end

    @test get(v1) === "a"
    @test get(v2) === "b"
    @test get(v3) === "c"
end

@testset "promote_type" begin
    @test promote_type(CategoricalValue{Int}, CategoricalValue{Float64}) ===
        CategoricalValue{Float64}
    @test promote_type(CategoricalValue{Int, UInt8}, CategoricalValue{Float64, UInt32}) ===
        CategoricalValue{Float64, UInt32}
    @test promote_type(CategoricalValue{Int, UInt8}, CategoricalValue{Float64}) ===
        CategoricalValue{Float64}
    # Tests that return Any before Julia 1.3 are due to JuliaLang/julia#29348
    if VERSION >= v"1.3.0-DEV"
        @test promote_type(CategoricalValue{Int},
                           Union{CategoricalValue{Float64}, Missing}) ===
            Union{Missing, Float64}
        @test promote_type(CategoricalValue{Int, UInt8},
                           Union{CategoricalValue{Float64, UInt32}, Missing}) ===
            Union{CategoricalValue{Float64, UInt32}, Missing}
        @test promote_type(Union{CategoricalValue{Int}, Missing},
                           CategoricalValue{Float64}) ===
            Union{Missing, Float64}
        @test promote_type(Union{CategoricalValue{Int, UInt8}, Missing},
                           CategoricalValue{Float64, UInt32}) ===
            Union{CategoricalValue{Float64, UInt32}, Missing}
        @test promote_type(Union{CategoricalValue{Int}, Missing},
                           Union{CategoricalValue{Float64}, Missing}) ===
            Union{Missing, Float64}
    else
        @test promote_type(CategoricalValue{Int},
                           Union{CategoricalValue{Float64}, Missing}) ===
            Any
        @test promote_type(CategoricalValue{Int, UInt8},
                           Union{CategoricalValue{Float64, UInt32}, Missing}) ===
            Union{CategoricalValue{Float64, UInt32}, Missing}
        @test promote_type(Union{CategoricalValue{Int}, Missing},
                           CategoricalValue{Float64}) ===
            Any
        @test promote_type(Union{CategoricalValue{Int, UInt8}, Missing},
                        CategoricalValue{Float64, UInt32}) ===
            Union{CategoricalValue{Float64, UInt32}, Missing}
        @test promote_type(Union{CategoricalValue{Int}, Missing},
                        Union{CategoricalValue{Float64}, Missing}) ===
            Any
    end
    @test promote_type(Union{CategoricalValue{Int, UInt8}, Missing},
                       Union{CategoricalValue{Float64, UInt32}, Missing}) ===
        Union{CategoricalValue{Float64, UInt32}, Missing}

    @test promote_type(CategoricalValue, Missing) === Union{CategoricalValue, Missing}
    @test promote_type(CategoricalValue{Int}, Missing) === Union{CategoricalValue{Int}, Missing}
    @test promote_type(CategoricalValue{Int, UInt32}, Missing) ===
        Union{CategoricalValue{Int, UInt32}, Missing}
    @test promote_type(CategoricalValue{Int, UInt32}, Any) === Any

    @test promote_type(CategoricalString{UInt8}, CategoricalString{UInt32}) ===
        CategoricalString{UInt32}
    @test promote_type(CategoricalString, CategoricalString{UInt32}) ===
        CategoricalString
    @test promote_type(Union{CategoricalString{UInt8}, Missing}, CategoricalString{UInt32}) ===
        Union{CategoricalString{UInt32}, Missing}
    @test promote_type(Union{CategoricalString{UInt8}, Missing}, Union{CategoricalString{UInt32}, Missing}) ===
        Union{CategoricalString{UInt32}, Missing}

    @test promote_type(CategoricalString{UInt8}, CategoricalValue{Int, UInt32}) ===
        CategoricalValue{Any, UInt32}
    @test promote_type(Union{CategoricalString{UInt8}, Missing}, CategoricalValue{Int, UInt32}) ===
        Union{CategoricalValue{Any, UInt32}, Missing}
end

@testset "convert() preserves `ordered`" begin
    pool = CategoricalPool([1, 2, 3], true)
    @test convert(CategoricalPool{Float64, UInt8}, pool).ordered === true
end

@testset "convert() with Union{T, Nothing}" begin
    pool = CategoricalPool([nothing, 2, 3])
    v1 = catvalue(1, pool)
    v2 = catvalue(2, pool)
    @test convert(Union{Int, Nothing}, v1) === nothing
    @test convert(Union{Int, Nothing}, v2) === 2
    @test convert(Union{Float64, Nothing}, v2) === 2.0
end

end


text:
# Math module
module Math

using Dierckx

# types for derivatives, etc
abstract FunctionType

type Derivative <: FunctionType end

type BernsteinPolynomial end

# misc function for comparison
function is_close{T <: Number}(x::T, y::T, n::Int = 42)
  if x == y
    return true
  end

  diff = abs(x - y)
  tol = n * eps(Float64) # machine epsilon

  if (x * y == 0.0) # x or y is 0
    return diff < (tol * tol)
  end

  return diff <= tol * abs(x) && diff <= tol * abs(y)
end

function close_enough{T <: Number}(x::T, y::T, n::Int = 42)
  if x == y
    return true
  end

  diff = abs(x -y)
  tol = n * eps()

  if x * y == 0
    return diff < (tol * tol)
  end

  return diff <= tol * abs(x) || diff <= tol * abs(y)
end

# misc functions - prob put in own file
function divide_array_by_self!{T, N <: Number}(a::Vector{T}, x::N)
  for i = 1:length(a)
    a[i] = a[i] / x
  end

  return a
end

function multiply_array_by_self!{T, N <: Number}(a::Vector{T}, x::N)
  for i = 1:length(a)
    a[i] = a[i] * x
  end

  return a
end

function get_factorial{I <: Integer}(i::I)
  if i > 20
    return Float64(factorial(BigInt(i)))
  else
    return factorial(i)
  end
end

function get_polynomial(::BernsteinPolynomial, i::Int, n::Int, x::Float64)
  coeff = get_factorial(n) / (get_factorial(n-1) * get_factorial(i))

  return coeff * (x ^ i) * (1.0 - x)^(n - i)
end

export is_close, close_enough, divide_array_by_self!, multiply_array_by_self!, get_factorial, get_polynomial, FunctionType, Derivative, BernsteinPolynomial

# Splines
type BSpline
  p::Integer
  n::Integer
  knots::Vector{Float64}
end

function spline_oper{I <: Integer}(spline::BSpline, i::I, x::Float64)
  i <= spline.n || error("i must not be greater than spline.n $i $(spline.n)")
  return N(spline, i, spline.p, x)
end

function N{I <: Integer}(spline::BSpline, i::I, p::I, x::Float64)
  if p == 0
    return (spline.knots[i] <= x && x < spline.knots[i + 1]) ? 1.0 : 0.0
  else
    return ((x - spline.knots[i]) / (spline.knots[i + p] - spline.knots[i])) * N(spline, i, p - 1, x) +
            ((spline.knots[i + p + 1] - x) / (spline.knots[i + p + 1] - spline.knots[i + 1])) * N(spline, i + 1, p - 1, x)
  end
end

export BSpline, spline_oper, N

# Constants
const EPS_VAL = eps()

# utilities.jl
include("utilities.jl")

# grid.jl
export bounded_log_grid
include("grid.jl")

# lmdif2.jl
export lmdif2!
include("lmdif2.jl")

# distribution.jl
export distribution_derivative, peizer_pratt_method_2_inversion
include("distributions.jl")

# integral.jl
export Integrator, IntegrationFunction, SegmentIntegral, operator, integrate, GaussLaguerreIntegration, get_order
include("integral.jl")

# svd.jl
export SVD
include("svd.jl")

# tridiagonal_operator.jl
export TridiagonalOperator, TridiagIdentity, set_mid_row!, set_last_row!, set_first_row!, apply_to, solve_for!
include("tridiagonal_operator.jl")

# matrix.jl
export NoneSalvagingAlgo, rank_reduced_sqrt, OrthogonalProjection, get_vector
include("matrix.jl")

# transformed_grid.jl
export TransformedGrid, LogGrid
include("transformed_grid.jl")

# rng.jl
export AbstractRandomSequenceGenerator, PseudoRandomRSG, InverseCumulativeRSG, SobolRSG, SobolInverseCumulativeRSG, next_sequence!, last_sequence, init_sequence_generator!
include("rng.jl")

# statistics.jl
export NonWeightedStatistics, GenericRiskStatistics, RiskStatistics, gen_RiskStatistics, GenericSequenceStats, reset!, add_sample!, adding_data!, error_estimate, stats_mean, stats_std_deviation,
stats_skewness, stats_kurtosis, stats_covariance
include("statistics.jl")

# sampled_curve.jl
export SampledCurve, get_size, set_log_grid!, set_grid!, sample!, value_at_center, first_derivative_at_center, second_derivative_at_center
include("sampled_curve.jl")

# general_linear_least_squares.jl
export GeneralLinearLeastSquares, get_coefficients
include("general_linear_least_squares.jl")

# interpolation.jl
export Interpolation, Spline, Lagrange, LogInterpolation, LogLinear, BicubicSpline, NaturalCubicSpline, CubicInterpolation, BackwardFlatInterpolation, update!, locate, initialize!, value, get_primitive,
derivative

include("interpolation/interpolation.jl")
include("interpolation/linear_interpolation.jl")
include("interpolation/log_interpolation.jl")
include("interpolation/backward_flat_interpolation.jl")
include("interpolation/spline_interpolation.jl")

# solvers.jl
export Solver1D, BrentSolver, NewtonSolver, FiniteDifferenceNewtonSafe, solve
include("solvers/solver.jl")
include("solvers/brent_solver.jl")
include("solvers/finite_difference.jl")
include("solvers/newton_solver.jl")

# optimization.jl
export lmdif!, Projection, CostFunction, Constraint, NoConstraint, PositiveConstraint, BoundaryConstraint, ProjectedConstraint, OptimizationMethod, LevenbergMarquardt, Simplex, Problem, EndCriteria,
project, minimize!, minimize_2!, include_params
include("optimization/lmdif.jl")
include("optimization/optimization.jl")
include("optimization/problem.jl")
include("optimization/simplex.jl")
include("optimization/levenberg_marquardt.jl")

end


text:
""
function constraint_power_balance_shunt_dispatch(pm::_PM.AbstractACPModel, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_gens, bus_storage, bus_shunts_var, bus_pd, bus_qd, bus_gs_const, bus_bs_const)
    vm   = var(pm, n, :vm, i)
    p    = get(var(pm, n),    :p, Dict()); _PM._check_var_keys(p, bus_arcs, "active power", "branch")
    q    = get(var(pm, n),    :q, Dict()); _PM._check_var_keys(q, bus_arcs, "reactive power", "branch")
    pg   = get(var(pm, n),   :pg, Dict()); _PM._check_var_keys(pg, bus_gens, "active power", "generator")
    qg   = get(var(pm, n),   :qg, Dict()); _PM._check_var_keys(qg, bus_gens, "reactive power", "generator")
    ps   = get(var(pm, n),   :ps, Dict()); _PM._check_var_keys(ps, bus_storage, "active power", "storage")
    qs   = get(var(pm, n),   :qs, Dict()); _PM._check_var_keys(qs, bus_storage, "reactive power", "storage")
    psw  = get(var(pm, n),  :psw, Dict()); _PM._check_var_keys(psw, bus_arcs_sw, "active power", "switch")
    qsw  = get(var(pm, n),  :qsw, Dict()); _PM._check_var_keys(qsw, bus_arcs_sw, "reactive power", "switch")
    p_dc = get(var(pm, n), :p_dc, Dict()); _PM._check_var_keys(p_dc, bus_arcs_dc, "active power", "dcline")
    q_dc = get(var(pm, n), :q_dc, Dict()); _PM._check_var_keys(q_dc, bus_arcs_dc, "reactive power", "dcline")

    bs = get(var(pm, n), :bs, Dict()); _PM._check_var_keys(bs, bus_shunts_var, "reactive power", "shunt")

    cstr_p = JuMP.@NLconstraint(pm.model, 0 == - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*vm^2)
    cstr_q = JuMP.@NLconstraint(pm.model, 0 == - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*vm^2 + sum(bs[s]*vm^2 for s in bus_shunts_var))

    if _IM.report_duals(pm)
        sol(pm, n, :bus, i)[:lam_kcl_r] = cstr_p
        sol(pm, n, :bus, i)[:lam_kcl_i] = cstr_q
    end
end

""
function constraint_power_balance_shunt_dispatch_soft(pm::_PM.AbstractACPModel, n::Int, i::Int, bus_arcs, bus_arcs_dc, bus_arcs_sw, bus_storage, bus_gens, bus_shunts_var, bus_pd, bus_qd, bus_gs_const, bus_bs_const)
    vm   = var(pm, n, :vm, i)
    p    = get(var(pm, n),    :p, Dict()); _PM._check_var_keys(p, bus_arcs, "active power", "branch")
    q    = get(var(pm, n),    :q, Dict()); _PM._check_var_keys(q, bus_arcs, "reactive power", "branch")
    pg   = get(var(pm, n),   :pg, Dict()); _PM._check_var_keys(pg, bus_gens, "active power", "generator")
    qg   = get(var(pm, n),   :qg, Dict()); _PM._check_var_keys(qg, bus_gens, "reactive power", "generator")
    ps   = get(var(pm, n),   :ps, Dict()); _PM._check_var_keys(ps, bus_storage, "active power", "storage")
    qs   = get(var(pm, n),   :qs, Dict()); _PM._check_var_keys(qs, bus_storage, "reactive power", "storage")
    psw  = get(var(pm, n),  :psw, Dict()); _PM._check_var_keys(psw, bus_arcs_sw, "active power", "switch")
    qsw  = get(var(pm, n),  :qsw, Dict()); _PM._check_var_keys(qsw, bus_arcs_sw, "reactive power", "switch")
    p_dc = get(var(pm, n), :p_dc, Dict()); _PM._check_var_keys(p_dc, bus_arcs_dc, "active power", "dcline")
    q_dc = get(var(pm, n), :q_dc, Dict()); _PM._check_var_keys(q_dc, bus_arcs_dc, "reactive power", "dcline")

    bs = get(var(pm, n), :bs, Dict()); _PM._check_var_keys(bs, bus_shunts_var, "reactive power", "shunt")

    p_delta_abs = var(pm, n, :p_delta_abs, i)
    q_delta_abs = var(pm, n, :q_delta_abs, i)

    p_delta = @NLexpression(pm.model, - sum(p[a] for a in bus_arcs) + sum(pg[g] for g in bus_gens) - sum(pd for pd in values(bus_pd)) - sum(gs for gs in values(bus_gs_const))*vm^2)
    q_delta = @NLexpression(pm.model, - sum(q[a] for a in bus_arcs) + sum(qg[g] for g in bus_gens) - sum(qd for qd in values(bus_qd)) + sum(bs for bs in values(bus_bs_const))*vm^2 + sum(bs[s]*vm^2 for s in bus_shunts_var))

    @NLconstraint(pm.model,  p_delta_abs >= p_delta)
    @NLconstraint(pm.model, -p_delta_abs <= p_delta)

    @NLconstraint(pm.model,  q_delta_abs >= q_delta)
    @NLconstraint(pm.model, -q_delta_abs <= q_delta)
end

""
function constraint_voltage_magnitude_link(pm::_PM.AbstractACPModel, n_1::Int, n_2::Int, i::Int)
    vm_1 = var(pm, n_1, :vm, i)
    vm_2 = var(pm, n_2, :vm, i)

    JuMP.@constraint(pm.model, vm_1 == vm_2)
end

""
function constraint_ohms_yt_from_goc(pm::_PM.AbstractACPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    p_fr  = var(pm, n,  :p, f_idx)
    q_fr  = var(pm, n,  :q, f_idx)
    vm_fr = var(pm, n, :vm, f_bus)
    vm_to = var(pm, n, :vm, t_bus)
    va_fr = var(pm, n, :va, f_bus)
    va_to = var(pm, n, :va, t_bus)

    JuMP.@NLconstraint(pm.model, p_fr ==  (g/tm^2+g_fr)*vm_fr^2 + (-g*tr+b*ti)/tm^2*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/tm^2*(vm_fr*vm_to*sin(va_fr-va_to)) )
    JuMP.@NLconstraint(pm.model, q_fr == -(b/tm^2+b_fr)*vm_fr^2 - (-b*tr-g*ti)/tm^2*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/tm^2*(vm_fr*vm_to*sin(va_fr-va_to)) )
end

""
function expression_bus_withdrawal(pm::_PM.AbstractACPModel, n::Int, i::Int, bus_storage, bus_pd, bus_qd, bus_gs, bus_bs)
    vm = var(pm, n, :vm, i)
    ps = get(var(pm, n), :ps, Dict()); _PM._check_var_keys(ps, bus_storage, "active power", "storage")
    qs = get(var(pm, n), :qs, Dict()); _PM._check_var_keys(ps, bus_storage, "reactive power", "storage")

    ps_total = 0.0
    qs_total = 0.0
    if length(bus_storage) > 0
        ps_total = sum(ps[s] for s in bus_storage)
        qs_total = sum(qs[s] for s in bus_storage)
    end

    pd_total = 0.0
    qd_total = 0.0
    if length(bus_pd) > 0
        pd_total = sum(pd for pd in values(bus_pd))
    end
        if length(bus_qd) > 0
        qd_total = sum(qd for qd in values(bus_qd))
    end

    gs_total = 0.0
    bs_total = 0.0
    if length(bus_gs) > 0
        gs_total = sum(gs for gs in values(bus_gs))*vm^2
    end
    if length(bus_bs) > 0
        bs_total = sum(bs for bs in values(bus_bs))*vm^2
    end

    var(pm, n, :bus_wdp)[i] = ps_total + pd_total + gs_total
    var(pm, n, :bus_wdq)[i] = qs_total + qd_total + bs_total
end

""
function expression_branch_power_ohms_yt_from_goc(pm::_PM.AbstractACPModel, n::Int, f_bus, t_bus, f_idx, t_idx, g, b, g_fr, b_fr, tr, ti, tm)
    vm_fr = var(pm, n, :vm, f_bus)
    vm_to = var(pm, n, :vm, t_bus)
    va_fr = var(pm, n, :va, f_bus)
    va_to = var(pm, n, :va, t_bus)

    var(pm, n, :p)[f_idx] = @NLexpression(pm.model,  (g/tm^2+g_fr)*vm_fr^2 + (-g*tr+b*ti)/tm^2*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/tm^2*(vm_fr*vm_to*sin(va_fr-va_to)) )
    var(pm, n, :q)[f_idx] = @NLexpression(pm.model, -(b/tm^2+b_fr)*vm_fr^2 - (-b*tr-g*ti)/tm^2*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/tm^2*(vm_fr*vm_to*sin(va_fr-va_to)) )
end


text:
using MPIMapReduce
using MPI

MPI.Init()

const comm = MPI.COMM_WORLD
const root = 0

y = MPIMapReduce.pmapgatherv(x -> x^2, vcat, 1:1)
y2 = mapreduce(x -> x^2, vcat, 1:1)

if MPI.Comm_rank(comm) == root
    @info "MPI value"
    show(stdout, MIME"text/plain"(), y)
    println()
    @info "Julia value"
    show(stdout, MIME"text/plain"(), y2)
    println()
end

y = MPIMapReduce.pmapgatherv(x -> x^2, vcat, 1:5)
y2 = mapreduce(x -> x^2, vcat, 1:5)

if MPI.Comm_rank(comm) == root
    @info "MPI value"
    show(stdout, MIME"text/plain"(), y)
    println()
    @info "Julia value"
	show(stdout, MIME"text/plain"(), y2)
	println()
end

y = MPIMapReduce.pmapgatherv(x -> ones(1,2) * x^2, vcat, 1:5)
y2 = mapreduce(x -> ones(1,2) * x^2, vcat, 1:5)

if MPI.Comm_rank(comm) == root
    @info "MPI value"
    show(stdout, MIME"text/plain"(), y)
    println()
    @info "Julia value"
	show(stdout, MIME"text/plain"(), y2)
	println()
end

y = MPIMapReduce.pmapgatherv(x -> ones(2,2,2) * x^2, vcat, 1:5)
y2 = mapreduce(x -> ones(2,2,2) * x^2, vcat, 1:5)

if MPI.Comm_rank(comm) == root
    @info "MPI value"
    show(stdout, MIME"text/plain"(), y)
    println()
    @info "Julia value"
	show(stdout, MIME"text/plain"(), y2)
	println()
end


text:
#this file contains functions/wrappers related to the solve function in FEniCS
#https://fenicsproject.org/olddocs/dolfin/1.3.0/python/programmers-reference/fem/solving/solve.html
using FEniCS

solve(A::Matrix,x, b::Matrix,solvers...) =fenics.solve(A.pyobject,x,b.pyobject,solvers...)
export solve

#lvsolve is the linear variational solver
lvsolve(a,L,u;solver_parameters::Dict=Dict("linear_solver"=>"default"),form_compiler_parameters::Dict=Dict("optimize"=>true))=fenics.solve(a.pyobject==L.pyobject, u.pyobject,solver_parameters=solver_parameters,form_compiler_parameters=form_compiler_parameters)

lvsolve(a,L,u,bcs=nothing;solver_parameters::Dict=Dict("linear_solver"=>"default"),form_compiler_parameters::Dict=Dict("optimize"=>true))=fenics.solve(a.pyobject==L.pyobject, u.pyobject, bcs=bcs.pyobject,solver_parameters=solver_parameters,form_compiler_parameters=form_compiler_parameters)
#allows BoundaryCondition to be provided in an AbstractArray (of type BoundaryCondition)
function lvsolve(a,L,u,bcs::AbstractArray;solver_parameters::Dict=Dict("linear_solver"=>"default"),form_compiler_parameters::Dict=Dict("optimize"=>true))
    bcs_py = [bc.pyobject for bc in bcs]
    fenics.solve(a.pyobject==L.pyobject, u.pyobject, bcs=bcs_py,solver_parameters=solver_parameters,form_compiler_parameters=form_compiler_parameters)
end

export lvsolve
#Dict("linear_solver"=>"default")
#Dict("optimize"=>true)
#nlvsolve is the non-linear variational solver
nlvsolve(F,u;J=nothing,solver_parameters::Dict=Dict("nonlinear_solver"=>"newton"),form_compiler_parameters::Dict=Dict("optimize"=>true))=fenics.solve(F.pyobject==0,u.pyobject,J=J,solver_parameters=solver_parameters,form_compiler_parameters=form_compiler_parameters)
nlvsolve(F,u,bcs=nothing;J=nothing,solver_parameters::Dict=Dict("nonlinear_solver"=>"newton"),form_compiler_parameters::Dict=Dict("optimize"=>true))=fenics.solve(F.pyobject==0,u.pyobject,J=J,solver_parameters=solver_parameters,form_compiler_parameters=form_compiler_parameters)

nlvsolve(F,u,bcs::BoundaryCondition;J=nothing,solver_parameters::Dict=Dict("nonlinear_solver"=>"newton"),form_compiler_parameters::Dict=Dict("optimize"=>true))=fenics.solve(F.pyobject==0,u.pyobject,bcs=bcs.pyobject,J=J,solver_parameters=solver_parameters,form_compiler_parameters=form_compiler_parameters)
#allows BoundaryCondition to be provided in an AbstractArray (of type BoundaryCondition)
function nlvsolve(F,u,bcs::AbstractArray;J=nothing,solver_parameters::Dict=Dict("nonlinear_solver"=>"newton"),form_compiler_parameters::Dict=Dict("optimize"=>true))
    bcs_py = [bc.pyobject for bc in bcs]
    fenics.solve(F.pyobject==0,u.pyobject,bcs=bcs,J=J,solver_parameters=solver_parameters,form_compiler_parameters=form_compiler_parameters)
end

export nlvsolve

#anlvsolve corresponds to the adapative non-linear solver.
anlvsolve(F,a,u,bcs,tol,M)=fenics.solve(F.pyobject==a.pyobject,u.pyobject,bcs=bcs.pyobject,tol=tol,M=M)
#this function hasnt been tested yet, so isnt exported

"""
errornorm is the function to calculate the error between our exact and calculated
solution. The norm kwarg defines the norm measure used (by default it is the L2 norm)
"""
errornorm(ans,sol;norm="L2") = fenics.errornorm(ans.pyobject,sol.pyobject,norm)
export errornorm

File(path::StringOrSymbol)=fenics.File(path) #used to store the solution in various formats

function File(path::StringOrSymbol,object::FeFunction)
    vtkfile = File(path)
    vtkfile << object.pyobject
end

function File(path::StringOrSymbol,object::FeFunction, time::Number)
    vtkfile = File(path)
    vtkfile << (object.pyobject,time)
end

function File(path::StringOrSymbol,object::Mesh)
    vtkfile = File(path)
    vtkfile << object.pyobject
end

function File(path::StringOrSymbol,object::Mesh, time::Number)
    vtkfile = File(path)
    vtkfile << (object.pyobject,time)
end

export File

XDMFFile(path::StringOrSymbol) = fenics.XDMFFile(path)
export XDMFFile

TimeSeries(path::StringOrSymbol) = fenics.TimeSeries(path)
retrieve(timeseries,placeholder,time) = timeseries.retrieve(placeholder,time)
export TimeSeries, retrieve

write(path,solution,time::Number) = path.write(solution.pyobject, time)
write(path,solution::PyObject,time::Number) = path.write(solution, time)

store(path,solution,time::Number) = path.store(solution.pyobject, time)
store(path,solution::PyObject,time::Number) = path.store(solution, time)

export write,store

array(matrix) = fenicspycall(matrix, :gather_on_zero)
vector(solution::FeFunction) = fenicspycall(solution,:vector) #
interpolate(ex, V::FunctionSpace) = FeFunction(fenics.interpolate(ex.pyobject, V.pyobject))

export vector, interpolate,array

"the following function is used to extract the array from a solution/form"
function get_array(form::Expression)
    assembled_form = assemble(form)
    return array(assembled_form)
end

function get_array(solution::FeFunction)
    generic_vector = vector(solution)
    instantiated_vector = fenics.Vector(generic_vector)
    return instantiated_vector.gather_on_zero()
end
"""
we return the array from an assembled form
"""
function get_array(assembled_form::Matrix)
    return array(assembled_form)
end

export get_array

"""
Return projection of given expression *v* onto the finite element space *V*
*Example of usage*
          v = Expression("sin(pi*x[0])")
          V = FunctionSpace(mesh, "Lagrange", 1)
          Pv = project(v, V)
"""
project(v::Union{FeFunction,Expression},V::FunctionSpace)  = FeFunction(fenics.project(v.pyobject,V.pyobject))
export project


text:
module Graphtypes

export
    Graph,
    SimpleGraph,
    StructureGraph,

    simple_graph,
    structure_graph,

    connected,
    count_connections,
    list_outgoing,
    total_edges,
    total_vertices,
    list_edges,
    connect!,
    disconnect!,
    reverse!

abstract Graph{directed, adjacent}

type SimpleGraph <: Graph
    adjacent::Matrix{Int64} #adjacency matrix
end

type StructureGraph <: Graph
    #undirected, no duplicate edges
    adjacent::Matrix{Bool}
end

function simple_graph (v_count::Int64)
    g = SimpleGraph(zeros(Int64, v_count, v_count))
    return g
end

function structure_graph (v_count::Int64)
    g = StructureGraph(zeros(Bool, v_count, v_count))
    return g
end

function list_edges(g::SimpleGraph)
    size = total_vertices(g)
    edges = Array(Int64, 0, 2)
    for i in 1:size
        for j in 1:size
            for k in 1:count_connections(g,i,j)
                edges = [edges, [i j]]
    end end end
    return edges
end

function list_edges(g::StructureGraph)
    size = total_vertices(g)
    edges = Array(Int64, 0, 2)
    for i in 1:size
        for j in i:size
            if g.adjacent[i,j]
                edges = [edges, [i j]]
    end end end
    return edges
end

function list_outgoing(g::SimpleGraph, from::Int64)
    return find( x->(x != 0), g.adjacent[from, :])
end

function list_outgoing(g::StructureGraph, from::Int64)
    return sort(unique([    find( x->(x != 0), g.adjacent[from, :]),
                    find( x->(x != 0), g.adjacent[:, from])]))
end

function connected(g::Graph, from::Int64, to::Int64)
    return g.adjacent[from, to] != 0
end

function count_connections(g::Graph, from::Int64, to::Int64)
    return g.adjacent[from, to]
end

function total_edges(g::SimpleGraph)
    return sum(g.adjacent)
end

function total_edges(g::StructureGraph)
    return sum(g.adjacent)/2
end

function total_vertices(g::Graph)
    return size(g.adjacent)[2]
end

function connect!(g::SimpleGraph, from::Int64, to::Int64)
    g.adjacent[from, to] += 1
end

function connect!(g::StructureGraph, from::Int64, to::Int64)
    g.adjacent[from, to] = true
    g.adjacent[to, from] = true
end

function disconnect!(g::SimpleGraph, from::Int64, to::Int64)
    if g.adjacent[from, to] > 0
        g.adjacent[from, to] -= 1
        return true
    else
        return false
    end
end

function disconnect!(g::StructureGraph, from::Int64, to::Int64)
    if g.adjacent[from, to]
        g.adjacent[from, to] = false
        return true
    else
        return false
    end
end
function reverse!(g::Graph, i1::Int64, i2::Int64)
    if disconnect!(g,i1,i2)
        connect!(g,i2,i1)
        return true
    end
    return false
end

end

max_stars_repo_name: AaltoScienceInst/pebblegame

max_issues_repo_name: AaltoScienceInst/pebblegame

max_forks_repo_name: AaltoScienceInst/pebblegame


text:
export sequenceDir, sequenceList, Sequence, triangle, hasEnableSequence,
       enableSequence,
       makeTriangleSequence, makeSineSequence, makeTriSweepSequence, makeTriSweepDeadSequence

#6.312341575312002
#13.705627903244254

sequenceDir() = @__DIR__

sequenceList() = String[ splitext(seq)[1] for seq in filter(a->contains(a,".toml"),readdir(sequenceDir()))]

mutable struct Sequence
  values::Matrix{Float64}
  numPeriodsPerPatch::Int64
  enable::Matrix{Bool}
end

Sequence(values, numPeriodsPerPatch) = Sequence(values, numPeriodsPerPatch, zeros(Bool,0,0))

MPIFiles.acqNumPeriodsPerFrame(s::Sequence) = size(s.values,2)*s.numPeriodsPerPatch
MPIFiles.acqNumPatches(s::Sequence) = size(s.values,2)
MPIFiles.acqNumPeriodsPerPatch(s::Sequence) = s.numPeriodsPerPatch

hasEnableSequence(seq::Sequence) = !isempty(seq.enable)
enableSequence(seq::Sequence) = seq.enable

function Sequence(name::String)
  filename = joinpath(sequenceDir(),name*".toml")
  p = TOML.parsefile(filename)
  values = reshape(p["values"], :, p["numPatches"])
  enable = haskey(p, "enable") ? reshape(p["enable"], :, p["numPatches"]) : zeros(Bool,0,0)
  return Sequence(values, p["numPeriodsPerPatch"], enable)
end

function saveSeq(name::String, seq::Sequence)
  p = Dict{String,Any}()
  p["numPeriodsPerPatch"] = seq.numPeriodsPerPatch
  p["values"] = vec(seq.values)
  p["numPatches"] = size(seq.values, 2)
  if hasEnableSequence(seq)
    p["enable"] = vec(seq.enable)
  end

  filename = joinpath(sequenceDir(),name*".toml")
  open(filename, "w") do f
    TOML.print(f, p)
  end
end

function makeSineSequence(name::String, minCurr, maxCurr, patches, periodsPerPatch)

  t = range(0,2*pi,length=patches+1)[1:end-1]

  A = ((maxCurr-minCurr)/2).*sin.(t).+(minCurr+maxCurr)/2
  B = -(((maxCurr-minCurr)/2).*sin.(t)).+(minCurr+maxCurr)/2

  C = cat(A',B', dims=1)

  saveSeq(name, Sequence(C, periodsPerPatch))
  return C
end

function triangle(t)
  t_ = mod(t, 2*pi)
  if 0 <= t_ < pi/2
    return t_/(pi/2)
  elseif pi/2 <= t_ < 3*pi/2
    return 2.0-t_/(pi/2)
  elseif 3*pi/2 <= t_ < 2*pi
    return t_/(pi/2) - 4.0
  end
end

function triangleNegSlope(t)
  t_ = mod(t, pi)
  if 0 <= t_ < pi
    return 1.0-t_/(pi/2)
  end
end
function trianglePosSlope(t)
  t_ = mod(t, pi)
  if 0 <= t_ < pi
    return t_/(pi/2)-1
  end
end
function triangleRampUp(t)
  t_ = mod(t, pi/2)
  if 0 <= t_ < pi/2
    return t_/(pi/2)
  end
end
function triangleRampDown(t)
  t_ = mod(t, pi/2)
  if 0 <= t_ < pi/2
    return 1.0-t_/(pi/2)
  end
end

function makeTriangleSequence(name::String, minCurr, maxCurr, patches, periodsPerPatch)

  t = range(0,2*pi,length=patches+1)[1:end-1]

  A = ((maxCurr-minCurr)/2).*triangle.(t).+(minCurr+maxCurr)/2
  B = -(((maxCurr-minCurr)/2).*triangle.(t)).+(minCurr+maxCurr)/2

  C = cat(A',B', dims=1)

  D = ones(Bool, patches)
  D[1:2:end] .= false
  E = ones(Bool, patches)

  F = cat(D',E', dims=1)

  saveSeq(name, Sequence(C, periodsPerPatch, F))
  return C
end

#e.g makeTriSweepSequence("test",-5.0:1.0:5.0,4,0:1.0:4.0,1)
function makeTriSweepSequence(name::String, XOffsetRange::StepRangeLen, periodsPerXOffset::Number, YOffsetRange::StepRangeLen, calibOffset::Bool)

  NumXOffsets = length(XOffsetRange)
  NumYOffsets = length(YOffsetRange)
  XminCurr = XOffsetRange[1]    #*calib factor
  XmaxCurr = XOffsetRange[end]  #*calib factor
  YminCurr = YOffsetRange[1]    #*calib factor
  YstepCurr = convert(Float64, YOffsetRange.step)    #*calib factor

  #toDo: conversion from current to mT for Y-channel?

  if calibOffset == true
    calibXoffsetCurr = 0.0
    calibYoffsetCurr = -0.31 #[A] offset, OPamp etc... ??
  else
    calibXoffsetCurr = 0.0
    calibYoffsetCurr = 0.0
  end

  #######################
  #crearte triangle with 2 datapoints in peak
  #which requires: slope%4 = 3
  #then every XOffset has both extrema as a datapoint
  if isodd(NumXOffsets)
    triangleSlope = 2*NumXOffsets+1
  elseif iseven(NumXOffsets)
    triangleSlope = 2*(NumXOffsets+1)+1
  end

  if mod(triangleSlope,4) != 3
    @error "Sequence is required to have odd amount of steps. Mod(len,4) should equal 3 to obtain triangle-plateau." triangleSlope%4
  end

  flagYZero = 0
  if iseven(NumYOffsets)
    NumYOffsets +=1
    flagYZero = 1
  end

  #########################
  t = range(0,2*pi,length=triangleSlope)[1:end-1]
  tri = ( ((XmaxCurr-XminCurr)/2).*triangle.(t).+(XminCurr+XmaxCurr)/2 )/maximum(triangle.(t))  #TODO: CAL offset due to plateau

  A = repeat(tri, ceil(Int,NumYOffsets/2) )
  B = YstepCurr .* zeros( ceil(Int,NumXOffsets/2)) .+YminCurr

  for i in 0:floor(Int,NumYOffsets)-1
    if flagYZero == 0
      B = vcat(B, YstepCurr*i .* ones(NumXOffsets) .+YminCurr)
    elseif flagYZero == 1 && i == floor(Int,NumYOffsets)-1
      B = vcat(B,zeros(NumXOffsets))   #back to zero in case of even No of Yoffsets
    else
      B = vcat(B, YstepCurr*i .* ones(NumXOffsets) .+YminCurr)
    end
  end

  B = vcat(B, zeros( floor(Int,NumXOffsets/2) )) .+ calibYoffsetCurr

  Seq = cat(A',B', dims=1)
  saveSeq(name, Sequence(Seq, periodsPerXOffset))
  return Seq
end

function makeTriSweepDeadSequence(name::String, XOffsetRange::StepRangeLen, periodsPerXOffset::Number, YOffsetRange::StepRangeLen, deadPatches::Number, calibOffset::Bool)

  NumXOffsets = length(XOffsetRange)
  NumYOffsets = length(YOffsetRange)
  XminCurr = XOffsetRange[1]    #*calib factor
  XmaxCurr = XOffsetRange[end]  #*calib factor
  YminCurr = YOffsetRange[1]    #*calib factor
  YstepCurr = convert(Float64, YOffsetRange.step)    #*calib factor

  #toDo: conversion from current to mT for Y-channel?

  if calibOffset == true
    calibXoffsetCurr = 0.0
    calibYoffsetCurr = -0.31 #[A] offset, OPamp etc... ??
  else
    calibXoffsetCurr = 0.0
    calibYoffsetCurr = 0.0
  end

  #######################
  #crearte triangle with 2 datapoints in peak
  #which requires: slope%4 = 3
  #then every XOffset has both extrema as a datapoint
  if isodd(NumXOffsets)
    triangleSlope = 2*NumXOffsets+1
  elseif iseven(NumXOffsets)
    triangleSlope = 2*(NumXOffsets+1)+1
  end

  if mod(triangleSlope,4) != 3
    @error "Sequence is required to have odd amount of steps. Mod(len,4) should equal 3 to obtain triangle-plateau." triangleSlope%4
  end
  #########################

  t_half = range(0,pi,length=NumXOffsets)[1:end-1]
  t_quart = range(0,pi/2,length=ceil(Int,NumXOffsets/2))[1:end-1]

  A = triangleRampUp.(t_quart)
  B = YstepCurr .* zeros(length(A)+1) .+YminCurr
  for i in 0:2:NumYOffsets  #if even, one offset more than necessary is created. throw away in data analysis.
    t1 = ones(deadPatches)
    t2 = triangleNegSlope.(t_half)
    B = vcat(B, YstepCurr*i .* ones(length(t1)+length(t2)) .+YminCurr)

    t3 = -ones(deadPatches)
    t4 = trianglePosSlope.(t_half)
    B = vcat(B, YstepCurr*(i+1) .* ones(length(t3)+length(t4)) .+YminCurr)

    A = vcat(A,t1,t2,t3,t4)
  end
  t5 = vcat(ones(1),triangleRampDown.(t_quart))

  A = ((XmaxCurr-XminCurr)/2).* vcat(A,t5) .+(XminCurr+XmaxCurr)/2
  B = vcat(B, zeros(length(t5)-1)) .+ calibYoffsetCurr

  Seq = cat(A',B', dims=1)
  saveSeq(name, Sequence(Seq, periodsPerXOffset))
  return Seq
end

#=
function makeSequence(name::String, minCurr, maxCurr, patches, periodsPerPatch)
  A=collect(linspace(minCurr,maxCurr,patches))

  C = cat(A',reverse(A',dims=2), dims=1)
  B=repeat(cat(C,reverse(C,dims=2), dims=2), inner=(1,periodsPerPatch))
  B = circshift(B,(0,div(size(B,2),4)))
  writedlm(joinpath(@__DIR__, name*".csv"), B, ',')
  return B
end

function makeSineSequence(name::String, minCurrCh1, maxCurrCh1,minCurrCh2, maxCurrCh2, patches, periodsPerPatch)

  t=range(0,2*pi,length=patches+1)[1:end-1]

  A=((maxCurrCh1-minCurrCh1)/2).*sin.(t).+(minCurrCh1+maxCurrCh1)/2
  B=-(((maxCurrCh2-minCurrCh2)/2).*sin.(t)).+(minCurrCh2+maxCurrCh2)/2

  C = cat(A',B', dims=1)
  D=repeat(C, inner=(1,periodsPerPatch))

  writedlm(joinpath(@__DIR__, name*".csv"), D, ',')
  return D
end
function makeSineSequenceTest(name::String, minCurr, maxCurr, patches, periodsPerPatch)

  t=range(0,2*pi,length=patches+1)[1:end-1]
  A=circshift(((maxCurr-minCurr)/2).*sin.(t).+(minCurr+maxCurr)/2,-6)
  B=circshift(-(((maxCurr-minCurr)/2).*sin.(t)).+(minCurr+maxCurr)/2,-6)

  C = cat(A',B', dims=1)
  D=repeat(C, inner=(1,periodsPerPatch))

  writedlm(joinpath(@__DIR__, name*".csv"), D, ',')
  return D
end

=#


text:
# script to parse CUDA headers and generate Julia wrappers

using Crayons

#
# Parsing
#

using Clang

function wrap(name, headers...; wrapped_headers=headers, library="lib$name", defines=[])
    include_dirs = map(dir->joinpath(dir, "include"), find_toolkit())
    filter!(isdir, include_dirs)

    clang_args = String[]
    append!(clang_args, map(dir->"-I$dir", include_dirs))
    for define in defines
        if isa(define, Pair)
            append!(clang_args, ["-D", "$(first(define))=$(last(define))"])
        else
            append!(clang_args, ["-D", "$define"])
        end
    end

    output_file = "lib$(name).jl"
    common_file = "lib$(name)_common.jl"

    context = init(;
                    headers = [headers...],
                    output_file = output_file,
                    common_file = common_file,
                    clang_includes = [include_dirs..., CLANG_INCLUDE],
                    clang_args = clang_args,
                    header_wrapped = (root, current)->any(header->endswith(current, header), wrapped_headers),
                    header_library = x->library,
                    clang_diagnostics = true,
                  )
    run(context)

    return output_file, common_file
end

#
# Pointer type database
#

using JSON, DataStructures

const db_path = joinpath(@__DIR__, "pointers.json")

const db = if isfile(db_path)
    JSON.parsefile(db_path; dicttype=DataStructures.OrderedDict)
else
    Dict{String, Any}()
end

function save_db()
    global db
    open(db_path, "w") do io
        JSON.print(io, db, 4)
    end
end

#
# Fixing-up
#

using CSTParser, Tokenize

## pass infrastructure

struct Edit{T}
    loc::T
    text::String
end

function pass(x, state, f = (x, state)->nothing)
    f(x, state)
    if x.args isa Vector
        for a in x.args
            pass(a, state, f)
        end
    else
        state.offset += x.fullspan
    end
    state
end

function apply(text, edit::Edit{Int})
    string(text[1:edit.loc], edit.text, text[nextind(text, edit.loc):end])
end
function apply(text, edit::Edit{UnitRange{Int}})
    # println("Rewriting '$(text[edit.loc])' to '$(edit.text)'")
    string(text[1:prevind(text, first(edit.loc))], edit.text, text[nextind(text, last(edit.loc)):end])
end

## rewrite passes

mutable struct State
    offset::Int
    edits::Vector{Edit}
end

# insert `@check` before each `ccall` when it returns a checked type,
# and make it a `@runtime_ccall`
const checked_types = [
    "cublasStatus_t",
    "cudnnStatus_t",
    "cufftResult",
    "curandStatus_t",
    "cusolverStatus_t",
    "cusparseStatus_t",
    "cutensorStatus_t",
]
function rewrite_ccall(x, state)
    if x isa CSTParser.EXPR && x.typ == CSTParser.Call && x.args[1].val == "ccall"
        # get the ccall return type
        rv = x.args[5]

        if rv.val in checked_types
            push!(state.edits, Edit(state.offset, "@check @runtime_"))
        else
            push!(state.edits, Edit(state.offset, "@runtime_"))
        end
    end
end

# change ::Ptr arguments to ::CuPtr / ::PtrOrCuPtr based on user input
function rewrite_pointers(x, state, headers)
    if x isa CSTParser.EXPR && x.typ == CSTParser.Call && x.args[1].val == "ccall"
        # get the ccall arguments, skipping comma's and parentheses
        handle = x.args[3]
        fn = handle.args[2].args[2].val
        rv = x.args[5]
        tt = x.args[7]
        args = x.args[9:2:end]

        # check for pointer types
        types = tt.args[2:2:end-1]
        is_pointer = Bool[x.typ == CSTParser.Curly && x.args[1].val == "Ptr" for x in types]
        offset = state.offset + sum(x->x.fullspan, x.args[1:6])
        if any(is_pointer)
            if haskey(db, fn)
                replacements = db[fn]

                # regenerate replacements with the new argument names
                changed = false
                old_replacements = collect(replacements)
                replacements = OrderedDict{String,Any}()
                for (i, arg) in enumerate(args)
                    if arg.val != old_replacements[i].first
                        changed = true
                    end
                    replacements[arg.val] = old_replacements[i].second
                end

                if changed
                    db[fn] = replacements
                    save_db()
                end
            else
                # print some context from the header (some mention host/device pointer)
                print(Crayon(foreground=:yellow))
                run(`awk "/\<$fn\>/,/;/" $headers`)
                println(Crayon(reset=true))

                # print pointer arguments and their types
                println(Crayon(foreground = :red), fn, Crayon(reset=true))
                for (i, arg) in enumerate(args)
                    if is_pointer[i]
                        println("- argument $i: $(arg.val)::$(Expr(types[i]))")
                    end
                end
                println()

                # prompt
                run(pipeline(`echo -n $fn`, `xclip -i -selection clipboard`));
                print("GPU pointers> ")
                gpu_pointers = parse.(Int, split(readline(stdin)))
                if gpu_pointers == [0]
                    # 0 is special match for all pointers
                    gpu_pointers = findall(is_pointer)
                elseif !isempty(gpu_pointers) && all(i->i<0, gpu_pointers)
                    # negative indicates all but these
                    gpu_pointers = map(i->-i, gpu_pointers)
                    @assert all(i->is_pointer[i], gpu_pointers) "You selected non-pointer arguments"
                    gpu_pointers = setdiff(findall(is_pointer), gpu_pointers)
                end
                print("Dual GPU/CPU pointers> ")
                dual_pointers = parse.(Int, split(readline(stdin)))
                @assert all(i->is_pointer[i], gpu_pointers ∪ dual_pointers) "You selected non-pointer arguments"

                # generate replacements
                replacements = OrderedDict{String,Any}()
                for (i, arg) in enumerate(args)
                    replacements[arg.val] = if is_pointer[i]
                        if i in gpu_pointers
                            "CuPtr"
                        elseif i in dual_pointers
                            "PtrOrCuPtr"
                        else
                            "Ptr"
                        end
                    else
                        nothing
                    end
                end

                db[fn] = replacements
                save_db()
            end

            # generate edits
            for (i, (_,replacement)) in enumerate(replacements)
                offset += tt.args[2*i-1].fullspan
                if replacement !== nothing && replacement != "Ptr"
                    ptr = types[i].args[1]
                    push!(state.edits, Edit(offset+1:offset+ptr.span, replacement))
                end
                offset += types[i].fullspan
            end

            println()
        end
    end
end

# replace handles from the CUDA runtime library with CUDA driver API equivalents
function rewrite_runtime(x, state)
    if x isa CSTParser.EXPR && x.typ == CSTParser.IDENTIFIER && x.val == "cudaStream_t"
        offset = state.offset
        push!(state.edits, Edit(offset+1:offset+x.span, "CUstream"))
    end
end

## indenting passes

mutable struct IndentState
    offset::Int
    lines
    edits::Vector{Edit}
end

function get_lines(text)
    lines = Tuple{Int,Int}[]
    pt = Tokens.EMPTY_TOKEN(Tokens.Token)
    for t in CSTParser.Tokenize.tokenize(text)
        if pt.endpos[1] != t.startpos[1]
            if t.kind == Tokens.WHITESPACE
                nl = findfirst("\n", t.val) != nothing
                if !nl
                    push!(lines, (length(t.val), 0))
                else
                end
            else
                push!(lines, (0, 0))
            end
        elseif t.startpos[1] != t.endpos[1] && t.kind == Tokens.TRIPLE_STRING
            nls = findall(x->x == '\n', t.val)
            for nl in nls
                push!(lines, (t.startpos[2] - 1, nl + t.startbyte))
            end
        elseif t.startpos[1] != t.endpos[1] && t.kind == Tokens.WHITESPACE
            push!(lines, (t.endpos[2], t.endbyte - t.endpos[2] + 1))
        end
        pt = t
    end
    lines
end

function wrap_at_comma(x, state, indent, offset, column)
    comma = nothing
    for y in x
        if column + y.fullspan > 92 && comma !== nothing
            column = indent
            push!(state.edits, Edit(comma, ",\n" * " "^column))
            column += offset - comma[1] - 1 # other stuff might have snuck between the comma and the current expr
            comma = nothing
        elseif y.typ == CSTParser.PUNCTUATION && y.kind == Tokens.COMMA
            comma = (offset+1):(offset+y.fullspan)
        end
        offset += y.fullspan
        column += y.fullspan
    end
end

function indent_ccall(x, state)
    if x isa CSTParser.EXPR && x.typ == CSTParser.MacroCall && x.args[1].args[2].val == "runtime_ccall"
        # figure out how much to indent by looking at where the expr starts
        line = findlast(y -> state.offset >= y[2], state.lines) # index, not the actual number
        line_indent, line_offset = state.lines[line]
        expr_indent = state.offset - line_offset
        indent = expr_indent + sum(x->x.fullspan, x.args[1:2])

        if length(x.args[7]) > 2    # non-empty tuple type
            # break before the tuple type
            offset = state.offset + sum(x->x.fullspan, x.args[1:6])
            push!(state.edits, Edit(offset:offset, "\n" * " "^indent))

            # wrap tuple type
            wrap_at_comma(x.args[7], state, indent+1, offset, indent+1)
        end

        if length(x.args) > 9
            # break before the arguments
            offset = state.offset + sum(x->x.fullspan, x.args[1:8])
            push!(state.edits, Edit(offset:offset, "\n" * " "^indent))

            # wrap arguments
            wrap_at_comma(x.args[9:end], state, indent, offset, indent)
        end
    end
end

function indent_definition(x, state)
    if x isa CSTParser.EXPR && x.typ == CSTParser.FunctionDef
        # figure out how much to indent by looking at where the expr starts
        line = findlast(y -> state.offset >= y[2], state.lines) # index, not the actual number
        line_indent, line_offset = state.lines[line]
        expr_indent = state.offset - line_offset
        indent = expr_indent + x.args[1].fullspan + sum(x->x.fullspan, x.args[2].args[1:2])

        if length(x.args[2]) > 2    # non-empty args
            offset = state.offset + x.args[1].fullspan + sum(x->x.fullspan, x.args[2].args[1:2])
            wrap_at_comma(x.args[2].args[3:end-1], state, indent, offset, indent)
        end
    end
end

#
# Main application
#

using CUDAapi

function process(name, headers...; kwargs...)
    output_file, common_file = wrap(name, headers...; kwargs...)

    for file in (output_file, common_file)
        text = read(file, String)

        ## rewriting passes

        state = State(0, Edit[])
        ast = CSTParser.parse(text, true)

        state.offset = 0
        pass(ast, state, rewrite_ccall)

        state.offset = 0
        pass(ast, state, (x,state)->rewrite_pointers(x,state,headers))

        state.offset = 0
        pass(ast, state, rewrite_runtime)

        # apply
        state.offset = 0
        sort!(state.edits, lt = (a,b) -> first(a.loc) < first(b.loc), rev = true)
        for i = 1:length(state.edits)
            text = apply(text, state.edits[i])
        end

        ## indentation passes

        lines = get_lines(text)
        state = IndentState(0, lines, Edit[])
        ast = CSTParser.parse(text, true)

        state.offset = 0
        pass(ast, state, indent_definition)

        state.offset = 0
        pass(ast, state, indent_ccall)

        # apply
        state.offset = 0
        sort!(state.edits, lt = (a,b) -> first(a.loc) < first(b.loc), rev = true)
        for i = 1:length(state.edits)
            text = apply(text, state.edits[i])
        end

        write(file, text)
    end

    ## manual patches

    patchdir = joinpath(@__DIR__, "patches")
    for entry in readdir(patchdir)
        if startswith(entry, name) && endswith(entry, ".patch")
            path = joinpath(patchdir, entry)
            run(`patch -p1 -i $path`)
        end
    end

    ## move to destination

    for src in (output_file, common_file)
        dst = joinpath(dirname(dirname(@__DIR__)), "src", name[3:end], src)
        cp(src, dst; force=true)
    end

    return
end

function main()
    # find the CUDA includes directory
    toolkit = find_toolkit()
    cuda_include = nothing
    for dir in toolkit
        cuda_include = joinpath(dir, "include")
        isdir(cuda_include) && break
    end
    isdir(cuda_include) || error("Could not find an 'include' directory in any of the toolkit directories ($(join(toolkit, ", ", " or ")))")
    @info "Found CUDA include directory at $cuda_include"

    function process_if_existing(name, headers...; kwargs...)
        env_var = uppercase(name) * "_INCLUDE"
        include = get(ENV, env_var, cuda_include)
        for header in headers
            if !isfile(joinpath(include, header))
                @error "Cannot wrap $name: could not find $header in $include (override the include directory with the $env_var environment variable)"
                return
            end
        end

        headers = map(header->joinpath(include, header), headers)
        process(name, headers...; kwargs...)
    end

    process_if_existing("cublas", "cublas_v2.h", "cublasXt.h";
                        wrapped_headers=["cublas_v2.h", "cublas_api.h", "cublasXt.h"],
                        defines=["CUBLASAPI"=>""])

    process_if_existing("cufft", "cufft.h")

    process_if_existing("curand", "curand.h")

    process_if_existing("cusparse", "cusparse.h")

    process_if_existing("cusolver", "cusolverDn.h", "cusolverSp.h";
                        wrapped_headers=["cusolver_common.h", "cusolverDn.h", "cusolverSp.h"])

    process_if_existing("cudnn", "cudnn.h")

    process_if_existing("cutensor", "cutensor.h";
                        wrapped_headers=["cutensor.h", "cutensor/types.h"])
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end


text:
module TestDeprecated
    using Compat, Compat.Test, DataFrames
    using DataFrames: Index

    i = Index()

    push!(i, :A)
    push!(i, :B)

    inds = Any[1.0, true, false,
               Any[1, missing], [1, missing],
               [true, missing], Any[true, missing],
               [:A, missing], Any[:A, missing],
               1.0:1.0, [1.0], Any[1.0]]
    for ind in inds
        if ind == :A || ndims(ind) == 0
            @test i[ind] == 1
        else
            @test (i[ind] == [1])
        end
    end

    # deprecation warning for automatically generating dup column name from indexing
    x = DataFrame(a = [1, 2, 3], b = [4, 5, 6])
    v = DataFrame(a = [5, 6, 7], b = [8, 9, 10])
    z = vcat(v, x)
    z2 = z[:, [1, 1, 2]]
    @test names(z2) == [:a, :a_1, :b]

    # TODO: uncomment the line below after deprecation, and move to dataframe.jl
    # @test_throws ArgumentError z[:, [1, 1, 2]]

    # old sort(df; cols=...) syntax
    df = DataFrame(a=[1, 3, 2], b=[6, 5, 4])
    @test sort(df; cols=[:a, :b]) == DataFrame(a=[1, 2, 3], b=[6, 4, 5])
    sort!(df; cols=[:a, :b])
    @test df == DataFrame(a=[1, 2, 3], b=[6, 4, 5])

    df = DataFrame(a=[1, 3, 2], b=[6, 5, 4])
    @test sort(df; cols=[:b, :a]) == DataFrame(a=[2, 3, 1], b=[4, 5, 6])
    sort!(df; cols=[:b, :a])
    @test df == DataFrame(a=[2, 3, 1], b=[4, 5, 6])
end


text:
## Pendulum
T = Float64

# model
model = TrajectoryOptimization.Dynamics.car
n = model.n; m = model.m
model_d = TO.discretize_model(model,:rk4)

# cost
Q = Array(1e-3*Diagonal(I,n))
R = Array(1e-3*Diagonal(I,m))
Qf = Array(Diagonal(I,n)*100.0)
x0 = zeros(n)
xf = [0.0;1.0;0.0]

# problem
N = 31
U = [ones(m) for k = 1:N-1]
obj = TrajectoryOptimization.LQRObjective(Q,R,Qf,xf,N)
dt = 0.15
prob = Problem(model_d,obj,U,dt=dt,x0=x0)
rollout!(prob)

## unconstrained
#bp
opts_ilqr = iLQRSolverOptions{T}()
ilqr_solver = iLQRSolver(prob,opts_ilqr)
jacobian!(prob,ilqr_solver)
TO.cost_expansion!(prob,ilqr_solver)
ΔV = TO._backwardpass!(prob,ilqr_solver)

#sqrt bp
opts_ilqr_sqrt = iLQRSolverOptions{T}(square_root=true)
ilqr_solver_sqrt = iLQRSolver(prob,opts_ilqr_sqrt)
jacobian!(prob,ilqr_solver_sqrt)
TO.cost_expansion!(prob,ilqr_solver_sqrt)
ΔV_sqrt = TO._backwardpass_sqrt!(prob,ilqr_solver_sqrt)
@test isapprox(ΔV_sqrt,ΔV)
@test all(isapprox(ilqr_solver.K,ilqr_solver_sqrt.K))
@test all(isapprox(ilqr_solver.d,ilqr_solver_sqrt.d))
for k = 1:N
    @test isapprox(ilqr_solver.S[k].xx,ilqr_solver_sqrt.S[k].xx'*ilqr_solver_sqrt.S[k].xx)
    @test isapprox(ilqr_solver.S[k].x,ilqr_solver_sqrt.S[k].x)
end

## constrained
opts_al = AugmentedLagrangianSolverOptions{T}(opts_uncon=opts_ilqr)
opts_al_sqrt = AugmentedLagrangianSolverOptions{T}(opts_uncon=opts_ilqr_sqrt)

# constraints
u_bnd = 5.
bnd = BoundConstraint(n,m,u_min=-u_bnd,u_max=u_bnd,trim=true)
bnd
goal_con = goal_constraint(xf)
con = [bnd, goal_con]
prob = update_problem(prob,constraints=Constraints(con,N))
rollout!(prob)

#bp
solver_ilqr = iLQRSolver(prob,opts_ilqr)
solver_al = AugmentedLagrangianSolver(prob,opts_al)
prob_al = TO.AugmentedLagrangianProblem(prob,solver_al)
TO.update_constraints!(prob_al.obj.C,prob_al.obj.constraints, prob.X, prob.U)
TO.update_active_set!(prob_al.obj)
jacobian!(prob_al,solver_ilqr)
TO.cost_expansion!(prob_al,solver_ilqr)
ΔV = TO._backwardpass!(prob_al,solver_ilqr)

#bp sqrt
solver_al_sqrt = AugmentedLagrangianSolver(prob,opts_al_sqrt)
solver_ilqr_sqrt = solver_al_sqrt.solver_uncon
prob_al_sqrt = TO.AugmentedLagrangianProblem(prob,solver_al_sqrt)
TO.update_constraints!(prob_al_sqrt.obj.C,prob_al_sqrt.obj.constraints, prob_al_sqrt.X, prob_al_sqrt.U)
TO.update_active_set!(prob_al_sqrt.obj)
jacobian!(prob_al_sqrt,solver_ilqr_sqrt)
TO.cost_expansion!(prob_al_sqrt,solver_ilqr_sqrt)
ΔV_sqrt = TO._backwardpass_sqrt!(prob_al_sqrt,solver_ilqr_sqrt)

@test isapprox(ΔV_sqrt,ΔV)
@test all(isapprox(solver_ilqr.K,solver_ilqr_sqrt.K))
@test all(isapprox(solver_ilqr.d,solver_ilqr_sqrt.d))
for k = 1:N
    @test isapprox(solver_ilqr.S[k].xx,solver_ilqr_sqrt.S[k].xx'*solver_ilqr_sqrt.S[k].xx)
    @test isapprox(solver_ilqr.S[k].x,solver_ilqr_sqrt.S[k].x)
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************

function add_boiler_tech_constraints(m, p; _n="")
    m[:TotalFuelCosts] += @expression(m,
        sum(m[:dvFuelUsage]["ExistingBoiler", ts] * p.s.existing_boiler.fuel_cost_series[ts] for ts in p.time_steps)
    )

    # Constraint (1e): Total Fuel burn for Boiler
    @constraint(m, [t in p.techs.boiler, ts in p.time_steps],
        m[:dvFuelUsage][t,ts] == p.hours_per_timestep * (
            m[Symbol("dvThermalProduction"*_n)][t,ts] / p.boiler_efficiency[t]
        )  # TODO removed p.production_factor[t,ts] * b/c all 1's for boiler; do we need it?
    )

    # Constraint (4f)-1: (Hot) Thermal production sent to storage must be less than technology's rated production
    # if !isempty(p.steam_techs)
    #     @constraint(m, [b in p.HotTES, t in p.techs.boiler, ts in p.time_steps],
    #         m[:dvProductionToStorage][b,t,ts] + m[:dvThermalToSteamTurbine][t,ts] <=
    #         p.production_factor[t,ts] * m[Symbol("dvThermalProduction"*_n)][t,ts]
    #     )
    # else
        # @constraint(m, [b in p.HotTES, t in p.techs.boiler, ts in p.time_steps],
        #     m[:dvProductionToStorage][b,t,ts] <= p.production_factor[t,ts] * m[Symbol("dvThermalProduction"*_n)][t,ts]
        # )
    # end

    # Constraint (7_heating_prod_size): Production limit based on size for boiler
    @constraint(m, [t in p.techs.boiler, ts in p.time_steps],
        m[Symbol("dvThermalProduction"*_n)][t,ts] <= m[Symbol("dvSize"*_n)][t]
    )
end

max_stars_repo_name: NREL/REoptLite

max_issues_repo_name: NREL/REoptLite

max_forks_repo_name: NREL/REoptLite


text:
#=
  From https://dtai.cs.kuleuven.be/problog/tutorial/basic/03_dice.html
  """
  We now consider yet another way to model dice, using fair ones only. This representation
  allows for convenient use of the results in arithmetic expressions, e.g., to add up the
  results from several dice. We query for the probabilities of the possible sums we can get
  from two dice given that the first number is even and the second odd.
  """

   Distributions of variable dice[1] (num:0)
   4.00000 =>    3357  (0.335700)
   2.00000 =>    3343  (0.334300)
   6.00000 =>    3300  (0.330000)
   Distributions of variable dice[2] (num:0)
   5.00000 =>    3569  (0.356900)
   1.00000 =>    3251  (0.325100)
   3.00000 =>    3180  (0.318000)
   Distributions of variable s (num:0)
   7.00000 =>    3402  (0.340200)
   9.00000 =>    2239  (0.223900)
   5.00000 =>    2114  (0.211400)
   11.00000 =>    1160  (0.116000)
   3.00000 =>    1085  (0.108500)

  See ~/blog/rolling_dice3.blog
      ~/psi/rolling_dice3.psi
=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function rolling_dice3()

    # d1 ~ DiscreteUniform(1,6)
    # d2 ~ DiscreteUniform(1,6)
    num_dice = 2
    dice = tzeros(2)
    for d in 1:num_dice
        dice[d] ~ DiscreteUniform(1,6)
    end

    s ~ Dirac(sum(dice))

    # Nice!
    function odd(d)
        return d % 2 == 1
    end

    # Nicer!
    function even(d)
        return !odd(d)
    end

    true ~ Dirac(even(dice[1]))
    true ~ Dirac(odd(dice[2]))

end

model = rolling_dice3()

num_chns = 2

# chns = sample(model, Prior(), MCMCThreads(), 10_000, num_chns)

# chns = sample(model, MH(), 10_000)
chns = sample(model, PG(15), 1_000)
# # chns = sample(model, IS(), 10_000)

chns = sample(model, SMC(), 10_000)

# chns = sample(model, HMC(0.1,10), 10_000) # nope

# chns = sample(model, NUTS(1000,0.65), MCMCThreads(), 40_000, num_chns)
# chns = sample(model, Gibbs(MH(:zlabels),NUTS(1000,0.65,:m,:b,:sigma)), MCMCThreads(), 40_000, num_chns)

display(chns)

show_var_dist_pct(chns,Symbol("dice[1]")) # Note: :"dice[1]" DON't work!
show_var_dist_pct(chns,Symbol("dice[2]"))
show_var_dist_pct(chns,:s)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
"""
    AbstractFilteringProblem{S1, S2, T1<:TimeType, T2<:TimeType}

Abstract type for a filtering problem for observations of type `S2` in `T2` and hidden states of type `S1` in `T1`.
"""
abstract type AbstractFilteringProblem{S1, S2, T1<:TimeType, T2<:TimeType} end

"""
    state_model(problem::AbstractFilteringProblem)

Return the hidden state model underlying `problem`.
"""
function state_model(problem::AbstractFilteringProblem) end

"""
    obs_model(problem::AbstractFilteringProblem)

Return the observation model underlying `problem`.
"""
function obs_model(problem::AbstractFilteringProblem) end

"""
    state_dim(problem::AbstractFilteringProblem)

Return the dimensionality of the hidden state in `problem`.
"""
state_dim(problem::AbstractFilteringProblem) = state_dim(state_model(problem))

"""
    obs_dim(problem::AbstractFilteringProblem)

Return the dimensionality of the observed state in `problem`.
"""
obs_dim(problem::AbstractFilteringProblem) = obs_dim(obs_model(problem))

state_type(problem::AbstractFilteringProblem{S1, S2, T1, T2}) where {S1, S2, T1, T2}       = S1
obs_type(problem::AbstractFilteringProblem{S1, S2, T1, T2}) where {S1, S2, T1, T2}         = S2
hidden_time_type(problem::AbstractFilteringProblem{S1, S2, T1, T2}) where {S1, S2, T1, T2} = T1
obs_time_type(problem::AbstractFilteringProblem{S1, S2, T1, T2}) where {S1, S2, T1, T2}    = T2

# How to add new filtering problem:
# * Add a struct which is a subtype of AbstractFilteringProblem{S1, S2, T1, T2}, where S1 and S2 are the types of hidden states and observations respectively, and T1 and T2 are the time types of the hidden and observed processes, respectively.
# * Implement methods for
#   - state_model
#   - obs_model

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
@doc raw"""
```
Shah(slant::Real = Inf, tilt::Real = Inf, iterations::Int = 200)
```
Creates a `Shah` fucntor which attempts to produce a heightmap from a grayscale
image using Shah's algorithm.

Under the assumption that the albedo is constant, and the surface is Lambertian,
the algorithm employs discrete approximations of p and q using finite
differences in order to linearize the reflectance map in terms of Z by taking
the first term of its Taylor series expansion.
# Output
Returns an M by N array (matching dimensions of original image) of Float `Z`
that represents the reconstructed height at the point and the gradients in
M by N arrays of Float 'p' and 'q'.
# Details
Given the reflectance map for the surface defined as bellow:
```math
R(x,y)=\dfrac{-i_xp-i_yq+i_z}{\sqrt{1+p^2+q^2}}=\dfrac{\cos\sigma+p\cos\tau\sin
\sigma+q\sin\tau\sin\sigma}{\sqrt{1+p^2+q^2}}
```
where ``i_x=\frac{I(1)}{I(3)}=\frac{\cos\tau\sin\sigma}{\cos\sigma}=\cos\tau
\tan\sigma`` and ``i_y=\frac{I(2)}{I(3)}=\frac{\sin\tau\sin\sigma}{\cos\sigma}
=\sin\tau\tan\sigma``.

and p and q are discretely approximated as:
```math
\begin{gathered}
p=Z(x,y)-Z(x-1,y)\\
q=Z(x,y)-Z(x,y-1)\\
\end{gathered}
```

Shah linearized the function ``f=E-R=0`` in terms of ``Z`` in the vicinity of
``Z^{k-1}`` by crating a system of linear equations which can be solved iteratively
using the Jacobi iterative scheme, simplifying the Taylor series expansion to the
first order to get the following:
```math
f(Z(x,y))=0\approx f(Z^{n-1}(x,y))+(Z(x,y)-Z^{n-1}(x,y))\dfrac{df(Z^{n-1}(x,y))}
{dZ(x,y)}
```
which by letting ``Z^n(x,y)=Z(x,y)`` gives:
```math
Z^n(x,y=Z^{n-1}(x,y)-\dfrac{f(Z^{n-1}(x,y))}{\dfrac{df(Z^{n-1}(x,y))}{dZ(x,y)}}
```
where,
```math
\dfrac{df(Z^{n-1}(x,y))}{dZ(x,y)}=\dfrac{(p+q)(pi_x+qi_y+1)}{\sqrt{(1+p^2+q^2)^3}
\sqrt{1+i_x+i_y}}-\dfrac{i_x+i_y}{\sqrt{1+p^2+q^2}\sqrt{1+i_x+i_y}}
```
which as ``Z^0(x,y)=0``, allows the algorithm to iteratively solve for ``Z(x,y)``.

# Arguments
The function arguments are described in more detail below.
## `iterations`
An `Int` that specifies the number of iterations the algorithm is to perform. If
left unspecified a default value of 200 is used.
## `slant`
A `Real` that specifies the slant value to be used by the algorithm. The `slant`
should be a value in the range [0,π/2]. If `slant` is specified to must the `tilt`.
Will defult to `Inf` which will trigger `estimate_img_properties` to be run.
## `tilt`
A `Real` that specifies the tilt value to be used by the algorithm. The `tilt`
should be a value in the range [0,2π]. If `tilt` is specified to must
the `slant`. Will defult to `Inf` which will trigger `estimate_img_properties` to
be run.
!!! note
    If `slant` and `tilt` are not defined (i.e. have defulted to `Inf`) they
    will be calculated at runtime using `estimate_img_properties`.
# Example
Compute the heightmap for a synthetic image generated by `generate_surface`.
```julia
using Images, Makie, ShapeFromShading

#generate synthetic image
img = generate_surface(SynthSphere(50), 1.0, [0.1,0,0.9])

#calculate the heightmap
shah = Shah(slant = 0.1106572211739, tilt = 0.0)
Z = shah(img)

#normalize to maximum of 1 (not necessary but makes displaying easier)
Z = Z./maximum(Z)

#display using Makie (Note: Makie can often take several minutes first time)
r = 0.0:0.1:2
surface(r, r, Z)
```
# Reference
1. T. Ping-Sing and M. Shah, "Shape from shading using linear approximation", Image and Vision Computing, vol. 12, no. 8, pp. 487-498, 1994. [doi:10.1016/0262-8856(94)90002-7](https://doi.org/10.1016/0262-8856(94)90002-7)
"""
function (algorithm::Shah)(img::AbstractArray)
    σ, τ = algorithm.slant, algorithm.tilt
    if σ == Inf || τ == Inf
        ρ, I, σ, τ = estimate_img_properties(img)
    end
    iterations = algorithm.iterations

    E = Float64.(img)
    E = E .* 255
    M, N = size(E)
    p = zeros(axes(E))
    q = zeros(axes(E))
    Z = zeros(axes(E))
    Zx = zeros(axes(E))
    Zy = zeros(axes(E))
    ix = cos(τ) * tan(σ)
    iy = sin(τ) * tan(σ)
    R = zeros(axes(E))
    δfδZ = zeros(axes(E))
    f = zeros(axes(E))
    @inbounds for i = 1:iterations
        #calculate reflectance map
        for i in CartesianIndices(R)
            R[i] = (cos(σ) + p[i] * cos(τ) * sin(σ) + q[i] * sin(τ) * sin(σ)) /
                sqrt(1 + p[i]^2 + q[i]^2)

            R[i] = max(0, R[i])
            f[i] = E[i] - R[i]
        end

        #calculate derivative of f in respect to Z and update Z
        for i in CartesianIndices(δfδZ)
            δfδZ[i] = (p[i] + q[i]) * (ix * p[i] + iy * q[i] + 1) / (sqrt((1
                + p[i]^2 + q[i]^2)^3) * sqrt(1 + ix^2 + iy^2)) - (ix + iy) /
                (sqrt(1 .+ p[i]^2 + q[i]^2) * sqrt(1 + ix^2 + iy^2))

            Z[i] = Z[i] - f[i] / (δfδZ[i] + eps())
        end

        #update surface normals
        for i = 2:M
            for j = 1:N
                Zx[i,j] = Z[i-1,j]
            end
        end
        for i = 2:N
            for j = 1:M
                Zy[j,i] = Z[j,i-1]
            end
        end
        for i in CartesianIndices(p)
            p[i] = Z[i] - Zx[i]
            q[i] = Z[i] - Zy[i]
        end
    end

    #smooth Z using a 21X21 median filter of the absolute values
    for i in CartesianIndices(Z)
        Z[i] = abs(Z[i])
    end
    Z = mapwindow(median!, Z, (21,21))
    return Z
end


text:
using SamIO
using Base.Test
using GZip

## BamReader

# test forward reads
reader = BamReader("data/small.bam", :forward, SamIO.assembly["GRCh38_alt"])
@test position(reader) == 10544
advance!(reader)
@test position(reader) == 11247
close(reader)

# test reverse reads
reader = BamReader("data/small.bam", :reverse, SamIO.assembly["GRCh38_alt"])
@test position(reader) == 10056
advance!(reader)
@test position(reader) == 10102
close(reader)

# test reading through a whole file
reader = BamReader("data/small.bam", :any, SamIO.assembly["GRCh38_alt"])
lastPos = -1
while !eof(reader)
	lastPos = position(reader)
	advance!(reader)
end
close(reader)

# test reading through a whole file with an iterator
reader = BamReader("data/small.bam", :reverse, SamIO.assembly["GRCh38_alt"])
for pos in eachposition(reader)
end
close(reader)

# test with a different assembly
reader = BamReader("data/small.bam", :forward, SamIO.assembly["GRCh38"])
@test position(reader) == 10544
advance!(reader)
@test position(reader) == 11247
close(reader)

# test with a different assembly
reader = BamReader("data/mm10.bam", :any, SamIO.assembly["mm10"])
@test position(reader) == 3007701
advance!(reader)
@test position(reader) == 3011073
close(reader)

## BedReader

# test first few positions
reader = BedReader("data/small.bed", SamIO.assembly["GRCh38_alt"])
@test position(reader) == 1001
advance!(reader)
@test position(reader) == 1002
advance!(reader)
@test position(reader) == 5001
close(reader)

# test reading through a whole file as a stream
f = open("data/small.bed")
reader = BedReader(f, SamIO.assembly["GRCh38_alt"])
lastPos = -1
while !eof(reader)
	lastPos = position(reader)
	advance!(reader)
end
close(reader)

# test reading through a whole file as a zipped stream
f = GZip.open("data/small.bed.gz")
reader = BedReader(f, SamIO.assembly["GRCh38_alt"])
lastPos = -1
while !eof(reader)
	lastPos = position(reader)

	advance!(reader)
end
close(reader)

# test reading through a whole file with an iterator
reader = BedReader("data/small.bed", SamIO.assembly["GRCh38_alt"])
for pos in eachposition(reader)
end
close(reader)

## BinningMap

reader = BinningMap(BamReader("data/small.bam", :any, SamIO.assembly["GRCh38_alt"]), 1000)
i = 1
while !eof(reader)
	if i == 1
		@test position(reader) == 11
		@test value(reader) == 4.0
	end
    advance!(reader)
	i += 1
end
close(reader)

reader = BinningMap(BedReader("data/small.bed", SamIO.assembly["GRCh38_alt"]), 1000)
i = 1
while !eof(reader)
	if i == 1
		@test position(reader) == 2
		@test value(reader) == 2.0
	end
    advance!(reader)
	i += 1
end
close(reader)