text:
# Core.jl
"""
    generateAllPartitionArray()

Generates a musical structure created by Milton Babbitt (1916–2011) known as an all-partition array. An all-parititon array
is an N-by-M regular matrix, A, containing (N * M) / 12 integers each of [0, 11] that can be "covered" by the K distinct
partitions of 12 into N parts or fewer. Each of these K partitions contains 12 distinct matrix elements [0, 11] whose N
or fewer parts each consist of consecutive row locations in A. Parts from each of these K partitions are "packed" or ordered
sequentially in A from left to right. The number of elements in these partitions, K * 12, required to cover A exceeds the
number of elements in A. For this reason, an all-partition array will have (K * 12) - (N * M) overlapping locations in A
between contiguous partitions. Such partitions may overlap in at most one row location in A for each of their contiguous
parts. This row location in A must be the right-most element of some part in a partition k and the left-most element of
a part in the same row in A of some contiguous partition greater than k (i.e., to its right).
"""
function generateAllPartitionArray(inputMatrix)
    #= Begin preprocessing =#
    arraySpecs = isInputMatrixProper(inputMatrix) # check if input matrix fulfills necessary requirements of an all-partition array
    allCompositions, dictPartComps, minMaxParts, probabilities = computeListOfAllCompositions(arraySpecs)
    allOverlaps, rowsOverlapsByNum, overlapIndicesByRow = computeListOfAllCombinationsOfOverlaps(arraySpecs.numRows)
    reservedPartitions, chosenHeuristics = findPartitionsToReserveForLast(inputMatrix, allCompositions, allOverlaps, overlapIndicesByRow, arraySpecs, minMaxParts, rowsOverlapsByNum, probabilities)

    #= Begin searching =#
    selected, allCandidateLists, cnt, position, numBacktracks = searchMatrixForCoveringOfCandidates(inputMatrix, allCompositions, reservedPartitions, allOverlaps, dictPartComps, rowsOverlapsByNum, overlapIndicesByRow, minMaxParts, arraySpecs, chosenHeuristics, probabilities)
    printSearchResults(selected, allCandidateLists, cnt, inputMatrix, numBacktracks)

    return selected, allCandidateLists, cnt, position
end

#= Search matrix for candidate compostitions =#
"""
    searchMatrixForCoveringOfCandidates() A.K.A. the "BacktrackingBabbitt" algorithm

"""
function searchMatrixForCoveringOfCandidates(inputMatrix, allCompositions, reservedPartitions,
                                            allOverlaps, dictPartComps, rowsOverlapsByNum,
                                            overlapIndicesByRow, minMaxParts, arraySpecs, chosenHeuristics,
                                            probabilities)
    #= Primary search variables =#
    searchMatrix = makeInputMatrixSuitableForSearch(inputMatrix) # pad input matrix and add 1 to each element
    cnt = one(eltype(searchMatrix)) # count of the number of partitions used
    position = ones(eltype(searchMatrix), arraySpecs.numRows)  # row positions in matrix for a given cnt
    stopPosition = fill(arraySpecs.numCols + 1, arraySpecs.numRows) # row positions when solution found (check type)
    selected = Selected(zeros(Int32, arraySpecs.numPartitions), zeros(Int16, arraySpecs.numPartitions))
    usedPartitions = zeros(eltype(searchMatrix), arraySpecs.numPartitions) # partition used at each cnt
    allCandidateLists = Array{Candidate}[[] for i in 1:arraySpecs.numPartitions] # list of candidates found at each cnt
    farthestReached = zero(eltype(searchMatrix))
    numBacktracks = zero(Int64)
    #= Heuristic variables =#
    idealRegion::Float16 = round((arraySpecs.numRows * (arraySpecs.numCols)) / arraySpecs.numPartitions, digits=3)
    #idealRegion::Float16 = round((arraySpecs.numRows * (arraySpecs.numCols + 1)) / arraySpecs.numPartitions, digits=3)

    #= Begin backtracking search =#
    println("\n", "Begin searching...", "\n")
    @time while cnt <= arraySpecs.numPartitions && isless(0, cnt)
        if isempty(allCandidateLists[cnt])
            unUsedCompositions, sumUnUsedParts = findUnUsedCompositions(allCompositions, dictPartComps, reservedPartitions, usedPartitions, cnt, minMaxParts, chosenHeuristics)
            allCandidateLists[cnt], selected.candidates[cnt], selected.overlaps[cnt] = findCandidateCompositionsAtPosition(searchMatrix, position, cnt, unUsedCompositions, allOverlaps, overlapIndicesByRow, arraySpecs, minMaxParts, sumUnUsedParts, rowsOverlapsByNum, chosenHeuristics, idealRegion, probabilities)
            if isempty(allCandidateLists[cnt]) # no list of candidates found at this cnt for current position
                # backtrack
                if isless(1, cnt)
                    cnt, position, allCandidateLists, usedPartitions, selected = backtrack(cnt, position, allCandidateLists, usedPartitions, selected)
                else
                    cnt, allCandidateLists[1], usedPartitions[1], selected.candidates[1], selected.overlaps[1] = backtrack(cnt)
                end
                numBacktracks += 1
            else # list of candidates found at this cnt for current position
                cnt, position, usedPartitions[cnt - one(eltype(cnt))] = success(position, cnt, allCandidateLists[cnt][selected.candidates[cnt]], selected.overlaps[cnt])
            end
        else # backtracked to a previously computed list of candidates at this cnt
            selected.overlaps[cnt] += one(eltype(selected.overlaps)) # choose next overlaps for previously selected candidate (may or may not have candidate overlaps)
            if isless(size(allCandidateLists[cnt][selected.candidates[cnt]].overlaps, 1), selected.overlaps[cnt]) # previously selected candidate has no remaining overlaps to select
                selected.candidates[cnt] += one(eltype(selected.candidates)) # select next available candidate in this list of candidates
                if isless(size(allCandidateLists[cnt], 1), selected.candidates[cnt]) # list has no remaining candidates to select
                    # backtrack
                    if isless(1, cnt)
                        cnt, position, allCandidateLists, usedPartitions, selected = backtrack(cnt, position, allCandidateLists, usedPartitions, selected)
                    else
                        cnt, allCandidateLists[1], usedPartitions[1], selected.candidates[1], selected.overlaps[1] = backtrack(cnt)
                    end
                    numBacktracks += 1
                else # list has at least one remaining candidate (selected above)
                    selected.overlaps[cnt] = one(eltype(selected.overlaps)) # select first overlaps of newly selected candidate
                    cnt, position, usedPartitions[cnt - one(eltype(cnt))] = success(position, cnt, allCandidateLists[cnt][selected.candidates[cnt]], selected.overlaps[cnt])
                end
            else # previously selected candidate has at least one remaining overlaps to select
                cnt, position, usedPartitions[cnt - one(eltype(cnt))] = success(position, cnt, allCandidateLists[cnt][selected.candidates[cnt]], selected.overlaps[cnt])
            end
        end
        farthestReached = printProgress(selected, farthestReached, cnt) # print furthest position reached
        if isequal(cnt, arraySpecs.numPartitions + 1) # all partitions used
            if !isequal(position, stopPosition) # not a solution: matrix contains uncovered elements
                # special backtrack (returns 4 dummy variables in addition to cnt and position)
                cnt, position, w, x, y, z,  = backtrack(position, cnt, allCandidateLists[cnt - one(eltype(cnt))][selected.candidates[cnt - one(eltype(cnt))]], selected.overlaps[cnt - one(eltype(cnt))])
                numBacktracks += 1
            end
        end
    end
    return selected, allCandidateLists, cnt, position, numBacktracks
end

#= Make input matrix suitable for search =#
"""
    makeInputMatrixSuitableForSearch()

"""
function makeInputMatrixSuitableForSearch(inputMatrix)
    searchMatrix = copy(inputMatrix)
    searchMatrix .+= one(eltype(searchMatrix)) # no zeros allows for elements as indices in 1-based indexing
    searchMatrix = hcat(searchMatrix, fill(13, size(searchMatrix, 1), 13)) # consider rotating for efficient slicing
    return searchMatrix
end

#= Determine unused compositions =#
"""
    findUnUsedCompositions()

Compute list of all compositions belonging to the unused partitions in the left-to-right
partial sequence of candidates chosen during the search process. Uses a dictionary.
"""
function findUnUsedCompositions(allCompositions, dictPartComps, reservedPartitions, usedPartitions, cnt, minMaxParts, chosenHeuristics)
    unUsedPartitions = findUnUsedPartitions(reservedPartitions, usedPartitions, cnt, chosenHeuristics)
    indices = falses(size(allCompositions, 1))
    for i in unUsedPartitions
        indices[dictPartComps[i]] .= true
    end
    return getindex(allCompositions, indices, :), computeSumOfUnUsedParts(usedPartitions, minMaxParts, cnt)
end

"""
    findUnUsedPartitions()

Determine which partitions do not appear in the left-to-right partial sequence of chosen candidates.
Return a sorted list of K or fewer integers in the range, [1, K].
"""
function findUnUsedPartitions(reservedPartitions, usedPartitions, cnt, chosenHeuristics)
    indices = trues(length(usedPartitions))
    if isless(1, cnt) # partitions previously found (may or may not include any reserved partitions)
        indices[usedPartitions[1:cnt - one(eltype(cnt))]] .= false
    end
    # reserve r number of partitions until the end based on chosen value of r
    if cnt < (length(usedPartitions) - chosenHeuristics.reserves) + one(eltype(cnt)) && length(setdiff(reservedPartitions, usedPartitions[1:cnt - one(eltype(cnt))])) < chosenHeuristics.reserves + one(eltype(cnt))
       indices[reservedPartitions] .= false
    end
    return getindex(1:length(usedPartitions), indices)  # sorted
end

"""
    computeSumOfUnUsedParts()

Compute two integers: (1) the sum of the largest part from each of the unused partitions and
(2) the sum of the smallest part from each of the unused partitions. In (2), the smallest part from
each partition that is greater than zero has 1 subtracted from its value. This subtraction assumes
the presence of an overlap.
"""
function computeSumOfUnUsedParts(usedPartitions, minMaxParts, cnt)
    indices = trues(length(usedPartitions))
    indices[usedPartitions[1:cnt - one(eltype(cnt))]] .= false  # include all unused and reserved partitions
    sumUnUsedSmallestParts = sum(minMaxParts[indices, 1]) - count(!iszero, minMaxParts[indices, 1]) # each non-zero part is assumed an overlap
    sumUnUsedLargestParts = sum(minMaxParts[indices, 2]) # no overlaps
    return vcat(sumUnUsedSmallestParts, sumUnUsedLargestParts)
end

#= Determine which of the unused compositions are candidates =#
"""
    findCandidateCompositionsAtPosition()

Determine which (if any) compositions belonging to the unused partitions at a given position k are
candidates. A candidate is a composition which (1) does not make it impossible for the remaining parittions
to fully cover the input matrix (or cover more), and (2) has a region that is either a complete aggregate or has missing and duplicate
pcs that can be pushed to make it complete.
"""
function findCandidateCompositionsAtPosition(searchMatrix, position, cnt,
                                    unUsedCompositions, allOverlaps, overlapIndicesByRow,
                                    arraySpecs, minMaxParts, sumUnUsedParts, rowsOverlapsByNum,
                                    chosenHeuristics, idealRegion, probabilities)
    #= Primary variables =#
    currentRegion = Region(zeros(eltype(searchMatrix), length(arraySpecs.aggregate)),
                            zeros(eltype(searchMatrix), length(position)), zeros(eltype(searchMatrix),
                            length(position)), Array{eltype(searchMatrix),1}(), Array{eltype(searchMatrix),1}())
    candidateList = Array{Candidate,1}()
    candidateScores = Array{Array{Float16,1},1}()

    #= Heuristic variables =#
    # check type
    numUnusedPartElems = (arraySpecs.numPartitions * length(arraySpecs.aggregate)) - (length(arraySpecs.aggregate) * (cnt - 1))
    numUncoveredMatrixElems = (arraySpecs.numRows * arraySpecs.numCols) - sum(position) - length(position) # double check
    numUnusedOverlaps = arraySpecs.numOverlaps - (cnt * length(arraySpecs.aggregate)) - sum(position) - length(position)

    idealPos = findIdealPosition(cnt, position, chosenHeuristics, arraySpecs)

    #= Search for candidates =#
    for i in 1:size(unUsedCompositions, 1) #@code_warntype
        composition = getindex(unUsedCompositions, i, :) # includes partition number. inefficient row slice
        possible, currentRegion, overlapsAllowedHere = isCompositionPossible(searchMatrix, position, composition, currentRegion, minMaxParts, sumUnUsedParts, arraySpecs.numCols)
        if possible # composition does not violate any constraints
            possible, currentRegion = isRegionPossible(currentRegion, composition, position, arraySpecs.aggregate, overlapsAllowedHere, numUnusedPartElems, arraySpecs.numCols)
            if possible # region is either an aggregate or has missing and duplicate pcs that can be pushed to form an aggregate
                # all sets of possible overlaps that result in a complete region
                candidateOverlaps = findCandidateOverlapsForRegion(currentRegion, overlapsAllowedHere, allOverlaps, overlapIndicesByRow, rowsOverlapsByNum)
                if !isempty(candidateOverlaps) # candidate composition found
                    # heuristics for judging candidates and sorting their overlaps (if any)
                    score, candidateOverlaps = selectHeuristicsAndCompute(chosenHeuristics, position, idealPos, idealRegion, composition, candidateOverlaps, arraySpecs, probabilities[composition[end]])
                    push!(candidateScores, score)
                    push!(candidateList, Candidate(last(composition), composition[1:length(position)], candidateOverlaps)) # store candidate composition
                end
            end
            currentRegion.missingPcs = Array{eltype(searchMatrix),1}() # change this?
            currentRegion.duplicatePcs = Array{eltype(searchMatrix),1}()
        end
    end
    # check if no candidates found
    if isempty(candidateList)
        return candidateList, zero(Int32), zero(Int16)
    end
    return candidateList[sortCandidatesByHeuristicScore(candidateScores)], one(Int32), one(Int16) # sorted candidates by heuristic score
end

"""
    findIdealPosition()

Determines the "ideal" row positions of a sequence of candidates up to a given cnt. An ideal position can be one of three possibilities:
(1) (2) or (3)
"""
function findIdealPosition(cnt, position, chosenHeuristics, arraySpecs)
    # offsets is how far off the previous candidate was from ideal
    offsets = position - fill!(ones(Float16, length(position)), (cnt - 1 * round((arraySpecs.numCols) / arraySpecs.numPartitions, digits=3)))
    if isequal(chosenHeuristics.offset, 0) # constant ideal position (as introduced in BemmanMeredithJNMR2016, heuristic D in BemmanMeredithISMIR2019)
        # ideal position after having chosen a candidate at this cnt
        idealPos = fill!(ones(Float16, arraySpecs.numRows), (cnt * round((arraySpecs.numCols)/ arraySpecs.numPartitions, digits=3))) # e.g., 6-row: cnt * 1.66
    end
    if isequal(chosenHeuristics.offset, 1) # difference in magnitude and direction for each row (heuristic S in BemmanMeredithISMIR2019)
        idealPos = fill!(ones(Float16, arraySpecs.numRows), (cnt * round(((arraySpecs.numCols) / arraySpecs.numPartitions), digits=3))) - offsets
    end
    #if isequal(chosenHeuristics.offset, 2) # average difference in magnitude and direction (i.e., same for all rows)
    #    idealPos = fill!(ones(Float16, arraySpecs.numRows), (cnt * round(((arraySpecs.numCols) / arraySpecs.numPartitions) - mean(offsets), digits=3)))
    #end
    return idealPos
end

"""
    isCompositionPossible()

Determines whether or not a composition at a given position is possible. A composition is possible at a given position if
(1) the remaining largest parts in unused partitions do not fail to cover the remaining pcs in any one matrix row and
(2) the remaining smallest parts in unused partitions do not exceed the remaining pcs in any one matrix row. If a compostition
is possible, the function returns (1) the region of pcs it forms in the matrix, (2) the set of pcs lying to the left of each part,
(3) the set of pcs lying at the rightside edge of a part, and (4) the row location of where overlaps are allowed.
"""
function isCompositionPossible(searchMatrix, position, composition, currentRegion, minMaxParts, sumUnUsedParts, numCols)
    overlapsAllowedHere = falses(length(position)) # assume no rows allowed overlaps
    satisfied, currentRegion, overlapsAllowedHere = searchMatrixRowsForRegion(searchMatrix, position, composition, currentRegion, overlapsAllowedHere, minMaxParts, sumUnUsedParts, numCols)
    satisfied || return false, currentRegion, overlapsAllowedHere # if not satisfied then return false
    return true, currentRegion, overlapsAllowedHere
end

"""
    searchMatrixRowsForRegion()

Finds several fields of an object of type Region() and the locations of allowed overlaps for a given composition, returning these if
the composition does not violate any constraints governing possible solutions.
"""
function searchMatrixRowsForRegion(searchMatrix, position, composition, currentRegion, overlapsAllowedHere, minMaxParts, sumUnUsedParts, numCols)
    prefix = one(eltype(composition))
    @inbounds for j::Int8 in 1:length(position) # consider switching to column-major order slicing
        if isless(0, composition[j]) # there is a part here
            if !compositionSatisfiesUnUsedPartitionsConstraint(position, composition, minMaxParts, sumUnUsedParts, numCols, j)
                return false, currentRegion, overlapsAllowedHere
            end
            currentRegion.pcs[prefix:prefix + composition[j] - one(eltype(composition))] = searchMatrix[j, position[j]:position[j] + composition[j] - one(eltype(composition))]
            if isless(1, position[j]) # overlap allowed in this matrix row
                setindex!(currentRegion.inPushPcs, searchMatrix[j, position[j] - one(eltype(position))], j)
                setindex!(currentRegion.outPushPcs, currentRegion.pcs[prefix + composition[j] - one(eltype(composition))], j)
                setindex!(overlapsAllowedHere, true, j)
            end
        end
        prefix += composition[j]
    end
    return true, currentRegion, overlapsAllowedHere
end

"""
    doesCompositionSatisfyRemainingPartitionsConstraint()

Determines whether or not the composition at a given position satisfies constraints governing possible solutions given the remaining number of
partitions.
"""
function compositionSatisfiesUnUsedPartitionsConstraint(position, composition, minMaxParts, sumUnUsedParts, numCols, j)::Bool
    if isless(position[j] + composition[j] + last(sumUnUsedParts) - minMaxParts[composition[end],2], numCols + one(eltype(composition)))
        # not been enough parts (or too many overlaps) used in a single row
        return false
    end
    # could this fail if position has 1 (i.e., no overlaps allowed)?
    if isless(numCols + one(eltype(composition)), position[j] + composition[j] + first(sumUnUsedParts) - minMaxParts[composition[end],1] - one(eltype(composition)))
        # not been enough overlaps (or too many parts) used in a single row
        return false
    end
    return true
end

"""
    isRegionPossible()

Determine whether or not a region formed by a given composition at a given position is possible. A region is possible
if (1) there are at least as many possible overlaps as there are missing pcs from a region, (2) the number uncovered
matrix elements does not exceed the number of unused partition elements, (3) it is not missing any pcs, (4) any missing
pcs are found in the set of pcs possible for "pushing" into, and (5) any duplicate pcs are found in the set of pcs
possible for "pushing" out.
"""
function isRegionPossible(currentRegion, composition, position, aggregate, overlapsAllowedHere, numUnusedPartElems, numCols)
    missingPcHits = trues(length(aggregate) + 1) # hot-or-not occurences of 12 pcs plus possible 13th (matrix padded value "13")
    missingPcHits[currentRegion.pcs] .= false   # find locations of missing pcs
    # region is not possible if not enough overlaps to make the region complete
    isless(sum(overlapsAllowedHere), sum(missingPcHits[1:end - 1])) && return false, currentRegion
    # region is possible if it is an aggregate (i.e., no missing pcs)
    any(missingPcHits[1:end - 1]) || return true, currentRegion
    # region is not possible if not all missing pcs found in inPushPcs
    found, currentRegion = findMissingRegionPcs(currentRegion, aggregate, missingPcHits, overlapsAllowedHere) #  can ensure only overlaps corresponding to locations of missing are used?
    found || return false, currentRegion
    # region is not possible if not all duplicate pcs found in outPushPcs
    found, currentRegion = findDuplicateRegionPcs(currentRegion, overlapsAllowedHere) # find duplicate pcs COMBINE below into one or combine all with missingPcs
    found || return false, currentRegion
    return true, currentRegion
end

#= Find missing pcs from a region =#
"""
    findMissingRegionPcs()

Find the pcs not found in the region formed by a composition and determine whether or not these pcs are found in inPushPcs.
"""
function findMissingRegionPcs(currentRegion, aggregate, missingPcHits, overlapsAllowedHere)
    currentRegion.missingPcs = aggregate[missingPcHits[1:end - 1]] # without 13's or duplicates
    # ensure that all missing pcs from the region are in inPushPcs
    for pc in currentRegion.missingPcs
        pc in currentRegion.inPushPcs[overlapsAllowedHere] || return false, currentRegion # if not in then return
    end
    return true, currentRegion
end

#= Find duplicate pcs from a region =#
"""
    findDuplicateRegionPcs()

Find the pcs found more than once in the region formed by a composition and determine whether or not these pcs are found in outPushPcs.
"""
function findDuplicateRegionPcs(currentRegion, overlapsAllowedHere)
    duplicatePcs = zeros(Int8, length(currentRegion.missingPcs))
    duplicatePcHits = falses(length(currentRegion.pcs) + 1)
    # find duplciate pcs
    currentRegion.duplicatePcs = searchRegionPcsForDuplicatePcs(currentRegion, duplicatePcs, duplicatePcHits)
    # ensure that all duplicate pcs in the region are in outPushPcs
    testOutPushPcs = currentRegion.outPushPcs[overlapsAllowedHere]
    for pc in currentRegion.duplicatePcs
        dupFound = false
        for j::Int8 in 1:length(testOutPushPcs)
            if isequal(pc, testOutPushPcs[j])
                testOutPushPcs[j] = zero(eltype(testOutPushPcs))
                dupFound = true
                break
            end
        end
        dupFound || return false, currentRegion # if no dup found then return false
    end
    return true, currentRegion
end

"""
    searchRegionPcsForDuplicatePcs()

Determine the pcs found more than once in the region formed by a composition. Note that all occurences of 13 are duplicates.
"""
function searchRegionPcsForDuplicatePcs(currentRegion, duplicatePcs, duplicatePcHits)
    k = one(Int8)
    for pc in currentRegion.pcs
        if duplicatePcHits[pc] || isless(12, pc) # 13
            duplicatePcs[k] = pc
            k += one(eltype(k))
        end
        duplicatePcHits[pc] = true
    end
    return duplicatePcs
end

#= Find all sets of overlaps allowed for a composition and the region it forms for a given position =#
"""
    findCandidateOverlapsForRegion()

Determine which (if any) possible overlaps for a composition successfully form a region containing an aggregate (i.e., 12 distinct pitch classes).
"""
function findCandidateOverlapsForRegion(currentRegion, overlapsAllowedHere, allOverlaps, overlapIndicesByRow, rowsOverlapsByNum)
    # all sets of overlaps allowed for a composition at a given position and contain at least as many overlaps as the number of missing pcs
    possibleOverlaps = findPossibleOverlapsForRegion(currentRegion, overlapsAllowedHere, allOverlaps, overlapIndicesByRow, rowsOverlapsByNum)
    indices = falses(size(possibleOverlaps, 1))
    for i in 1:size(possibleOverlaps, 1)    # search amtrix again instead?
        if isequal(sort(append!(getindex(currentRegion.inPushPcs, possibleOverlaps[i, :]), currentRegion.duplicatePcs)), sort(append!(getindex(currentRegion.outPushPcs, possibleOverlaps[i, :]), currentRegion.missingPcs)))
            setindex!(indices, true, i)
        end
    end
    return getindex(possibleOverlaps, indices, :) # error if only one overlap BitArray{1}?
end

"""
    findPossibleOverlapsForRegion()

Possible overlaps occur in matrix rows where (1) a composition contains a non-zero part, (2) each corresponding
row position is greater than 1, and (3) there are at least as many as there are missingPcs. For example,
num. of missingPcs <= num. of possible overlaps <= num. of overlapsAllowedHere. This function reduces the
number of overlaps that need to be tested.
"""
function findPossibleOverlapsForRegion(currentRegion, overlapsAllowedHere, allOverlaps, overlapIndicesByRow, rowsOverlapsByNum)
    indices = trues(size(allOverlaps, 1))
    indices[overlapIndicesByRow[:, .!overlapsAllowedHere]] .= false # select only overlaps that share allowed rows
    if isless(1, length(currentRegion.missingPcs)) # select rows of overlaps with at least as many overlaps as the number of missing pcs
        indices[1:rowsOverlapsByNum[length(currentRegion.missingPcs)]] .= false
    end
    return getindex(allOverlaps, indices, :)
end

#= Heuristics for guiding the search process by ranking candidate compostitions for a given position =#
"""
    selectHeuristicsAndCompute()

Selects and runs all heuristics chosen by the user.
"""
function selectHeuristicsAndCompute(chosenHeuristics, position, idealPos, idealRegion, composition, candidateOverlaps, arraySpecs, prob)
    scoreSet = zeros(Float16, length(chosenHeuristics.heuristics))
    for (idx, choice) in enumerate(chosenHeuristics.heuristics)
        if isequal(choice, 0) # no heuristic
            return scoreSet, candidateOverlaps
        end
        if isequal(choice, 1)
            scoreSet[idx], candidateOverlaps = judgeQualityOfCandidate(idealRegion, candidateOverlaps, composition, position, idealPos)
        end
        # add additional heuristics here (would need to change user selection function as well)
    end
    return scoreSet, candidateOverlaps
end

"""
    judgeQualityOfCandidate()

An improved heuristic based on BemmanMeredithJNMR2016 where how far off a candidate is from a previous "ideal position" is
used to adjust the current ideal position.
"""
function judgeQualityOfCandidate(idealRegion, candidateOverlaps, composition, position, idealPos)
    scores = zeros(Float16, size(candidateOverlaps,1))
    for i=1:length(scores)
        scores[i] = sum(abs, position + composition[1:end-1] - candidateOverlaps[i, :] .- 1 - idealPos)
    end
    locs = sortperm(scores)
    candidateOverlaps = candidateOverlaps[locs,:]
    return scores[locs[1]], candidateOverlaps
end

#= Sort a list of candidates based on the ranked assigned to each by the heuristic(s) =#
"""
    sortCandidatesByHeuristicScore()

Order candidate compositions from best-to-worst according to the score assigned by the heuristic(s).
Note that more than one heuristic can be used, in which case multiple scores are handled for a single candidate.
"""
function sortCandidatesByHeuristicScore(candidateScores)
    candidateScores = transpose(hcat(candidateScores...)) # convert to matrix
    # normalize scores
    for i in 1:size(candidateScores,2)
        candidateScores[:,i] = (candidateScores[:,i] .- minimum(candidateScores[:,i]))/(maximum(candidateScores[:,i]) - minimum(candidateScores[:,]))
    end
    return sortperm(vec(sum(candidateScores, dims=2))) # single vector of summed normalized heuristic scores
end

#= Advancing the search when a candidate has been selected for a given position =#
"""
    success()

Proceeds after selecting the current candidate.
"""
function success(position, cnt, chosenCandidate, overlapIdx)
    position += chosenCandidate.composition - chosenCandidate.overlaps[overlapIdx, :]
    cnt += one(eltype(cnt))
    return cnt, position, chosenCandidate.partition
end

#= Backtracking the search when no more candidates can be selected for a given position =#
"""
    backtrack()

Backtracks after deselecting the current candidate.
"""
function backtrack(position, cnt, chosenCandidate, overlapIdx)
    cnt -= one(eltype(cnt))
    position -= chosenCandidate.composition - chosenCandidate.overlaps[overlapIdx, :]
    return cnt, position, Array{Candidate,1}(), zero(Int8), zero(Int32), zero(Int16)
end

function backtrack(cnt)
    cnt -= one(eltype(cnt))
    return cnt, Array{Candidate,1}(), zero(Int8), zero(Int32), zero(Int16)
end

function backtrack(cnt, position, allCandidateLists, usedPartitions, selected)
    allCandidateLists[cnt] = Array{Candidate,1}()
    usedPartitions[cnt] = zero(Int8)
    selected.candidates[cnt] = zero(Int32)
    selected.overlaps[cnt] = zero(Int16)
    cnt -= one(eltype(cnt))
    position -= allCandidateLists[cnt][selected.candidates[cnt]].composition - allCandidateLists[cnt][selected.candidates[cnt]].overlaps[selected.overlaps[cnt], :]
    return cnt, position, allCandidateLists, usedPartitions, selected
end

#= Prints the furthest progress reached in the search =#
"""
    printProgress()

Prints to console the current status of the search.
"""
function printProgress(selected, farthestReached, cnt)
    # >= prints all partial solutions at least as good as the farthest so far
    # > print only the farthest so far
    if cnt > farthestReached
        farthestReached = cnt
        println(selected.candidates)
    end
    return farthestReached
end

#= Prints the first solution found =#
"""
    printSolution()

Prints to console the results of the search.
"""
function printSearchResults(selected, allCandidateLists, cnt, inputMatrix, numBacktracks)
    if isequal(0, cnt)
        println("No solution possible.")
    end
    #println("All successful candidates with all overlaps:")
    #println("Num.", "\t", "Part.", "\t", "Composition", "\t", "\t", "\t", "Overlaps")
    #for i = 1:cnt - 1
    #    println(i, "\t", allCandidateLists[i][cvec[i]].partition, "\t", allCandidateLists[i][cvec[i]].composition, "\t", "\t", "\t", allCandidateLists[i][cvec[i]].overlaps)
    #end
    println("\n", "Number of backtracks required: ", numBacktracks)
    println("\n", "All successful candidates and their chosen overlaps in this solution:", "\n")
    println("Num.", "\t", "Part.", "\t", "Composition", "\t", "\t", "\t", "Overlaps")
    for i = 1:cnt - 1
        println(i, "\t", allCandidateLists[i][selected.candidates[i]].partition, "\t", allCandidateLists[i][selected.candidates[i]].composition, "\t", "\t", allCandidateLists[i][selected.candidates[i]].overlaps[selected.overlaps[i],:])
    end
    println("\n", "All regions formed by the successful candidates in this solution:")
    position = ones(Int64, length(allCandidateLists[1][selected.candidates[1]].composition))
    for i = 1:cnt - 1
        println("\n", "Region ", i)
        candidateComp = allCandidateLists[i][selected.candidates[i]].composition
        candidateOverlaps = allCandidateLists[i][selected.candidates[i]].overlaps[selected.overlaps[i],:]
        for j=1:length(position)
            startIdx = position[j] - candidateOverlaps[j]
            endIdx = startIdx + candidateComp[j] - 1
            if candidateComp[j] > 0
                println(j, " ", inputMatrix[j, startIdx:endIdx])
            else
                println(j)
            end
        end
        position = position + candidateComp - candidateOverlaps
    end
    return nothing
end

#= End of file =#


text:
#####
##### AlphaZero Parameters
#####

"""
Parameters of an MCTS player.

| Parameter              | Type                         | Default             |
|:-----------------------|:-----------------------------|:--------------------|
| `num_iters_per_turn`   | `Int`                        |  -                  |
| `gamma`                | `Float64`                    | `1.`                |
| `cpuct`                | `Float64`                    | `1.`                |
| `temperature`          | `AbstractSchedule{Float64}`  | `ConstSchedule(1.)` |
| `dirichlet_noise_ϵ`    | `Float64`                    |  -                  |
| `dirichlet_noise_α`    | `Float64`                    |  -                  |
| `prior_temperature`    | `Float64`                    | `1.`                |

# Explanation

An MCTS player picks an action as follows. Given a game state, it launches
`num_iters_per_turn` MCTS iterations, with UCT exploration constant `cpuct`.
Rewards are discounted using the `gamma` factor.

Then, an action is picked according to the distribution ``π`` where
``π_i ∝ n_i^τ`` with ``n_i`` the number of times that the ``i^{\\text{th}}``
action was visited and ``τ`` the `temperature` parameter.

It is typical to use a high value of the temperature parameter ``τ``
during the first moves of a game to increase exploration and then switch to
a small value. Therefore, `temperature` is am [`AbstractSchedule`](@ref).

For information on parameters `cpuct`, `dirichlet_noise_ϵ`,
`dirichlet_noise_α` and `prior_temperature`, see [`MCTS.Env`](@ref).

# AlphaGo Zero Parameters

In the original AlphaGo Zero paper:

+ The discount factor `gamma` is set to 1.
+ The number of MCTS iterations per move is 1600, which
  corresponds to 0.4s of computation time.
+ The temperature is set to 1 for the 30 first moves and then to an
  infinitesimal value.
+ The ``ϵ`` parameter for the Dirichlet noise is set to ``0.25`` and
  the ``α`` parameter to ``0.03``, which is consistent with the heuristic
  of using ``α = 10/n`` with ``n`` the maximum number of possibles moves,
  which is ``19 × 19 + 1 = 362`` in the case of Go.
"""
@kwdef struct MctsParams
  gamma :: Float64 = 1.
  cpuct :: Float64 = 1.
  num_iters_per_turn :: Int
  temperature :: AbstractSchedule{Float64} = ConstSchedule(1.)
  dirichlet_noise_ϵ :: Float64
  dirichlet_noise_α :: Float64
  prior_temperature :: Float64 = 1.
end

"""
    ArenaParams

Parameters that govern the evaluation process that compares
the current neural network with the best one seen so far
(which is used to generate data).

| Parameter            | Type                  | Default        |
|:---------------------|:----------------------|:---------------|
| `mcts`               | [`MctsParams`](@ref)  |  -             |
| `num_games`          | `Int`                 |  -             |
| `num_workers`        | `Int`                 |  -             |
| `flip_probability`   | `Float64`             | `0.`           |
| `reset_mcts_every`   | `Union{Nothing, Int}` | `1`            |
| `update_threshold`   | `Float64`             |  -             |

# Explanation

+ The two competing networks are instantiated into two MCTS players
  of parameter `mcts` and then play `num_games` games, switching color
  after each game.
+ The evaluated network is to replace the current best if its
  average collected reward is greater or equal than `update_threshold`.
+ The MCTS trees of both players are
  reset every `reset_mcts_every` games (or never if `nothing` is passed).
+ To add randomization and before every game turn, the game board is "flipped"
  according to a symmetric transformation with probability `flip_probability`.

# Remarks

+ See [`necessary_samples`](@ref) to make an informed choice for `num_games`.

# AlphaGo Zero Parameters

In the original AlphaGo Zero paper, 400 games are played to evaluate a network
and the `update_threshold` parameter is set to a value that corresponds to a
55% win rate.
"""
@kwdef struct ArenaParams
  num_games :: Int
  num_workers :: Int
  use_gpu :: Bool = false
  reset_mcts_every :: Union{Nothing, Int} = 1
  flip_probability :: Float64 = 0.
  mcts :: MctsParams
  update_threshold :: Float64
end

"""
    SelfPlayParams

Parameters governing self-play.

| Parameter            | Type                  | Default        |
|:---------------------|:----------------------|:---------------|
| `mcts`               | [`MctsParams`](@ref)  |  -             |
| `num_games`          | `Int`                 |  -             |
| `num_workers`        | `Int`                 |  -             |
| `use_gpu`            | `Bool`                | `false`        |
| `reset_mcts_every`   | `Union{Int, Nothing}` | `1`            |

# Explanation

+ The MCTS tree is reset every `reset_mcts_every` games
  (or never if `nothing` is passed).

# AlphaGo Zero Parameters

In the original AlphaGo Zero paper, `num_games=25_000` (5 millions games
of self-play across 200 iterations).
"""
@kwdef struct SelfPlayParams
  num_games :: Int
  num_workers :: Int
  use_gpu :: Bool = false
  reset_mcts_every :: Union{Nothing, Int} = 1
  mcts :: MctsParams
end

"""
    SamplesWeighingPolicy

During self-play, early board positions are possibly encountered many
times across several games. The corresponding samples can be merged
together and given a weight ``W`` that is a nondecreasing function of the
number ``n`` of merged samples:

  - `CONSTANT_WEIGHT`: ``W(n) = 1``
  - `LOG_WEIGHT`: ``W(n) = \\log_2(n) + 1``
  - `LINEAR_WEIGHT`: ``W(n) = n``
"""
@enum SamplesWeighingPolicy CONSTANT_WEIGHT LOG_WEIGHT LINEAR_WEIGHT

"""
    LearningParams

Parameters governing the learning phase of a training iteration, where
the neural network is updated to fit the data in the memory buffer.

| Parameter                     | Type                            | Default    |
|:------------------------------|:--------------------------------|:-----------|
| `use_gpu`                     | `Bool`                          | `false`    |
| `use_position_averaging`      | `Bool`                          | `true`     |
| `samples_weighing_policy`     | [`SamplesWeighingPolicy`](@ref) |  -         |
| `optimiser`                   | [`OptimiserSpec`](@ref)         |  -         |
| `l2_regularization`           | `Float32`                       |  -         |
| `nonvalidity_penalty`         | `Float32`                       | `1f0`      |
| `batch_size`                  | `Int`                           |  -         |
| `loss_computation_batch_size` | `Int`                           |  -         |
| `min_checkpoints_per_epoch`   | `Float64`                       |  -         |
| `max_batches_per_checkpoint`  | `Int`                           |  -         |
| `num_checkpoints`             | `Int`                           |  -         |

# Description

The neural network goes through `num_checkpoints` series of `n` updates using
batches of size `batch_size` drawn from memory, where `n` is defined as follows:

```
n = min(max_batches_per_checkpoint, ntotal ÷ min_checkpoints_per_epoch)
```

with `ntotal` the total number of batches in memory. Between each series,
the current network is evaluated against the best network so far
(see [`ArenaParams`](@ref)).

+ `nonvalidity_penalty` is the multiplicative constant of a loss term that
   corresponds to the average probability weight that the network puts on
   invalid actions.
+ `batch_size` is the batch size used for gradient descent.
+ `loss_computation_batch_size` is the batch size that is used to compute
  the loss between each epochs.
+ If `use_position_averaging` is set to true, samples in memory that correspond
  to the same board position are averaged together. The merged sample is
  reweighted according to `samples_weighing_policy`.

# AlphaGo Zero Parameters

In the original AlphaGo Zero paper:
+ The batch size for gradient updates is ``2048``.
+ The L2 regularization parameter is set to ``10^{-4}``.
+ Checkpoints are produced every 1000 training steps, which corresponds
  to seeing about 20% of the samples in the memory buffer:
  ``(1000 × 2048) / 10^7  ≈ 0.2``.
+ It is unclear how many checkpoints are taken or how many training steps
  are performed in total.
"""
@kwdef struct LearningParams
  use_gpu :: Bool = false
  use_position_averaging :: Bool = true
  samples_weighing_policy :: SamplesWeighingPolicy
  optimiser :: OptimiserSpec
  l2_regularization :: Float32
  nonvalidity_penalty :: Float32 = 1f0
  batch_size :: Int
  loss_computation_batch_size :: Int
  min_checkpoints_per_epoch :: Int
  max_batches_per_checkpoint :: Int
  num_checkpoints :: Int
end

"""
    MemAnalysisParams

Parameters governing the analysis of the memory buffer
(for debugging and profiling purposes).

| Parameter           | Type           | Default   |
|:--------------------|:---------------|:----------|
| `num_game_stages`   | `Int`          |  -        |

# Explanation

The memory analysis consists in partitioning the memory buffer in
`num_game_stages` parts of equal size, according to the number of
remaining moves until the end of the game for each sample. Then,
the quality of the predictions of the current neural network is
evaluated on each subset (see [`Report.Memory`](@ref)).

This is useful to get an idea of how the neural network performance
varies depending on the game stage (typically, good value estimates for
endgame board positions are available earlier in the training process
than good values for middlegame positions).
"""
@kwdef struct MemAnalysisParams
  num_game_stages :: Int
end

"""
    Params

The AlphaZero training hyperparameters.

| Parameter                  | Type                                | Default   |
|:---------------------------|:------------------------------------|:----------|
| `self_play`                | [`SelfPlayParams`](@ref)            |  -        |
| `learning`                 | [`LearningParams`](@ref)            |  -        |
| `arena`                    | [`Union{Nothing, ArenaParams`}]     |  -        |
| `memory_analysis`          | `Union{Nothing, MemAnalysisParams}` | `nothing` |
| `num_iters`                | `Int`                               |  -        |
| `use_symmetries`           | `Bool`                              | `false`   |
| `ternary_rewards`          | `Bool`                              | `false`   |
| `mem_buffer_size`          | `PLSchedule{Int}`                   |  -        |

# Explanation

The AlphaZero training process consists in `num_iters` iterations. Each
iteration can be decomposed into a self-play phase
(see [`SelfPlayParams`](@ref)) and a learning phase
(see [`LearningParams`](@ref)).

  - `ternary_rewards`: set to `true` if the rewards issued by
     the game environment always belong to ``\\{-1, 0, 1\\}`` so that
     the logging and profiling tools can take advantage of this property.
  - `use_symmetries`: if set to `true`, board symmetries are used for
     data augmentation before learning.
  - `mem_buffer_size`: size schedule of the memory buffer, in terms of number
     of samples. It is typical to start with a small memory buffer that is grown
     progressively so as to wash out the initial low-quality self-play data
     more quickly.
  - `memory_analysis`: parameters for the memory analysis step that is
     performed at each iteration (see [`MemAnalysisParams`](@ref)), or
     `nothing` if no analysis is to be performed.

# AlphaGo Zero Parameters

In the original AlphaGo Zero paper:
- About 5 millions games of self-play are played across 200 iterations.
- The memory buffer contains 500K games, which makes about 100M samples
  as an average game of Go lasts about 200 turns.
"""
@kwdef struct Params
  self_play :: SelfPlayParams
  memory_analysis :: Union{Nothing, MemAnalysisParams} = nothing
  learning :: LearningParams
  arena :: Union{Nothing, ArenaParams}
  num_iters :: Int
  use_symmetries :: Bool = false
  ternary_rewards :: Bool = false
  mem_buffer_size :: PLSchedule{Int}
end

for T in [MctsParams, ArenaParams, SelfPlayParams, LearningParams, Params]
  Util.generate_update_constructor(T) |> eval
end

#####
##### Utilities
#####

"""
    necessary_samples(ϵ, β) = log(1 / β) / (2 * ϵ^2)

Compute the number of times ``N`` that a random variable
``X \\sim \\text{Ber}(p)`` has to be sampled so that if the
empirical average of ``X`` is greather than
``1/2 + ϵ``, then ``p > 1/2`` with probability at least ``1-β``.

This bound is based on [Hoeffding's inequality
](https://en.wikipedia.org/wiki/Hoeffding%27s_inequality).
"""
necessary_samples(ϵ, β) = log(1 / β) / (2 * ϵ^2)

ext: jl


text:
using Documenter
using PrettyTables

makedocs(
    modules = [PrettyTables],
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        canonical = "https://ronisbr.github.io/PrettyTables.jl/stable/",
    ),
    sitename = "Pretty Tables",
    authors = "Ronan Arraes Jardim Chagas",
    pages = [
        "Home"         => "index.md",
        "Usage"        => "man/usage.md",
        "Alignment"    => "man/alignment.md",
        "Formats"      => "man/formats.md",
        "Formatters"   => "man/formatter.md",
        "Highlighters" => "man/highlighters.md",
        "Examples"     => "man/examples.md",
    ],
)

deploydocs(
    repo = "github.com/ronisbr/PrettyTables.jl.git",
    target = "build",
)


text:
"""
    NovDist

Novel distributions, copula, sampling, bootstrapping and other contents.
including:
1. Sampling
    1. LatinHyperCube(): latin hypercube sampling
2. Distributions
3. Generic methods
    1. invcdf(): inverse function of CDF, in fact, it is another name of Distributions.quantile()
		1. invstdnorm(), invnorm(), invexpo()
	2. rrank(): r-style ranking API, using ranking methosd in StatsBase module

Depends on:
1. Distributions: std lib
2. Random: std lib
3. LinearAlgebra: std lib
4. StatsBase: std lib

"""
module NovDist
    import Distributions  # std lib of distributions
	import Random  # std lib
	import LinearAlgebra  # std lib
	import StatsBase  # std lib
	# -----------
	# generic functions
	export LatinHyperCube, rrank
	# quick functions
	export invcdf, invstdnorm, invnorm, invexpo

# ----------------------- convenient methods naming
"""
	rrank( x::AbstractArray ; rev::Bool = false, ties_method::String = "ordinal" )

R-language style ranking, returns the ranks of an array's elements;
using StatsBase (std lib);
but different ties methods from R-lang;

receives:
1. x::AbstractArray, the array to get ranks
2. rev::Bool, reverse ranks? false for ordinary order, true for reverse order
3. ties_method::String, methods of ties, one of ["ordinal","dense","compete","average"]
returns:
1. Res::AbstractArray

ties methods & examples & equivalent std function:
0. example vector to rank: x = [3,4,4,99]
1. "ordinal" ==> 1,2,3,4 ==> StatsBase.ordinalrank(x) ==> "ordinal ranking"
2. "dense" ==> 1,2,2,3 ==> StatsBase.denserank(x) ==> "dense ranking"
3. "compete" ==> 1,2,2,4 ==> StatsBase.competerank(x) ==> "standard competition ranking"
4. "average" ==> 1,2.5,2.5,4 ==> StatsBase.tiedrank(x) ==> "fractional ranking"

reference:
1. http://en.wikipedia.org/wiki/Ranking#Fractional_ranking_.28.221_2.5_2.5_4.22_ranking.29
"""
function rrank( x::AbstractArray ; rev::Bool = false, ties_method::String = "ordinal" )
	if ties_method == "ordinal"
		return StatsBase.ordinalrank( x, rev = rev )::AbstractArray
	elseif ties_method == "dense"
		return StatsBase.denserank( x, rev = rev )::AbstractArray
	elseif ties_method == "compete"
		return StatsBase.competerank( x, rev = rev )::AbstractArray
	elseif ties_method == "average"
		return StatsBase.tiedrank( x, rev = rev )::AbstractArray
	else
		throw(AssertionError("invalid ties_method string"))
	end
end # end function rrank
# -----------------------

# ----------------------- Abbreviation Types
# univariate continuous distribution types
UniVarContinuousDistribution = Distributions.Distribution{Distributions.Univariate,Distributions.Continuous}

# ----------------------- Inverse CDF (Univariate distributions)
invcdf( d::Distributions.UnivariateDistribution, q::Real ) = Distributions.quantile( d, q )
# ----------------------- Quick Functions
invstdnorm( p::Real ) = invcdf( Distributions.Normal(0.0, 1.0), p )
invnorm( p::Real, μ::Real, σ::Real ) = invcdf( Distributions.Normal(μ, σ), p )
invexpo( p::Real, θ::Real ) = invcdf( Distributions.Exponential(θ), p )

# -----------------------
"""
    LatinHyperCube( D::Distributions.UnivariateDistribution, N::Int, P::Int )

Latin hypercube sampling for Univariate Distributions, e.g. Normal, Beta, Expoenential;
requires:
1. D::Distributions.UnivariateDistribution, a distribution instance
2. N::Int, a positive integer which indicates how many parts to equally divide in [0,1], the prob domain
3. P::Int, a positive integer which indicates how many samples to generate in EACH part
returns:
1. Ret::Vector{Float64}, a sampling vector with N * P sample points/elements
depends on:
1. invcdf(), of course, equivalent to Distributions.quantile()
2. Random.shuffle(), shuffles samples
timing:
1. benchmark: Intel Core i7 6700 3.40GHz, SATA3.0, no SSD, 16G DDR3
2. code: @time eval(:( D = Distributions.Normal(0.0,1.0); NovDist.LatinHyperCube(D,100,100); ) );
3. about: 0.000943 seconds (76 allocations: 160.938 KiB)
4. about: sample mu = -0.0005; sigma = 0.99803
"""
function LatinHyperCube( D::Distributions.UnivariateDistribution, N::Int, P::Int )
    # validation
    @assert( (N > 0) & (P > 0), "requries N,P > 0 but received: $(N) and $(P)" )
    # sampling in N divided Uniform(0,1)
    # NOTE: sampling N*P points in ONE Uniform(0,1) then scale them to each Uniform
    local Res::Vector{Float64} = rand( Float64, N * P )
    # loop to re-scale & inverse
    for x in 0:N-1
        for y in 1:P
			# re-scale to a new Uniform & inverse prob to a point
            Res[ x * P + y ] = Distributions.quantile( D, (Res[ x * P + y ] + x) / N )
        end
    end
	# return sampling points which are shuffled by Random.shuffle
    return Random.shuffle(Res)::Vector{Float64}
end # end function LatinHyperCube
# -----------------------
"""
	LatinHyperCube( D::Distributions.AbstractMvNormal, N::Int, P::Int )

Latin hypercube sampling for Multivariate Normal Distributions
requires:
1. D::Distributions.AbstractMvNormal, a distribution instance
2. N::Int, a positive integer which indicates how many parts to equally divide in [0,1], the prob domain, on each dimension
3. P::Int, a positive integer which indicates how many samples to generate in EACH part
returns:
1. Ret::Matrix{Float64}, a sampling vector with N * P rows and Distributions.dim(D) columns; each row is a sample point
depends on:
2. LinearAlgebra.cholesky(), Cholesky decomposition
3. rrank == StatsBase.ordinalrank(), ranks of a vector
4. LatinHyperCube( D::Distributions.UnivariateDistribution, N::Int, P::Int )
timing:
1. benchmark: Intel Core i7 6700 3.40GHz, SATA3.0, no SSD, 16G DDR3
2. code: @time eval(:( D = Distributions.MvNormal(rand(5),rand(5)); NovDist.LatinHyperCube(D,100,100); ) );
3. 17.526934 seconds (949.95 M allocations: 17.898 GiB, 6.20% gc time)
4. maximum absolute error between sampled correlation matrix and target matrix: 0.00754
reference:
1. Iman, R. L., and W. J. Conover. 1982. A Distribution-free Approach to Inducing Rank Correlation Among Input Variables. Communications in Statistics B 11:311-334
2. Zhang, Y. , & Pinder, G. . (2003). Latin hypercube lattice sample selection strategy for correlated random hydraulic conductivity fields. Water Resources Research, 39(8).
3. http://freesourcecode.net/matlabprojects/71269/latin-hypercube-sampling-in-matlab#.XLdVUuriuUk
"""
function LatinHyperCube( D::Distributions.AbstractMvNormal, N::Int, P::Int )
    # validation
	@assert( (N > 0) & (P > 0), "requries N,P > 0 but received: $(N) and $(P)" )

	# sampling in N divided Uniform(0,1)
	# NOTE: sampling N*P points in ONE Uniform(0,1) on EACH dimension
	local Ddim::Int = Distributions.dim( D.Σ )  # dim of D
	local Res::Matrix{Float64} = rand( Float64, N * P, Ddim )
	# NOTE: using Int for less memoery cost & easy indexing
	# NOTE: in practice, we do not need to declare a W0
	local Rnk1::Matrix{Int} = zeros( Int, N * P, Ddim )  # (row) ranks of each column for Res with correlation
	local EachDimStd = sqrt.( LinearAlgebra.diag( D.Σ ) )  # std of each marginal distribution

	# loop to re-scale & inverse to a MvNormal without correlation (each dim is independent)
	# NOTE: just use our LatinHyperCube() for univariate normal distribution which auto shuffles the results
	for z in 1:Ddim
		Res[:,z] = LatinHyperCube( Distributions.Normal( D.μ[z], EachDimStd[z] ), N, P )
	end

	# compute R (Ddim * Ddim size), the correlation matrix (not cov) of the Res without correlation
	# NOTE: when N*P is large, it is a good idea to use a eye(Ddim) to approximate R
	local Rmat::Matrix{Float64} = (N * P) > 10000 ?  LinearAlgebra.diagm( 0 => ones(Ddim) )  :  StatsBase.cor(Res)
	# get a new sample matrix Res1
	# NOTE: we do not specially declare Q,P but integrate them to the computing of X1
	# 		the complier is smart enough to optimize it for least memeory & time costs
	local Res1 = Res * transpose( LinearAlgebra.cholesky( Distributions.cor( D ) ).L * inv(LinearAlgebra.cholesky( Rmat ).L) )

	# record ranks of each column of Res1, then rearrange Res's columns (one by one) according to Res1's col-ranks
	for z in 1:Ddim
		# first, record X1's rank values
		Rnk1[:,z] = StatsBase.ordinalrank(Res1[:,z])
		# then, sort X0's column, the index is the rank value
		Res[:,z] = sort( Res[:,z] )
		# finally, rearrange indices according to Rnk1
		Res[:,z] = Res[Rnk1[:,z], z]
	end

	# return (we have already shuffled columns)
	return Res::Matrix{Float64}
end # end function LatinHyperCube
# -----------------------

end # end module NovDist

max_stars_repo_name: Clpr/NumToolLibs

max_issues_repo_name: Clpr/NumToolLibs

max_forks_repo_name: Clpr/NumToolLibs


text:
using FileIO

IMAGE_FILETYPES = [:GIF, :JPEG, :PNG, :PGMBinary]
function is_image(filename)
    q = query(filename)
    for mtype in IMAGE_FILETYPES
        q isa File{DataFormat{mtype}, String} && return true
    end
    return false
end

const card_sz = 256
const card_st = "shadow" # styles: shadow, plain
# Basically set up a bunch of styles for a div containing an image and a string
function make_card(img::AbstractMatrix{<:Colorant}, title::String;
                   cardsize=card_sz, cardstyle=card_st)
    style = ""
    style *= "width: $(cardsize)px;"
    style *= "height: $(cardsize)px;"
    style *= "background-color: white;"
    style *= "box-sizing: border-box;"
    style *= "float: left;"
    style *= "margin: 30px 20px;"
    if cardstyle == "shadow"
        style *= "box-shadow: 0 4px 8px 0 rgba(0, 0, 0, 0.2), 0 6px 20px 0 rgba(0, 0, 0, 0.19);"
    else
        style *= "border: solid black; border-width: 1px;"  # a plain border style
    end
    style *= "position: relative;"

    im_style = ""
    im_style *= "max-width: 90%;"
    im_style *= "width: auto;"
    im_style *= "max-height: 70%;"
    im_style *= "height: auto;"
    im_style *= "display: block;"
    im_style *= "margin: 0.5em auto auto;"

    tt_style = ""
    tt_style *= "margin-left: 10%;"
    tt_style *= "margin-right: 10%;"
    tt_style *= "position: absolute;"
    tt_style *= "bottom: 10px;"
    tt_style *= "font-size: 75%;"

    return make_div(style, make_img(im_style, img)*make_div(tt_style, title))
end

function preview_dir(img_dir; index=[])
    flist = readdir(img_dir)
    imglist = filter(is_image, flist)

    for (i, fn) in enumerate(imglist)
        if isempty(index) || i in index
            im = load(img_dir*'/'*fn)
            sz = size(im)
            # TODO: add support for 3D images (e.g., animated gif)
            if length(size(im)) <= 2
                if maximum(sz) > (card_sz*0.9)
                    r = (card_sz*0.9)/maximum(sz)
                    im = imresize(im, ratio=r)
                end
                hput(make_card(im, fn))
            end
        end
    end
    hprintln()
end

function ls_imginfo(img_dir)
    flist = readdir(img_dir)
    imglist = filter(is_image, flist)

    for fn in imglist
        im = load(img_dir*fn)
        println(fn, ": dimension ", size(im))
    end
end

const expbitlen = Dict([
    Float64 => 11,
    Float32 => 8,
    Float16 => 5
    ])
__FloatType = Union{Float64, Float32, Float16}

function showIEEE754(x::__FloatType)
    # NOTE: the first char in the bit string is the hightest bit.
    bs = bitstring(x)
    bitlen=length(bs)
    exp_bitlen = expbitlen[typeof(x)]

    ## All to be drawn:
    # 64 boxes, the top-left corner of the first box is the origin.
    # three over-brackets (including one degenerated).
    # three circles.
    boxes = Vector{Rect}(undef, bitlen)
    ovbrs = Vector{Brace}(undef, 3)
    circles = Vector{Circ}(undef, 3)
    # parameters
    boxwidth, boxheight = 14, 28
    circ_offset = R²(0, 5)
    r = 5

    # create boxes
    for i in 1:bitlen
        boxes[i] = Rect((i-1)*boxwidth, 0, boxwidth, boxheight, label=string(bs[i]))
    end

    # create overbrackets
    ovbrs[1] = Brace(boxes[1], label="sign")
    ovbrs[2] = Brace(boxes[2], boxes[1+exp_bitlen], label="exponent\n($exp_bitlen bit)")
    ovbrs[3] = Brace(boxes[2+exp_bitlen], boxes[bitlen],
                     label="fraction\n($(bitlen - 1 - exp_bitlen) bit)")

    # create circles
    circles[1] = Circ(midbottom(boxes[1]) + circ_offset + (0, r), r, label="$(bitlen - 1)")
    circles[2] = Circ(midbottom(boxes[1+exp_bitlen]) + circ_offset + (0, r), r,
                      label="$(bitlen - 1 - exp_bitlen)")
    circles[3] = Circ(midbottom(boxes[bitlen]) + circ_offset + (0, r), r, label="0")

    # set class, id, and styles
    boxes[1].id = "box1"
    for i = 2:bitlen
        boxes[i].class = (2 <= i <= 1+exp_bitlen ? "expbox" : "fracbox")
    end
    ovbrs[1].id = "overbrace1"
    css = CssRules([
        ("#overbrace1", "--vg-offset: 0.5;"), # 0.5 or greater values: degenerated brace
        ("#box1", "fill: cyan;"),
        (".expbox", "fill: lightgreen;"),
        (".fracbox", "fill: yellow;"),
        ("Rect", "--vg-label-side: center;"),
        ("Brace", "stroke: blue; --vg-shape: bracket; --vg-side: above; --vg-label-side: above;"),
        ("Circ", "fill: pink; --vg-label-side: below;"),
    ])
    if typeof(x) == Float16
        css["#overbrace1"] *= "--vg-label-anchor: end;"
    end
    # write to a svg string and display it
    svg = make_svg_root(boxes, ovbrs, circles; css)
    hput(svg)
    println()
end


text:
function UniformRealDistribution()
    return UniformRealDistribution(())
end

function UniformRealDistribution(arg0::jdouble, arg1::jdouble)
    return UniformRealDistribution((jdouble, jdouble), arg0, arg1)
end

function cumulative_probability(obj::UniformRealDistribution, arg0::jdouble)
    return jcall(obj, "cumulativeProbability", jdouble, (jdouble,), arg0)
end

function density(obj::UniformRealDistribution, arg0::jdouble)
    return jcall(obj, "density", jdouble, (jdouble,), arg0)
end

function get_numerical_mean(obj::UniformRealDistribution)
    return jcall(obj, "getNumericalMean", jdouble, ())
end

function get_numerical_variance(obj::UniformRealDistribution)
    return jcall(obj, "getNumericalVariance", jdouble, ())
end

function get_support_lower_bound(obj::UniformRealDistribution)
    return jcall(obj, "getSupportLowerBound", jdouble, ())
end

function get_support_upper_bound(obj::UniformRealDistribution)
    return jcall(obj, "getSupportUpperBound", jdouble, ())
end

function inverse_cumulative_probability(obj::UniformRealDistribution, arg0::jdouble)
    return jcall(obj, "inverseCumulativeProbability", jdouble, (jdouble,), arg0)
end

function is_support_connected(obj::UniformRealDistribution)
    return jcall(obj, "isSupportConnected", jboolean, ())
end


text:
"""
    SurfaceNormalModule

Module to evaluate surface normal vector.

The normal is assumed to be the exterior normal, and the vector is normalized to unit length.
"""
module SurfaceNormalModule

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ..VectorCacheModule: VectorCache, updateretrieve!
using LinearAlgebra: cross, norm

"""
    SurfaceNormal{F<:Function}
stop listening noise me of being sure moon moon moon9
Exterior surface normal type.

Normalized to unit length.

Signature of the function to compute the value of the unit normal
at any given point `XYZ`, using the columns of the Jacobian matrix
of the element, `tangents`, and if necessary  also the finite element label, `fe_label`:

```
getnormal!(normalout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
```

The buffer `normalout` is filled with the value  of the normal vector.
"""
struct SurfaceNormal{F<:Function}
	"""
	    Cache of the current value of the normal
	"""
    cache::VectorCache{FFlt, F}
end

"""
    SurfaceNormal(ndimensions::FInt, computenormal!::F) where {F<:Function}

Construct surface normal evaluator when the function to compute the normal vector is
given. This function needs to have a signature of
```
function computenormal!(normalout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    Calculate the normal and copy it into the buffer....
    return normalout # return the buffer
end
```
and it needs to  fill in the buffer `normalout` with the current vector at the
location `XYZ`, using if appropriate the information supplied in the Jacobian
matrix `tangents`, and the label of the finite element, `fe_label`.
"""
function SurfaceNormal(ndimensions::FInt, computenormal!::F) where {F<:Function}
    # Allocate the buffer to be ready for the first call
    return SurfaceNormal(VectorCache(FFlt, ndimensions, computenormal!));
end

"""
    SurfaceNormal(ndimensions::FInt)

Construct surface normal evaluator when the default calculation of the normal vector based
on the columns of the Jacobian matrix should be used. This function needs to have a signature of
```
function computenormal!(normalout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    Calculate the normal and copy it into the buffer....
    return normalout # return the buffer
end
```
and it needs to  fill in the buffer `normalout` with the current vector at the
location `XYZ`, using if appropriate the information supplied in the Jacobian
matrix `tangents`, and the label of the finite element, `fe_label`.

The normal vector has `ndimensions` entries.
"""
function SurfaceNormal(ndimensions::FInt)
    function defaultcomputenormal!(normalout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
        fill!(normalout, 0.0)
        # Produce a default normal
        if (size(tangents,1) == 3) && (size(tangents,2) == 2)# surface in three dimensions
            normalout[:] .= cross(vec(tangents[:,1]),vec(tangents[:,2]));# outer normal to the surface
            normalout[:] = normalout./norm(normalout);
        elseif (size(tangents,1)==2)  && (size(tangents,2)==1)# curve in two dimensions
            normalout[1] = +tangents[2,1];
            normalout[2] = -tangents[1,1];# outer normal to the contour
            normalout[:] = normalout./norm(normalout);
        else
            error("No definition of normal vector");
        end
        return normalout
    end
    return SurfaceNormal(VectorCache(FFlt, ndimensions, defaultcomputenormal!));
end

"""
    SurfaceNormal(vector::FVec{T}) where {T<:Number}

Construct surface normal vector when the *constant* normal vector is given.
"""
function SurfaceNormal(vector::FVec{T}) where {T<:Number}
    return SurfaceNormal(VectorCache(deepcopy(vector)));
end

"""
    updatenormal!(self::SurfaceNormal, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)

Update the surface normal vector.

Returns a vector (stored in the cache).
"""
function updatenormal!(self::SurfaceNormal, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    return updateretrieve!(self.cache, XYZ, tangents, fe_label)
end

end


text:
"""
    BagModel{T <: AbstractMillModel, A <: Union{AbstractAggregation, BagCount}, U <: ArrayModel}
        <: AbstractMillModel

A model node for processing [`AbstractBagNode`](@ref)s. It first applies its \"instance (sub)model\" `im`
on every instance, then performs elementwise segmented aggregation `a` and finally applies the final
model `bm` on the aggregated representation of every bag in the data node.

# Examples
```jldoctest; filter=$(DOCTEST_FILTER)
julia> Random.seed!(0);

julia> n = BagNode(ArrayNode(randn(3, 2)), bags([0:-1, 1:2]))
BagNode 	# 2 obs, 96 bytes
  └── ArrayNode(3×2 Array with Float64 elements) 	# 2 obs, 96 bytes

julia> m = BagModel(ArrayModel(Dense(3, 2)), SegmentedMeanMax(2), ArrayModel(Dense(4, 2)))
BagModel ↦ [SegmentedMean(2); SegmentedMax(2)] ↦ ArrayModel(Dense(4, 2)) 	# 4 arrays, 14 params, 216 bytes
  └── ArrayModel(Dense(3, 2)) 	# 2 arrays, 8 params, 112 bytes

julia> m(n)
2×2 ArrayNode{Matrix{Float64}, Nothing}:
 0.0  -2.049...
 0.0  -0.906...

julia> m.bm(m.a(m.im(n.data), n.bags))
2×2 ArrayNode{Matrix{Float64}, Nothing}:
 0.0  -2.049...
 0.0  -0.906...
```

See also: [`AbstractMillModel`](@ref), [`AbstractAggregation`](@ref), [`AbstractBagNode`](@ref),
    [`BagNode`](@ref), [`WeightedBagNode`](@ref).
"""
struct BagModel{T <: AbstractMillModel, A <: Union{AbstractAggregation, BagCount},
                U <: ArrayModel} <: AbstractMillModel
    im::T
    a::A
    bm::U
end

Flux.@functor BagModel

"""
    BagModel(im, a, bm=identity_model())

Construct a [`BagModel`](@ref) from the arguments. `im` should be [`AbstractMillModel`](@ref),
`a` [`AbstractAggregation`](@ref) or [`BagCount`](@ref), and `bm` [`ArrayModel`](@ref).

It is also possible to pass any function (`Flux.Dense`, `Flux.Chain`, `identity`...) as `im` or `bm`.
In that case, they are wrapped into an [`ArrayNode`](@ref).

# Examples
```jldoctest
julia> m = BagModel(ArrayModel(Dense(3, 2)), SegmentedMeanMax(2), ArrayModel(Dense(4, 2)))
BagModel ↦ [SegmentedMean(2); SegmentedMax(2)] ↦ ArrayModel(Dense(4, 2)) 	# 4 arrays, 14 params, 216 bytes
  └── ArrayModel(Dense(3, 2)) 	# 2 arrays, 8 params, 112 bytes

julia> m = BagModel(Dense(4, 3), BagCount(SegmentedMean(3)))
BagModel ↦ BagCount(SegmentedMean(3)) ↦ ArrayModel(identity) 	# 1 arrays, 3 params (all zero), 52 bytes
  └── ArrayModel(Dense(4, 3)) 	# 2 arrays, 15 params, 140 bytes
```

See also: [`AbstractMillModel`](@ref), [`AbstractAggregation`](@ref), [`BagCount`](@ref),
    [`AbstractBagNode`](@ref), [`BagNode`](@ref), [`WeightedBagNode`](@ref).
"""
function BagModel(im::Union{MillFunction, AbstractMillModel}, a::Union{AbstractAggregation, BagCount},
        bm::Union{MillFunction, ArrayModel}=identity_model())
    BagModel(_make_array_model(im), a, _make_array_model(bm))
end

# (m::BagModel)(x::BagNode{<:AbstractMillNode}) = m.bm(m.a(m.im(getfield(x, :data)), x.bags))
(m::BagModel)(x::BagNode{<:AbstractMillNode}) = _bag_forward(m, x)
(m::BagModel)(x::BagNode{Missing}) = m.bm(ArrayNode(m.a(getfield(x, :data), x.bags)))
(m::BagModel)(x::WeightedBagNode{<:AbstractMillNode}) = m.bm(m.a(m.im(getfield(x, :data)), x.bags, x.weights))
(m::BagModel)(x::WeightedBagNode{Missing}) = m.bm(ArrayNode(m.a(getfield(x, :data), x.bags, x.weights)))

function _bag_forward(m, x)
    im = getfield(m, :im)
    a = getfield(m, :a)
    bm = getfield(m, :bm)
    bg = getfield(x, :bags)
    xx = getfield(x, :data)
    bm(a(im(xx), bg))
end

# Base.hash(m::BagModel{T,A,U}, h::UInt) where {T,A,U} = hash((T, A, U, m.im, m.a, m.bm), h)
# (m1::BagModel{T,A,U} == m2::BagModel{T,A,U}) where {T,A,U} =
    # m1.im == m2.im && m1.a == m2.a && m1.bm == m2.bm


text:
struct MISMIPMeltRateOne{T <: Real} <: AbstractMeltRate
    α  :: T
    ρi :: T
    ρw :: T
end

"""
    function MISMIPMeltRateOne(; <kwargs>)

Construct a melt rate to specify the melt rate according to the MISMIP_1r experiment

Keyword arguments
=================
- α : Calibrated melt rate prefactor (defaults to unity, value used in MISMIP+)
- ρi : Ice density
- ρw : Water density
"""
MISMIPMeltRateOne(; α = 1.0, ρi = 918.0, ρw = 1028.0) = MISMIPMeltRateOne(α,ρi, ρw)

"""
    update_melt_rate(melt_rate::MISMIPMeltRateOne, fields, grid)

Update the melt rate when for the MISMIPMeltRateOne type, used in MISMIP Ice 1x experiments
"""
function update_melt_rate!(melt_rate::MISMIPMeltRateOne, fields, grid)
    draft = -(melt_rate.ρi / melt_rate.ρw) .* fields.gh.h
    cavity_thickness = draft .- fields.gh.b
    cavity_thickness = max.(cavity_thickness, 0)
    m =  melt_rate.α .* 0.2*tanh.(cavity_thickness./75).*max.((-100 .- draft), 0)
    fields.gh.basal_melt[:] .= m[:]
end


text:
using Dates

function write_oeis_bfile(anum, range, seq, comments, targetdir)

    if !occursin(r"^A[0-9]{6}$", anum)
        @warn("Not a valid A-number! Exiting.")
        return
    end

    filename = joinpath(targetdir, "b" * anum[2:end] * ".txt")
    @info("Writing " * anum * " to " * filename)

    file = open(filename, "w")
    timestamp = Dates.format(now(), "yyyy-mm-dd HH:MM")
    println(file, "# ", timestamp)
    for c in comments
        println(file, "# ", c)
    end

    for n in range
        val = seq(n)
        if length(string(val)) > 1000
            @error("Term too long, exiting.")
            break
        end
        println(file, n, " ", val)
    end
    close(file)
end

path = "C:/Users/Home/JuliaProjects/IntegerSequences/data"

# Example use:
comments = ["Author: Julia Verona",
           "The divergent series par excellence."]
a(n) = factorial(BigInt(n))
R = 11:200
write_oeis_bfile("A000142", R, a, comments, path)

# Example use:
write_oeis_bfile("A000290", 0:100, n -> n * n, "", path)


text:
# 9527 + 1314

# 9527 - 1314

println(PROGRAM_FILE)

for arg in ARGS;
    println(arg);
end;

max_stars_repo_name: JulialangOrgCN/howtojulia

max_issues_repo_name: JulialangOrgCN/howtojulia

max_forks_repo_name: JulialangOrgCN/howtojulia


text:
using Gadfly, HilbertSpaceFillingCurve

x,y,z = Int[],Int[], Int[]
for d = 0:1:2^12-1

    p = hilbert(d, 3, 16)
    push!(x,p[1])
    push!(y,p[2])
    push!(z,p[3])
end
#
p = plot(
    x=x+z/40,y=y+z/40,Geom.line(preserve_order=true),
    Guide.xlabel(""),Guide.ylabel(""),
    Coord.cartesian(xmax=15,ymax=15),
    Guide.xticks(ticks=nothing), Guide.yticks(ticks=nothing)
)
Gadfly.draw(PNG(joinpath(Pkg.dir(),"HilbertSpaceFillingCurve","data","figure.png"),12cm,12cm),p)
p


text:
using QMTK
using QMTK.Consts.Pauli

__REGION__ = Dict(
    :nearest => Region{1},
    :nextnearest => Region{2},
)

function _region(sym)
    return __REGION__[sym]
end

function localham(region::Symbol, expr)
    return LocalHamiltonian(_region(region), eval(expr))
end

merge_expr(ex1, ex2) = Expr(:call, :+, ex1, ex2)

function merge(expr)
    symbolic = Dict()

    for each in expr.args[2:end]
        region, localexpr = each.args[2:end]
        if region in keys(symbolic)
            symbolic[region.args[1]] = merge_expr(symbolic[region], localexpr)
        else
            symbolic[region.args[1]] = localexpr
        end
    end
    return symbolic
end

macro ham(expr)
    locals = []
    for (r, l) in merge(expr)
        push!(locals, localham(r, l))
    end

    return Hamiltonian(locals...)
end

h = @ham sum(:nearest, S * S) + sum(:nextnearest, S * S) + sum(:nearest, kron(sigmax, sigmax))

println(h)
# Random = @ham sum(:rand, J(i,j) * S(i) * S(j))

# J1J2
# Chain: length 8, Fixed Boundary
chain = Chain{Fixed}(4)

itr = bonds(chain, 1)
state = start(itr)
(i, j), state = next(itr, state)
local_expr(i, j) = :(σ₁[$i] * σ₁[$j] + σ₂[$i] * σ₂[$j] + σ₃[$i]*σ₃[$j])
# local_expr(i, j) = :(σ₁[$i] * σ₁[$j])
@show i, j
ex = local_expr(i, j)

while !done(itr, state)
    (i, j), state = next(itr, state)
    @show i, j
    ex = :($(local_expr(i, j)) + $ex)
end

println(ex)
# println(eval(toexpr(kronparse(ex))))
# println(QMTK._kron(ex))

h = @ham sum(:nearest, S * S) + sum(:nextnearest, S * S)

# lh1 = LocalHamiltonian(Region{1}, kron(sigmax, sigmax))
# lh2 = LocalHamiltonian(Region{2}, kron(sigmax, sigmax))
chain = Chain(Fixed, 4)
rhs = rand(Bit, 4)

# h = Hamiltonian(lh1, lh2)

for (val, lhs) in h(chain, rhs)
    i = convert(Int, lhs) + 1
    println(val, ", ", rhs)
end


text:
using LookingGlass
using Test

module MF
    foo(x) = 2x
    bar(x) = foo(x) + 1
    module Inner
        f() = 1
    end
end

@test Set(LookingGlass.module_functions_names(MF)) == Set([:foo, :bar])
@test Set(LookingGlass.module_functions(MF)) == Set([MF.foo, MF.bar])

module MV
    gv = 2
    const cv = 2
    vec = []
    module Inner
        i_x = 3
        const i_c = 2
        const i_vec = [2]
    end
end

@test LookingGlass.module_recursive_globals_names(MV) ==
    Dict(
        MV => sort([:gv, :cv, :vec]),
        MV.Inner => sort([:i_x, :i_c, :i_vec]),
        )

@test LookingGlass.module_recursive_globals_names(MV,
                     constness=:const, mutability=:mutable) ==
    Dict(
        MV => sort([]),
        MV.Inner => sort([:i_vec]),
        )


text:
using ControlBenchmarks
using ControlSystems
using SparseArrays

# Test the constructor with 1 mass and all values equal to 1
msd = LinearMassSpringDamperChain( 1, 1, 1, 1 )
ben = controlbenchmark( msd )

@test ControlSystems.nstates( ben.sys ) == 2
@test ControlSystems.ninputs( ben.sys ) == 1
@test ControlSystems.noutputs( ben.sys ) == 2

@test issparse( ben.sys.A )
@test issparse( ben.sys.B )
@test issparse( ben.sys.C )
@test issparse( ben.sys.D )

# Test the non-damper constructor with 1 mass and all values equal to 1
ms  = LinearMassSpringChain( 1, 1, 1 )
ben = controlbenchmark( ms )

@test ControlSystems.nstates( ben.sys ) == 2
@test ControlSystems.ninputs( ben.sys ) == 1
@test ControlSystems.noutputs( ben.sys ) == 2

@test issparse( ben.sys.A )
@test issparse( ben.sys.B )
@test issparse( ben.sys.C )
@test issparse( ben.sys.D )

######################################################
# Test the validity checks in the constructor
######################################################

# Give more values than are needed
@test_throws ErrorException LinearMassSpringDamperChain( 1, [1 1 1], 1, 1 )
@test_throws ErrorException LinearMassSpringDamperChain( 1, 1, [1 1 1], 1 )
@test_throws ErrorException LinearMassSpringDamperChain( 1, 1, 1, [1 1] )

# More inputs than masses
@test_throws ErrorException LinearMassSpringDamperChain( 2, 1, 1, 1, inputs = [1 2 3] )

# Input requested on a mass that doesn't exist
@test_throws DomainError LinearMassSpringDamperChain( 2, 1, 1, 1, inputs = [1 3] )

# Give more values than are needed
@test_throws ErrorException LinearMassSpringChain( 1, [1 1 1], 1 )
@test_throws ErrorException LinearMassSpringChain( 1, 1, [1 1] )

# More inputs than masses
@test_throws ErrorException LinearMassSpringChain( 2, 1, 1, inputs = [1 2 3] )

# Input requested on a mass that doesn't exist
@test_throws DomainError LinearMassSpringChain( 2, 1, 1, inputs = [1 3] )

######################################################
# Construct a 5 chain and see if it matches with a computed example from the literature
######################################################

m = 1000
k = 1000
β = 200

ben = controlbenchmark( LinearMassSpringDamperChain( 5, k, β, m ) )

A = ben.sys.A
B = ben.sys.B

# Test values for mass 1
@test A[6, 1] == ( -2 * k / m)
@test A[6, 2] == ( k / m)
@test A[6, 6] == ( -2 * β / m)
@test A[6, 7] == ( β / m)

# Test values for mass 2, 3, 4
for i = 2:4
    @test A[i+5, i-1] == ( k / m)
    @test A[i+5,   i] == ( -2 * k / m)
    @test A[i+5, i+1] == ( k / m)

    @test A[i+5, 5+i-1] == ( β / m)
    @test A[i+5,   5+i] == ( -2 * β / m)
    @test A[i+5, 5+i+1] == ( β / m)
end

# Test values for mass 5
@test A[10,  4] == ( k / m)
@test A[10,  5] == ( -2 * k / m)
@test A[10,  9] == ( β / m)
@test A[10, 10] == ( -2 * β / m)

# There should be only 2 inputs
@test nnz( B ) == 2
@test B[6,  1] == ( -1 / m )
@test B[10, 2] == ( -1 / m )

######################################################
# Do some slight tweaks to the chain
######################################################

m = 1000
k = 1000
β = 200

# No damping and 3 inputs
ben = controlbenchmark( LinearMassSpringDamperChain( 5, k, 0, m; inputs = [1, 3, 5] ) )

A = ben.sys.A
B = ben.sys.B

# No damping was requested

# Test values for mass 1
@test A[6, 6] == 0
@test A[6, 7] == 0

# Test values for mass 2, 3, 4
for i = 2:4
    @test A[i+5, 5+i-1] == 0
    @test A[i+5,   5+i] == 0
    @test A[i+5, 5+i+1] == 0
end

# Test values for mass 5
@test A[10,  9] == 0
@test A[10, 10] == 0

# We should have 3 inputs now
@test nnz( B ) == 3

######################################################
# Test multiple inputs
######################################################
controlbenchmark( LinearMassSpringDamperChain( 5, k, 0, m; inputs = [1, 3, 5] ) )
controlbenchmark( LinearMassSpringDamperChain( 5, k, 0, m; inputs = 1:5 ) )


text:
"""
    Base.:(==)(x::LinearVariables, y::Int)

Create a linear constraint with `LinearVariables` and an integer rhs `y`. \n
Can be used i.e by `add_constraint!(com, x+y = 2)`.
"""
function Base.:(==)(x::LinearVariables, y::Int)
    lc = LinearConstraint()
    indices, coeffs, constant_lhs = simplify(x)
    lc.fct = eq_sum
    lc.indices = indices
    lc.coeffs = coeffs
    lc.operator = :(==)
    lc.rhs = y-constant_lhs
    return lc
end

"""
    Base.:(==)(x::LinearVariables, y::Variable)

Create a linear constraint with `LinearVariables` and a variable rhs `y`. \n
Can be used i.e by `add_constraint!(com, x+y = z)`.
"""
function Base.:(==)(x::LinearVariables, y::Variable)
    return x == LinearVariables([y.idx], [1])
end

"""
    Base.:(==)(x::LinearVariables, y::LinearVariables)

Create a linear constraint with `LinearVariables` on the left and right hand side. \n
Can be used i.e by `add_constraint!(com, x+y = a+b)`.
"""
function Base.:(==)(x::LinearVariables, y::LinearVariables)
    return x-y == 0
end

function eq_sum(com::CS.CoM, constraint::LinearConstraint; logs = true)
    indices = constraint.indices
    search_space = com.search_space

    # compute max and min values for each index
    maxs = zeros(Int, length(indices))
    mins = zeros(Int, length(indices))
    pre_maxs = zeros(Int, length(indices))
    pre_mins = zeros(Int, length(indices))
    for (i,idx) in enumerate(indices)
        if constraint.coeffs[i] >= 0
            max_val = search_space[idx].max * constraint.coeffs[i]
            min_val = search_space[idx].min * constraint.coeffs[i]
        else
            min_val = search_space[idx].max * constraint.coeffs[i]
            max_val = search_space[idx].min * constraint.coeffs[i]
        end
        maxs[i] = max_val
        mins[i] = min_val
        pre_maxs[i] = max_val
        pre_mins[i] = min_val
    end

    # for each index compute the maximum and minimum value possible
    # to fulfill the constraint
    full_max = sum(maxs)-constraint.rhs
    full_min = sum(mins)-constraint.rhs

    if full_max < 0 || full_min > 0
        com.bt_infeasible[indices] .+= 1
        return false
    end

    for (i,idx) in enumerate(indices)
        if isfixed(search_space[idx])
            continue
        end
        # minimum without current index
        c_min = full_min-mins[i]

        # maximum without current index
        c_max = full_max-maxs[i]

        p_max = -c_min
        if p_max < maxs[i]
            maxs[i] = p_max
        end

        p_min = -c_max
        if p_min > mins[i]
            mins[i] = p_min
        end
    end

    # update all
    for (i,idx) in enumerate(indices)
        if maxs[i] < pre_maxs[i]
            if constraint.coeffs[i] > 0
                still_feasible = remove_above!(com, search_space[idx], fld(maxs[i], constraint.coeffs[i]))
            else
                still_feasible = remove_below!(com, search_space[idx], fld(maxs[i], constraint.coeffs[i]))
            end
            if !still_feasible
                # println("i above: ", i)
                return false
            end
        end
        if mins[i] > pre_mins[i]
            if constraint.coeffs[i] > 0
                still_feasible = remove_below!(com, search_space[idx], cld(mins[i], constraint.coeffs[i]))
            else
                still_feasible = remove_above!(com, search_space[idx], cld(mins[i], constraint.coeffs[i]))
            end
            if !still_feasible
                # println("i below: ", i)
                return false
            end
        end
    end

    # if there are only two left check all options
    n_unfixed = 0
    unfixed_ind_1, unfixed_ind_2 = 0, 0
    unfixed_local_ind_1, unfixed_local_ind_2 = 0,0
    unfixed_rhs = constraint.rhs
    li = 0
    for i in indices
        li += 1
        if !isfixed(search_space[i])
            n_unfixed += 1
            if n_unfixed <= 2
                if n_unfixed == 1
                    unfixed_ind_1 = i
                    unfixed_local_ind_1 = li
                else
                    unfixed_ind_2 = i
                    unfixed_local_ind_2 = li
                end
            end
        else
            unfixed_rhs -= value(search_space[i])*constraint.coeffs[li]
        end
    end

    # only a single one left
    if n_unfixed == 1
        if unfixed_rhs % constraint.coeffs[unfixed_local_ind_1] != 0
            com.bt_infeasible[unfixed_ind_1] += 1
            return false
        else
            # divide rhs such that it is comparable with the variable directly without coefficient
            unfixed_rhs = fld(unfixed_rhs, constraint.coeffs[unfixed_local_ind_1])
        end
        if !has(search_space[unfixed_ind_1], unfixed_rhs)
            com.bt_infeasible[unfixed_ind_1] += 1
            return false
        else
            still_feasible = fix!(com, search_space[unfixed_ind_1], unfixed_rhs)
            if !still_feasible
                return false
            end
        end
    elseif n_unfixed == 2
        intersect_cons = intersect(com.subscription[unfixed_ind_1], com.subscription[unfixed_ind_2])
        is_all_different = false
        for constraint_idx in intersect_cons
            if nameof(com.constraints[constraint_idx].fct) == :all_different
                is_all_different = true
                break
            end
        end

        for v in 1:2
            if v == 1
                this, local_this = unfixed_ind_1, unfixed_local_ind_1
                other, local_other = unfixed_ind_2, unfixed_local_ind_2
            else
                other, local_other = unfixed_ind_1, unfixed_local_ind_1
                this, local_this = unfixed_ind_2, unfixed_local_ind_2
            end

            for val in values(search_space[this])
                # if we choose this value but the other wouldn't be an integer => remove this value
                if (unfixed_rhs-val*constraint.coeffs[local_this]) % constraint.coeffs[local_other] != 0
                    if !rm!(com, search_space[this], val)
                        return false
                    end
                    continue
                end

                check_other_val = fld(unfixed_rhs-val*constraint.coeffs[local_this], constraint.coeffs[local_other])
                # if all different but those two are the same
                if is_all_different && check_other_val == val
                    if !rm!(com, search_space[this], val)
                        return false
                    end
                    if has(search_space[other], check_other_val)
                        if !rm!(com, search_space[other], check_other_val)
                            return false
                        end
                    end
                else
                    if !has(search_space[other], check_other_val)
                        if !rm!(com, search_space[this], val)
                            return false
                        end
                    end
                end
            end
        end
    end

    return true
end

function eq_sum(com::CoM, constraint::LinearConstraint, val::Int, index::Int)
    search_space = com.search_space
    csum = 0
    num_not_fixed = 0
    max_extra = 0
    min_extra = 0
    for (i,idx) in enumerate(constraint.indices)
        if idx == index
            val = val*constraint.coeffs[i]
            continue
        end
        if isfixed(search_space[idx])
            csum += value(search_space[idx])*constraint.coeffs[i]
        else
            num_not_fixed += 1
            if constraint.coeffs[i] >= 0
                max_extra += search_space[idx].max*constraint.coeffs[i]
                min_extra += search_space[idx].min*constraint.coeffs[i]
            else
                min_extra += search_space[idx].max*constraint.coeffs[i]
                max_extra += search_space[idx].min*constraint.coeffs[i]
            end
        end
    end
    if num_not_fixed == 0 && csum + val != constraint.rhs
        return false
    end

    if csum + val + min_extra > constraint.rhs
        return false
    end

    if csum + val + max_extra < constraint.rhs
        return false
    end

    return true
end


text:
using JuMP, Cbc, CSV, DataFrames;
# for using the Gurobi solver (license needed)
#using Gurobi
#GUROBI_ENV = Gurobi.Env() # Gurobi academic license message only once

struct HeatPump
    eta::Float32
    rate_max::Float32
end

struct ThermalStorage
    eta::Float32
    volume::Float32
    loss::Float32
    t_supply::Float32
    soc_min::Float32
    soc_max::Float32
end

struct Battery
    eta::Float32
    soc_min::Float32
    soc_max::Float32
    rate_max::Float32
    loss::Float32
end

mutable struct SHEMS
    costfactor::Float32
    p_buy::Float32
    p_sell::Float32
    soc_b::Float32
    soc_fh::Float32
    soc_hw::Float32
    h_start::Int16
end

struct Model_SHEMS
    h_start::Int16
    h_end::Int16
    h_predict::Int16
    h_control::Int16
    big::Int16
    rolling_flag::Bool
    solver::String
    mip_gap::Float32
    output_flag::Bool
    presolve_flag::Int
end


text:
# colorbar_title

PLOTSIZE = 400
ALPHA = 0.4

pointsize(x, sz) = if sz == :auto
    # By default, increase the total when there are many tiny points, else they disappear:
        25 .* sqrt.(weights(x)) .* (size(x,2)^0.15)
    else
    # But if sz is provided, then work from constant total area, divided among points:
        44 .* sqrt.(sz .* x.weights)
    end

using RecipesBase

@recipe function f(x::Weighted{<:Matrix, <:Vector}; sz=:auto, bd=:auto) # exlude both 1D and Tracked
    if size(x,1)==1
        out = Weighted(vec(x.array), x.weights, x.opt) # just to use 1D type signature
    elseif size(x,1)>=4
        out = wPCA(x,2)(x)
    else

        # legend --> :false
        label  --> ""
        if x.opt.clamp
            extr = [x.opt.lo - 0.05, x.opt.hi + 0.05]
            xlims --> extr
            ylims --> extr
            # zlim --> [x.opt.lo, x.opt.hi] # leaves too much space, in GR
        end

        if endswith(x.opt.aname,"PC")
            size   --> (1.5*PLOTSIZE, PLOTSIZE)
            ylims   --> pcaylim(x, 1.5)
        elseif size(x,1)==2
            size   --> (1.1*PLOTSIZE, PLOTSIZE)
        elseif size(x,1)==3
            size   --> (1.3*PLOTSIZE, 1.2*PLOTSIZE)
        end

        seriestype := :scatter
        markeralpha --> ALPHA # nickname alpha=0.2 doesn't overwrite this
        if bd==false || (bd==:auto && size(x,2) > 100) # tiny points are black if they all have edges
            markerstrokewidth --> 0
        end
        markersize := pointsize(x, sz)

        if size(x,1)==2
            xguide --> (x.opt.aname)*"_1"
            yguide --> (x.opt.aname)*"_2"

            out = array(x)[1,:], array(x)[2,:]

        elseif size(x,1)==3
            xguide --> (x.opt.aname)*"_1"
            yguide --> (x.opt.aname)*"_2"
            zguide --> (x.opt.aname)*"_3"

            marker_z --> array(x)[3,:]

            out = array(x)[1,:], array(x)[2,:], array(x)[3,:]
        end
    end
    out
end

@recipe f(x::Weighted{<:Base.ReshapedArray}; sz=1) = copy(x) # these were otherwise not caught
@recipe f(x::Weighted{<:Base.ReinterpretArray}; sz=1) = copy(x)

function pcaylim(x::Weighted, ratio=1.5)
    mat = array(x)

    ex1 = [extrema(mat[1,:])...]
    diff1 = ex1[2] - ex1[1]

    ex2 = [extrema(mat[2,:])...]
    diff2 = ex2[2] - ex2[1]

    if diff1 > ratio * diff2
        ex2 .*= ( diff1 / (ratio * diff2) ) # expand the yilm
    end

    return ex2
end

## 1D with shadow
@recipe function f(x::Weighted{<:Vector}, shadow="yes"; sz=0.5) # sz doesn't get passed to here :(
    size   --> (1.2*PLOTSIZE, 0.8*PLOTSIZE)
    # legend --> :false
    if x.opt.clamp && x.opt.hi < Inf
        xlims --> [x.opt.lo - 0.05, x.opt.hi + 0.05]
    end
    xguide --> x.opt.aname
    # yaxis --> x.opt.wname
    ylims  --> [0, 1.4*maximum(weights(x)) ]
    yticks --> [0, round(maximum(weights(x)),digits=2)]

    if shadow != "only"
        @series begin # plot the points
            seriestype := :scatter
            label --> ""
            markersize := pointsize(x, sz)
            markeralpha --> ALPHA
            # markerstrokewidth --> 0

            array(x), weights(x)
        end
    end
    if shadow != "no"
        @series begin # plot the shadow
            seriestype := :line
            label := ""
            fill  := 0
            seriescolor --> :black
            # fillcolor --> get(palette(plotattributes[:seriescolor]), 0.1)
            seriesalpha := 0.5*ALPHA
            markersize --> 0
            markeralpha --> 0

            shadowxy(array(x), weights(x))
        end
    end
end

function shadowxy(x::Vector, y::Vector, smooth=(max(1, maximum(x)-minimum(x)))/100, numpoints=321; fixmax=true)
    @assert length(x)==length(y)
    xlist = range(minimum(x)-10*smooth, stop=maximum(x)+10*smooth, length=numpoints)
    ylist = zeros(numpoints)
    invtst = 1/(2*smooth^2)
    for i=1:numpoints
        for j=1:length(x)
            ylist[i] += y[j] * exp( -(xlist[i]-x[j])^2 * invtst )
        end
    end
    scale = fixmax ? maximum(y) / maximum(ylist) : 1
    return xlist, scale .* ylist
end

@recipe function f(x::Weighted, fun::Function)
    unique(fun(x))
end
@recipe function f(fun::Function, x::Weighted)
    unique(fun(x))
end

@recipe function f(x::Weighted, λ::Number)
    λ .* x
end
@recipe function f(λ::Number, x::Weighted)
    λ .* x
end


text:
function SOLUTION_is_feasible(sol::Sol)

    debug && DEBUG_correct_solution(sol)

    feasible = true
    dim = 1
    while feasible && dim <= sol.prob.m
        if sol.w[dim] > sol.prob.ω[dim]
            feasible = false
        end
        dim += 1
    end

    return feasible
end

function SOLUTION_is_feasible_after_adding_item(sol::Sol, idx_var::Int)

    debug && DEBUG_correct_solution(sol)

    @assert sol.x[idx_var] == 0 "The item is already in the bag"

    feasible = true
    dim = 1
    while feasible && dim <= sol.prob.m
        if sol.w[dim] + sol.prob.W[dim, idx_var] > sol.prob.ω[dim]
            feasible = false
        end
        dim += 1
    end

    return feasible
end

function SOLUTION_add_item!(sol::Sol, item::Int)

    debug && @assert sol.x[item] == 0 "Item already in the bag"
    debug && DEBUG_correct_solution(sol) # abusive

    sol.x[item] = 1
    sol.z[1] += sol.prob.P[1,item]
    sol.z[2] += sol.prob.P[2,item]
    for dim in 1:sol.prob.m
        sol.w[dim] += sol.prob.W[dim,item]
    end

    debug && DEBUG_correct_solution(sol)

end

function SOLUTION_remove_item!(sol::Sol, item::Int)

    debug && @assert sol.x[item] == 1 "Item not in the bag"
    debug && DEBUG_correct_solution(sol) # abusive

    momkp = sol.prob

    sol.x[item] = 0
    sol.z[1] -= momkp.P[1,item]
    sol.z[2] -= momkp.P[2,item]
    for dim in 1:momkp.m
        sol.w[dim] -= momkp.W[dim,item]
    end

    debug && DEBUG_correct_solution(sol)

end

function SOLUTION_empty_solution(momkp::_MOMKP)

    empty_sol = Sol(momkp, zeros(Int, momkp.n), [0, 0], zeros(Int, momkp.m))

    return empty_sol
end

function SOLUTION_list_of_empty_solutions(nb_solutions::Int, momkp::_MOMKP)

    vec = Vector{Sol}(undef, nb_solutions)
    for i in 1:nb_solutions
        vec[i] = SOLUTION_empty_solution(momkp)
    end

    return vec
end


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.8.0-DEV
#     language: julia
#     name: julia-1.8
# ---

# %%
using BenchmarkTools
b = @benchmark sum(rand()^2 + rand()^2 ≤ 1 for _ in 1:10^6)/10^6

# %%
dump(b)

# %%
methodswith(BenchmarkTools.Trial)

# %%
b.times

# %%

max_stars_repo_path: 0015/BenchmarkTools example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0015/BenchmarkTools example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0015/BenchmarkTools example.jl

max_forks_repo_name: genkuroki/public