text:
using StatsBase
using Random
using Distributions

#Simulation process
#Big function that takes dict of isoforms => [counts] for each gene as input
#TODO

################################################################################
#1.Rank Isoforms. See
#https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1005761#sec011
################################################################################

function FindRank(InputDict)

    arr = []

    #find sums of values in InputDict
    for (key, value) in InputDict
        push!(arr, sum(value))
    end

    return ordinalrank(arr, rev=true)
end

function FindRankbyNoCells(InputDict)

    arr = []

    for (key, value) in InputDict
        push!(arr, sum(value .> 0))
    end

    return ordinalrank(arr, rev = true)
end

function FindStd(InputDict)

    arr = []

    for (key, value) in InputDict
        push!(arr, sum(value .> 0))
    end
    @show arr
    return std(arr)
end

function FindPObs(InputDict)

    arr = []

    for (key, value) in InputDict
        pObs = sum(x -> x>0, value) / length(value)
        push!(arr, pObs)
    end

    return arr
end

function FindHarmonicMean(numIsoforms)

    if numIsoforms ==1
        throw(MethodError("It doesn't make sense to find the Harmonic Mean for
        genes with only 1 isoform!"))
    end

    harmonicMean = 0

    for i in 1:numIsoforms
        harmonicMean += (1/i) * exp(-(1 + i/numIsoforms)^2)
    end

    #@show harmonicMean
    return harmonicMean
end

function FindMedianFrequencies(numIsoforms, harmonicMean)

    arr = Float64[]

    for i in 1:numIsoforms
        top = exp(-(1 + i/numIsoforms)^2)
        bottom = i * harmonicMean
        push!(arr, top/bottom)
    end

    return arr
end

function ScaleFrequencies(frequencies)

    total = sum(frequencies)
    frequencies = frequencies ./ total
end

function FindCumulativeFrequencies(scaledFreqs)

    arr = Float64[]
    scaledFreqs = reverse(sort(scaledFreqs))

    for i in 1:length(scaledFreqs)
        push!(arr, sum(scaledFreqs[1:i]))
    end

    return arr
end

#Function to pick isoforms based on cumulative frequencies and a random number
function PickRank(cumulativeFreqs, random)

    for i in 1:length(cumulativeFreqs)
        if random < cumulativeFreqs[i]
            return i
            break
        end
    end
end

#Function to rescale frequencies after one isoform has been picked
function RescaleFreqs(cumulativeFreqs, rank)

    if rank!=length(cumulativeFreqs)
        for i in (rank + 1):(length(cumulativeFreqs))
            cumulativeFreqs[i] -= cumulativeFreqs[rank]
        end

    end

    deleteat!(cumulativeFreqs, rank)
    cumulativeFreqs = ScaleFrequencies(cumulativeFreqs)
    cumulativeFreqs = FindCumulativeFrequencies(cumulativeFreqs)
end

#Function to rescale rankings after one isoform has been picked
function RescaleRankings(rankings, rank, index)

    deleteat!(rankings, index)

    #small internal rescaling function
    function Rescale(element)
        if element > rank
            return element - 1
        end
        return element
    end

    return Rescale.(rankings)
end

#Wrapper function - only have to do this once for each numIsoforms
function GetIsoformChoiceProbabilities(numIsoforms)

    if numIsoforms > 1
        harmonicMean = FindHarmonicMean(numIsoforms)
        frequencies = FindMedianFrequencies(numIsoforms, harmonicMean)
        frequencies = ScaleFrequencies(frequencies)
        #@show frequencies
        #frequencies = FindCumulativeFrequencies(frequencies)
    else
        frequencies = [1.0]
    end
    return frequencies
end

#Wrapper function - only have to do this once for each numIsoforms
function GetIsoformChoiceProbabilitiesGene(numIsoforms)

    if numIsoforms > 1
        harmonicMean = FindHarmonicMean(numIsoforms)
        frequencies = FindMedianFrequencies(numIsoforms, harmonicMean)
        frequencies = ScaleFrequencies(frequencies)
        frequencies = FindCumulativeFrequencies(frequencies)
    else
        frequencies = [1.0]
    end
    return frequencies
end

function RandomIsoformChoiceProbabilities(numIsoforms, numGenes)

    rand_probs = rand(numGenes, numIsoforms)
    output = Array{Float64}(undef, numGenes, numIsoforms)

    for i in 1:numGenes
        rollingSum = sum(rand_probs[i, :])

        for j in 1:numIsoforms
            output[i,j] = rand_probs[i,j] / rollingSum
        end
    end
    return output
end

################################################################################
#2.Find P(dropout) for each isoform
################################################################################
function FindPDropout(InputDict, Km)

    arr = []

    #find PDropout for each value in InputDict
    for (key, value) in InputDict
        S = mean(value)
        PDropout = 1 - (S / (Km + S))
        push!(arr, PDropout)
    end

    return arr
end

function FindMeanExpr(InputDict)
    arr = []

    for (key, value) in InputDict
        S = mean(value)
        push!(arr, S)
    end

    return arr
end

#Introduce quantification errors for expressed isoforms using global pFN and pFP
function QuantErrorsForExpressed(BooleanValue::Bool, randomNumber::Float64)

    if BooleanValue == true
        if randomNumber < pFN
            return false
        end
    else
        if randomNumber < pFP
            return true
        end
    end

    return BooleanValue
end

function RandomQuantErrorsForExpressed(BooleanValue::Bool, randomNumber::Float64)

    if BooleanValue == true
        if randomNumber < rand(Uniform(0,0.5))
            return false
        end
    else
        if randomNumber < rand(Uniform(0,0.5))
            return true
        end
    end

    return BooleanValue
end

function RandomQuantErrorsForUnexpressed(randomNumber::Float64)

    if randomNumber < rand(Uniform(0,0.5))
        return true
    end
    return false
end

################################################################################
# Simulation function
################################################################################

function Simulate(cumulativeFreqs::Array{Float64, 1}, rankings::Array{Int64, 1},
    pDropouts::Array{Float64, 1}, numCells::Int64, IsoformsPerGenePerCell::Int64,
    NumIsoformsForGene::Int64, test::Bool)

    #Initialise some temporary variables
    tmpFreqs = deepcopy(cumulativeFreqs)
    tmpRankings = deepcopy(rankings)
    tmpDropouts = deepcopy(pDropouts)
    staticRankings = deepcopy(rankings)
    ObservedIsoformsPerCell = Float64[]

    #Use a seed when testing so results are reproducible and thus testable
    if test == true
        Random.seed!(1234)
    end

    #Iterate over cells
    for i in 1:numCells

        #Re-initialise some variables
        ObservedExpression = fill(false, NumIsoformsForGene)

        #Necessary to reset temporary variables if IsoformsPerGenePerCell > 1
        if tmpFreqs != cumulativeFreqs
            tmpFreqs = deepcopy(cumulativeFreqs)
            tmpRankings = deepcopy(rankings)
            tmpDropouts = deepcopy(pDropouts)
            staticRankings = deepcopy(rankings)
        end

        #Iterate over the no. expressed isoforms per cell
        for isoform in 1:IsoformsPerGenePerCell

            #Pick isoforms
            choice = PickRank(tmpFreqs, rand())

            #Make two indexes. index is for use with rescaled (and shrinking)
            #temporary variable arrays used in this loop
            tmp = findall(x -> x==choice, tmpRankings)

            index = findall(x -> x==choice, tmpRankings)[1]

            #ObsExIndex is used for indexing into ObservedExpression and needs
            #to not shrink as the temporary variable arrays shrink
            ObsExIndex = staticRankings[index]

            #See whether random number < Pdropout
            if rand() > tmpDropouts[index] #Note if this is not true, ObservedExpression[index] stays as false, which is the default
                ObservedExpression[ObsExIndex] = true
            end

            #Rescaling necessary before next iteration
            if isoform != IsoformsPerGenePerCell
                tmpFreqs = RescaleFreqs(tmpFreqs, choice)
                tmpRankings = RescaleRankings(tmpRankings, choice, index)
                tmpDropouts = deleteat!(tmpDropouts, index)
                staticRankings = deleteat!(staticRankings, index)
            end

        end

        #Simulate quant errors
        for obs in 1:length(ObservedExpression)
            if ObservedExpression[obs] == false
                if rand() < pFP
                     ObservedExpression[obs] = true
                end
            else
                if rand() < pFN
                    ObservedExpression[obs] = false
                end
            end
        end

        #Save the observed num isoforms per cell
        push!(ObservedIsoformsPerCell, sum(ObservedExpression))
    end

    #Return the mean
    return(mean(ObservedIsoformsPerCell))
end

#For GPU version should use curand from CuArrays.CURAND to make random number array before starting
#Or maybe implement a GPU version of rand

################################################################################
#globalArraySimulation() helper functions

#Function that returns 4D array of pDropouts corresponding to isoform choice
function pickRandIsoforms( IsoChoiceProbabilities, NumSimulations,
    genesList, pDropoutDict, numCells, NumIsoformsToSimulate)

    numGenes = length(genesList) #convenient binding
    IsoChoiceArr = Array{Float64}(undef, NumSimulations, numGenes,numCells,NumIsoformsToSimulate)

    for sim in 1:NumSimulations

        for i in 1:length(genesList)

            #get dropout info
            dropoutArr::Array{Float64, 1} = pDropoutDict[genesList[i]]

            #for cell in cells
            for j in 1:numCells

                choices = sample(collect(1:length(IsoChoiceProbabilities[i, :])),
                ProbabilityWeights(IsoChoiceProbabilities[i, :]),
                NumIsoformsToSimulate, replace = false)

                for k in 1:length(choices)
                    IsoChoiceArr[sim,i,j,k] = dropoutArr[choices[k]]
                end
            end
        end
    end
    return IsoChoiceArr
end

function FindMeanUnifObsProbabilities(numIsoforms, genesList, pObsDict,
    pDropoutDict, rankingDict)

    numGenes = length(genesList)

    df = DataFrame(Rank=Int64[], pDropout = Float64[], pChoice = Float64[])

    for i in 1:numGenes

        #get dropout info and get pObs info
        dropoutArr::Array{Float64, 1} = pDropoutDict[genesList[i]]
        pObsArr::Array{Float64} = pObsDict[genesList[i]]

        #Get pChoices
        pChoiceArr::Array{Float64} = []
        for iso in 1:length(pObsArr)
            pChoice = (pObsArr[iso] - pFP)/ (((1-dropoutArr[iso]) * (1 - pFN)) - pFP)

            #if pChoice >1, set to 1 as pChoice is a probability
            #I think this sometimes happens for stochastic reasons or
            #because we have estimated parameters poorly.
            if pChoice > 1
                pChoice = 1.0
           end
            push!(pChoiceArr, pChoice)
        end

        newRow = DataFrame(Rank = rankingDict[genesList[i]],
        pDropout = dropoutArr,
        pChoice = pChoiceArr)

        append!(df, newRow)

    end

    return df
end

function PickUnifObsIsoforms(numIsoforms, genesList, pObsDict,
     pDropoutDict, numCells, NumSimulations)

     numGenes = length(genesList)
     output = Array{Float64}(undef, NumSimulations, numGenes,numCells,numIsoforms)

     for sim in 1:NumSimulations

         for i in 1:numGenes

             #get dropout info and get pObs info
             dropoutArr::Array{Float64, 1} = pDropoutDict[genesList[i]]
             pObsArr::Array{Float64} = pObsDict[genesList[i]]

             #Get pChoices
             pChoiceArr::Array{Float64} = []
             for iso in 1:length(pObsArr)
                 pChoice = abs((pObsArr[iso] - pFP))/ abs((((1-dropoutArr[iso]) * (1 - pFN)) - pFP))

                 #if pChoice >1, set to 1 as pChoice is a probability
                 #I think this sometimes happens for stochastic reasons or
                 #because we have estimated parameters poorly.
                 if pChoice > 1
                     pChoice = 1.0
                end
                 push!(pChoiceArr, pChoice)
             end

             for j in 1:numCells

                 choices = sample(collect(1:length(pChoiceArr)),
                 ProbabilityWeights(pChoiceArr),
                 numIsoforms, replace = false)

                 for k in 1:length(choices)
                     output[sim,i,j,k] = dropoutArr[choices[k]]
                end
            end
        end
    end
    return output
end

function PickCellVarIsoforms(numIsoforms, genesList, pObsDict,
     pDropoutDict, numCells, NumSimulations)

     numGenes = length(genesList)
     output = Array{Float64}(undef, NumSimulations, numGenes,numCells,numIsoforms)
     sigma = 0.002

     for sim in 1:NumSimulations

         for i in 1:numGenes

             #get dropout info and get pObs info
             dropoutArr::Array{Float64, 1} = pDropoutDict[genesList[i]]
             pObsArr::Array{Float64} = pObsDict[genesList[i]]

             #Get pChoices
             pChoiceArr::Array{Float64} = []
             for iso in 1:length(pObsArr)
                 #I wonder if this will ever be negative? Might need to find absolute value
                 pChoice = abs((pObsArr[iso] - pFP))/ abs((((1-dropoutArr[iso]) * (1 - pFN)) - pFP))

                 #if pChoice >1, set to 1 as pChoice is a probability
                 #I think this sometimes happens for stochastic reasons or
                 #because we have estimated parameters poorly.
                 if pChoice > 1
                     pChoice = 1.0
                end
                 push!(pChoiceArr, pChoice)
             end

             for j in 1:numCells

                 cellpChoiceArr::Array{Float64} = []

                 for iso in 1:length(pChoiceArr)

                     mu = pChoiceArr[iso]

                     #This condition ensures alpha and beta are positive.
                     #If it is not fulfilled the beta distrubution is not
                     #appropriate, so we just use the original probability.
                     if mu > 0 && mu < 0.99999

                         #find alpha and beta
                         #@show mu
                         alpha = (((1 - mu)/(sigma^2)) - (1/mu)) * mu^2
                         beta = alpha * ((1/mu) - 1)
                         #@show alpha
                         #@show beta

                         #find p(choice) for this cell by pulling from beta distribution
                         CellChoice = rand(Beta(alpha, beta), 1)
                         push!(cellpChoiceArr, CellChoice[1])
                    else
                        push!(cellpChoiceArr, mu)
                    end
                end

                #@show pObsArr
                #@show pChoiceArr
                #@show cellpChoiceArr
                #@show numIsoforms

                 choices = sample(collect(1:length(cellpChoiceArr)),
                 ProbabilityWeights(cellpChoiceArr),
                 numIsoforms, replace = false)

                 for k in 1:length(choices)
                     output[sim,i,j,k] = dropoutArr[choices[k]]
                end
            end
        end
    end
    return output
end

function PickEqualProbabilities(numIsoforms, genesList, pObsDict,
     pDropoutDict, numCells, NumSimulations, nameDistribution)

     numGenes = length(genesList)
     output = Array{Float64}(undef, NumSimulations, numGenes,numCells,numIsoforms)
     sigma = 0.002

     for sim in 1:NumSimulations

         for i in 1:numGenes

             #get dropout info and get pObs info
             dropoutArr::Array{Float64, 1} = pDropoutDict[genesList[i]]

             #Get pChoices
             pChoiceArr::Array{Float64} = [0.25,0.25,0.25,0.25] #Equal probabilities

             for j in 1:numCells

                 cellpChoiceArr::Array{Float64} = []

                 for iso in 1:length(pChoiceArr)

                     mu = pChoiceArr[iso]

                     #This condition ensures alpha and beta are positive.
                     #If it is not fulfilled the beta distrubution is not
                     #appropriate, so we just use the original probability.
                     if nameDistribution == "Bernoulli"

                         #find p(choice) for this cell by pulling from beta distribution
                         CellChoice = rand(Bernoulli(0.25))
                         push!(cellpChoiceArr, Float64(CellChoice))
                    end

                    if nameDistribution == "Normal"
                        push!(cellpChoiceArr, rand(TruncatedNormal(mu, 0.06, 0.0, 1.0))) #sd chosen to get a large range but >0
                    end
                    if nameDistribution == "Uniform"
                        push!(cellpChoiceArr, mu)
                    end
                end

                #@show pObsArr
                #@show pChoiceArr
                #@show cellpChoiceArr
                #@show numIsoforms

                choices::Array{Int64} = []

                if nameDistribution == "Bernoulli" && sum(cellpChoiceArr) == 1
                    choices = findall(x -> x==1, cellpChoiceArr)
                elseif nameDistribution == "Bernoulli" && sum(cellpChoiceArr) == 0
                    cellpChoiceArr = [0.25,0.25,0.25,0.25]  #collapse back to uniform
                    choices = sample(collect(1:length(cellpChoiceArr)),
                    ProbabilityWeights(cellpChoiceArr),
                    numIsoforms, replace = false)
                else
                    choices = sample(collect(1:length(cellpChoiceArr)),
                    ProbabilityWeights(cellpChoiceArr),
                    numIsoforms, replace = false)
                end

                 for k in 1:length(choices)
                     output[sim,i,j,k] = dropoutArr[choices[k]]
                end
            end
        end
    end
    return output
end

#Function that returns 4D array of pDropouts corresponding to isoform choice
function pickIsoforms(IsoChoiceProbabilities, NumSimulations,
    genesList, rankingDict, pDropoutDict, numCells,
    NumIsoformsToSimulate)

    numGenes = length(genesList) #convenient binding
    IsoChoiceArr = Array{Float64}(undef, NumSimulations, numGenes,numCells,NumIsoformsToSimulate)

    #for gene in genes
    for sim in 1:NumSimulations

        for i in 1:length(genesList)

            #get ranking info
            rankArr::Array{Int64, 1} = rankingDict[genesList[i]]
            dropoutArr::Array{Float64, 1} = pDropoutDict[genesList[i]]

            #@show rankArr
            #@show dropoutArr

            #for cell in cells
            for j in 1:numCells

                #Sample NumIsoformsToSimulate from IsoChoiceProbabilities
                #@show IsoChoiceProbabilities
                ranks = sample(collect(1:length(IsoChoiceProbabilities)),
                ProbabilityWeights(IsoChoiceProbabilities),
                NumIsoformsToSimulate, replace = false)

                #Find the index of the isoform at the chosen ranks, then write
                #p(dropout) to the appropriate location in the array
                for k in 1:length(ranks)
                    index = findall(x -> x==ranks[k], rankArr)[1]
                    IsoChoiceArr[sim,i,j,k] = dropoutArr[index]
                end
            end
            #@show index
        end
    end

    return IsoChoiceArr
end

function pickIsoformsReportRank(IsoChoiceProbabilities, NumSimulations,
    genesList, rankingDict, pDropoutDict, numCells,
    NumIsoformsToSimulate)

    numGenes = length(genesList) #convenient binding
    IsoChoiceArr = Array{Float64}(undef, NumSimulations, numGenes,numCells,NumIsoformsToSimulate)
    RecordedRankArr = IsoChoiceArr = Array{Float64}(undef, NumSimulations, numGenes,numCells,NumIsoformsToSimulate)

    #for gene in genes
    for sim in 1:NumSimulations

        for i in 1:length(genesList)

            #get ranking info
            rankArr::Array{Int64, 1} = rankingDict[genesList[i]]
            dropoutArr::Array{Float64, 1} = pDropoutDict[genesList[i]]

            #@show rankArr
            #@show dropoutArr

            #for cell in cells
            for j in 1:numCells

                #Sample NumIsoformsToSimulate from IsoChoiceProbabilities
                #@show IsoChoiceProbabilities
                ranks = sample(collect(1:length(IsoChoiceProbabilities)),
                ProbabilityWeights(IsoChoiceProbabilities),
                NumIsoformsToSimulate, replace = false)

                #Find the index of the isoform at the chosen ranks, then write
                #p(dropout) to the appropriate location in the array
                for k in 1:length(ranks)
                    index = findall(x -> x==ranks[k], rankArr)[1]
                    IsoChoiceArr[sim,i,j,k] = dropoutArr[index]
                    RecordedRankArr[sim, i, j, k] = ranks[k]
                end
            end
            #@show index
        end
    end

    return (IsoChoiceArr, RecordedRankArr)
end

#Make a modified version of this for IsoformDependence model
function pickUnexpressedErrors(MaxUnexpressedDetectedIsoforms, NumSimulations,
    genesList, numCells, unexpr_dict, NumIsoformsToSimulate, filteringThreshold)

    numGenes = length(genesList) #convenient binding
    UnexpressedDetectedIsoforms = Array{Int64}(undef, NumSimulations, numGenes,
    numCells)

    #Extract correct number of isoforms per gene
    for sim in 1:NumSimulations

        for i in 1:length(genesList)
            numUnexpr = unexpr_dict[genesList[i]] +
            (filteringThreshold - NumIsoformsToSimulate)

            for j in 1:numCells

                numDetected=0

                for k in 1:numUnexpr
                    numDetected += MaxUnexpressedDetectedIsoforms[sim,i,j,k]
                end

                UnexpressedDetectedIsoforms[sim,i,j] = numDetected
            end
        end
    end

    return UnexpressedDetectedIsoforms
end

#pickUnexpressedErrors for IsoformDependence model
function UnexpressedErrorsIsoformDependence(unexprRandArr, NumSimulations,
    genesList, numCells, unexpr_dict, NumIsoformsToSimulate, filteringThreshold)

    numGenes = length(genesList) #convenient binding
    UnexpressedDetectedIsoforms = Array{Int64}(undef, NumSimulations, numGenes,
    numCells)

    for sim in 1:NumSimulations

        for i in 1:length(genesList)

            #set error rate
            numUnexpr = unexpr_dict[genesList[i]] +
            (filteringThreshold - NumIsoformsToSimulate)
            pFP = ((numUnexpr + NumIsoformsToSimulate) * 0.01)

            for j in 1:numCells

                numDetected = 0

                for k in 1:numUnexpr
                    isDetected = unexprRandArr[sim,i,j,k] < pFP
                    numDetected += isDetected
                end

                UnexpressedDetectedIsoforms[sim,i,j] = numDetected
            end
        end
    end

    return UnexpressedDetectedIsoforms
end

function ExpressedErrorsIsoformDependence(readCapturedIsoforms, randIsoChoiceArr,
    NumSimulations, genesList, numCells, unexpr_dict, NumIsoformsToSimulate,
    filteringThreshold)

    numGenes = length(genesList) #convenient binding
    expressedDetectedIsoforms = Array{Bool}(undef, NumSimulations, numGenes,
    numCells, NumIsoformsToSimulate)
    pFN = 0.04

    for sim in 1:NumSimulations

        for i in 1:numGenes

            #set error rate
            numUnexpr = unexpr_dict[genesList[i]] +
            (filteringThreshold - NumIsoformsToSimulate)
            pFP = ((numUnexpr + NumIsoformsToSimulate) * 0.01)

            for j in 1:numCells
                for k in 1:NumIsoformsToSimulate
                    if readCapturedIsoforms[sim,i,j,k] == true
                        if randIsoChoiceArr[sim,i,j,k] < pFN
                            expressedDetectedIsoforms[sim,i,j,k] = false
                        else
                            expressedDetectedIsoforms[sim,i,j,k] = true
                        end
                    else
                        if randIsoChoiceArr[sim,i,j,k] < pFP
                            expressedDetectedIsoforms[sim,i,j,k] = true
                        else
                            expressedDetectedIsoforms[sim,i,j,k] = false
                        end
                    end
                end
            end
        end
    end
    return expressedDetectedIsoforms
end

#Negative control simulations
function negControlSimulation(NumSimulations, genesList, rankingDict, pObsDict,
    numCells, NumIsoformsToSimulate, filteringThreshold, pDropoutDict,
    MaxNumIsos)

    #make a 3D array of dims numGenes by numCells by NumIsoforms and fill with
    #random numbers
    numGenes = length(genesList) #convenient binding
    randIsoChoiceArr = rand(NumSimulations, numGenes,numCells,NumIsoformsToSimulate)

    #Do Weibull model of isoform choice
    local WeibullIsoChoiceProbabilities = Array{Float64}(undef, NumIsoformsToSimulate)
    WeibullIsoChoiceProbabilities = GetIsoformChoiceProbabilities(filteringThreshold)

    #Pick isoforms
    IsoChoiceArr = pickIsoforms(WeibullIsoChoiceProbabilities,
    NumSimulations, genesList, rankingDict, pDropoutDict, numCells,
    NumIsoformsToSimulate)

    #Do a more than comparison to eliminate dropouts (should be none as pdropout=0)
    #This line returns false when the random number is less than the probability of dropout
    readCapturedIsoforms::Array{Bool} = randIsoChoiceArr .>= IsoChoiceArr

    #Find sum of detected isoforms
    totalExpressedDetected = reshape(sum(readCapturedIsoforms, dims = 4),
    NumSimulations, numGenes, numCells)

    #Find mean no. isoforms detected per gene per cell
    MeanIsoformsPerGenePerCell = reshape(mean(totalExpressedDetected, dims = 3),
    NumSimulations * numGenes)

    return MeanIsoformsPerGenePerCell

end

function FindOverlap(totalExpressedDetected, NumIsoformsToSimulate,
    NumSimulations, numGenes)

    overlap = totalExpressedDetected ./ NumIsoformsToSimulate
    meanOverlap = reshape(mean(overlap, dims = 3),
    NumSimulations * numGenes)

    return meanOverlap
end

#Global simulation function using arrays (hopefully faster)
function globalArraySimulation(NumSimulations, genesList, rankingDict, pObsDict,
    pDropoutDict, numCells, NumIsoformsToSimulate, filteringThreshold,
    unexpr_dict, MaxNumIsos, error_model, choice_model)

    #@show choice_model

    #make a 3D array of dims numGenes by numCells by NumIsoforms and fill with
    #random numbers
    numGenes = length(genesList) #convenient binding
    randIsoChoiceArr = rand(NumSimulations, numGenes,numCells,NumIsoformsToSimulate)
    #@show size(randIsoChoiceArr)

    local WeibullIsoChoiceProbabilities = Array{Float64}(undef, NumIsoformsToSimulate)
    local IsoChoiceProbabilities = Array{Float64}(undef, numGenes, NumIsoformsToSimulate)
    local IsoChoiceArr = Array{Float64}(undef, NumSimulations, numGenes,
    numCells, NumIsoformsToSimulate)

    if choice_model == "Weibull"
        #Get general isoform choice probabilities
        WeibullIsoChoiceProbabilities = GetIsoformChoiceProbabilities(filteringThreshold)

        #Pick isoforms
        IsoChoiceArr = pickIsoforms(WeibullIsoChoiceProbabilities,
        NumSimulations, genesList, rankingDict, pDropoutDict, numCells,
        NumIsoformsToSimulate)

    elseif choice_model =="Random"

        #Get probabilities and pick isoforms
        IsoChoiceProbabilities = RandomIsoformChoiceProbabilities(
        filteringThreshold, length(genesList))

        IsoChoiceArr = pickRandIsoforms(IsoChoiceProbabilities,
        NumSimulations, genesList, pDropoutDict, numCells, NumIsoformsToSimulate)

    elseif choice_model == "UniformObserved"

        IsoChoiceArr = PickUnifObsIsoforms(NumIsoformsToSimulate, genesList,
        pObsDict, pDropoutDict, numCells, NumSimulations)

    elseif choice_model == "CellVariable"

        IsoChoiceArr = PickCellVarIsoforms(NumIsoformsToSimulate, genesList,
        pObsDict, pDropoutDict, numCells, NumSimulations)

    elseif choice_model == "Uniform"

        IsoChoiceArr = PickEqualProbabilities(NumIsoformsToSimulate, genesList,
        pObsDict, pDropoutDict, numCells, NumSimulations, "Uniform")

    elseif choice_model == "Normal"

        IsoChoiceArr = PickEqualProbabilities(NumIsoformsToSimulate, genesList,
        pObsDict, pDropoutDict, numCells, NumSimulations, "Normal")

    elseif choice_model == "Bernoulli"

        IsoChoiceArr = PickEqualProbabilities(NumIsoformsToSimulate, genesList,
        pObsDict, pDropoutDict, numCells, NumSimulations, "Bernoulli")

    else
        throw(ArgumentError("choice_model not recognised."))
    end

    #Do a more than comparison to eliminate dropouts
    #This line returns false when the random number is less than the probability of dropout
    readCapturedIsoforms::Array{Bool} = randIsoChoiceArr .>= IsoChoiceArr
    pDropoutArr = reshape(IsoChoiceArr, NumSimulations * numGenes *
    numCells * NumIsoformsToSimulate)

    #Fill with new random numbers
    randIsoChoiceArr = rand(NumSimulations, numGenes,numCells,NumIsoformsToSimulate)

    local expressedDetectedIsoforms::Array{Bool}

    #Quant errors for expressed isoforms
    if error_model == "none"
        expressedDetectedIsoforms = map(QuantErrorsForExpressed,
        readCapturedIsoforms, randIsoChoiceArr)
    elseif error_model == "random"
        expressedDetectedIsoforms = map(RandomQuantErrorsForExpressed,
        readCapturedIsoforms, randIsoChoiceArr)
    elseif error_model == "IsoformDependence"
        expressedDetectedIsoforms = ExpressedErrorsIsoformDependence(
        readCapturedIsoforms, randIsoChoiceArr, NumSimulations, genesList,
        numCells, unexpr_dict, NumIsoformsToSimulate, filteringThreshold)
    else
        throw(ArgumentError("Unrecognised error_model"))
    end

    #Find number of expressed detected isoforms per gene per cell
    totalExpressedDetected = reshape(sum(expressedDetectedIsoforms, dims = 4),
    NumSimulations, numGenes, numCells)

    overlapArr = FindOverlap(totalExpressedDetected, NumIsoformsToSimulate,
    NumSimulations, numGenes)

    #Now we need to do quant errors for unexpressed isoforms
    #First, make the biggest array of random numbers we could possibly need
    unexprRandArr = rand(NumSimulations, numGenes,numCells,
    MaxNumIsos - NumIsoformsToSimulate) #Do something about hard coded numbers

    local MaxUnexpressedDetectedIsoforms::Array{Bool}

    #Apply quant errors
    if error_model == "none"
        MaxUnexpressedDetectedIsoforms = map(x -> x < pFP,
        unexprRandArr)
        UnexpressedDetectedIsoforms = pickUnexpressedErrors(
        MaxUnexpressedDetectedIsoforms, NumSimulations, genesList, numCells,
        unexpr_dict, NumIsoformsToSimulate, filteringThreshold)
    elseif error_model == "random"
        MaxUnexpressedDetectedIsoforms = map(
        RandomQuantErrorsForUnexpressed, unexprRandArr)
        UnexpressedDetectedIsoforms = pickUnexpressedErrors(
        MaxUnexpressedDetectedIsoforms, NumSimulations, genesList, numCells,
        unexpr_dict, NumIsoformsToSimulate, filteringThreshold)
    elseif error_model == "IsoformDependence"
        UnexpressedDetectedIsoforms = UnexpressedErrorsIsoformDependence(
        unexprRandArr, NumSimulations, genesList, numCells, unexpr_dict,
        NumIsoformsToSimulate, filteringThreshold)
    else
        throw(ArgumentError("Unrecognised error_model"))
    end

    #Find total num detected isoforms and convert to array of means
    DetectedIsoforms = map(+, totalExpressedDetected, UnexpressedDetectedIsoforms)

    MeanIsoformsPerGenePerCell = reshape(mean(DetectedIsoforms, dims = 3),
    NumSimulations * numGenes)

    return (MeanIsoformsPerGenePerCell, overlapArr, pDropoutArr)
end

#############################################
# reviewer 2 ranking
#############################################

function countRanks(detectedIsos, rankArr,
    numGenes, numCells, NumIsoformsToSimulate)

    spearmans_rho_arr = []

    sim = 1 # gonna assume 1 simulation
    for gene in 1:numGenes

        ones = 0
        twos = 0
        threes = 0
        fours = 0

        for cell in 1:numCells
            for iso in 1:NumIsoformsToSimulate
                if rankArr[sim,gene,cell,iso] == 1 &&
                    detectedIsos[sim,gene,cell,iso] == true

                    ones += 1

                elseif rankArr[sim,gene,cell,iso] == 2 &&
                    detectedIsos[sim,gene,cell,iso] == true

                    twos += 1

                elseif rankArr[sim,gene,cell,iso] == 3 &&
                    detectedIsos[sim,gene,cell,iso] == true

                    threes += 1

                elseif rankArr[sim,gene,cell,iso] == 4 &&
                    detectedIsos[sim,gene,cell,iso] == true

                    fours += 1
                else
                    continue
                end
            end
        end

        if ones == twos == threes == fours
            continue
        end

        spearmans_rho = corspearman([1,2,3,4], [ones, twos, threes, fours])
        push!(spearmans_rho_arr, spearmans_rho)
    end

    return spearmans_rho_arr
end

#Global simulation function using arrays (hopefully faster)
function rankingSimulation(NumSimulations, genesList, rankingDict, pObsDict,
    pDropoutDict, numCells, NumIsoformsToSimulate, filteringThreshold,
    unexpr_dict, MaxNumIsos, error_model, choice_model)

    #@show choice_model

    #make a 3D array of dims numGenes by numCells by NumIsoforms and fill with
    #random numbers
    numGenes = length(genesList) #convenient binding
    randIsoChoiceArr = rand(NumSimulations, numGenes,numCells,NumIsoformsToSimulate)
    #@show size(randIsoChoiceArr)

    local WeibullIsoChoiceProbabilities = Array{Float64}(undef, NumIsoformsToSimulate)
    local IsoChoiceProbabilities = Array{Float64}(undef, numGenes, NumIsoformsToSimulate)
    local IsoChoiceArr = Array{Float64}(undef, NumSimulations, numGenes,
    numCells, NumIsoformsToSimulate)

    if choice_model == "Weibull"
        #Get general isoform choice probabilities
        WeibullIsoChoiceProbabilities = GetIsoformChoiceProbabilities(filteringThreshold)

        #Pick isoforms
        out = pickIsoformsReportRank(WeibullIsoChoiceProbabilities,
        NumSimulations, genesList, rankingDict, pDropoutDict, numCells,
        NumIsoformsToSimulate)
        IsoChoiceArr = out[1]
        #@show IsoChoiceArr
        RecordedRankArr = out[2]
        #@show RecordedRankArr
    else
        throw(ArgumentError("choice_model not recognised."))
    end

    #Do a more than comparison to eliminate dropouts
    #This line returns false when the random number is less than the probability of dropout
    readCapturedIsoforms::Array{Bool} = randIsoChoiceArr .>= IsoChoiceArr
    pDropoutArr = reshape(IsoChoiceArr, NumSimulations * numGenes *
    numCells * NumIsoformsToSimulate)

    #Fill with new random numbers
    randIsoChoiceArr = rand(NumSimulations, numGenes,numCells,NumIsoformsToSimulate)

    local expressedDetectedIsoforms::Array{Bool}

    #Quant errors for expressed isoforms
    if error_model == "none"
        expressedDetectedIsoforms = map(QuantErrorsForExpressed,
        readCapturedIsoforms, randIsoChoiceArr)
    elseif error_model == "random"
        expressedDetectedIsoforms = map(RandomQuantErrorsForExpressed,
        readCapturedIsoforms, randIsoChoiceArr)
    elseif error_model == "IsoformDependence"
        expressedDetectedIsoforms = ExpressedErrorsIsoformDependence(
        readCapturedIsoforms, randIsoChoiceArr, NumSimulations, genesList,
        numCells, unexpr_dict, NumIsoformsToSimulate, filteringThreshold)
    else
        throw(ArgumentError("Unrecognised error_model"))
    end

    #Find number of expressed detected isoforms per gene per cell
    #noCellsExpressing = reshape(sum(expressedDetectedIsoforms, dims = 3),
    #numGenes, NumIsoformsToSimulate) #NumSimulations must always be 1

    spearmans_rho = countRanks(expressedDetectedIsoforms, RecordedRankArr,
        numGenes, numCells, NumIsoformsToSimulate)

    return spearmans_rho

end

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
# the base function searchsorted does exactly what the exercise asks for
const binarysearch = searchsorted

ext: jl


text:
########## Parameters for scripts ##########

######## moving window estimators ########

###### bspline-scalogram ######
cscl = Int(Npd*sfct)  # central scale
wscl = 10  # half scale window width, e.g. <=10
minscl, maxscl = max(1, cscl-wscl), cscl+wscl
# maxscl = Npd - 0*Npd÷4
# # minscl = maxscl ÷ 2
# minscl = maxscl - Npd÷2

wfct = 40  # factor of moving window

parms = (
    estimator = "bspline_scalogram",
    method = :optim,  # method of optimizer

    # Parameters for B-Spline wavelet transform
    scalefactor = sfct,
    sclrng = 2*(minscl:maxscl),  # index range of inertial scales
    vm = 2,  # vanishing moment of the wavelet
    pows = 1.0:0.1:2.0, # multifractal moments
    # pows = [2.0],
    reweight = false,  # reweighted regression, warning: heavy computation if turned on
    fmode = :causal,  # mode of filtration

    # Parameters for moving window
    windowfactor = wfct,
    wsize = cscl*wfct,  # size of moving window
    pov = Npd÷1,  # step of moving window
    boundary = :soft,  # boundary condition
)

# multi-fractal estimator
mono_estimator = (x,p) -> FracFin.bspline_scalogram_estim(x, parms.sclrng, parms.vm; pow=p, method=parms.method, reweight=parms.reweight)
base_estimator = x -> FracFin.multifractal_estim(mono_estimator, x, parms.pows)

# mono-fractal estimator
# base_estimator = x -> FracFin.bspline_scalogram_estim(x, parms.sclrng, parms.vm; pow=parms.pow, method=parms.method, reweight=parms.reweight)

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
#
# Halley's method
#
# Reference:
# Y. Nakatsukasa, Z. Bai and F. Gygi, Optimizing Halley's iteration
# for computing the matrix polar decomposition, SIAM, J. Mat. Anal.
# Vol. 31, Num 5 (2010) pp. 2700-2720
#

type HalleyAlg{T} <: PolarAlg
    maxiter::Int
    verbose::Bool
    tol::T

    function HalleyAlg(;maxiter::Integer=100,
                       verbose::Bool=false,
                       tol::Real = cbrt(eps(T)))
        maxiter > 1 || error("maxiter must be greater than 1.")
        tol > 0 || error("tol must be positive.")

        new(maxiter,
            verbose,
            tol)
    end
end

function solve!{T}(alg::HalleyAlg,
                X::Matrix{T}, U::Matrix{T}, H::Matrix{T})
    common_iter!(HalleyUpdater(), X, U, H, alg.maxiter, alg.verbose, alg.tol)
end

immutable HalleyUpdater <: PolarUpdater end

function update_U!{T}(upd::HalleyUpdater, U::Matrix{T})
    UtU = Array(T, size(U))
    At_mul_B!(UtU, U, U)
    copy!(U, U * (3*I + UtU)* inv(I + 3*UtU))
end

#
# QR-based Dynamically Weighted Halley (QDWH) algorithm
#
# Reference: Optimizing Halley's iteration for computing the matrix
#            polar decomposition, Yuji Nakatsukasa, Zhaojun Bai and
#            Francois Gygi, SIAM, J. Mat. Anal. Vol. 31, Num 5 (2010)
#            pp. 2700-2720
#
# Limitations: 1. the QDWH method should support m > n matrix itself in the future.
#              2. the computing of alpha and L in solve! can be improved by using
#                 norm and condition number estimate, which are currently
#                 not available in Julia.
#
type QDWHAlg{T} <: PolarAlg
    maxiter::Int
    verbose::Bool
    piv::Bool       # whether to pivot
    tol::T

    function QDWHAlg(;maxiter::Integer=100,
                     verbose::Bool=false,
                     piv::Bool=true,
                     tol::Real=cbrt(eps(T)))
        maxiter > 1 || error("maxiter must be greater than 1.")
        tol > 0 || error("tol must be positive.")

        new(maxiter,
            verbose,
            piv,
            tol)
    end
end

function solve!{T}(alg::QDWHAlg,
                X::Matrix{T}, U::Matrix{T}, H::Matrix{T})
    # alpha is an estimate of the largest singular value of the
    # original matrix
    X_temp = Array(T, size(X))
    copy!(X_temp, X)

    n = size(X_temp, 1)
    alpha = vecnorm(X_temp)
    for i in length(X_temp)
        X_temp[i] /= alpha # form X0
    end

    # L is a lower bound for the smallest singular value of X0
    smin_est = norm(X_temp, 1)/cond(X_temp, 1)
    L  = smin_est/convert(T, sqrt(n))

    common_iter!(QDWHUpdater(alg.piv, L), X, U, H, alg.maxiter, alg.verbose, alg.tol)
end

type QDWHUpdater{T} <: PolarUpdater
    piv::Bool   # whether to pivot QR factorization
    L::T  # a lower bound for the smallest singluar value of each update matrix U
end

function update_U!{T}(upd::QDWHUpdater, U::Matrix{T})
    piv = upd.piv
    L = upd.L
    m, n = size(U)
    B = Array(T, m+n, n)
    Q1 = Array(T, n, n)
    Q2 = Array(T, n, n)
    # Compute paramters L, a, b, c
    L2 = L^2
    dd = try
        (4 * (1 - L2)/L2^2)^(1/3)
        catch
        (complex(4 * (1 -L2)/L2^2, 0))^(1/3)
    end
    sqd = sqrt(1+dd)
    a = sqd + 0.5 * sqrt(8 - 4 * dd + 8 * (2 - L2)/(L2 * sqd))
    a = real(a)
    b = (a - 1)^2 / 4
    c = a + b - 1

    # update L
    upd.L = L * (a + b * L2)/(1 + c * L2)

    copy!(B, [sqrt(c)*U; eye(n)])
    if piv
        VERSION < v"0.4.0-dev+1827" ? F = qrfact(B, pivot = true) : F = qrfact(B, Val{true})
    else
        F = qrfact(B)
    end
    copy!(Q1, full(F[:Q])[1:m, :])
    copy!(Q2, full(F[:Q])[m+1:end, :])
    copy!(U, b / c * U + (a - b / c) / sqrt(c) * Q1 * Q2')

end


text:
snd(x::String) = append!(played, registers[x])
set(x::String, y::String) = registers[x] = registers[y]
set(x::String, y::Int) = registers[x] = y
add(x::String, y::String) = registers[x] += registers[y]
add(x::String, y::Int) = registers[x] += y
mul(x::String, y::String) = registers[x] *= registers[y]
mul(x::String, y::Int) = registers[x] *= y
mod(x::String, y::String) = registers[x] = rem(registers[x], registers[y])
mod(x::String, y::Int) = registers[x] = rem(registers[x], y)
rcv(x::String, y::String) = (registers[x] != 0) ? registers[x] = last(played) : nothing
# rcv(x::String, y::Int) =

function puzzle18(path::String=joinpath(@__DIR__, "..", "data/18.txt"))

    instructions = readlines(path)

    # Initialize registers...
    registers = Dict();
    [registers[x] = 0 for x in string.('a':'z')];

    played = [];
    idx = 1;
    while idx < length(instructions)
        # Parse instruction...
        instruction = instructions[idx]
        command = instruction[1:3]
        x = instruction[5:5]
        y = instruction[7:end]
        if command == "snd"
            append!(played, registers[x])
            idx += 1
        elseif command == "set"
            try
                y = parse(Int, y)
            catch
                y = registers[string(y)]
            end
            registers[x] = y
            idx += 1
        elseif command == "add"
            y = parse(Int, instruction[7:end])
            if typeof(y) == Symbol
                y = registers[string(y)]
            end
            registers[x] += y
            idx += 1
        elseif command == "mul"
            y = parse(Int, instruction[7:end])
            if typeof(y) == Symbol
                y = registers[string(y)]
            end
            registers[x] *= y
            idx += 1
        elseif command == "mod"
            try
                y = parse(Int, instruction[7:end])
            catch
                y = Symbol(instruction[7:end])
            end
            if typeof(y) == Symbol
                y = registers[string(y)]
            end
            registers[x] = rem(registers[x], y)
            idx += 1
        elseif command == "rcv"
            if registers[x] != 0
                registers[x] = last(played)
                break
            end
            idx += 1
        elseif command == "jgz"
            try
                x = parse(Int, x)
            catch
                x = registers[x]
            end
            if x > 0
                idx += parse(Int, instruction[7:end])
            else
                idx += 1
            end
        end
    end

    return [last(played), nothing]
end


text:
abstract type Branch <: Device end
abstract type ACBranch <: Branch end
abstract type DCBranch <: Branch end

get_from_bus(b::T) where {T <: Branch} = b.arc.from
get_to_bus(b::T) where {T <: Branch} = b.arc.to


text:
module Etcetera

export longest_repeated_substrings
include("longest_repeated_substrings.jl")

include("longest_repeated.jl")
include("longest_repeated_no_overlap.jl")

module PrefixTree
include("PrefixTree.jl")
end
using .PrefixTree

end


text:
function czt(A::Complex128,W::Complex128,k::Int,xn::Vector{Complex128})
    #m = length of xn
    #k = desired output length
    #L : nextpow2(m-1+k))
    #yn: xn[n]*A^(-n)*W^(n^2/2), running from n = 0:(m-1), padded to L
    #vn: the inner filter W^(-n^2/2), running from n = (-m+1):(k-1), padded to L
    #vk: the outer product W^(m^2), running from m = 0:(k-1), padded to L
    #t : temporary variable
    m = length(xn)
    L = nextpow2(m+k-1)
    M = maximum([(m-1) (k-1)])
    yn = zeros(Complex128,L)
    vn = zeros(Complex128,L)
    vk = zeros(Complex128,L)
    W2 = sqrt(W)
    for i=1:(M+1)
        p = i-1
        t = W2^(p*p)
        yn[i] = (p<m) ? xn[i]*t*A^(-p) : zero(Complex128)
        vk[i] = (p<k) ? t : zero(Complex128)
        l1 = (i+m-1)
        if (l1<L)
            vn[i+m-1] = (p<k) ? one(Complex128)/t : zero(Complex128)
        end
        l2 = m-i+1
        if l2>0
            vn[l2] = one(Complex128)/t
        end
    end
    gn = ifft(fft(vn).*fft(yn))
    return [(gn[m+i-1]*vk[i])::Complex128 for i=1:k]
end

function generateU(ul::Vector{Complex128},xn::Vector{Complex128})
    L = length(ul) #Num of basis functions
    N = length(xn) #signal length
    p = 1 #calculate U matrices up to this num (1) for U0,U1
    K = ifloor((N-1-p)/2)
    NPOW = 3

    ul_inv = [1/ul[i] for i=1:L]
    ul_invK = [Base.power_by_squaring(ul_inv[i],K) for i=1:L]

    #initialize memory
    g0 = zeros(Complex128,L)
    g0_K = zeros(Complex128,L)
    g0_NK = zeros(Complex128,L)
    D0 = zeros(Complex128,L)
    U0 = zeros(Complex128,(L,L))
    U1 = zeros(Complex128,(L,L))
    ul_invk = ones(Complex128,L)
    for i=0:(K)
        for j=0:(L-1)
            g0[j+1] += ul_invk[j+1]*xn[i+1]
            g0_K[j+1] += ul_invk[j+1]*xn[K+1+i+1]
            g0_NK[j+1] += ul_invk[j+1]*xn[N-K+i]
            D0[j+1] += (i + 1)*xn[i+1]*ul_invk[j+1] + (K - i)*xn[i + K + 1+1]*ul_invk[j+1]*ul_inv[j+1]*ul_invK[j+1]
            ul_invk[j+1] = ((i%NPOW)==(NPOW-1)) ? Base.power_by_squaring(ul_inv[j+1],i+1) : ul_invk[j+1]*ul_inv[j+1]
        end
    end

    for i=1:L
        for j=1:i
            U0[i,j] = (1/(ul[i]-ul[j]))*(ul[i]*g0[j] - ul[j]*g0[i]
                    + ul_invK[j]*g0_K[i] - ul_invK[i]*g0_K[j])
            U0[j,i] = U0[i,j]
            U1[i,j] = (1/(ul[i]-ul[j]))*(ul[i]*ul[j]*g0[j] - ul[j]*ul[i]*g0[i]
                    + ul_invK[j]*ul[i]*g0_K[i] - ul_invK[i]*ul[j]*g0_K[j])
            U1[j,i] = U1[i,j]
        end
        U0[i,i] = D0[i]
        U1[i,i] = D0[i]*ul[i]-ul[i]*g0[i]+ul_invK[i]*g0_K[i]
    end
    return (U0,U1,g0,g0_K, g0_NK)
end

function xczt(xn::Vector{Complex128},A::Complex128,W::Complex128,K::Int,q::Float64)
    N = length(xn)
    Q = N*sum(abs(xn))
    q = Q*q
    Nul = minimum((N,100))
    ul = [exp(im*2*pi/Nul)^(i-1) for i=1:Nul] #fourier basis for expanding Krylov operators U0,U1
    Z = [A*W^(-(i-1)) for i=1:K] #this is the spiral the CZT is to be evaluated on
    M = ifloor((N-2)/2)+1
    (U0,U1,g0,g0_K,g0_NK) = generateU(ul,xn)
    genR(z,U0,U1) = U0-U1/z
    U0dagU0 = ctranspose(U0)*U0
    U1dagU1 = ctranspose(U1)*U1
    U1dagU0 = ctranspose(U1)*U0
    U0dagU1 = ctranspose(U0)*U1
    genRdagG0(R,g0) = ctranspose(R)*g0
    genRdagR(z,U0dagU0,U1dagU1,U1dagU0,U0dagU1) = U0dagU0 + U1dagU1 -z*U1dagU0 - U0dagU1/z
    Spec = zeros(Complex128,length(Z))
    for k=1:length(Z)
        R  = genR(Z[k],U0,U1)
        RC2 = genRdagG0(R,g0_NK*Z[k]^(-N+M))
        RR = genRdagR(Z[k],U0dagU0,U1dagU1,U1dagU0,U0dagU1)
        Spec[k] = (transpose(g0_K*Z[k]^(-M))*(\((RR+q^2),RC2)))[1]
    end
    return Spec+czt(A,W,K,xn)
end

function rczt(xn::Vector{Complex128},A::Complex128,W::Complex128,K::Int,q::Float64)
    N = length(xn)
    Q = N*sum(abs(xn))
    q = Q*q
    Nul = minimum((N,100))
    ul = [exp(im*2*pi/Nul)^(i-1) for i=1:Nul] #fourier basis for expanding Krylov operators U0,U1
    Z = [A*W^(-(i-1)) for i=1:K] #this is the spiral the CZT is to be evaluated on
    M = ifloor((N-2)/2)+1
    (U0,U1,g0,g0_K,g0_NK) = generateU(ul,xn)
    genR(z,U0,U1) = U0-U1/z
    U0dagU0 = ctranspose(U0)*U0
    U1dagU1 = ctranspose(U1)*U1
    U1dagU0 = ctranspose(U1)*U0
    U0dagU1 = ctranspose(U0)*U1
    genRdagG0(R,g0) = ctranspose(R)*g0
    genRdagR(z,U0dagU0,U1dagU1,U1dagU0,U0dagU1) = U0dagU0 + U1dagU1 -z*U1dagU0 - U0dagU1/z
    Spec = zeros(Complex128,length(Z))
    for k=1:length(Z)
        R  = genR(Z[k],U0,U1)
        RC = genRdagG0(R,g0)
        RR = genRdagR(Z[k],U0dagU0,U1dagU1,U1dagU0,U0dagU1)
        Spec[k] = (transpose(g0)*(\((RR+q^2),RC)))[1]
    end
    return Spec
end

function zTransformHorner(z::Complex128,xn::Vector{Complex128})
    #uses Horner's method
    x = one(Complex128)/z
    r = xn[end]
    for k=2:length(xn)
        r = r*x+xn[end-k+1]
    end
    return r*x
end

function zTransform(z::Complex128,xn::Vector{Complex128})
    x = one(Complex128)/z
    r = zero(Complex128)
    for k=1:length(xn)
        Z = x^k
        r += xn[k]*Z
    end
    return r
end

function alphaj(a::Complex128)
    return sqrt(1 - abs(a)^2)
end

function dictionaryElement(z::Complex128,a::Complex128)
    return alphaj(a)/(z-a)
end

function blaschkeElement(z::Complex128,a::Complex128)
    return (1-conj(a)*z)/(z-a)
end

function invBlaschkeElement(z::Complex128,a::Complex128)
    return (z-a)/(1-conj(a)*z)
end

function blaschkeProd(z::Complex128,aj::Vector{Complex128})
    r = one(z)
    for i=1:length(aj)
        r *= blaschkeElement(z,aj[i])
    end
    return r
end

function blaschkeProdExceptQ(z::Complex128,aj::Vector{Complex128},q::Int)
    r = one(z)
    for i=1:length(aj)
        if i!=q
            r *= blaschkeElement(z,aj[i])
        end
    end
    return r
end

function innerProduct(g::Complex128,a::Complex128)
    #a must be INTERIOR to the unit disk
    return alphaj(a)*(1/conj(a))*g
end

function innerProduct(G::Vector{Complex128},Z::Vector{Complex128})
    #Z is a collection poles INTERIOR to the unit disk
    L = length(Z)
    if length(G)!=L
        error("Length of G must match that of Z")
    end
    OverlapOld = zero(Float64)
    OverlapNew = zero(Float64)
    I = zero(Complex128)
    maxind = zero(Int64)
    f(On,Oo,Oin,Oio,Ivn,Ivo) = On > Oo ? (On,Oin,Ivn) : (Oo,Oio,Ivo)
    for i=1:L
        t = innerProduct(G[i],Z[i])
        OverlapNew = abs2(t)
        (OverlapOld,maxind,I) = f(OverlapNew,OverlapOld,i,maxind,t,I)
    end
    return (Z[maxind],I,maxind)
end

function OverlapVector(G::Vector{Complex128},Z::Vector{Complex128})
    #Z is a collection poles INTERIOR to the unit disk
    L = length(Z)
    if length(G)!=L
        error("Length of G must match that of Z")
    end
    O = zeros(Float64,L)
    for i=1:L
        O[i] = abs2(innerProduct(G[i],Z[i]))
    end
    return O
end

function g1!(G::Vector{Complex128},f::Function,Q::Vector{Complex128})
    #Q is a collection of poles EXTERIOR to the unit disk
    L = length(G)
    if length(Q)!=L
        error("input G vector must have same length as input Q vector")
    end
    for i=1:L
        G[i] = f(Q[i])
    end
end

function g1czt!(G::Vector{Complex128},xn::Vector{Complex128},Q::Vector{Complex128})
    zxn = vcat(zero(Complex128),xn)
    #Q is a collection of poles EXTERIOR to the unit disk
    lcol(i,l) = ((i-1)*l+i):((i-1)*l+((l+i-1)-2*(i-1)))
    rcol(i,l) = (((l-i)*l+((l+i-1)-2*(i-1)):-1:((l-i)*l+i)))
    brow(i,l) = (i*l+((l+i-1)-2*(i-1))):l:((l-i-1)*l+((l+i-1)-2*(i-1)))
    trow(i,l) = ((l-i-1)*l+i):-l:(i*l+i)
    L = length(G)
    l = int(sqrt(L)) #L should be a perfect square, and l is a power of 2 by construction
    if length(Q)!=L
        error("input G vector must have same length as input Q vector")
    end
    ulCorner(i,l) = ((i-1)*l+i) #upper left corner
    for i=1:int(l/2)
        side_len = l-2*(i-1)
        A = Q[ulCorner(i,l)]
        K = (4*(l-2*(i-1))-4)
        W = exp(-2*pi*im/K)
        temp = czt(A,W,K,zxn)
        G[lcol(i,l)]=temp[1:side_len]
        G[brow(i,l)]=temp[(side_len+1):(2*(side_len-1))]
        G[rcol(i,l)]=temp[(2*side_len-1):(3*side_len-2)]
        G[trow(i,l)]=temp[(3*side_len-1):(4*(side_len-1))]
    end
end

function dictBlaschke(z::Complex128,a::Complex128)
    return alphaj(a)/(1-conj(a)*z)
end

function gNext!(Gn::Vector{Complex128},Gnm1::Vector{Complex128},Anm1::Complex128,Onm1::Complex128,Q::Vector{Complex128})
    #this function modifies Gn, so returns nothing.
    #Q is a collection of poles EXTERIOR to the unit disk
    L = length(Q)
    if length(Gnm1)!=L
        error("input gnm1 vector must have same length as input z vector")
    end
    if length(Gn)!=L
        error("input gn vector must have same length as input z vector")
    end
    for i=1:L
        Gn[i] = (Gnm1[i]*invBlaschkeElement(Q[i],Anm1) - Onm1*dictBlaschke(Q[i],Anm1))
        #Gn[i] = ((Gnm1[i] - Onm1*dictionaryElement(Q[i],Anm1))*invBlaschkeElement(Q[i],Anm1))
    end
end

function impulseTakenakaMalmquist(n,aj::Vector{Complex128})
    if n<=zero(n)
        return zero(Complex128)
    end
    temp1 = zero(Complex128)
    temp2 = zero(Complex128)
    for i=1:(length(aj)-1)
        temp1 += aj[i]^(n-1)*dictionaryElement(aj[i],aj[end])*blaschkeProdExceptQ(aj[i],aj[1:(length(aj)-1)],i)*(1-abs(aj[i])^2)
    end
    temp2 += aj[end]^(n-1)*alphaj(aj[end])*blaschkeProd(aj[end],aj[1:(length(aj)-1)])
    return temp1 + temp2
end

function unitDiskGrid(n::Int,delta::Float64)
    N = nextpow2(n)
    M = N>>1
    v = linspace(delta,one(Float64)-delta,M)
    Zjk = zeros(Complex128,N,N)
    for j=1:(M)
             #bottom, running counter-clockwise
            Zjk[(M+j),((M-(j-1)):(M+1+(j-1)))] = v[j]*exp(im*linspace(5*pi/4,7*pi/4,2*j))
            #top, running clockwise
            Zjk[(M-j+1),((M-(j-1)):(M+1+(j-1)))] = v[j]*exp(im*linspace(3*pi/4,1*pi/4,2*j))
    end
    for j=1:(M)
            #right, running clockwise
            Zjk[((M-(j-1)):(M+1+(j-1))),M+j] = v[j]*exp(im*linspace(pi/4,-pi/4,2*j)) #right
            #left, running counter-clockwise
            Zjk[((M-(j-1)):(M+1+(j-1))),M-j+1] = v[j]*exp(im*linspace(3*pi/4,5*pi/4,2*j)) #left
    end
    Qjk = one(Complex128)./conj(Zjk)
    return (reshape(Zjk,N^2), reshape(Qjk,N^2))
end

function healG!(G::Vector{Complex128},x::Vector{Int},L::Int)
    for k=1:L
        G[x[k]] = zero(Complex128)
    end
end

function reconstructTimeDomain(poles::Vector{Complex128},amps::Vector{Complex128},t)
    if length(poles)!=length(amps)
        error("Input pole vector must have same length as amplitude vector")
    end
    r = zero(Complex128)
    for k=1:length(amps)
        p = k > 1 ? poles[1:k] : [poles[1]]
        r += amps[k]*impulseTakenakaMalmquist(t,p)
    end
    return r
end

function residualSumOfSquares(model::Vector{Complex128},data::Vector{Complex128})
    if length(model)!=length(data)
        error("Input model vector must have same length as the data vector!")
    end
    r = zero(Float64)
    for k=1:length(data)
        r += abs2(data[k]-model[k])
    end
    return r
end

function akaikeInformationCriterion(RSS,numParams,dataLen)
    #This gives the sample size corrected Akaike Information Criteria
    #commonly known as AICc in the literature
    K = numParams+1
    return dataLen*log(RSS/dataLen)+2*K+2*K*((K+1)/(dataLen-K-1))
end

function afd(xn::Vector{Complex128},Nmax::Int,Nz::Int)
    Nz = int(floor(sqrt(Nz)))
    (Z,Q) = unitDiskGrid(Nz,0.005)
    L = length(Z)
    rNz = int(sqrt(L))
    Gnm1 = zeros(Complex128,L)
    Gn   = zeros(Complex128,L)
    polesList = zeros(Complex128,Nmax)
    polesAmps = zeros(Complex128,Nmax)
    poleIndList = zeros(Int64,Nmax)
    g1czt!(Gnm1,xn,Q)
    (poleNm1,poleAmpNm1,poleIndNm1) = innerProduct(Gnm1,Z)
    polesList[1] = poleNm1
    polesAmps[1] = poleAmpNm1
    poleIndList[1] = poleIndNm1
    for i=2:Nmax
        gNext!(Gn,Gnm1,poleNm1,poleAmpNm1,Q)
        healG!(Gn,poleIndList,i-1)
        (poleNm1,poleAmpNm1,poleIndNm1) = innerProduct(Gn,Z)
        polesList[i]=poleNm1
        polesAmps[i]=poleAmpNm1
        poleIndList[i] = poleIndNm1
        Gnm1 = Gn
    end
    return polesList, polesAmps
end

function xafd(xn::Vector{Complex128},Nmax::Int,Nz::Int,xN::Int,q::Float64)
    Nz = int(floor(sqrt(Nz)))
    (Z,Q) = unitDiskGrid(Nz,0.005)
    #extrapolate the time domain using xczt to xN points
    xhat = ifft(xczt(xn,one(Complex128),exp(-im*2*pi/xN),xN,q))
    L = length(Z)
    rNz = int(sqrt(L))
    Gnm1 = zeros(Complex128,L)
    Gn   = zeros(Complex128,L)
    polesList = zeros(Complex128,Nmax);
    polesAmps = zeros(Complex128,Nmax);
    poleIndList = zeros(Int64,Nmax);
    g1czt!(Gnm1,xhat,Q)
    (poleNm1,poleAmpNm1,poleIndNm1) = innerProduct(Gnm1,Z)
    polesList[1] = poleNm1
    polesAmps[1] = poleAmpNm1
    poleIndList[1] = poleIndNm1
    for i=2:Nmax
        gNext!(Gn,Gnm1,poleNm1,poleAmpNm1,Q)
        healG!(Gn,poleIndList,i-1)
        (poleNm1,poleAmpNm1,poleIndNm1) = innerProduct(Gn,Z)
        polesList[i]=poleNm1
        polesAmps[i]=poleAmpNm1
        poleIndList[i] = poleIndNm1
        Gnm1 = Gn
    end
    return polesList, polesAmps
end

function optimalAfdTruncation(poles::Vector{Complex128},amps::Vector{Complex128},xn::Vector{Complex128})
    L = length(poles)
    Nt = length(xn)
    if length(amps)!=L
        error("# of poles should be the same as number of amplitudes")
    end
    infVec = zeros(Float64,L)
    for k=1:L
        model = [reconstructTimeDomain(poles[1:k],amps[1:k],i)::Complex128 for i=linspace(1,Nt,Nt)]
        rssval = residualSumOfSquares(model,xn)
        # one complex value = 2 parameters, so each pole and amp is 4 parameters, each complex valued input
        # is worth 2 data points of input in this thinking, thus:
        infVec[k] = akaikeInformationCriterion(rssval,4*k,2*length(xn))
    end
    infMinInd = indmin(infVec)
    aicVec = exp((-infVec+infVec[infMinInd])/2)
    return indmax(aicVec)
end


text:
module TestScoring

import Wordle
import Test: @test, @test_throws

trials = Wordle.Trial[
    Wordle.Trial(
        Wordle.Candidate("alone"),
        Wordle.Result("⬜⬜🟩⬜⬜")
    ),
]

never_seen_chars = Wordle.never_seen(trials)
@test !never_seen_chars[Wordle.to_index('a')]
@test !never_seen_chars[Wordle.to_index('l')]
@test !never_seen_chars[Wordle.to_index('o')]
@test !never_seen_chars[Wordle.to_index('n')]
@test !never_seen_chars[Wordle.to_index('e')]
@test sum(never_seen_chars) === 26 - 5

trials = Wordle.Trial[
    Wordle.Trial(
        Wordle.Candidate("alone"),
        Wordle.Result("⬜⬜🟩⬜⬜")
    ),
    Wordle.Trial(
        Wordle.Candidate("trust"),
        Wordle.Result("⬜⬜⬜⬜⬜")
    ),
]

never_seen_chars = Wordle.never_seen(trials)
@test !never_seen_chars[Wordle.to_index('a')]
@test !never_seen_chars[Wordle.to_index('l')]
@test !never_seen_chars[Wordle.to_index('o')]
@test !never_seen_chars[Wordle.to_index('n')]
@test !never_seen_chars[Wordle.to_index('e')]
@test !never_seen_chars[Wordle.to_index('t')]
@test !never_seen_chars[Wordle.to_index('r')]
@test !never_seen_chars[Wordle.to_index('u')]
@test !never_seen_chars[Wordle.to_index('s')]
@test !never_seen_chars[Wordle.to_index('t')]

@test sum(never_seen_chars) === 26 - 5 - 4

end # module


text:
using stoosc
rc("font", family="arial")

# Model and its parameters
model_name = "model01"
parameter_names = ["A", "τ", "I"]
parameter_values = [0.1, 3.66, 0.3]

# 3 days in DD to see the effect of the intital conditions and then LD cycle
events = create_events([(:DD, 3), (:LD, 3, 0.5, 0.5)])

# Set up ODE model
ode_model = load_model(model_name, "ode")
ode_model = set_parameter_value(ode_model, parameter_names, parameter_values)
ode_model = set_initial_conditions(ode_model, [0.0, 0.0, 0.1])
ode_model = set_timespan(ode_model, (0.0, 6.0))
ode_model = set_input(ode_model, events)

# Simulate ODE model
ode_solution = simulate_population(ode_model)
t_ode = ode_solution.time
x_ode = ode_solution.mean[:, 3]

# Set up SDE model
sde_model = load_model(model_name, "sde")
sde_model = set_parameter_value(sde_model, parameter_names, parameter_values)
sde_model = set_parameter_value(sde_model, "σ", 0)
sde_model = set_initial_conditions(sde_model, [0.0, 0.0, 0.1])
sde_model = set_input(sde_model, events)
sde_model = set_timespan(sde_model, (0.0, 6.0))

# Simulate SDE model and create a figure for varying integration step dt
fig, ax_arr = subplots(4, 1, figsize=(7, 7))
for (i, dt) = enumerate([0.1, 0.01, 0.001, 0.0001])

    sde_model_2 = set_solver(sde_model, dt=dt)
    sde_solution = simulate_population(sde_model_2)
    t = sde_solution.time
    x = sde_solution.mean[:, 3]

    ax_arr[i].plot(t_ode, x_ode, color="black", label="ODE")
    ax_arr[i].plot(t, x, color="blue", label="SDE (σ = 0)")
    plot_events(events, ylims=[0.04, 0.16], ax=ax_arr[i])
    ax_arr[i].set_ylim(0.04, 0.16)
    ax_arr[i].set_xlabel("Time")
    ax_arr[i].set_ylabel("Z")
    ax_arr[i].set_title("dt = $dt")
    ax_arr[i].legend(ncol=2, edgecolor="black", framealpha=1.0)
    ax_arr[i].set_xlim(0, 6)
end
fig.tight_layout(pad=0.2)

# Save the figure
save_figure(fig, "varying_dt.svg")

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
include("util.jl")

lines = getlines("inputs/8.txt")
line = lines[1]
nums = map(x->parse(Int64,x), split(line, " "))

mutable struct Node
    children::Array{Node}
    metadata::Array{Int64}
end

function getnode(xs::Array{Int64})::Tuple{Node,Array{Int64}}
    numnodes = xs[1]
    nummetadatas = xs[2]
    thisnode = Node([],[])
    newxs = xs[3:end]
    for n in 1:numnodes
        node, newxs = getnode(newxs)
        push!(thisnode.children, node)
    end
    for m in 1:nummetadatas
        push!(thisnode.metadata, newxs[1])
        newxs = newxs[2:end]
    end
    return thisnode, newxs
end

function sumtree(x::Node)::Int64
    this_sum = sum(x.metadata)
    for c in x.children
        this_sum += sumtree(c)
    end
    return this_sum
end

function part_2_sumtree(x::Node)::Int64
    if length(x.children) == 0
        return sumtree(x)
    else
        this_sum = 0
        for m in x.metadata
            try
                cnode = x.children[m]
                this_sum += part_2_sumtree(cnode)
            catch
                this_sum += 0
            end
        end
        return this_sum
    end
end

tree = getnode(nums)[1]

part_1 = sumtree(tree)
@time part_1 = sumtree(tree)
part_2 = part_2_sumtree(tree)
@time part_2 = part_2_sumtree(tree)


text:
using NeighbourLists
using JuLIP
using BenchmarkTools

include("test_aux.jl")

print("N-body energy and forces benchmark: ")
println("# Threads = ", Base.Threads.nthreads())

at = bulk(:Cu, cubic=true) * 5
println("bulk :Cu with nat = $(length(at))")
r0 = rnn(:Cu)
rcut = 2.1 * r0
X = positions(at)
C = cell(at)
f, f_d = gen_fnbody(rcut, r0)
nlist = PairList(X, rcut, C, (false, false, false), sorted = true)

for M in [2, 3, 4, 5]
   println("M = $M")
   print("  energy: ")
   @btime n_body($f, $M, $nlist)
   print("  forces: ")
   @btime grad_n_body($f_d, $M, $nlist)
end

lj = LennardJones(r0, 1.0) * C1Shift(cutoff)
@btime forces($lj, $at)


text:
using DelimitedFiles

data= readdlm("input_day2", '\n')

valid_pass = 0

for line in data

   i1 = findfirst(isequal('-'),line)
   i2 = findfirst(isequal(' '),line)

   #get minimum number of occurrences
   min_o = parse(Int, line[1:i1-1])
   #get maximum number of occurrences
   max_o = parse(Int, line[i1+1:i2-1])

   letter = line[i2+1]
   password = line[i2+4:end]

   rx = Regex(string(letter))
   m = eachmatch(rx,password)
   ocurrences = size(collect(m))[1]

   if  min_o <= ocurrences <= max_o
     global  valid_pass += 1
   end

end

println(valid_pass)


text:
# LPA_wb_plus.jl
# Label propagation algorithm for weighted bipartite networks that finds modularity.
# Contains the LPAwb+ and the DIRTLPAwb+ algorithms
# Author :  Stephen Beckett ( https://github.com/sjbeckett/weighted-modularity-LPAwbPLUS )
# MIT License

###############################################

function BarbersMatrix(MATRIX)
	return MATRIX - (sum(MATRIX,2)*sum(MATRIX,1))/sum(MATRIX);
end

###############################################

function WEIGHTEDMODULARITY(BMatrix,Matsum,redlabels::Array{Int64,1},bluelabels::Array{Float64,2})
	#see equation 8
	holdsum = 0;

	for rr = 1:length(redlabels)
	    for cc = 1:length(bluelabels)
	        kroneckerdelta = redlabels[rr] == bluelabels[cc];
	        holdsum += BMatrix[rr,cc] * kroneckerdelta;
	    end
	end
	return holdsum/Matsum;
end

###############################################
function TRACE(MATRIX)
	return sum(diag(MATRIX))
end

###############################################

function WEIGHTEDMODULARITY2(BMatrix,Matsum,redlabels,bluelabels)
	#see equation 9
	UNIred = unique(redlabels)
	Lred = length(redlabels)
	UNIblu = unique(bluelabels)
	Lblu = length(bluelabels)
	LABELMAT1 = zeros(length(UNIred),Lred)
	LABELMAT2 = zeros(Lblu,length(UNIblu))

	for aa = 1:Lred
		LABELMAT1[findin(UNIred,redlabels[aa]),aa] = 1
	end

	for aa = 1:Lblu
		LABELMAT2[aa,findin(UNIblu,bluelabels[aa])] = 1
	end

	return TRACE(LABELMAT1 * BMatrix * LABELMAT2)/Matsum

end

###############################################
function DIVISION(redlabels::Array{Int64,1},bluelabels::Array{Float64,2})
	divisionsFound= Int64[];

	for aa = 1:min(maximum(redlabels),maximum(bluelabels))
		if( (isempty(findin(redlabels,aa))==false) && (isempty(findin(bluelabels,aa))==false) )
			push!(divisionsFound,aa)
		end
	end

	return divisionsFound;
end

###############################################

function StageTwo_LPAwbdash(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels, Qb_now)

	divisionsFound = DIVISION(redlabels,bluelabels);
	NUMdiv = length(divisionsFound);
	IterateFlag = 1;
	while IterateFlag == 1
		CombinedDivisionsThisTime = 0;
		if NUMdiv > 1
			for div1check = 1:NUMdiv-1
				Mod1 = divisionsFound[div1check]
				for div2check = div1check+1:NUMdiv
					CHECK_RED = copy(redlabels);
					CHECK_RED[findin(redlabels,Mod1)] = divisionsFound[div2check];
					CHECK_BLUE = copy(bluelabels);
					CHECK_BLUE[findin(bluelabels,Mod1)] = divisionsFound[div2check];
					QQ = WEIGHTEDMODULARITY2(BMatrix,Matsum,CHECKRED,CHECKBLUE);
					if QQ > Qb_now #If a division will improve modularity - find the best way to do this
                    				FoundBetter = 0;
                    				for aa = 1:NUMdiv
                        				CHECK_RED2 = copy(redlabels);
                        				CHECK_RED2[findin(redlabels,divisionsFound[aa])] = Mod1;
                        				CHECK_BLUE2 = copy(bluelabels);
                        				CHECK_BLUE2[findin(bluelabels,divisionsFound[aa])] = Mod1;
                        				if WEIGHTEDMODULARITY2(BMatrix,Matsum,CHECK_RED2,CHECK_BLUE2) > QQ
                            					FoundBetter = 1;
                        				end
                        				CHECK_RED2 = copy(redlabels);
                        				CHECK_RED2[findin(redlabels,divisionsFound[aa])] = divisionsFound[div2check];
                        				CHECK_BLUE2 = copy(bluelabels);
                        				CHECK_BLUE2[findin(bluelabels,divisionsFound[aa])] = divisionsFound[div2check];
                        				if WEIGHTEDMODULARITY2(BMatrix,Matsum,CHECK_RED2,CHECK_BLUE2) > QQ
                            					FoundBetter = 1;
                        				end
                    				end
                    				if FoundBetter == 0 #If no better configuration found - JOIN.
                        				redlabels = copy(CHECK_RED);
                        				bluelabels = copy(CHECK_BLUE);
                        				CombinedDivisionsThisTime += 1;
                    				end
                			end
				end
			end
			if CombinedDivisionsThisTime == 0 #If no divisions were joined move on
		            IterateFlag = 0;
		        end
		else
		IterateFlag = 0;
		end

		(redlabels,bluelabels,Qb_now) = StageOne_LPAwbdash(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels);
		divisionsFound = DIVISION(redlabels,bluelabels);
		NUMdiv = length(divisionsFound);
	end

return redlabels, bluelabels, Qb_now;

end

###############################################

function StageOne_LPAwbdash(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels)
	#Create storage containers for total marginals attached to each red(row)
	#label and blue(column) label

	BLUELABELLENGTH = length(bluelabels);
	REDLABELLENGTH = length(redlabels);
	TotalRedDegrees = NaN*ones(maximum(redlabels),1);

	TotalBlueDegrees = NaN*ones(max(BLUELABELLENGTH,REDLABELLENGTH),1);

	#Fill up these containers according to current labels
	#Red
	for aa = 1:REDLABELLENGTH
	    if isnan(TotalRedDegrees[redlabels[aa]])
        	TotalRedDegrees[redlabels[aa]] = row_marginals[aa];
	    else
        	TotalRedDegrees[redlabels[aa]] += row_marginals[aa];
	    end
	end

	#Blue
	if sum(isnan(bluelabels)) != BLUELABELLENGTH #occurs first time through as blue nodes unlabelled
	    for bb = 1:BLUELABELLENGTH
        	if isnan(TotalBlueDegrees[bluelabels[bb]])
	            TotalBlueDegrees[bluelabels[bb]] = col_marginals[bb];
        	else
        	    TotalBlueDegrees[bluelabels[bb]] += col_marginals[bb];
        	end
	    end
	else
	    TotalBlueDegrees = zeros(max(BLUELABELLENGTH,REDLABELLENGTH),1);
	end

	#locally maximise modularity!
	(redlabels,bluelabels,Qb_now) = LOCALMAXIMISATION(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels,TotalRedDegrees,TotalBlueDegrees);

	return redlabels, bluelabels, Qb_now;

end

###############################################

function LOCALMAXIMISATION(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels,TotalRedDegrees,TotalBlueDegrees)

	#Find score for current partition
	QbAfter = WEIGHTEDMODULARITY2(BMatrix,Matsum,redlabels,bluelabels);

	IterateFlag = 1;
	while IterateFlag == 1
		#Save old information
		QbBefore = QbAfter;
		old_redlabels = copy(redlabels);
		old_bluelabels = copy(bluelabels);
		old_TRD = copy(TotalRedDegrees);
		old_TBD = copy(TotalBlueDegrees); #13

		#Update Blue Nodes using red node information (see equation 10)
		bluelabelchoices = 1:maximum(redlabels);

		for bb = 1:length(bluelabels)
			if isnan(bluelabels[bb]) == false
				TotalBlueDegrees[bluelabels[bb]] -= col_marginals[bb];
			end
			changebluelabeltest = NaN*ones(1,length(bluelabelchoices));

			for ww = 1:length(bluelabelchoices)
				changebluelabeltest[ww] = sum( (redlabels .== bluelabelchoices[ww]) .* MATRIX[:,bb])  -  col_marginals[bb]*TotalRedDegrees[bluelabelchoices[ww]]/Matsum ;
			end

			#assign new label based on maximisation of above condition

			labels = findin(changebluelabeltest, maximum(changebluelabeltest));
			newlabelindex = labels[rand(1:length(labels),1)];
			bluelabels[bb] = bluelabelchoices[newlabelindex[1]];

			if bluelabels[bb] > length(TotalBlueDegrees)
	                	TotalBlueDegrees[bluelabels[bb]] = 0;
	        	end

        		#Update total marginals on new labelling
            		TotalBlueDegrees[bluelabels[bb]] += col_marginals[bb];
        	end#40

		#Now update red node labels based on blue node information (see equation 10)
		redlabelchoices = 1:maximum(bluelabels);

		for aa = 1:length(redlabels)
			TotalRedDegrees[redlabels[aa]] -= row_marginals[aa];
			changeredlabeltest = NaN*ones(1,length(redlabelchoices));

			for ww = 1:length(redlabelchoices)
				changeredlabeltest[ww] = sum( (bluelabels .== redlabelchoices[ww]) .* MATRIX[aa,:])  -  row_marginals[aa]*TotalBlueDegrees[redlabelchoices[ww]]/Matsum ;
			end

			#assign new label based on maximisation of above condition
			labels = findin(changeredlabeltest,maximum(changeredlabeltest));
			newlabelindex = labels[rand(1:length(labels),1)];
 			redlabels[aa] = redlabelchoices[newlabelindex[1]];

			if redlabels[aa] > length(TotalRedDegrees)
				TotalRedDegrees[redlabels[aa]] = 0;
			end
			TotalRedDegrees[redlabels[aa]] +=  row_marginals[aa];
		end

		#Find the new modularity score based on node label updates.
        	QbAfter = WEIGHTEDMODULARITY(BMatrix,Matsum,redlabels,bluelabels);

        	#If this modularity is not as good as previous stop iterating and
        	#use that previous best information
        	if QbAfter <= QbBefore
            		redlabels = copy(old_redlabels);
	            	bluelabels = copy(old_bluelabels);
        	    	TotalRedDegrees = copy(old_TRD);
        	    	TotalBlueDegrees = copy(old_TBD);
        	    	IterateFlag = 0;
        	end

	end

	Qb_now = QbAfter;

	return redlabels, bluelabels, Qb_now;
end

###############################################

function LPA_wb_plus(MATRIX, initialmoduleguess = NaN)

	#Make sure the smallest matrix dimension represent the red labels by making them the rows (If matrix is transposed here, will be transposed back at the end)
	flipped = 0;
	if size(MATRIX,1) > size(MATRIX,2)
		MATRIX = MATRIX';
		flipped = 1;
	end

	Matsum = sum(MATRIX);
	col_marginals = sum(MATRIX,1);
	row_marginals = sum(MATRIX,2);
	BMatrix = BarbersMatrix(MATRIX);

	#initiliase labels
	bluelabels = [NaN*ones(1,size(MATRIX,2))];
	if isnan(initialmoduleguess)
		redlabels = [1:size(MATRIX,1)];
	else
		redlabels = rand(1:(initialmoduleguess+1),size(MATRIX,1));
	end

	#Run Phase 1: Locally update lables to maximise Qb
	(redlabels,bluelabels,Qb_now) = StageOne_LPAwbdash(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels);

	#Run Phase 2: Connect divisions from top-down if it improves Qb, then run
	#phase 1 again. Repeat till Qb cannot be improved.
	(redlabels,bluelabels,Qb_now) = StageTwo_LPAwbdash(row_marginals,col_marginals,MATRIX,BMatrix,Matsum,redlabels,bluelabels,Qb_now);

	if flipped==1 #If matrix was flipped, swap row and column labels
    		holder = copy(redlabels);
		redlabels = bluelabels;
		bluelabels = holder;
	end

	return redlabels, bluelabels, Qb_now;
end

###############################################

function DIRT_LPA_wb_plus(MATRIX, mini = 4, reps = 10)

	A = LPA_wb_plus(MATRIX);

	mods = length(unique(A[1]));

	if (mods-mini) > 0
		for aa = mini:mods
			for bb = 1:reps
				B = LPA_wb_plus(MATRIX,aa)
				if B[3] > A[3]
					A = B
				end
			end
		end
	end

	return A[1], A[2], A[3];
end

max_stars_repo_name: sjbeckett/weighted-modularity-LPAwbPLUS

max_issues_repo_name: sjbeckett/weighted-modularity-LPAwbPLUS

max_forks_repo_name: sjbeckett/weighted-modularity-LPAwbPLUS