text:
# Preprocess.jl
"""
    isInputMatrixProper()

Determine whether or not the input matrix can fulfill the minimum non-musical requirements of an all-partition array.
"""
function isInputMatrixProper(inputMatrix)
    print("\n", "Begin preprocessing...")
    numRows = size(inputMatrix, 1)
    numCols = size(inputMatrix, 2)
    if isless(numCols, numRows)
        throw(DimensionMismatch("input matrix must have fewer rows than columns"))
    end
    if numRows > 12 || numRows < 4
        throw(DimensionMismatch("input matrix must have more than 3 rows and fewer than 13"))
    end
    if numCols > 96 || numCols < 72
        throw(DimensionMismatch("input matrix must have more than 71 columns and fewer than 97"))
    end
    aggregate = sort(unique(inputMatrix))
    numPcs = length(aggregate)
    if numPcs != 12
        throw(ErrorException("number of distinct elements in input matrix must be equal to 12"))
    end
    if !isequal(aggregate, [0,1,2,3,4,5,6,7,8,9,10,11])
        throw(ErrorException("elements of input matrix must be in the range [0, 11]"))
    end
    counts = zeros(numPcs)
    for i in 1:numPcs
        counts[i] = count(isequal(aggregate[i]), inputMatrix)
    end
    if !all(counts[1] .== counts)
        throw(ErrorException("input matrix must contain an equal number of each element in the range [0, 11]"))
    end
    aggregate .+= 1 # add one to match matrix adjustment done later
    numElements::Int16 = length(inputMatrix)
    numPartitions = findNumOfPartitions(numRows, numPcs)
    numOverlaps::Int16 = (numPartitions * numPcs) - (numRows * numCols)
    if isless((convert(Int16, numPartitions) * convert(Int16, numPcs)), numElements) # FIX
        throw(ErrorException("input matrix must contain a number of elements equal to or fewer than the number required of its partitions"))
    end
    return APArray(numRows, numCols, numElements, numOverlaps, numPartitions, aggregate)
end

"""
    findNumOfPartitions()

Find the number of required partitions given the number of rows in the input matrix and the size of the aggregate.
"""
function findNumOfPartitions(numRows, numPcs)
    numPartitions = zero(eltype(numRows))
    for i in 1:numRows
        numPartitions += length(collect(partitions(numPcs, i)))
    end
    return numPartitions
end

"""
    requestUserInputToSelectHeuristics()

Halts the continuation of the program until the user provides appropriate input to either continue or stop the program.
If the user opts to continue, they can select whether or not to use searhc heuristics. If they choose to use heuristics,
they can select one of two options.
"""
function requestUserInputToSelectHeuristics(arraySpecs, reservedPartitions, allCompositions)
    userChoiceOffset = zero(Int8)
    userChoiceHeuristic = Array{Int8,1}()
    userChoiceReserves = zero(Int8)
    println("The generateAllPartitionArray() algorithm will search through the input matrix and attempt to generate an all-partition array with the following specifications:","\n")
    printInputMatrixSummary(arraySpecs, reservedPartitions, allCompositions) # print details of input matrix and its an all-partition array
    # ask for user input to begin preprocessing and searching for a solution
    userInput = ""
    while (userInput == "")
        println("")
        userInput = getUserInput("This process can include the use of heuristics that work globally to help guide the search. To select from these heuristics and begin searching, press 'y' then RETURN. To terminate the program now, press 'n' then RETURN.")
        #userInput = parse(AbstractString, getUserInput("This process can include several heuristics that work globally to help guide the search. To select from these heuristcs and begin searching, press 'y' then RETURN. To terminate the program now, press 'n' then RETURN."))
        if isequal(userInput, "n")
            return exit()
        elseif !isequal(userInput, "y")
            println("Incorrect key entered", "\n") # use ErrorException?
            userInput = ""
        end
    end
    # ask for user input for selecting a heuristic or not
    userInput = ""
    while (userInput == "")
        println("")
        # ask for user input for selecting offset
        userInput = getUserInput("Input ONE of the following numbers: '0' for no search heuristic (not recommended), or '1' for use of a search heuristic. Afterwards, press RETURN.")
        userChoiceHeuristic = [parse(Int8,i) for i in userInput]
        if length(userChoiceHeuristic) == 0
            println("Incorrect key entered", "\n")
            userInput = ""
        end
    end
    if userChoiceHeuristic[1] != 0
        println("\n", "Select the options below if you are using any of the following input matrix types:","\n")
        println("Babbitt(6,96):","\t","1 search heuristic and ", 13, " (i.e., maximum)", " reserved partitions","\n")
        println("Smalley(6,96):","\t","0 search heuristic and ", 22, " (i.e., maximum)", " reserved partitions","\n")
        println("Babbitt(4,96):","\t","1 search heuristic and ", 1, " reserved partitions","\n")

        # ask for user input for selecting the search heuristic
        userInput = ""
        while (userInput == "")
            println("")
            userInput = getUserInput("Input ONE of the following numbers for judging candidates with respect to an 'ideal position': '0' for only the magnitude of difference, or '1' for the individual magnitude and direction of difference. Afterwards, press RETURN.")
            userChoiceOffset = parse(Int8, userInput)
            if length(userChoiceOffset) != 1
                println("Input must be a single number", "\n")
                userInput = ""
            end
            if length(userChoiceOffset) == 0
                println("Incorrect key entered", "\n")
                userInput = ""
            end
        end
    end
    # ask for user input for selecting reserved partitions
    userInput = ""
    while (userInput == "")
        println("")
        userInput = getUserInput("Input the number of final region partitions––from 0 to the number of possible reserved partitions (above) to reserve in the search. Afterwards, press RETURN.")
        userChoiceReserves = parse(Int8, userInput)
        if userChoiceReserves > length(reservedPartitions) || userChoiceReserves < 0
            println("Input must be a number between 0 and ", length(reservedPartitions),"\n")
            userInput = ""
        end
    end
    return Heuristic(userChoiceOffset, userChoiceHeuristic, userChoiceReserves)
end

"""
    printInputMatrixSummary()

"""
function printInputMatrixSummary(arraySpecs, reservedPartitions, allCompositions) # print more info
    println("Input matrix size:", "\t", "\t", "\t", "\t", arraySpecs.numRows, " x ", arraySpecs.numCols)
    println("Aggregate size:", "\t", "\t", "\t", "\t", "\t", length(arraySpecs.aggregate))
    println("Num. of required partitions:", "\t", "\t", "\t", arraySpecs.numPartitions, "\t", "(from ", size(allCompositions, 1), " compositions)")
    println("Num. of possible reserved partitions:", "\t", "\t", length(reservedPartitions))
    println("Num. of require overlaps:", "\t", "\t", "\t", arraySpecs.numOverlaps)
    #println("Summary: Input matrix is proper.", "\n")
    return nothing
end

"""
    getUserInput()

Prompt user for input and return user input as a string.
"""
function getUserInput(prompt::AbstractString="")
    println(prompt * " ")
    convert(String, chomp(readline()))
end

"""
    computeListOfAllCompositions()

Computes a list of all compositions needed for the given input matrix.
Also finds the sum of the largest and smallest parts of all partitions and creates a helper dictionary variable for searching all compositions.
"""
function computeListOfAllCompositions(arraySpecs)
    lexCompositions = computeLexicalOrderCompostitions(length(arraySpecs.aggregate), arraySpecs.numRows)
    allCompositions, dictPartComps, minMaxParts, probabilities = groupCompositionsByPartition(lexCompositions, arraySpecs.numRows, arraySpecs.numPartitions)
    return allCompositions, dictPartComps, minMaxParts, probabilities
end

"""
    computeLexicalOrderCompostitions()

Computes a lexicographically ordered list of all compositions with n parts or fewer that sum to k, where n is the number
of input matrix rows and k is the number of pitch classes in an aggregate.
"""
function computeLexicalOrderCompostitions(numPcs, numRows)
    # k: sum, n: number of non-negative integers (i.e., n sum k)
    m = binomial(numPcs + numRows - 1, numRows - 1)
    d1 = zeros(eltype(m), m, 1)
    d2 = collect(combinations(collect(1:numPcs + numRows - 1), numRows - 1))
    d2 = transpose(hcat(d2...)) #
    d3 = fill!(ones(eltype(m), m, 1), numPcs + numRows)
    dividers = hcat(d1, d2, d3)
    return convert(Array{Int8,2}, diff(dividers, dims=2) .- 1)
end

"""
    groupCompositionsByPartition()

Groups compositions by partition ordered lexicographically. Returns this list in addition to some helper variables.
"""
function groupCompositionsByPartition(lexCompositions, numRows, numPartitions)
    allCompositions = zeros(eltype(lexCompositions), size(lexCompositions, 1), size(lexCompositions, 2)+1)
    dictPartComps = Dict{eltype(lexCompositions), UnitRange}()
    minMaxParts = zeros(Int16, numPartitions, 2)
    uniqueMatrix = unique(sort(lexCompositions, dims=2), dims=1)
    probabilities = zeros(Float16, numPartitions)
    start = 1
    for i in 1:numPartitions
        testComp = uniqueMatrix[i, :]
        indices = Array{Int,1}()
        for j in 1:size(lexCompositions, 1)
            if testComp == sort(lexCompositions[j, :])
                push!(indices, j)
            end
        end
        allCompositions[start:start+length(indices)-1, 1:numRows] = lexCompositions[indices, 1:numRows] # store composition
        allCompositions[start:start+length(indices)-1, end] .= i # store partition number
        dictPartComps[i] = start:start+length(indices)-1 # create dictionary entry of key: partition number, value: composition index range
        minMaxParts[i, 1] = first(testComp) # store smallest part
        minMaxParts[i, 2] = last(testComp) # store largest part
        # Compute probabilities here
        #probabilities[i] = length(indices) / size(allCompositions, 1) # number of compositions
        probabilities[i] = length(testComp) - count(iszero, testComp)  # number of zeros in composition
        start = start + length(indices)
    end
    # normalize probabilities
    probabilities = (probabilities .- minimum(probabilities))/(maximum(probabilities) - minimum(probabilities))
    return allCompositions, dictPartComps, minMaxParts, probabilities
end

"""
    computeListOfAllCombinationsOfOverlaps()

Find all possible overlaps for a composition with at most a number of parts equal to the number of rows
in the input matrix. The number of possible overlaps is equal to all combinations of the number of matrix rows.
"""
function computeListOfAllCombinationsOfOverlaps(numRows)
    indices = collect(Array{Int16,1}, combinations(1:numRows))
    allOverlaps = falses(length(indices) + 1, numRows)
    # rowsOverlapsByNum = [1, 7, 22, 42, 57, 63] # e.g., greater than x[1] is 1 or more overlaps
    rowsOverlapsByNum = Array{Int16,1}()
    setSize = 0
    for i::Int16 in 2:length(indices) + 1
        allOverlaps[i, indices[i - 1]] .= true
        if length(indices[i - 1]) != setSize
            push!(rowsOverlapsByNum, i - one(eltype(rowsOverlapsByNum)))
            setSize = length(indices[i - 1])
        end
    end
    overlapIndicesByRow = findOverlapRowIndices(allOverlaps)
    return allOverlaps, rowsOverlapsByNum, overlapIndicesByRow
end

"""
    findOverlapRowIndices()

Find all rows of overlaps that share one overlap and return a list of the indices of these rows in
the set of all combinations of overlaps (one list for each input matrix row).
"""
function findOverlapRowIndices(allOverlaps)
    # group row indices of overlaps by shared rows (12-row input matrix: max 4095)
    overlapIndicesByRow = zeros(Int16, Int(size(allOverlaps, 1) / 2), size(allOverlaps, 2))
    for i in 1:size(allOverlaps, 2)
        overlapIndicesByRow[:, i] = findall(allOverlaps[:, i])
    end
    return overlapIndicesByRow
end

"""
    findPartitionsToReserveForLast()

Finds all candidate partitions at the righthand side of the imput matrix (from right to left) corresponding
to the possible partitions for the final position in the primary left-to-right search. These partitions can
be excluded from the primary search so that at least one remains available for selection for the final position.
"""
function findPartitionsToReserveForLast(inputMatrix, allCompositions, allOverlaps, overlapIndicesByRow, arraySpecs, minMaxParts, rowsOverlapsByNum, probabilities)
    matrix = reverse(inputMatrix, dims=2) # work backwards from matrix
    matrix .+= one(eltype(inputMatrix))
    position = ones(eltype(inputMatrix), arraySpecs.numRows)   # matrix columns for each row
    cnt = one(eltype(inputMatrix))
    usedPartitions = zeros(Int, arraySpecs.numPartitions)
    sumUnUsedParts = computeSumOfUnUsedParts(usedPartitions, minMaxParts, cnt)
    idealRegion::Float16 = round((arraySpecs.numRows * (arraySpecs.numCols + 1)) / arraySpecs.numPartitions, digits=3)
    chosenHeuristics = Heuristic(0, Array{Int8,1}([3]), 0)

    # find final position compositions
    # two dummy variables (x, y)
    candidateList, x, y = findCandidateCompositionsAtPosition(matrix, position, cnt, allCompositions, allOverlaps, overlapIndicesByRow, arraySpecs, minMaxParts, sumUnUsedParts, rowsOverlapsByNum, chosenHeuristics, idealRegion, probabilities)

    # find partitions from the final position compositions
    reservedLocs = falses(arraySpecs.numPartitions)
    for i in 1:length(candidateList)
        reservedLocs[candidateList[i].partition] = true
    end
    reservedPartitions::Array{eltype(inputMatrix),1} = 1:arraySpecs.numPartitions
    reservedPartitions = reservedPartitions[reservedLocs]

    #= Request user input to select heuristics used in searching =#
    print("[complete]")
    println("\n")
    chosenHeuristics = requestUserInputToSelectHeuristics(arraySpecs, reservedPartitions, allCompositions) # user-selected heuristics

    return reservedPartitions, chosenHeuristics
end

#= End of file =#


text:
using NaNMath
using Base.Test

@test isnan(NaNMath.log(-10))
@test isnan(NaNMath.log1p(-100))
@test isnan(NaNMath.pow(-1.5,2.3))
@test isnan(NaNMath.sqrt(-5))
@test NaNMath.sqrt(5) == Base.sqrt(5)
@test NaNMath.sum([1., 2., NaN]) == 3.0
@test isnan(NaNMath.sum([NaN, NaN]))
@test NaNMath.sum(Float64[]) == 0.0
@test NaNMath.sum([1f0, 2f0, NaN32]) === 3.0f0
@test NaNMath.maximum([1., 2., NaN]) == 2.0
@test NaNMath.minimum([1., 2., NaN]) == 1.0
@test NaNMath.mean([1., 2., NaN]) == 1.5
@test NaNMath.var([1., 2., NaN]) == 0.5
@test NaNMath.std([1., 2., NaN]) == 0.7071067811865476

ext: jl


text:
import ..UncertainValues.AbstractUncertainValue
import ..SamplingConstraints:
    SamplingConstraint,
	NoConstraint,
	TruncateLowerQuantile,
	TruncateUpperQuantile,
	TruncateQuantiles,
	TruncateMinimum,
	TruncateMaximum,
	TruncateRange,
    TruncateStd
import Distributions: truncated

########################################################################
# Resampling with constraints
########################################################################

"""
Resample an uncertain value by first truncating its furnishing distribution with the
provided sampling `constraint`.
"""
function resample(uv::AbstractUncertainValue, constraint::SamplingConstraint) end

"""
    resample(uv::AbstractUncertainValue, constraint::NoConstraint)

Resample an uncertain value without contraints (use the full furnishing distribution).

## Example

```julia
uncertainval = UncertainValue(0, 0.2, Normal)

# Resample the uncertain value by resampling the full distribution once.
resample(uncertainval, NoConstraint())
```
"""
resample(uv::AbstractUncertainValue, constraint::NoConstraint) = resample(uv)

"""
    resample(uv::AbstractUncertainValue, constraint::NoConstraint, n::Int)

Resample an uncertain value without contraints (use the full furnishing distribution).

## Example

```julia
uncertainval = UncertainValue(0, 0.2, Normal)

# Resample the uncertain value by resampling the full distribution 1000 times.
resample(uncertainval, NoConstraint(), 1000)
```
"""
resample(uv::AbstractUncertainValue, constraint::NoConstraint, n::Int) = resample(uv, n)

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateUpperQuantile)

Resample by first truncating the distribution representing the value at a lower
quantile, then performing the resampling.

## Example

```julia
uncertainval = UncertainValue(0, 0.2, Normal)
constraint = TruncateLowerQuantile(0.16)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution once.
resample(uncertainval, constraint)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateLowerQuantile)
    # Apply (another level of) truncation, then sample
    lower_bound = quantile(uv.distribution, constraint.lower_quantile)
    upper_bound = maximum(uv)
    rand(truncated(uv.distribution, lower_bound, upper_bound))
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateLowerQuantile, n::Int)

Resample by first truncating the distribution representing the value at a lower
quantile, then performing the resampling.

## Example

```julia
uncertainval = UncertainValue(0, 0.2, Normal)
constraint = TruncateLowerQuantile(0.16)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 1000 times.
resample(uncertainval, constraint, 1000)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateLowerQuantile, n::Int)
    # Apply (another level of) truncation, then sample
    lower_bound = quantile(uv.distribution, constraint.lower_quantile)
    upper_bound = maximum(uv)
    rand(truncated(uv.distribution, lower_bound, upper_bound), n)
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateUpperQuantile)

Resample by first truncating the distribution representing the value at an upper
quantile, then performing the resampling.

## Example

```julia
uncertainval = UncertainValue(0, 0.2, Normal)
constraint = TruncateUpperQuantile(0.8)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution once.
resample(uncertainval, constraint)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateUpperQuantile)
    # Apply (another level of) truncation, then sample
    lower_bound = minimum(uv)
    upper_bound = quantile(uv.distribution, constraint.upper_quantile)
    rand(truncated(uv.distribution, lower_bound, upper_bound))
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateUpperQuantile, n::Int)

Resample by first truncating the distribution representing the value at an upper
quantile, then performing the resampling.

## Example

```julia
uncertainval = UncertainValue(0, 0.2, Normal)
constraint = TruncateUpperQuantile(0.8)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 1000 times.
resample(uncertainval, constraint, 1000)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateUpperQuantile, n::Int)
    # Apply (another level of) truncation, then sample
    lower_bound = minimum(uv)
    upper_bound = quantile(uv.distribution, constraint.upper_quantile)
    rand(truncated(uv.distribution, lower_bound, upper_bound), n)
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateQuantiles)

Resample by first truncating the distribution representing the value at a set of
qunatiles, then performing the resampling.

## Example

```julia
uncertainval = UncertainValue(0, 1, Uniform)
constraint = TruncateLowerQuantile(0.2)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution once.
resample(uncertainval, constraint)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateQuantiles)
    # Apply (another level of) truncation, then sample
    lower_bound = quantile(uv.distribution, constraint.lower_quantile)
    upper_bound = quantile(uv.distribution, constraint.upper_quantile)
    rand(truncated(uv.distribution, lower_bound, upper_bound))
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateQuantiles, n::Int)

Resample by first truncating the distribution representing the value at a set of
qunatiles, then performing the resampling.

## Example

```julia
uncertainval = UncertainValue(0, 1, Uniform)
constraint = TruncateLowerQuantile(0.2)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 1000 times.
resample(uncertainval, constraint, 1000)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateQuantiles, n::Int)
    # Apply (another level of) truncation, then sample
    lower_bound = quantile(uv.distribution, constraint.lower_quantile)
    upper_bound = quantile(uv.distribution, constraint.upper_quantile)
    rand(truncated(uv.distribution, lower_bound, upper_bound), n)
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateStd, n::Int;
        n_draws::Int = 10000)

Resample by first truncating the distribution representing the value to some multiple
of its standard deviation around the mean.

## Example

```julia
uncertainval = UncertainValue(0, 0.8, Normal)
constraint = TruncateStd(1.1) # accept values only in range 1.1*stdev around the mean

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 1000 times.
resample(uncertainval, constraint, 1000)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateStd, n::Int;
        n_draws::Int = 10000)
    # Apply (another level of) truncation, then sample
    stdev = std(resample(uv, n_draws))
    m = mean(resample(uv, n_draws))
    lower_bound = m - stdev
    upper_bound = m + stdev

    rand(truncated(uv.distribution, lower_bound, upper_bound), n)
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateMinimum)

Resample by first truncating the distribution representing the value at some
minimum value, then performing the resampling.

## Example

```julia
uncertainval = UncertainValue(0, 0.2, Normal)
constraint = TruncateMinimum(-0.5) # accept no values less than -0.5

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution once.
resample(uncertainval, constraint)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateMinimum)
    # Apply (another level of) truncation, then sample
    lower_bound = constraint.min
    upper_bound = maximum(uv)
    lower_bound > upper_bound ? error("lower bound > upper_bound") : nothing
    rand(truncated(uv.distribution, lower_bound, upper_bound))
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateMinimum, n::Int)

Resample by first truncating the distribution representing the value at some
minimum value, then performing the resampling.

## Example

```julia
uncertainval = UncertainValue(0, 0.2, Normal)
constraint = TruncateMinimum(-0.5)

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 1000 times.
resample(uncertainval, constraint, 1000)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateMinimum, n::Int)
    # Apply (another level of) truncation, then sample
    lower_bound = constraint.min
    upper_bound = maximum(uv)
    lower_bound > upper_bound ? error("lower bound > upper_bound") : nothing
    rand(truncated(uv.distribution, lower_bound, upper_bound), n)
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateMaximum)

Resample by first truncating the distribution representing the value at some
minimum value, then performing the resampling.

## Example

```julia
uncertainval = UncertainValue(0, 0.8, Normal)
constraint = TruncateMaximum(1.1) # accept no values larger than 1.1

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution once.
resample(uncertainval, constraint)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateMaximum)
    # Apply (another level of) truncation, then sample
    upper_bound = constraint.max
    lower_bound = minimum(uv)
    lower_bound > upper_bound ? error("lower bound > upper_bound") : nothing
    rand(truncated(uv.distribution, lower_bound, upper_bound))
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateMaximum, n::Int)

Resample by first truncating the distribution representing the value at some
minimum value, then performing the resampling.

## Example

```julia
uncertainval = UncertainValue(0, 0.8, Normal)
constraint = TruncateMaximum(1.1) # accept no values larger than 1.1

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 1000 times.
resample(uncertainval, constraint, 1000)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateMaximum, n::Int)
    # Apply (another level of) truncation, then sample
    lower_bound = minimum(uv)
    upper_bound = constraint.max
    lower_bound > upper_bound ? error("lower bound > upper_bound") : nothing
    rand(truncated(uv.distribution, lower_bound, upper_bound), n)
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateRange)

Resample by first truncating the distribution representing the value at some
minimum value, then performing the resampling.

## Example

```julia
uncertainval = UncertainValue(0, 0.8, Normal)
constraint = TruncateRange(-0.7, 1.1) # accept values only in range [-0.7, 1.1]

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution once.
resample(uncertainval, constraint)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateRange)
    # Apply (another level of) truncation, then sample
    lower_bound = constraint.min
    upper_bound = constraint.max
    lower_bound > upper_bound ? error("lower bound > upper_bound") : nothing

    rand(truncated(uv.distribution, lower_bound, upper_bound))
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateRange, n::Int)

Resample by first truncating the distribution representing the value at some
minimum value, then performing the resampling.

## Example

```julia
uncertainval = UncertainValue(0, 0.8, Normal)
constraint = TruncateRange(-0.7, 1.1) # accept values only in range [-0.7, 1.1]

# Resample the uncertain value `1000` times by truncating the distribution furnishing it,
# then resampling the new distribution `1000` times
resample(uncertainval, constraint, 1000)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateRange, n::Int)
    # Apply (another level of) truncation, then sample
    lower_bound = constraint.min
    upper_bound = constraint.max
    lower_bound > upper_bound ? error("lower bound > upper_bound") : nothing

    rand(truncated(uv.distribution, lower_bound, upper_bound), n)
end

"""
    resample(uv::AbstractUncertainValue, constraint::TruncateStd; n_draws::Int = 10000)

Resample by first truncating the distribution representing the value to some multiple
of its standard deviation around the mean.

## Example

```julia
uncertainval = UncertainValue(0, 0.8, Normal)
constraint = TruncateStd(1.1) # accept values only in range 1.1*stdev around the mean

# Resample the uncertain value by truncating the distribution furnishing it,
# then resampling the new distribution 1000 times.
resample(uncertainval, constraint, 1000)
```
"""
function resample(uv::AbstractUncertainValue, constraint::TruncateStd; n_draws::Int = 10000)
    # Apply (another level of) truncation, then sample
    stdev = std(resample(uv, n_draws))
    m = mean(resample(uv, n_draws))
    lower_bound = m - stdev
    upper_bound = m + stdev
    lower_bound > upper_bound ? error("lower bound > upper_bound") : nothing
    rand(truncated(uv.distribution, lower_bound, upper_bound))
end

export resample


text:
using Jupyter2Pluto
using Test
using JSON

const line_markdown = "The diameter of a"
const multi_markdown = """The diameter of a pizza is often stated on a menu so let's define a **formula** to calculate the area of a pizza given the diameter **d**.

We do this by writing a formula like this: `area(d) = pi * (d/2)^2`

Let's write that below:"""

const line_markdown_wrap = "md\"$line_markdown\""
const multiline_markdown_wrap = """md\"\"\"$multi_markdown\"\"\"
"""

const code_line = """
area(d) = pi * (d / 2)^2
"""
function test_pluto_parser()
    multiline_markdown_str = "2ed4bb92-d45c-11ea-0b31-2d8e32ce7b44\n$multiline_markdown_wrap"
    line_markdown_str = "03664f5c-d45c-11ea-21b6-91cd647a07aa\n$line_markdown_wrap"

    pmark_line_cell = Jupyter2Pluto.parse_pluto_cell(line_markdown_str)
    pmark_multiline_cell = Jupyter2Pluto.parse_pluto_cell(multiline_markdown_str)
    @testset "parse pluto line markdown cell" begin
        @test pmark_line_cell.cell_id == Base.UUID("03664f5c-d45c-11ea-21b6-91cd647a07aa")
        @test pmark_line_cell.content == line_markdown
    end
    @testset "parse pluto multiline markdown cell" begin
        @test pmark_multiline_cell.cell_id == Base.UUID("2ed4bb92-d45c-11ea-0b31-2d8e32ce7b44")
        @test pmark_multiline_cell.content == multi_markdown
    end
    pcode_cell = Jupyter2Pluto.parse_pluto_cell("""14158eb0-d45c-11ea-088f-330e45412320
$code_line""")

    @testset "parse pluto code cell" begin
        @test pcode_cell.cell_id == Base.UUID("14158eb0-d45c-11ea-088f-330e45412320")
        @test pcode_cell.content == code_line
    end
    orders = Jupyter2Pluto.parse_pluto_end("""Cell order:
# ╟─03664f5c-d45c-11ea-21b6-91cd647a07aa
# ╠═14158eb0-d45c-11ea-088f-330e45412320
# ╠═2ed4bb92-d45c-11ea-0b31-2d8e32ce7b44
# ╟─e5e0a0da-d45c-11ea-1042-e9b5d0654d4f
""")
    @testset "Pluto cell order" begin
        @test orders == ["03664f5c-d45c-11ea-21b6-91cd647a07aa",
                        "14158eb0-d45c-11ea-088f-330e45412320",
                        "2ed4bb92-d45c-11ea-0b31-2d8e32ce7b44",
                         "e5e0a0da-d45c-11ea-1042-e9b5d0654d4f"]
         @test Jupyter2Pluto.order(Jupyter2Pluto.PlutoCell[pmark_line_cell, pmark_multiline_cell, pcode_cell], orders) == [pmark_line_cell, pcode_cell, pmark_multiline_cell]
        @test Jupyter2Pluto.order(Jupyter2Pluto.PlutoCell[pmark_line_cell, pmark_multiline_cell, pcode_cell], reverse(orders)) == [pmark_multiline_cell, pcode_cell, pmark_line_cell]
    end
    jmark_line_cell = Jupyter2Pluto.JupyterCell(pmark_line_cell)
    jmark_multiline_cell = Jupyter2Pluto.JupyterCell(pmark_multiline_cell)
    jcodecell = Jupyter2Pluto.JupyterCell(pcode_cell)
    @testset "Pluto to jupyter conversion"  begin
        @test jmark_line_cell.content == line_markdown
        @test jmark_multiline_cell.content == multi_markdown
        @test jcodecell.content == split(code_line, "\n")
        @test jcodecell.execution_count != 0
    end
    @testset "Jupyter to dict"  begin
        jmark_line_dict = Dict(jmark_line_cell)
        @test jmark_line_dict["cell_type"] == "markdown"
        @test isempty(jmark_line_dict["metadata"])
        @test jmark_line_dict["source"] == [line_markdown*"\n"]

        jmark_multiline_dict = Dict(jmark_multiline_cell)
        @test jmark_multiline_dict["cell_type"] == "markdown"
        @test isempty(jmark_multiline_dict["metadata"])
        @test jmark_multiline_dict["source"] == string.(split(multi_markdown, "\n")).*"\n"

        jcode_dict = Dict(jcodecell)
        @test jcode_dict["cell_type"] == "code"
        @test isempty(jcode_dict["metadata"])
        @test jcode_dict["execution_count"] != 0
        code_lines = string.(split(code_line, "\n"))
        for (index, line) in enumerate(jcode_dict["source"])
            if (line == "\n")
                @test false
            end
            if !isempty(line)
                @test code_lines[index]*"\n" == line
            else
                @test code_lines[index] == line
            end
        end
    end
end

function test_pluto2jupyter()
    test_notebook = joinpath(@__DIR__, "Basic mathematics.jl")
    pluto2jupyter(test_notebook)
    notebook = JSON.parsefile(test_notebook*".ipynb")
    @testset "check converted " begin
        @test notebook["nbformat"] == 4
        @test notebook["nbformat_minor"] == 2
        @test notebook["metadata"]["kernelspec"]["display_name"] == "Julia $(VERSION)"
        @test notebook["metadata"]["kernelspec"]["language"] == "julia"
        @test notebook["metadata"]["kernelspec"]["name"] == "julia-$(VERSION.major).$(VERSION.minor)"
        @test notebook["metadata"]["language_info"]["file_extension"] == ".jl"
        @test notebook["metadata"]["language_info"]["mimetype"] == "application/julia"
        @test notebook["metadata"]["language_info"]["name"] == "julia"
        @test string(notebook["metadata"]["language_info"]["version"]) == string(VERSION)
        @test !isempty(notebook["cells"])
   end

end

function main()
    test_pluto_parser()
    test_pluto2jupyter()
end
main()


text:
#=
--------------------------------------------------------------------------------------------------------------------------------
Types and constructors
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    VARSettings(...)

Define an immutable structure used to initialise the estimation routine for VAR(q) models.

# Arguments
- `Y`: observed measurements (`nxT`)
- `n`: Number of series
- `T`: Number of observations
- `q`: Order of the autoregressive polynomial
- `nq`: n*q
- `m`: n*q
- `λ`: overall shrinkage hyper-parameter for the elastic-net penalty
- `α`: weight associated to the LASSO component of the elastic-net penalty
- `β`: additional shrinkage for distant lags (p>1)
- `Γ`: Diagonal matrix used to input the hyperparameters in the estimation - see Pellegrino (2022) for details
- `ε`: Small number (default: 1e-4)
- `tol`: tolerance used to check convergence (default: 1e-4)
- `max_iter`: maximum number of iterations for the estimation algorithm (default: 1000)
- `prerun`: number of iterations prior the actual estimation algorithm (default: 2)
- `check_quantile`: check the quantile of the relative change of the parameters for convergence purposes (default: false)
- `verb`: Verbose output (default: true)
"""
struct VARSettings <: EstimSettings
    Y::Union{FloatMatrix, JMatrix{Float64}}
    n::Int64
    T::Int64
    q::Int64
    nq::Int64
    m::Int64
    λ::Float64
    α::Float64
    β::Float64
    Γ::DiagMatrix
    ε::Float64
    tol::Float64
    max_iter::Int64
    prerun::Int64
    check_quantile::Bool
    verb::Bool
end

# VARSettings constructor
function VARSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, q::Int64, λ::Float64, α::Float64, β::Float64; ε::Float64=1e-4, tol::Float64=1e-4, max_iter::Int64=1000, prerun::Int64=2, check_quantile::Bool=false, verb::Bool=true)

    # Compute missing inputs
    n, T = size(Y);
    m = n*q;
    Γ = build_Γ(n, q, λ, β);

    # Return VARSettings
    return VARSettings(Y, n, T, q, n*q, m, λ, α, β, Γ, ε, tol, max_iter, prerun, check_quantile, verb);
end

#=
--------------------------------------------------------------------------------------------------------------------------------
Initialisation
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    check_model_bounds(estim::VARSettings)

Check whether the arguments for the estimation of a VAR(q) model are correct.
"""
function check_model_bounds(estim::VARSettings)
    check_bounds(estim.q, 0);
    check_bounds(estim.λ, 0);
    check_bounds(estim.α, 0, 1);
    check_bounds(estim.β, 1);
    check_bounds(estim.max_iter, 3);
    check_bounds(estim.max_iter, estim.prerun);
    check_bounds(estim.n, 2); # It supports only multivariate models (for now ...)
end

"""
    initialise(estim::VARSettings)

Initialise the ECM algorithm for VAR(q) models.

# Arguments
- `estim`: settings used for the estimation

# References
Pellegrino (2022)
"""
function initialise(estim::VARSettings)

    # Interpolate data
    data = interpolate_series(estim.Y, estim.n, estim.T);

    # VAR(q) data
    Y, Y_lagged = lag(data, estim.q);

    # Estimate ridge VAR(q) coefficients
    C_star_init = Y*Y_lagged'/Symmetric(Y_lagged*Y_lagged' + estim.Γ);
    enforce_causality_and_invertibility!(C_star_init);

    # Estimate var-cov matrix of the residuals
    U = Y - C_star_init*Y_lagged;
    Q_init = Symmetric((U*U')/(estim.T - estim.q))::SymMatrix;

    # Initialise B_star
    B_star_init = zeros(estim.n, estim.m-estim.n);

    sspace = KalmanSettings(estim.Y, sspace_representation(estim, B_star_init, C_star_init, Q_init)..., compute_loglik=false);

    # Coordinates
    coordinates_transition_current = collect(1:estim.n);
    coordinates_transition_lagged = collect(1:estim.m);
    coordinates_transition_PPs = coordinates_transition_lagged .+ estim.n;
    coordinates_transition_P0 = LinearIndices(sspace.P0)[:];
    coordinates_free_params_C = CartesianIndices(sspace.C)[1:estim.n, 1:estim.m][:];

    # Useful views
    C_star = @view sspace.C[1:estim.n, 1:estim.m]; # in the case of the VAR this is equivalent to the cm view
    Q_view = @view sspace.DQD.data[1:sspace.Y.n, 1:sspace.Y.n];

    # Return output
    return sspace, nothing, C_star, Q_view, nothing, coordinates_transition_current, coordinates_transition_lagged, coordinates_transition_PPs, coordinates_transition_P0, nothing, coordinates_free_params_C;
end

#=
--------------------------------------------------------------------------------------------------------------------------------
CM-step
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    cm_step!(estim::VARSettings, sspace::KalmanSettings, B_star::Nothing, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, F_raw::FloatMatrix, G::FloatMatrix, H_raw::FloatMatrix, M::Nothing, N::Nothing, O::Nothing, coordinates_free_params_B::Nothing, coordinates_free_params_C::CoordinatesVector)

Run the CM-step for elastic-net VAR(q) models.
"""
function cm_step!(estim::VARSettings, sspace::KalmanSettings, B_star::Nothing, C_star::SubArray{Float64}, Q_view::SubArray{Float64}, F_raw::FloatMatrix, G::FloatMatrix, H_raw::FloatMatrix, M::Nothing, N::Nothing, O::Nothing, coordinates_free_params_B::Nothing, coordinates_free_params_C::CoordinatesVector)

    # Make sure `F_raw` and `H_raw` are symmetric
    F = Symmetric(F_raw)::SymMatrix;
    H = Symmetric(H_raw)::SymMatrix;

    # Companion forms
    C_companion_old = sspace.C[1:estim.m, 1:estim.m];
    C_companion_new = @view sspace.C[1:estim.m, 1:estim.m];

    # Compute inverse of Q
    inv_Q = inv(Symmetric(Q_view));

    # Loop over the free parameters
    for ij in coordinates_free_params_C

        # Coordinates
        i,j = ij.I;

        # Compute new value for sspace.C[ij]
        C_star[ij] = soft_thresholding(turbo_dot(inv_Q[:,i], G[:,j] - C_star*H[:,j]) + inv_Q[i,i]*C_star[ij]*H[j,j], 0.5*estim.α*estim.Γ[j,j]);
        C_star[ij] /= inv_Q[i,i]*H[j,j] + (1-estim.α)*estim.Γ[j,j];
    end

    # Adjust coefficients (if needed)
    enforce_causality_and_invertibility!(C_companion_new, C_companion_old);

    # CM-step for Q
    Q_cm_step = (F-G*C_star'-C_star*G'+C_star*H*C_star')/estim.T;
    for ij in eachindex(sspace.Q) # computing eachindex on sspace.Q is convenient since it gives Cartesian indices aligned with sspace.DQD
        Q_view[ij] = Q_cm_step[ij];
    end
end

#=
--------------------------------------------------------------------------------------------------------------------------------
Model validation
--------------------------------------------------------------------------------------------------------------------------------
=#

"""
    VARSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, model_args::Tuple, model_kwargs::Nothing, p::Int64, λ::Float64, α::Float64, β::Float64)
    VARSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, model_args::Tuple, model_kwargs::NamedTuple, p::Int64, λ::Float64, α::Float64, β::Float64)

Generate the relevant VARSettings structure by using the candidate hyperparameters `p`, `λ`, `α`, `β`.
"""
VARSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, model_args::Tuple, model_kwargs::Nothing, p::Int64, λ::Float64, α::Float64, β::Float64) = VARSettings(Y, p, λ, α, β, verb=false);
VARSettings(Y::Union{FloatMatrix, JMatrix{Float64}}, model_args::Tuple, model_kwargs::NamedTuple, p::Int64, λ::Float64, α::Float64, β::Float64) = VARSettings(Y, p, λ, α, β; model_kwargs...);


text:
function TGEAR(THTL)
      if (THTL<=0.77)
        return tgearv = 64.94*THTL
      else
        return tgearv = 217.38*THTL-117.38
      end
end


text:
module ScHoLP

using Arpack
using Base.Threads
using Combinatorics
using DataStructures
using DelimitedFiles
using Krylov
using LinearAlgebra
using LinearOperators
using Printf
using Random
using SparseArrays
using Statistics
using StatsBase

# This needs to go first
include("util.jl")

# Computations of various dataset statistics
include("simplicial_closure_probs.jl")
include("summary_statistics.jl")
include("temporal_asynchronicity.jl")

# Triangle prediction
include("walk_scores.jl")
include("local_scores.jl")
include("simplicial_ppr_scores.jl")

# Lifecycle analysis
include("lifecycles.jl")

end # module


text:
# Benchmark on column-wise distance evaluation

using Distances

macro bench_colwise_dist(repeat, dist, x, y)
    quote
        println("bench ", typeof($dist))

        # warming up
        r1 = evaluate($dist, ($x)[:,1], ($y)[:,1])
        colwise($dist, $x, $y)

        # timing

        t0 = @elapsed for k = 1 : $repeat
            n = size($x, 2)
            r = Array(typeof(r1), n)
            for j = 1 : n
                r[j] = evaluate($dist, ($x)[:,j], ($y)[:,j])
            end
        end
        @printf "    loop:     t = %9.6fs\n" (t0 / $repeat)

        t1 = @elapsed for k = 1 : $repeat
            r = colwise($dist, $x, $y)
        end
        @printf "    colwise:  t = %9.6fs  |  gain = %7.4fx\n" (t1 / $repeat) (t0 / t1)
        println()
    end
end

m = 200
n = 10000

x = rand(m, n)
y = rand(m, n)
w = rand(m)

Q = rand(m, m)
Q = Q' * Q

@bench_colwise_dist 20 SqEuclidean() x y
@bench_colwise_dist 20 Euclidean() x y
@bench_colwise_dist 20 Cityblock() x y
@bench_colwise_dist 20 Chebyshev() x y
@bench_colwise_dist  5 Minkowski(3.0) x y
@bench_colwise_dist 20 Hamming() x y

@bench_colwise_dist 20 CosineDist() x y
@bench_colwise_dist 10 CorrDist() x y
@bench_colwise_dist 20 ChiSqDist() x y
@bench_colwise_dist 10 KLDivergence() x y
@bench_colwise_dist  5 JSDivergence() x y

@bench_colwise_dist 10 BhattacharyyaDist() x y
@bench_colwise_dist 10 HellingerDist() x y

@bench_colwise_dist 20 WeightedSqEuclidean(w) x y
@bench_colwise_dist 20 WeightedEuclidean(w) x y
@bench_colwise_dist 20 WeightedCityblock(w) x y
@bench_colwise_dist  5 WeightedMinkowski(w, 3.0) x y
@bench_colwise_dist 20 WeightedHamming(w) x y

@bench_colwise_dist 10 SqMahalanobis(Q) x y
@bench_colwise_dist 10 Mahalanobis(Q) x y


text:
function medianSD(x)

    x = vec(x)

    if mod.(length(x), 2) == 0

        #do find mid location
        midLocation = (length(x) + 1) / 2
        # evaluate mean of the upper and lower index
        meanSD([sort(x)[Int64(floor(midLocation))], sort(x)[Int64(ceil(midLocation))]])

    elseif mod.(length(x), 2) != 0

        midLocation = (length(x) + 1) / 2
        sort(x)[Int64(midLocation)]

    else

        println("No median")

    end

end


text:
module JuliaBLAS

include("lowlevel.jl")
include("gemm.jl")

export addmul!

end


text:
"""
RenoiseOSC.jl is a collection of wrappers around the [Renoise OSC API functions](https://tutorials.renoise.com/wiki/Open_Sound_Control)
"""
module RenoiseOSC

export luaeval,
    tempo,
    editmode,
    octave,
    patternfollow,
    editstep,
    setmacroparam,
    monophonic,
    monophonicglide,
    phraseplayback,
    phraseprogram,
    quantizationmode,
    scalekey,
    scalemode,
    transpose,
    volume,
    volumedb,
    linesperbeat,
    metronome,
    metronomeprecount,
    quantization,
    quantizationstep,
    scheduleadd,
    scheduleset,
    slotmute,
    slotunmute,
    trigger,
    ticksperline,
    bypass,
    setparam,
    mute,
    unmute,
    solo,
    outputdelay,
    postfxpanning,
    postfxvolume,
    postfxvolumedb,
    prefxpanning,
    prefxvolume,
    prefxvolumedb,
    prefxwidth,
    loopblock,
    loopblockmovebackwards,
    loopblockmoveforwards,
    looppattern,
    start,
    stop,
    resume,
    midievent,
    noteon,
    noteoff,
    sethost!,
    setport!

using OpenSoundControl
using Sockets

const settings = Dict{String,Union{Sockets.InetAddr,UDPSocket}}("server" => Sockets.InetAddr(ip"127.0.0.1", 8000))

"""
    sethost!(host::Union{AbstractString,Sockets.IPAddr})

Set the host of the Renoise OSC server. Default value is `ip"127.0.0.1"`.

!!! compat "Julia 1.3"
    Setting the host by passing an AbstractString requires at least Julia 1.3.

# Example

```jldoctest; setup=:(using Sockets: @ip_str; using RenoiseOSC: sethost!)
julia> sethost!(ip"127.0.0.1")
Sockets.InetAddr{Sockets.IPv4}(ip"127.0.0.1", 8000)
```
"""
sethost!(host::Union{AbstractString,IPAddr}) = settings["server"] = Sockets.InetAddr(host, settings["server"].port)

"""
    setport!(port::Integer)

Set the port of the Renoise OSC server. Default value is `8000`.
"""
setport!(port::Integer) = settings["server"] = Sockets.InetAddr(settings["server"].host, port)

"""
    setaddress!(host::Union{AbstractString,Sockets.IPAddr}, port::Integer)

Set the host and the port of the OSC Server.
"""
setaddress!(host::Union{AbstractString,IPAddr}, port::Integer) = settings["server"] = Sockets.InetAddr(host, port)

function postmessage(path::AbstractString, argtypes::AbstractString="", args...)
    send(
        get!(UDPSocket, settings, "socket"),
        settings["server"].host,
        settings["server"].port,
        OpenSoundControl.message("/renoise/" * path, argtypes, args...).data
    )
end

"""
    luaeval(expr::AbstractString)

Evaluate a Lua expression inside Renoise's scripting environment.

# Example

```julia
luaeval("renoise.song().transport.bpm = 132")
```
"""
luaeval(expr::AbstractString) = postmessage("evaluate", "s", expr)

"""
    tempo(bpm::Integer)

Set the song's current bpm `[32:999]`
"""
function tempo(bpm::Integer)
    if 32 <= bpm <= 999
        postmessage("song/bpm", "i", Int32(bpm))
    else
        @warn "Can't set bpm to $(bpm), bpm must be in 32:999"
    end
end

"""
    editmode(on::Bool)

Set the song's global edit mode on or off
"""
editmode(on::Bool) = postmessage("song/edit/mode", on ? "T" : "F")

"""
    octave(oct::Integer)

Sets the song's current octave [0:8]
"""
function octave(oct::Integer)
    if oct in 0:8
        postmessage("song/edit/octave", "i", Int32(oct))
    else
        @warn "Can't set octave to $(oct), octave must be in 0:8"
    end
end

"""
    patternfollow(on::Bool)

Enable or disable global pattern follow mode
"""
patternfollow(on::Bool) = postmessage("song/edit/pattern_follow", on ? "T" : "F")

"""
    editstep(step::Integer)

Set the song's current edit step `[0:8]`
"""
function editstep(step::Integer)
    if step in 0:8
        postmessage("song/edit/step", "i", Int32(step))
    else
        @warn "Can't set edit step to $(step), edit step must be in 0:8"
    end
end

"""
    setmacroparam(param::Integer, value::Real; instrument::Integer = -1)

Set `instrument`'s macro parameter value `[0.0:1.0]`.
Default to the currently selected instrument.
"""
function setmacroparam(param::Integer, value::Real; instrument::Integer=-1)
    postmessage("song/instrument/$(instrument)/$(param)", "d", Float64(value))
end

"""
    monophonic(mono::Bool; instrument::Integer=-1)

Enable or disable `instrument`'s mono mode.
Default to the currently selected instrument.
"""
monophonic(mono::Bool; instrument::Integer=-1) = postmessage("song/instrument/$(instrument)/monophonic", mono ? "T" : "F")

"""
    monophonicglide(glide::Integer; instrument::Integer=-1)

Set `instrument`'s glide amount `[0:255]`.
Default to the currently selected instrument.
"""
function monophonicglide(glide::Integer; instrument::Integer=-1)
    if glide in 0:255
        postmessage("song/instrument/$(instrument)/monophonic_glide", "i", Int32(glide))
    else
        @warn "Can't set monophonic glide to $(glide), glide must be in 0:255"
    end
end

"""
    phraseplayback(mode::AbstractString; instrument::Integer=-1)

Set `instrument`'s phraseplayback mode `["Off", "Program", "Keymap"]`
Default to the currently selected instrument.
"""
function phraseplayback(mode::AbstractString; instrument::Integer=-1)
    if mode in Set(["Off", "Program", "Keymap"])
        postmessage("song/instrument/$(instrument)/phrase_playback", "s", mode)
    else
        @warn "Can't set phrase playback mode to \"$(mode)\", mode must be one of `[\"Off\", \"Program\", \"Keymap\"]`"
    end
end

"""
    phraseprogram(program::Integer; instrument::Integer=-1)

Set `instrument`'s phrase program number [0:127]
Default to the currently selected instrument.
"""
function phraseprogram(program::Integer; instrument::Integer=-1)
    if program in 0:127
        postmessage("song/instrument/$(instrument)/phrase_playback", "i", Int32(program))
    else
        @warn "Can't set phrase program to $(program), program must be in 0:127"
    end
end

"""
    quantizationmode(mode::AbstractString; instrument::Integer=-1)

Set `instrument`'s quantization method, one of `"None"`, `"Line"`, `"Beat"`, `"Bar"`.
Default to the currently selected instrument.
"""
function quantizationmode(mode::AbstractString; instrument::Integer=-1)
    if mode in Set(["None", "Line", "Beat", "Bar"])
        postmessage("song/instrument/$(instrument)/quantize", "s", mode)
    else
        @warn "Can't set quantization mode to \"$(mode)\", mode must be one of `[\"None\", \"Line\", \"Beat\", \"Bar\"]`"
    end
end

"""
    scalekey(key::AbstractString; instrument::Integer=-1)

Set `instrument`'s note scaling key.
Default to the currently selected instrument.
"""
function scalekey(key::AbstractString; instrument::Integer=-1)
    if key in Set(["C", "C#", "D", "D#", "E", "F", "F#", "G", "G#", "A", "A#", "B"])
        postmessage("song/instrument/$(instrument)/scale_key", "s", key)
    else
        @warn "Can't set scale key to \"$(key)\", key must be one of [\"C\", \"C#\", \"D\", \"D#\", \"E\", \"F\", \"F#\", \"G\", \"G#\", \"A\", \"A#\", \"B\"]"
    end
end

"""
    scalemode(mode::AbstractString; instrument::Integer=-1)

Set `instrument`'s note scaling mode.
Default to the currently selected instrument.
"""
function scalemode(mode::AbstractString; instrument::Integer=-1)
    postmessage("song/instrument/$(instrument)/scale_mode", "s", mode)
end

"""
    transpose(pitch::Integer; instrument::Integer=-1)

Set `instrument`'s global pitch transpose `[-120:120]``.
Default to the currently selected instrument.
"""
function transpose(pitch::Integer; instrument::Integer=-1)
    if pitch in -120:120
        postmessage("song/instrument/$(instrument)/transpose", "i", pitch)
    else
        @warn "Can't transpose by $(pitch), `pitch` must be in `-120:120`"
    end
end

"""
    volume(level::Real; instrument::Integer=-1)

Set `instrument`'s global volume to `level` `[0.0:db2lin(6)]≈[0.0:2.0]`.
Default to the currently selected instrument.
"""
function volume(level::Real; instrument::Integer=-1)
    postmessage("song/instrument/$(instrument)/volume", "d", clamp(Float64(level), 0.0, 1.9952623149688795))
end

"""
    volumedb(level::Real; instrument::Integer=-1)

Set `instrument`'s global volume to `level` in decibels `[0.0:6.0]`.
Default to the currently selected instrument.
"""
function volumedb(level::Real; instrument::Integer=-1)
    postmessage("song/instrument/$(instrument)/volume_db", "d", clamp(Float64(level), 0.0, 6.0))
end

"""
    linesperbeat(lpb::Integer)

Set the song's current lines per beat [1:255]
"""
function linesperbeat(lpb::Integer)
    if lpb in 1:255
        postmessage("song/lpb", "i", Int32(lpb))
    else
        @warn "Can't set lines per beat to $(lpb), `lpb` must be in `1:255`"
    end
end

"""
    metronome(on::Bool)

Enable or disable the global metronome
"""
metronome(on::Bool) = postmessage("song/record/metronome", on ? "T" : "F")

"""
    metronomeprecount(on::Bool)

Enable or disable the global metronome precount
"""
metronomeprecount(on::Bool) = postmessage("song/record/metronome_precount", on ? "T" : "F")

"""
    quantization(on::Bool)

Enable or disable the global record quantization
"""
quantization(on::Bool) = postmessage("song/record/quantization", on ? "T" : "F")

"""
    quantizationstep(step::Integer)

Set the global record quantization step [1:32]
"""
function quantizationstep(step::Integer)
    if step in 1:32
        postmessage("song/record/quantization_step", b ? "T" : "F")
    else
        @warn "Can't set quantization step to $(step), `step` must be in `1:32`"
    end
end

"""
    scheduleadd(sequence::Integer)

Add a scheduled sequence playback pos
"""
scheduleadd(sequence::Integer) = postmessage("song/sequence/schedule_add", "i", Int32(sequence))

"""
    scheduleset(sequence::Integer)

Replace the current sequence playback pos
"""
scheduleset(sequence::Integer) = postmessage("song/sequence/schedule_set", "i", Int32(sequence))

"""
    slotmute(track::Integer, sequence::Integer)
    slotmute(mute::Bool, track::Integer, sequence::Integer)

Mute the given track, sequence slot in the matrix.
"""
slotmute(track::Integer, sequence::Integer) = postmessage("song/sequence/slot_mute", "ii", Int32(track), Int32(sequence))
slotmute(mute::Bool, track::Integer, sequence::Integer) = postmessage("song/sequence/slot_$(mute ? "" : "un")mute", "ii", Int32(track), Int32(sequence))

"""
    slotunmute(track::Integer, sequence::Integer)

Mute the given track, sequence slot in the matrix.
"""
slotunmute(track::Integer, sequence::Integer) = postmessage("song/sequence/slot_unmute", "ii", Int32(track), Int32(sequence))

"""
    trigger(sequence::Integer)

Set playback pos to the specified sequence pos
"""
trigger(sequence::Integer) = postmessage("song/sequence/trigger", "i", Int32(sequence))

"""
    ticksperline(tpl::Integer)

Set the song's current ticks per line [1:16]
"""
function ticksperline(tpl::Integer)
    if tpl in 1:16
        postmessage("song/tpl", "i", Int32(tpl))
    else
        @warn "Can't set ticks per line to $(tpl), ticks per line must be in `1:16`"
    end
end

"""
    bypass(bypassed::Bool; track::Integer=-1, device::Integer=-1)

Set the bypass status of a device, set `true` to bypass the device.
When unspecified, `track` and `device` default to the currently selected track and device.
"""
function bypass(bypassed::Bool; track::Integer=-1, device::Integer=-1)
    postmessage("song/track/$(track)/device/$(device)/bypass", bypassed ? "T" : "F")
end

"""
    setparam(key::Union{AbstractString,Integer}, value::Real; track::Integer=-1, device::Integer=-1)

Set the parameter of any device by its name or index. `value` is clamped between `0.0` and `1.0`.
When unspecified, track and device default to the currently selected track and device.
"""
function setparam(key::Integer, value::Real; track::Integer=-1, device::Integer=-1)
    postmessage("song/track/$(track)/device/$(device)/set_parameter_by_index", "id", Int32(key), Float64(clamp(value, 0.0, 1.0)))
end
function setparam(key::AbstractString, value::Real; track::Integer=-1, device::Integer=-1)
    postmessage("song/track/$(track)/device/$(device)/set_parameter_by_name", "sd", key, Float64(clamp(value, 0.0, 1.0)))
end

"""
    mute(muted::Bool = true; track::Integer = -1)

Mute or unmute the given track. Default to the currently selected track.
"""
mute(muted::Bool=true; track::Integer=-1) = postmessage("song/track/$(track)/$(muted ? "" : "un")mute")

"""
    unmute(; track::Integer=-1)
    mute(false; track::Integer=-1)

Unmute the given track. Default to the currently selected track.
"""
unmute(; track::Integer=-1) = postmessage("song/track/$(track)/unmute")

"""
    solo(; track::Integer)

Toggle solo for the given track. Default to the currently selected track.
"""
solo(; track::Integer=-1) = postmessage("song/track/$(track)/solo")

"""
    outputdelay(δt::Real; track::Integer)

Set the given track's output delay in milliseconds [-100:100]. Default to the currently selected track.
"""
function outputdelay(δt::Real; track::Integer=-1)
    if -100.0 <= δt <= 100.0
        postmessage("song/track/$(track)/output_delay", "d", Float64(δt))
    else
        @warn "Can't set output delay to `$(δt)`. The output delay time must be between `-100.0` and `100.0`"
    end
end

"""
    postfxpanning(pan::Integer; track::Integer=-1)

Set `track`'s post-FX panning, `[-50:50]` left to right.
Default to the currently selected track.
"""
function postfxpanning(pan::Integer; track::Integer=-1)
    if -50 <= pan <= 50
        postmessage("song/track/$(track)/postfx_panning", "i", Int32(pan))
    else
        @warn "Can't set post-FX panning to $(pan), `pan` must be in `-50:50`"
    end
end

"""
    postfxvolume(level::Real; track::Integer=-1)

Set `track`s post-FX volume, `[0.0:db2lin(6.0)]`, -Inf:+6dB.
Default to the currently selected track.
"""
function postfxvolume(level::Real; track::Integer=-1)
    postmessage("song/track/$(track)/postfx_volume", "d", clamp(Float64(level), 0.0, 1.9952623149688795))
end

"""
    postfxvolumedb(level::Real; track::Integer=-1)

Set `track`s post-FX volume in decibels, `[-200.0:3.0]`.
Default to the currently selected track.
"""
function postfxvolumedb(level::Real; track::Integer=-1)
    postmessage("song/track/$(track)/postfx_volume_db", "d", clamp(Float64(level), -200.0, 3.0))
end

"""
    prefxpanning(pan::Integer; track::Integer=-1)

Set `track`'s pre-FX panning, `[-50:50]` left to right.
Default to the currently selected track.
"""
function prefxpanning(pan::Integer; track::Integer=-1)
    if -50 <= pan <= 50
        postmessage("song/track/$(track)/prefx_panning", "i", Int32(pan))
    else
        @warn "Can't set pre-FX panning to $(pan), `pan` must be in `-50:50`"
    end
end

"""
    prefxvolume(level::Real; track::Integer=-1)

Set `track`s pre-FX volume, `[0.0:db2lin(6.0)]`, -Inf:+6dB.
Default to the currently selected track.
"""
function prefxvolume(level::Real; track::Integer=-1)
    postmessage("song/track/$(track)/prefx_volume", "d", clamp(Float64(level), 0.0, 1.9952623149688795))
end

"""
    prefxvolumedb(level::Real; track::Integer=-1)

Set `track`s pre-FX volume in decibels, `[-200.0:3.0]`.
Default to the currently selected track.
"""
prefxvolumedb(level::Real; track::Integer=-1) = postmessage("song/track/$(track)/prefx_volume_db", "d", clamp(Float64(level), -200.0, 3.0))

"""
    prefxwidth(width::Real; track::Integer=-1)

Set `track`'s pre-FX width `[0.0:1.0]`.
Default to the currently selected track.
"""
prefxwidth(width::Real; track::Integer=-1) = postmessage("song/track/$(track)/prefx_width", "d", clamp(Float64(width), 0.0, 1.0))

"""
    loopblock(loop::Bool)

Enable or disable pattern block loop.
"""
loopblock(loop::Bool) = postmessage("/transport/loop/block", loop ? "T" : "F")

"""
    loopblockmovebackwards()

Move the block loop one segment backwards
"""
loopblockmovebackwards() = postmessage("transport/loop/block_move_backwards")

"""
    loopblockmoveforwards()

Move the block loop one segment forwards
"""
loopblockmoveforwards() = postmessage("transport/loop/block_move_forwards")

"""
    looppattern(loop::Bool)

Enable or disable looping the current pattern.
"""
looppattern(loop::Bool) = postmessage("transport/loop/pattern", loop ? "T" : "F")

"""
    start()

Start playback or restart playing the current pattern
"""
start() = postmessage("transport/start")

"""
    stop()

Stop playback
"""
stop() = postmessage("transport/stop")

"""
    panic()

Stop playback and reset all playing instruments and DSPs
"""
panic() = postmessage("transport/panic")

"""
    resume()

Continue playback
"""
resume() = postmessage("transport/continue")

"""
    midievent(portid::UInt8, status::UInt8, data1::UInt8, data2::UInt8)
    midievent(event::Union{UInt32,UInt64,Array{<:Integer},NTuple{4,<:Integer}})

Fire a raw MIDI event.
"""
midievent(portid::UInt8, status::UInt8, data1::UInt8, data2::UInt8) = postmessage("trigger/midi", "m", [portid, status, data1, data2])
function midievent(event::Union{Array{<:Integer},NTuple{4,<:Integer}})
    if all(<=(255), event)
        postmessage("trigger/midi", "m", UInt8.(event))
    else
        @warn "Can't send MIDI event"
    end
end
midievent(portid::Integer, status::Integer, data1::Integer, data2::Integer) = midievent([portid, status, data1, data2])
midievent(event::UInt64) = postmessage("trigger/midi", "h", event)
midievent(event::UInt32) = postmessage("trigger/midi", "i", event)

"""
    noteon(pitch::Integer, velocity::Integer; instrument::Integer=-1, track::Integer=-1)

Turn on the note with the velocity.
Default to the currently selected instrument and track.
"""
function noteon(pitch::Integer, velocity::Integer; instrument::Integer=-1, track::Integer=-1)
    postmessage("trigger/note_on", "iiii", Int32.((instrument, track, pitch, velocity))...)
end

"""
    noteoff(pitch::Integer, velocity::Integer; instrument::Integer=-1, track::Integer=-1)

Turn off the note
"""
function noteoff(pitch::Integer; instrument::Integer=-1, track::Integer=-1)
    postmessage("trigger/note_off", "iii", Int32.((instrument, track, pitch))...)
end

end # module

max_stars_repo_name: stellartux/RenoiseOSC

max_issues_repo_name: stellartux/RenoiseOSC

max_forks_repo_name: stellartux/RenoiseOSC


text:
"""

```
writePDB(atoms::Vector{Atom}, filename, selection)
```

Function that writes a PDB file given the vector of atoms, with
optional definition of a selection to be print.

"""
function writePDB(atoms::AbstractVector{Atom}, filename, selection)
    query = parse_query(selection)
    writePDB(atoms, filename, only = atom -> apply_query(query, atom))
end

function writePDB(atoms::AbstractVector{Atom}, filename; only = all)
    file = open(filename, "w")
    for atom in atoms
        if only(atom)
            println(file, write_atom(atom))
        end
    end
    close(file)
end

max_stars_repo_name: mcubeg/PDBTools

max_issues_repo_name: mcubeg/PDBTools

max_forks_repo_name: mcubeg/PDBTools


text:
using Pkg
Pkg.activate(".")
Pkg.instantiate()

using RadMod2D
using GLMakie

include("./models2D.jl")
include("./plot2D.jl")

#### calculation
elemsize = 0.01
n = 50

m = model_line_to_line_with_obstacles(1.0, 1.0, 0.5, 6, elemsize)

# tiles
dx, dy = create_tiles(m, n)
t_occ = check_tile_occupation(m, dx, dy, n, blockparts = 3:m.npar)
tile_occ_analysis(t_occ, printit = true)

println("elem size / tile size: ", round(elemsize/((dx+dy)/2), digits = 2), " %")

# classic algo
vfmat = zeros(Float64, m.nelem, m.nelem)
existing_vf!(m, vfmat)
@time blocking_vf_with_tiles!(m, vfmat, dx, dy, n, t_occ)
calculating_vf!(m, vfmat, normit = false)
vfmatp = vfmat_to_parts(m, vfmat, normit = false)

# simplified algo
vfmat_simp = zeros(Float64, m.nelem, m.nelem)
existing_vf!(m, vfmat_simp)
@time blocking_vf_with_tiles_simplified!(m, vfmat_simp, dx, dy, n, t_occ, elem_in_t = 1, skipped_t = 1)
calculating_vf!(m, vfmat_simp, normit = false)
vfmatp_simp = vfmat_to_parts(m, vfmat_simp, normit = false)

#### 2D plot
fig = Figure(resolution = (900, 900))
ax = fig[1, 1] = Axis(fig)
ax.xlabel = "X in m"
ax.ylabel = "Y in m"
ax.aspect = DataAspect()
plot_empty_tiles(fig, ax, dx, dy, n)
plot_occupied_tiles(fig, ax, dx, dy, n, t_occ)
plot_model(fig, ax, m, shownvec = true, shownodes = false, showcom = false, showleg = true)
display(fig)

#### compare results of classic algo with simplified algo

vf12_classic = vfmatp[1,2]
vf12_simp = vfmatp_simp[1,2]

diff = vf12_classic - vf12_simp

println("classic algo: ", vf12_classic)
println("simplified algo: ", vf12_simp)
println("diff: ", diff)

println("error rate: ", round((1-vf12_simp/vf12_classic)*100, digits = 2), " %")

# include("./test/run_algo_vs_simpalgo2.jl")


text:
using Literate: Literate
using Pkg: Pkg

# Print `@debug` statements (https://github.com/JuliaDocs/Documenter.jl/issues/955)
if haskey(ENV, "GITHUB_ACTIONS")
    ENV["JULIA_DEBUG"] = "Documenter"
end

const EXPERIMENTS = ("ols", "synthetic", "friedman")
const INPUT = joinpath(@__DIR__, "..", "experiments")
const OUTPUT = joinpath(@__DIR__, "src", "generated")

ispath(OUTPUT) && rm(OUTPUT; recursive=true)
mkpath(OUTPUT)

# Link existing data to avoid expensive computations
if isdir(joinpath(INPUT, "data"))
    symlink(joinpath(INPUT, "data"), joinpath(OUTPUT, "data"))
end

# Add links and explanations below the first heading of level 1
function preprocess(content)
    sub = SubstitutionString(
        """
\\0
#
#md # [![](https://img.shields.io/badge/show-nbviewer-579ACA.svg)](@__NBVIEWER_ROOT_URL__/generated/@__NAME__.ipynb)
#
# You are seeing the
#md # HTML output generated by [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl) and
#nb # notebook output generated by
# [Literate.jl](https://github.com/fredrikekre/Literate.jl) from the
# [Julia source file](@__REPO_ROOT_URL__/experiments/src/@__NAME__/script.jl).
# The corresponding
#md # notebook can be viewed in [nbviewer](@__NBVIEWER_ROOT_URL__/generated/@__NAME__.ipynb),
#nb # HTML output can be viewed [here](https://devmotion.github.io/Calibration_ICLR2021/dev/generated/@__NAME__/),
# and the plain script output can be found [here](./@__NAME__.jl).
#
#note # !!!
#note # If you want to run the experiments, make sure you have an identical environment.
#note # Please use Julia $(VERSION) and activate and instantiate the environment using
#note # [this Project.toml file](@__REPO_ROOT_URL__/experiments/src/@__NAME__/Project.toml)
#note # and [this Manifest.toml file](@__REPO_ROOT_URL__/experiments/src/@__NAME__/Manifest.toml).
#note #
#note # [The Github repository](https://github.com/devmotion/Calibration_ICLR2021/) contains
#note # [more detailed instructions](@__REPO_ROOT_URL__/experiments/README.md) and a
#note # `nix` project environment with a pinned Julia binary for improved reproducibility.
#md #
#md # ```@setup @__NAME__
#md # using Pkg: Pkg
#md # Pkg.activate("$(INPUT)/src/@__NAME__")
#md # Pkg.instantiate()
#md # ```
#
        """,
    )
    return replace(content, r"^# # .*$"m => sub; count=1)
end

# Add admonitions
function md_note(content)
    return replace(
        replace(content, r"^#note # !!!$"m => "#md # !!! note"),
        r"^#note #(.*)$"m => s"#md #    \1",
    )
end
function nb_note(content)
    return replace(
        replace(content, r"^#note # !!!$"m => "#nb # > **Note**"),
        r"^#note #(.*)$"m => s"#nb # >\1",
    )
end

for name in EXPERIMENTS
    file = joinpath(INPUT, "src", name, "script.jl")

    # Activate project environment
    Pkg.activate(dirname(file)) do
        Pkg.instantiate()

        # Export and execute output formats
        Literate.markdown(
            file,
            OUTPUT;
            name=name,
            documenter=true,
            preprocess=md_note ∘ preprocess,
            repo_root_url="https://github.com/devmotion/Calibration_ICLR2021/blob/main",
        )
        Literate.notebook(
            file,
            OUTPUT;
            name=name,
            documenter=true,
            preprocess=nb_note ∘ preprocess,
            repo_root_url="https://github.com/devmotion/Calibration_ICLR2021/blob/main",
        )
        Literate.script(file, OUTPUT; name=name, documenter=true)
    end
end

using Documenter

makedocs(;
    authors="David Widmann <david.widmann@it.uu.se>",
    sitename="Calibration tests beyond classification",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://devmotion.github.io/Calibration_ICLR2021",
        assets=String[],
    ),
    pages=[
        "index.md",
        "software.md",
        "Experiments" => [joinpath("generated", "$(name).md") for name in EXPERIMENTS],
    ],
)

# Remove link to existing (would be broken when deployed)
if islink(joinpath(@__DIR__, "build", "generated", "data"))
    rm(joinpath(@__DIR__, "build", "generated", "data"))
end

deploydocs(;
    repo="github.com/devmotion/Calibration_ICLR2021.git",
    push_preview=true,
    devbranch="main",
)


text:
using Test

using DescentMethods
using DescentMethods: argumentvec, gradientvec, step_origin, mcholesky!

using  LinearAlgebra

using Random

const test_rng = MersenneTwister(123581321)

OPT_TYPES = (
    SteepestDescent,
    HyperGradDescent,
    FixedRateDescent,
    MomentumDescent,
    NesterovMomentum,
    CGDescent,
    BFGS,
    CholBFGS,
)

@testset "Testing package DescentMethods" begin
# Test that all methods have a constructor
include("constructors.jl")

# Test that all methods have interface functions
include("interface.jl")

# Test linesearch on quadratic functions
include("linesearch.jl")

# Test modified Cholesky decomposition
include("cholesky.jl")

# Test convergence on the Rosenbrock function
include("convergence.jl")
end


text:
import Base: min, max, sum, count

export max, imax, min, imin, sum, count

immutable Reduction <: AFImpl
    max::Ptr{Void}
    maxAll::Ptr{Void}
    imax::Ptr{Void}
    imaxAll::Ptr{Void}
    min::Ptr{Void}
    minAll::Ptr{Void}
    imin::Ptr{Void}
    iminAll::Ptr{Void}
    sum::Ptr{Void}
    sumNan::Ptr{Void}
    sumAll::Ptr{Void}
    sumNanAll::Ptr{Void}
    count::Ptr{Void}
    countAll::Ptr{Void}

    function Reduction(ptr)
        new(
            Libdl.dlsym(ptr, :af_max),
            Libdl.dlsym(ptr, :af_max_all),
            Libdl.dlsym(ptr, :af_imax),
            Libdl.dlsym(ptr, :af_imax_all),
            Libdl.dlsym(ptr, :af_min),
            Libdl.dlsym(ptr, :af_min_all),
            Libdl.dlsym(ptr, :af_imin),
            Libdl.dlsym(ptr, :af_imin_all),
            Libdl.dlsym(ptr, :af_sum),
            Libdl.dlsym(ptr, :af_sum_nan),
            Libdl.dlsym(ptr, :af_sum_all),
            Libdl.dlsym(ptr, :af_sum_nan_all),
            Libdl.dlsym(ptr, :af_count),
            Libdl.dlsym(ptr, :af_count_all)
        )
    end
end

macro minMaxRed(regular, all, indexed, indexedAll)
    quote
        function $(esc(regular))(arr::AFArray, dim = -1)
            verifyAccess(arr)
            af = arr.af
            result = af.results.ptr
            err = ccall(af.reduction.$regular,
                Cint, (Ptr{Ptr{Void}}, Ptr{Void}, Int32),
                result, arr.ptr, Int32(dim < 0 ? firstDim(af, arr.ptr) - 1 : dim))
            assertErr(err)
            array(af, result[])
        end

        function $(esc(regular)){T<:Real}(::Type{T}, arr::AFArray)
            verifyAccess(arr)
            af = arr.af
            real = af.results.double
            imag = af.results.double2
            err = ccall(af.reduction.$all,
                Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Void}),
                real, imag, arr.ptr)
            assertErr(err)
            T(real[])
        end

        function $(esc(indexed))(arr::AFArray, dim = -1)
            verifyAccess(arr)
            af = arr.af
            result = af.results.ptr
            idx = af.results.ptr2
            err = ccall(af.reduction.$indexed,
                Cint, (Ptr{Ptr{Void}}, Ptr{Ptr{Void}}, Ptr{Void}, Int32),
                result, idx, arr.ptr, Int32(dim < 0 ? firstDim(af, arr.ptr) - 1 : dim))
            assertErr(err)
            array(af, result[]), array(af, idx[])
        end

        function $(esc(indexed)){T<:Real}(::Type{T}, arr::AFArray)
            verifyAccess(arr)
            af = arr.af
            real = af.results.double
            imag = af.results.double2
            idx = af.results.uint
            err = ccall(af.reduction.$indexedAll,
                Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{UInt32}, Ptr{Void}),
                real, imag, idx, arr.ptr)
            assertErr(err)
            T(real[]), idx[]
        end
    end
end

@minMaxRed(max, maxAll, imax, imaxAll)

@minMaxRed(min, minAll, imin, imaxAll)

function sum(arr::AFArray, dim::Int = -1)
    verifyAccess(arr)
    af = arr.af
    result = af.results.ptr
    err = ccall(af.reduction.sum,
        Cint, (Ptr{Ptr{Void}}, Ptr{Void}, Int32),
        result, arr.ptr, Int32(dim < 0 ? firstDim(af, arr.ptr) - 1 : dim))
    assertErr(err)
    array(af, result[])
end

function sum{T<:Union{Float32,Float64}}(arr::AFArray, dim::Int, nanVal::T)
    verifyAccess(arr)
    af = arr.af
    result = af.results.ptr
    err = ccall(af.reduction.sumNan,
        Cint, (Ptr{Ptr{Void}}, Ptr{Void}, Int32, Float64),
        result, arr.ptr, Int32(dim < 0 ? firstDim(af, arr.ptr) - 1 : dim), nanVal)
    assertErr(err)
    array(af, result[])
end

function sum{T<:Real}(::Type{T}, arr::AFArray)
    verifyAccess(arr)
    af = arr.af
    real = af.results.double
    imag = af.results.double2
    err = ccall(af.reduction.sumAll,
        Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Void}),
        real, imag, arr.ptr)
    assertErr(err)
    T(real[])
end

function sum{T<:Union{Float32,Float64}}(arr::AFArray, nanVal::T)
    verifyAccess(arr)
    af = arr.af
    real = af.results.double
    imag = af.results.double2
    err = ccall(af.reduction.sumNanAll,
        Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Void}, Float64),
        real, imag, arr.ptr, nanVal)
    assertErr(err)
    real[]
end

function count(arr::AFArray, dim::Int = -1)
    verifyAccess(arr)
    af = arr.af
    result = af.results.ptr
    err = ccall(af.reduction.count,
        Cint, (Ptr{Ptr{Void}}, Ptr{Void}, Int32),
        result, arr.ptr, Int32(dim < 0 ? firstDim(af, arr.ptr) - 1 : dim))
    assertErr(err)
    array(af, result[])
end

function count{T<:Real}(::Type{T}, arr::AFArray)
    verifyAccess(arr)
    af = arr.af
    real = af.results.double
    imag = af.results.double2
    err = ccall(af.reduction.countAll,
        Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{Void}),
        real, imag, arr.ptr)
    assertErr(err)
    T(real[])
end

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdradiogroup

"""
    'feffery'_antdradiogroup(;kwargs...)

An AntdRadioGroup component.

Keyword arguments:
- `id` (String; optional)
- `buttonStyle` (String; optional)
- `className` (String; optional)
- `defaultValue` (String; optional)
- `disabled` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `optionType` (String; optional)
- `options` (optional): . options has the following type: Array of lists containing elements 'label', 'value', 'disabled'.
Those elements have the following types:
  - `label` (String; required)
  - `value` (String; required)
  - `disabled` (Bool; optional)s
- `size` (String; optional)
- `style` (Dict; optional)
- `value` (String; optional)
"""
function 'feffery'_antdradiogroup(; kwargs...)
        available_props = Symbol[:id, :buttonStyle, :className, :defaultValue, :disabled, :loading_state, :optionType, :options, :size, :style, :value]
        wild_props = Symbol[]
        return Component("'feffery'_antdradiogroup", "AntdRadioGroup", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: Beddow/feffery-antd-components

max_issues_repo_name: Beddow/feffery-antd-components

max_forks_repo_name: Beddow/feffery-antd-components


text:
export jacobian

# 1-polytopes
jacobian(l::LineSegment, r) = l[2] - l[1]
jacobian(q::QuadraticSegment, r) = (4r - 3)*(q[1] - q[3]) +
                                   (4r - 1)*(q[2] - q[3])

# 2-polytopes
function jacobian(tri::Triangle, r, s)
    ∂ᵣF = tri[2] - tri[1]
    ∂ₛF = tri[3] - tri[1]
    return hcat(∂ᵣF, ∂ₛF)
end

function jacobian(quad::Quadrilateral, r, s)
    ∂ᵣF = (1 - s)*(quad[2] - quad[1]) + s*(quad[3] - quad[4])
    ∂ₛF = (1 - r)*(quad[4] - quad[1]) + r*(quad[3] - quad[2])
    return hcat(∂ᵣF, ∂ₛF)
end

function jacobian(tri6::QuadraticTriangle, r, s)
    ∂ᵣF = (4r + 4s - 3)*(tri6[1] - tri6[4]) +
               (4r - 1)*(tri6[2] - tri6[4]) +
                   (4s)*(tri6[5] - tri6[6])

    ∂ₛF = (4r + 4s - 3)*(tri6[1] - tri6[6]) +
               (4s - 1)*(tri6[3] - tri6[6]) +
                   (4r)*(tri6[5] - tri6[4])
    return hcat(∂ᵣF, ∂ₛF)
end

function jacobian(quad8::QuadraticQuadrilateral, r, s)
    # Chain rule
    # ∂F   ∂F ∂ξ     ∂F      ∂F   ∂F ∂η     ∂F
    # -- = -- -- = 2 -- ,    -- = -- -- = 2 --
    # ∂r   ∂ξ ∂r     ∂ξ      ∂s   ∂η ∂s     ∂η
    ξ = 2r - 1; η = 2s - 1
    ∂F_∂ξ = (η*(1 - η)/2)*(quad8[1] - quad8[2]) +
            (η*(1 + η)/2)*(quad8[3] - quad8[4]) +
            (ξ*(1 - η)  )*(
                          (quad8[1] - quad8[5]) +
                          (quad8[2] - quad8[5])
                          ) +
            (ξ*(1 + η)  )*(
                          (quad8[3] - quad8[7]) +
                          (quad8[4] - quad8[7])
                          ) +
        ((1 + η)*(1 - η))*(quad8[6] - quad8[8])

    ∂F_∂η = (ξ*(1 - ξ)/2)*(quad8[1] - quad8[4]) +
            (ξ*(1 + ξ)/2)*(quad8[3] - quad8[2]) +
            (η*(1 - ξ)  )*(
                          (quad8[1] - quad8[8]) +
                          (quad8[4] - quad8[8])
                          ) +
            (η*(1 + ξ)  )*(
                          (quad8[3] - quad8[6]) +
                          (quad8[2] - quad8[6])
                          ) +
        ((1 + ξ)*(1 - ξ))*(quad8[7] - quad8[5])

    return hcat(∂F_∂ξ, ∂F_∂η)
end

# 3-polytopes
function jacobian(tet::Tetrahedron, r, s, t)
    ∂ᵣF = tet[2] - tet[1]
    ∂ₛF = tet[3] - tet[1]
    ∂ₜF = tet[4] - tet[1]
    return hcat(∂ᵣF, ∂ₛF, ∂ₜF)
end

function jacobian(hex::Hexahedron, r, s, t)
    ∂ᵣF = ((1 - s)*(1 - t))*(hex[2] - hex[1]) +
          ((    s)*(1 - t))*(hex[3] - hex[4]) +
          ((1 - s)*(    t))*(hex[6] - hex[5]) +
          ((    s)*(    t))*(hex[7] - hex[8])

    ∂ₛF = ((1 - r)*(1 - t))*(hex[4] - hex[1]) +
          ((    r)*(1 - t))*(hex[3] - hex[2]) +
          ((1 - r)*(    t))*(hex[8] - hex[5]) +
          ((    r)*(    t))*(hex[7] - hex[6])

    ∂ₜF = ((1 - r)*(1 - s))*(hex[5] - hex[1]) +
          ((    r)*(1 - s))*(hex[6] - hex[2]) +
          ((1 - r)*(    s))*(hex[8] - hex[4]) +
          ((    r)*(    s))*(hex[7] - hex[3])
    return hcat(∂ᵣF, ∂ₛF, ∂ₜF)
end

function jacobian(tet::QuadraticTetrahedron, r, s, t)
    u = 1 - r - s - t
    ∂ᵣF =      (tet[ 1] - tet[2]) +
          (4u)*(tet[ 5] - tet[1]) +
          (4r)*(tet[ 2] - tet[5]) +
          (4s)*(tet[ 6] - tet[7]) +
          (4t)*(tet[ 9] - tet[8])

    ∂ₛF =      (tet[ 1] - tet[3]) +
          (4u)*(tet[ 7] - tet[1]) +
          (4r)*(tet[ 6] - tet[5]) +
          (4s)*(tet[ 3] - tet[7]) +
          (4t)*(tet[10] - tet[8])

    ∂ₜF =      (tet[ 1] - tet[4]) +
          (4u)*(tet[ 8] - tet[1]) +
          (4r)*(tet[ 9] - tet[5]) +
          (4s)*(tet[10] - tet[7]) +
          (4t)*(tet[ 4] - tet[8])
    return hcat(∂ᵣF, ∂ₛF, ∂ₜF)
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
visualize(x; customizations...) = visualize(string(x); customizations...)
visualize(x::Function; customizations...) = visualize(methods(x); customizations...)

#################################################################################################################################
#Text Rendering:
TEXT_DEFAULTS = @compat Dict(
:Default => @compat Dict(
  :start            => Vec3(0.0),
  :offset           => Vec2(1.0, 1.5), #Multiplicator for advance, newline
  :color            => rgba(248.0/255.0, 248.0/255.0,242.0/255.0, 1.0),
  :backgroundcolor  => rgba(0,0,0,0),
  :model            => eye(Mat4),
  :newline          => -Vec3(0, getfont().props[1][2], 0),
  :advance          => Vec3(getfont().props[1][1], 0, 0),
  :screen           => ROOT_SCREEN,
  :font             => getfont(),
  :stride           => 1024
))

# High Level text rendering for one line or multi line text, which is decided by searching for the occurence of '\n' in text
visualize(text::String,                           style::Style=Style(:Default); customization...) = visualize(style, text, mergedefault!(style, TEXT_DEFAULTS, customization))
# Low level text rendering for multiple line text
visualize(text::Texture{GLGlyph{Uint16}, 4, 2},   style::Style=Style(:Default); customization...) = visualize(style, text, mergedefault!(style, TEXT_DEFAULTS, customization))
visualize(text::Vector{GLGlyph{Uint16}},          style::Style=Style(:Default); customization...) = visualize(style, text, mergedefault!(style, TEXT_DEFAULTS, customization))
visualize(text::Matrix{GLGlyph{Uint16}},          style::Style=Style(:Default); customization...) = visualize(style, text, mergedefault!(style, TEXT_DEFAULTS, customization))

# END Text Rendering
#################################################################################################################################

#################################################################################################################################
#Surface Rendering:

SURFACE_DEFAULTS = @compat Dict(
:Default => @compat Dict(
    :primitive      => SURFACE(),     # can also be CUBES(), CIRCLES(), POINT()
    :x              => (-1,1),        # can also be a matrix
    :y              => (-1,1),        # can also be a matrix
    :color          => rgba(1,0,0,1), # can also be Array/Texture{RGB/RGBA, 1/2}, with "/" meaning OR.
                                      # A 1D Array of color values is assumed to be a colormap.
                                      # A 2D Array can have higher or lower resolution, and will be automatically mapped on the data points.
    :light_position => Vec3(20, 20, -20),

    :screen         => ROOT_SCREEN,
    :model          => eye(Mat4),
))

begin
local PointType = Union(AbstractFixedVector, Real)
# Visualizes a matrix of 1D Values as a surface, whereas the values get interpreted as z-values
visualize{T <: PointType}(zpoints::Matrix{T},        attribute = :z,  style::Style=Style(:Default); customization...) = visualize(style, zpoints, attribute, mergedefault!(style, SURFACE_DEFAULTS, customization))
visualize{T <: PointType}(zpoints::Texture{T, 1, 2}, attribute = :z,  style::Style=Style(:Default); customization...) = visualize(style, zpoints, attribute, mergedefault!(style, SURFACE_DEFAULTS, customization))
visualize{T <: PointType}(x::Matrix{T}, y::Matrix{T}, z::Matrix{T},   style::Style=Style(:Default); customization...) = visualize(style, zpoints, attribute, mergedefault!(style, SURFACE_DEFAULTS, customization))
end
# END Surface Rendering
#################################################################################################################################

#################################################################################################################################
# Image Rendering:
IMAGE_DEFAULTS = @compat(Dict(
:Default => @compat(Dict(
    :normrange      => Vec2(0,1),   # stretch the value 0-1 to normrange: normrange.x + (color * (normrange.y - normrange.x))
    :kernel         => 1f0,         # kernel can be a matrix or a float, whereas the float gets interpreted as a multiplicator
    :model          => eye(Mat4),
    :screen         => ROOT_SCREEN,
    :model          => eye(Mat4),
)),
:GaussFiltered => @compat(Dict(
    :normrange      => Vec2(0,1),   # stretch the value 0-1 to normrange: normrange.x + (color * (normrange.y - normrange.x))
    :kernel         => Float32[1 2 1; 2 4 2; 1 2 1] / 16f0,         # kernel can be a matrix or a float, whereas the float gets interpreted as a multiplicator
    :model          => eye(Mat4),
    :screen         => ROOT_SCREEN,
    :model          => eye(Mat4),
)),
:LaPlace => @compat(Dict(
    :normrange      => Vec2(0,1),   # stretch the value 0-1 to normrange: normrange.x + (color * (normrange.y - normrange.x))
    :kernel         => Float32[-1 -1 -1; -1 9 -1; -1 -1 -1],         # kernel can be a matrix or a float, whereas the float gets interpreted as a multiplicator
    :model          => eye(Mat4),
    :screen         => ROOT_SCREEN,
    :model          => eye(Mat4),
))))
begin
local PixelType = Union(ColorValue, AbstractAlphaColorValue, Images.ColorTypes.AlphaColor)
visualize{T <: PixelType, CDim}(image::Texture{T, CDim, 2}, style::Style=Style(:Default); customization...) = visualize(style, image, mergedefault!(style, IMAGE_DEFAULTS, customization))
end

# END Image Rendering
#################################################################################################################################

#################################################################################################################################
# Volume Rendering:
VOLUME_DEFAULTS = @compat(Dict(
:Default => @compat(Dict(
  :spacing        => [1f0, 1f0, 1f0],
  :stepsize       => 0.001f0,
  :isovalue       => 0.5f0,
  :algorithm      => 1f0,
  :color          => Vec3(0,0,1),
  :light_position => Vec3(2, 2, -2),
  :model          => eye(Mat4),
  :screen         => ROOT_SCREEN
))
))
begin
local PointType = Union(RGB, Real, RGBA)
visualize{T <: PointType}(intensities::Array{T, 3},         style::Style=Style(:Default); customization...) = visualize(style, intensities, mergedefault!(style, VOLUME_DEFAULTS, customization))
visualize{T <: PointType}(intensities::Image{T, 3},         style::Style=Style(:Default); customization...) = visualize(style, intensities, mergedefault!(style, VOLUME_DEFAULTS, customization))
visualize{T <: PointType}(intensities::Texture{T, 1, 3},    style::Style=Style(:Default); customization...) = visualize(style, intensities, mergedefault!(style, VOLUME_DEFAULTS, customization))
end
# END Volume Rendering
#################################################################################################################################

#################################################################################################################################
# Color Rendering:
COLOR_DEFAULTS = @compat(Dict(
:Default => @compat(Dict(
  :screen                   => ROOT_SCREEN,
  :middle                   => Vec2(0.5),

  :swatchsize               => 0.1f0,
  :border_color             => rgba(1, 1, 0.99, 1),
  :border_size              => 0.02f0,

  :hover                    => Input(false),
  :hue_saturation           => Input(false),
  :brightness_transparency  => Input(false),
  :antialiasing_value       => 0.01f0,
  :model                    => scalematrix(Vec3(200,200,1))
))
))
visualize(color::AbstractAlphaColorValue, style::Style=Style(:Default); customization...) = visualize(style, color, mergedefault!(style, COLOR_DEFAULTS, customization))

# END Color Rendering
#################################################################################################################################