text:
module iPredict
include("TableProcessing.jl")
include("FindKm.jl")
include("Simulate.jl")
include("Plots.jl")
include("pValues.jl")
import Cairo
import Fontconfig
import DataFrames
using Distributions
using Statistics
using DelimitedFiles

function InvestigateIsoformChoice(countsMatrixPath::String, GeneIsoTablePath::String,
     NumIsoformsToSimulate::Int64,
     filteringThreshold::Int64)

     global pFN = 0.04
     global pFP = 0.01

     #Read in counts matrix
     println("Data wrangling commencing!")
     println("Reading in counts matrix.") #Quick
     countsMatrix = ReadCountsMatrix(countsMatrixPath)

     #Read in gene-isoform table
     println("Reading in gene-isoform relationships")
     GeneIsoformRelationships = ReadGeneIsoformRelationships(GeneIsoTablePath)

     #Make sure gene-isoform table matches counts matrix
     if CheckGenesMatchIsoforms(countsMatrix, GeneIsoformRelationships) == false
         throw(ArgumentError("Your counts matrix doesn't match your list of genes
         and isoforms - there is something wrong with your input files."))
     end

     #Add gene information to isoforms table
     println("Adding Gene information to countsMatrix")
     countsMatrixWithGenes = AddGeneInfo(countsMatrix, GeneIsoformRelationships)

     #Clean data for finding Km
     println("Cleaning data for finding Km")
     MMdf = CleanDataForMM(countsMatrix)

     #Find Km
     println("Calculating Km....")
     MMParams = FindKm(MMdf)
     Km = MMParams[1]
     sigma = MMParams[2]

     #Extract genes with filteringThreshold expressed isoforms
     println("Converting counts to dict")
     GenesWithNIsoforms = KeepGenesWithNIsoforms(countsMatrixWithGenes,
     filteringThreshold)
     workingDict = GenesWithNIsoforms[1]
     unexpr_dict = GenesWithNIsoforms[2]
     MaxNumIsos = GenesWithNIsoforms[3]

     #Get ranking and dropout info for each gene
     rankingDict = Dict()
     pObsDict = Dict()
     pDropoutDict = Dict()
     genesList = []

     for (key, value) in workingDict
         rankingDict[key] = FindRank(value)
         pObsDict[key] = FindPObs(value)
         pDropoutDict[key] = FindPDropout(value, Km)
         push!(genesList, key)
     end

     df = FindMeanUnifObsProbabilities(NumIsoformsToSimulate,
     genesList, pObsDict, pDropoutDict, rankingDict)

     return df
end

#Function that performs global level predictions
function globalPredict(countsMatrixPath::String, GeneIsoTablePath::String,
     NumSimulations::Int64, output::String, NumIsoformsToSimulate::Int64,
     filteringThreshold::Int64, test::Bool, alpha, beta, pFNeg::Float64,
     pFPos::Float64, error_model::String, choice_model::String, path::String)

    #Set probability(FP|N) and probability(FN|P)
    if error_model == "none"
        global pFN = pFNeg
        global pFP = pFPos
    end

    #Read in counts matrix
    println("Data wrangling commencing!")
    println("Reading in counts matrix.") #Quick
    countsMatrix = ReadCountsMatrix(countsMatrixPath)

    #Read in gene-isoform table
    println("Reading in gene-isoform relationships")
    GeneIsoformRelationships = ReadGeneIsoformRelationships(GeneIsoTablePath)

    #Make sure gene-isoform table matches counts matrix
    if CheckGenesMatchIsoforms(countsMatrix, GeneIsoformRelationships) == false
        throw(ArgumentError("Your counts matrix doesn't match your list of genes
        and isoforms - there is something wrong with your input files."))
    end

    #Add gene information to isoforms table
    println("Adding Gene information to countsMatrix")
    countsMatrixWithGenes = AddGeneInfo(countsMatrix, GeneIsoformRelationships)

    #Only find Km if no Beta distribution parameters provided
    if isnothing(alpha) || isnothing(beta)

        #Clean data for finding Km
        println("Cleaning data for finding Km")
        MMdf = CleanDataForMM(countsMatrix)

        #Find Km
        println("Calculating Km....")
        MMParams = FindKm(MMdf)
        Km = MMParams[1]
        sigma = MMParams[2]
    end

    #Extract genes with filteringThreshold expressed isoforms
    println("Converting counts to dict")
    GenesWithNIsoforms = KeepGenesWithNIsoforms(countsMatrixWithGenes,
    filteringThreshold)
    workingDict = GenesWithNIsoforms[1]
    unexpr_dict = GenesWithNIsoforms[2]
    MaxNumIsos = GenesWithNIsoforms[3]

    #Get ranking and dropout info for each gene
    rankingDict = Dict()
    pObsDict = Dict()
    pDropoutDict = Dict()
    genesList = []

    #How we calculate dropouts depends on whether or not beta parameters are provided
    if isnothing(alpha) || isnothing(beta)
        for (key, value) in workingDict
            rankingDict[key] = FindRank(value)
            pObsDict[key] = FindPObs(value)
            pDropoutDict[key] = FindPDropout(value, Km)
            push!(genesList, key)
        end
    else
        for (key, value) in workingDict
            rankingDict[key] = FindRank(value)
            pObsDict[key] = FindPObs(value)
            pDropoutDict[key] = rand(Beta(alpha, beta), length(value))
            push!(genesList, key)
        end
    end

    plot_arr = []
    log_plot_arr = []
    mean_arr = []
    median_arr = []
    std_arr = []

    println("Beginning simulations!")

    for i in 1:NumIsoformsToSimulate

        #I pre-allocate arr in the hope it will make my code faster
        arr = Array{Float64}(undef, NumSimulations * length(genesList))
        arr_index = 1

        overlapArr = Array{Float64}(undef, NumSimulations * length(genesList))
        overlap_index = 1

        pDropoutArr = Array{Float64}(undef, NumSimulations * length(genesList) *
        (ncol(countsMatrix) - 1) * i)
        drop_index = 1

        #Need to add a for loop with num simulations here
        for j in 1:NumSimulations

            #@show arr_index

            #Get an array of mean isoforms per gene per cell
            simulation_out = globalArraySimulation(1, genesList, rankingDict, pObsDict,
            pDropoutDict, ncol(countsMatrix) - 1, i, filteringThreshold,
            unexpr_dict, MaxNumIsos, error_model, choice_model)

            arr[arr_index:(arr_index + length(genesList) - 1)] = simulation_out[1]
            overlapArr[overlap_index:(overlap_index + length(genesList) - 1)] = simulation_out[2]
            pDropoutArr[drop_index:(drop_index + length(simulation_out[3]) - 1)] = simulation_out[3]

            #Increment arr_index
            arr_index += length(genesList)
            overlap_index += length(genesList)
            drop_index += length(simulation_out[3])
        end

        #Save histogram of data
        if i == 1
            p1 = makeNormalisedHistogram(arr, string(i) * " Isoform")
        else
            p1 = makeNormalisedHistogram(arr, string(i) * " Isoforms")
        end
        push!(plot_arr, p1[1])
        push!(log_plot_arr, p1[2])
        push!(mean_arr, mean(arr))
        push!(std_arr, std(arr))
        push!(median_arr, median(arr))

        #Write to file if required
        if path != "none"
            writedlm(path * "_" * string(i) * "_isoforms_array.txt", arr)
            writedlm(path * "_" * string(i) * "_isoforms_overlap_array.txt", overlapArr)
            writedlm(path * "_" * string(i) * "_isoforms_pdrop_array.txt", pDropoutArr)
        end
    end

    #Get an array and a dict of mean isoforms per gene per cell in real data
    println("Finding means in real data.")
    findRealMeans = MeanNumIsoformsPerGenePerCellReal(workingDict)
    realMeansArr = findRealMeans[1]

    #Write to file if required
    if path != "none"
        writedlm(path * "_" * "real_isoforms_array.txt", realMeansArr)
    end

    #Plot real data
    real_plot = makeNormalisedHistogramReal(realMeansArr, "Real")
    push!(plot_arr, real_plot[1])
    push!(log_plot_arr, real_plot[2])

    #Arrange plots and write to output
    graph = plot(plot_arr..., layout=(length(plot_arr),1),legend=false, size = (300, 600), bottom_margin = 1.0mm)
    savefig(graph, output)
    graph = plot(log_plot_arr..., layout=(length(plot_arr),1),legend=false, size = (300, 600))
    savefig(graph, output * "_log.png")

    return (mean_arr, std_arr, median_arr)

end

#Function that makes negative control model
function NegativeControl(countsMatrixPath::String, GeneIsoTablePath::String,
     NumSimulations::Int64, output::String, NumIsoformsToSimulate::Int64,
     filteringThreshold::Int64, path::String)

     #Read in counts matrix
     println("Data wrangling commencing!")
     println("Reading in counts matrix.") #Quick
     countsMatrix = ReadCountsMatrix(countsMatrixPath)

     #Read in gene-isoform table
     println("Reading in gene-isoform relationships")
     GeneIsoformRelationships = ReadGeneIsoformRelationships(GeneIsoTablePath)

     #Make sure gene-isoform table matches counts matrix
     if CheckGenesMatchIsoforms(countsMatrix, GeneIsoformRelationships) == false
         throw(ArgumentError("Your counts matrix doesn't match your list of genes
         and isoforms - there is something wrong with your input files."))
     end

     #Add gene information to isoforms table
     println("Adding Gene information to countsMatrix")
     countsMatrixWithGenes = AddGeneInfo(countsMatrix, GeneIsoformRelationships)

     println("Converting counts to dict")
     GenesWithNIsoforms = KeepGenesWithNIsoforms(countsMatrixWithGenes,
     filteringThreshold)
     workingDict = GenesWithNIsoforms[1]
     unexpr_dict = GenesWithNIsoforms[2]
     MaxNumIsos = GenesWithNIsoforms[3]

     #Get ranking and dropout info for each gene
     rankingDict = Dict()
     pObsDict = Dict()
     pDropoutDict = Dict()
     genesList = []

     #Fill arrays
     for (key, value) in workingDict
         rankingDict[key] = FindRank(value)
         pObsDict[key] = FindPObs(value)

         push!(genesList, key)
     end

     println("Beginning simulations!")

     for i in 1:NumIsoformsToSimulate

         #I pre-allocate arr in the hope it will make my code faster
         arr = Array{Float64}(undef, NumSimulations * length(genesList))
         arr_index = 1

         #Need to add a for loop with num simulations here
         for j in 1:NumSimulations

             #@show arr_index

             #Get an array of mean isoforms per gene per cell
             arr[arr_index:(arr_index + length(genesList) - 1)] =
             negControlSimulation(1, genesList, rankingDict, pObsDict,
             ncol(countsMatrix) - 1, i, filteringThreshold, pDropoutDict,
             MaxNumIsos)

             #Increment arr_index
             arr_index += length(genesList)
         end

         #Write to file if required
         if path != "none"
             writedlm(path * "_" * string(i) * "_isoforms_array.txt", arr)
         end
     end
     return true
end

#Function that performs global level predictions
function NegativeControlDrop(countsMatrixPath::String, GeneIsoTablePath::String,
     NumSimulations::Int64, output::String, NumIsoformsToSimulate::Int64,
     filteringThreshold::Int64, test::Bool, pFNeg::Float64,
     pFPos::Float64, error_model::String, choice_model::String, path::String)

    #Set probability(FP|N) and probability(FN|P)
    if error_model == "none"
        global pFN = pFNeg
        global pFP = pFPos
    end

    #Read in counts matrix
    println("Data wrangling commencing!")
    println("Reading in counts matrix.") #Quick
    countsMatrix = ReadCountsMatrix(countsMatrixPath)

    #Read in gene-isoform table
    println("Reading in gene-isoform relationships")
    GeneIsoformRelationships = ReadGeneIsoformRelationships(GeneIsoTablePath)

    #Make sure gene-isoform table matches counts matrix
    if CheckGenesMatchIsoforms(countsMatrix, GeneIsoformRelationships) == false
        throw(ArgumentError("Your counts matrix doesn't match your list of genes
        and isoforms - there is something wrong with your input files."))
    end

    #Add gene information to isoforms table
    println("Adding Gene information to countsMatrix")
    countsMatrixWithGenes = AddGeneInfo(countsMatrix, GeneIsoformRelationships)

    #Extract genes with filteringThreshold expressed isoforms
    println("Converting counts to dict")
    GenesWithNIsoforms = KeepGenesWithNIsoforms(countsMatrixWithGenes,
    filteringThreshold)
    workingDict = GenesWithNIsoforms[1]
    unexpr_dict = GenesWithNIsoforms[2]
    MaxNumIsos = GenesWithNIsoforms[3]

    #Get ranking and dropout info for each gene
    rankingDict = Dict()
    pObsDict = Dict()
    pDropoutDict = Dict()
    genesList = []

    for (key, value) in workingDict
        rankingDict[key] = FindRank(value)
        pObsDict[key] = FindPObs(value)
        pDropoutDict[key] = Float64[0.0, 0.0, 0.0, 0.0]
        push!(genesList, key)
    end

    plot_arr = []
    log_plot_arr = []
    mean_arr = []
    median_arr = []
    std_arr = []

    println("Beginning simulations!")

    for i in 1:NumIsoformsToSimulate

        #I pre-allocate arr in the hope it will make my code faster
        arr = Array{Float64}(undef, NumSimulations * length(genesList))
        arr_index = 1

        overlapArr = Array{Float64}(undef, NumSimulations * length(genesList))
        overlap_index = 1

        pDropoutArr = Array{Float64}(undef, NumSimulations * length(genesList) *
        (ncol(countsMatrix) - 1) * i)
        drop_index = 1

        #Need to add a for loop with num simulations here
        for j in 1:NumSimulations

            #@show arr_index

            #Get an array of mean isoforms per gene per cell
            simulation_out = globalArraySimulation(1, genesList, rankingDict, pObsDict,
            pDropoutDict, ncol(countsMatrix) - 1, i, filteringThreshold,
            unexpr_dict, MaxNumIsos, error_model, choice_model)

            arr[arr_index:(arr_index + length(genesList) - 1)] = simulation_out[1]
            overlapArr[overlap_index:(overlap_index + length(genesList) - 1)] = simulation_out[2]
            pDropoutArr[drop_index:(drop_index + length(simulation_out[3]) - 1)] = simulation_out[3]

            #Increment arr_index
            arr_index += length(genesList)
            overlap_index += length(genesList)
            drop_index += length(simulation_out[3])
        end

        #Save histogram of data
        if i == 1
            p1 = makeNormalisedHistogram(arr, string(i) * " Isoform")
        else
            p1 = makeNormalisedHistogram(arr, string(i) * " Isoforms")
        end
        push!(plot_arr, p1[1])
        push!(log_plot_arr, p1[2])
        push!(mean_arr, mean(arr))
        push!(std_arr, std(arr))
        push!(median_arr, median(arr))

        #Write to file if required
        if path != "none"
            writedlm(path * "_" * string(i) * "_isoforms_array.txt", arr)
            writedlm(path * "_" * string(i) * "_isoforms_overlap_array.txt", overlapArr)
            writedlm(path * "_" * string(i) * "_isoforms_pdrop_array.txt", pDropoutArr)
        end
    end

    #Get an array and a dict of mean isoforms per gene per cell in real data
    println("Finding means in real data.")
    findRealMeans = MeanNumIsoformsPerGenePerCellReal(workingDict)
    realMeansArr = findRealMeans[1]

    #Write to file if required
    if path != "none"
        writedlm(path * "_" * "real_isoforms_array.txt", realMeansArr)
    end

    #Plot real data
    real_plot = makeNormalisedHistogramReal(realMeansArr, "Real")
    push!(plot_arr, real_plot[1])
    push!(log_plot_arr, real_plot[2])

    #Arrange plots and write to output
    graph = plot(plot_arr..., layout=(length(plot_arr),1),legend=false, size = (300, 600), bottom_margin = 1.0mm)
    savefig(graph, output)
    graph = plot(log_plot_arr..., layout=(length(plot_arr),1),legend=false, size = (300, 600))
    savefig(graph, output * "_log.png")

    return (mean_arr, std_arr, median_arr)

end

function oneGenePrediction(countsMatrixPath::String, GeneIsoTablePath::String,
    gene::String, NumSimulations::Int64, output::String, test::Bool)

    global pFN = 0.04
    global pFP = 0.01

    #Read in counts matrix
    println("Data wrangling commencing!")
    println("Reading in counts matrix.")
    countsMatrix = ReadCountsMatrix(countsMatrixPath)

    #Read in gene-isoform table
    println("Reading in gene-isoform relationships")
    GeneIsoformRelationships = ReadGeneIsoformRelationships(GeneIsoTablePath)

    #Make sure gene-isoform table matches counts matrix
    if CheckGenesMatchIsoforms(countsMatrix, GeneIsoformRelationships) == false
        throw(ArgumentError("Your counts matrix doesn't match your list of genes
        and isoforms - there is something wrong with your input files."))
    end

    #Add gene information to isoforms table
    println("Adding Gene information to countsMatrix")
    countsMatrixWithGenes = AddGeneInfo(countsMatrix, GeneIsoformRelationships)

    #Make a dictionary just for the expressed isoforms in gene - no point
    #simulating unexpressed isoforms.
    tmp = filter(r -> any(occursin.([gene], r.Genes)), countsMatrixWithGenes)
    NumIsos = nrow(tmp)  #Get number of isoforms (incl. unexpressed) for gene of interest
    countsDict = ConvertExpressedIsosToDict(tmp)

    #Clean data for finding Km
    println("Cleaning data for finding Km")
    MMdf = CleanDataForMM(countsMatrix)

    #Find Km
    println("Calculating Km....")
    MMParams = FindKm(MMdf)
    Km = MMParams[1]
    sigma = MMParams[2]

    println("Preparing for simulations...")

    #Find ranking info
    cumulativeFreqs = GetIsoformChoiceProbabilitiesGene(length(countsDict[gene]))

    meansDict = Dict()

    #For each gene - bit unnecessary since just 1 gene, but convenient way of
    #getting value.
    for (key, value) in countsDict

        #Rank isoforms
        rankArr::Array{Int64, 1} = FindRank(value)

        #Find p(dropouts)
        dropoutArr::Array{Float64, 1} = FindPDropout(value, Km)

        for i in 1:length(countsDict[gene])

            meansArr = []

            for j in 1:NumSimulations
                mean = Simulate(cumulativeFreqs, rankArr, dropoutArr,
                ncol(countsMatrix) - 1, i, NumIsos, test)
                push!(meansArr, mean)
            end

            meansDict[i] = meansArr

        end
    end

    #Make a plot and write it to file
    OneGeneHistPlot(meansDict, tmp, output)

    #Find and return p-values
    mean = MeanIsoPerGenePerCell(tmp)
    return FindpValues(meansDict, mean, 40) #40 is an arbitrary number of required unique values

end

#Function that performs global level predictions
function rankingPredict(countsMatrixPath::String, GeneIsoTablePath::String,
     NumSimulations::Int64, output::String, NumIsoformsToSimulate::Int64,
     filteringThreshold::Int64, test::Bool, alpha, beta, pFNeg::Float64,
     pFPos::Float64, error_model::String, choice_model::String, path::String)

    #Set probability(FP|N) and probability(FN|P)
    if error_model == "none"
        global pFN = pFNeg
        global pFP = pFPos
    end

    #Read in counts matrix
    println("Data wrangling commencing!")
    println("Reading in counts matrix.") #Quick
    countsMatrix = ReadCountsMatrix(countsMatrixPath)

    #Read in gene-isoform table
    println("Reading in gene-isoform relationships")
    GeneIsoformRelationships = ReadGeneIsoformRelationships(GeneIsoTablePath)

    #Make sure gene-isoform table matches counts matrix
    if CheckGenesMatchIsoforms(countsMatrix, GeneIsoformRelationships) == false
        throw(ArgumentError("Your counts matrix doesn't match your list of genes
        and isoforms - there is something wrong with your input files."))
    end

    #Add gene information to isoforms table
    println("Adding Gene information to countsMatrix")
    countsMatrixWithGenes = AddGeneInfo(countsMatrix, GeneIsoformRelationships)

    #Only find Km if no Beta distribution parameters provided
    if isnothing(alpha) || isnothing(beta)

        #Clean data for finding Km
        println("Cleaning data for finding Km")
        MMdf = CleanDataForMM(countsMatrix)

        #Find Km
        println("Calculating Km....")
        MMParams = FindKm(MMdf)
        Km = MMParams[1]
        sigma = MMParams[2]
    end

    #Extract genes with filteringThreshold expressed isoforms
    println("Converting counts to dict")
    GenesWithNIsoforms = KeepGenesWithNIsoforms(countsMatrixWithGenes,
    filteringThreshold)
    workingDict = GenesWithNIsoforms[1]
    unexpr_dict = GenesWithNIsoforms[2]
    MaxNumIsos = GenesWithNIsoforms[3]

    #Get ranking and dropout info for each gene
    rankingDict = Dict()
    pObsDict = Dict()
    pDropoutDict = Dict()
    genesList = []

    #How we calculate dropouts depends on whether or not beta parameters are provided
    if isnothing(alpha) || isnothing(beta)
        for (key, value) in workingDict
            rankingDict[key] = FindRank(value)
            pObsDict[key] = FindPObs(value)
            pDropoutDict[key] = FindPDropout(value, Km)
            push!(genesList, key)
        end
    else
        for (key, value) in workingDict
            rankingDict[key] = FindRank(value)
            pObsDict[key] = FindPObs(value)
            pDropoutDict[key] = rand(Beta(alpha, beta), length(value))
            push!(genesList, key)
        end
    end

    plot_arr = []
    log_plot_arr = []
    mean_arr = []
    median_arr = []
    std_arr = []

    println("Beginning simulations!")

    for i in 4:NumIsoformsToSimulate #Only interested in 4 isoform sims

        #I pre-allocate arr in the hope it will make my code faster
        arr = []
        arr_index = 1

        std_arr = []
        std_arr_index = 1

        #Need to add a for loop with num simulations here
        for j in 1:NumSimulations

            #@show arr_index

            #Get an array of mean isoforms per gene per cell
            simulation_out = rankingSimulation(1, genesList, rankingDict, pObsDict,
            pDropoutDict, ncol(countsMatrix) - 1, i, filteringThreshold,
            unexpr_dict, MaxNumIsos, error_model, choice_model)

            #@show simulation_out

            append!(arr, simulation_out)
            #append!(std_arr, simulation_out[2])

            #arr[arr_index:(arr_index + length(simulation_out[1]) - 1)] = simulation_out[1]
            #Increment arr_index
            #arr_index += length(genesList)

            #std_arr[std_arr_index:(std_arr_index + length(simulation_out[2]) - 1)] = simulation_out[2]
            #Increment arr_index
            #std_arr_index += length(genesList)

        end

        if path != "none"
            writedlm(path * "_spearmans_array.txt", arr)
            writedlm(path * "_std_array.txt", std_arr)
        end
    end

    return ("Simulations complete")

end

#Function that performs global level predictions
function rankingOriginal(countsMatrixPath::String, GeneIsoTablePath::String,
     NumSimulations::Int64, output::String, NumIsoformsToSimulate::Int64,
     filteringThreshold::Int64, test::Bool, alpha, beta, pFNeg::Float64,
     pFPos::Float64, error_model::String, choice_model::String, path::String)

    #Set probability(FP|N) and probability(FN|P)
    if error_model == "none"
        global pFN = pFNeg
        global pFP = pFPos
    end

    #Read in counts matrix
    println("Data wrangling commencing!")
    println("Reading in counts matrix.") #Quick
    countsMatrix = ReadCountsMatrix(countsMatrixPath)

    #Read in gene-isoform table
    println("Reading in gene-isoform relationships")
    GeneIsoformRelationships = ReadGeneIsoformRelationships(GeneIsoTablePath)

    #Make sure gene-isoform table matches counts matrix
    if CheckGenesMatchIsoforms(countsMatrix, GeneIsoformRelationships) == false
        throw(ArgumentError("Your counts matrix doesn't match your list of genes
        and isoforms - there is something wrong with your input files."))
    end

    #Add gene information to isoforms table
    println("Adding Gene information to countsMatrix")
    countsMatrixWithGenes = AddGeneInfo(countsMatrix, GeneIsoformRelationships)

    #Only find Km if no Beta distribution parameters provided
    if isnothing(alpha) || isnothing(beta)

        #Clean data for finding Km
        println("Cleaning data for finding Km")
        MMdf = CleanDataForMM(countsMatrix)

        #Find Km
        println("Calculating Km....")
        MMParams = FindKm(MMdf)
        Km = MMParams[1]
        sigma = MMParams[2]
    end

    #Extract genes with filteringThreshold expressed isoforms
    println("Converting counts to dict")
    GenesWithNIsoforms = KeepGenesWithNIsoforms(countsMatrixWithGenes,
    filteringThreshold)
    workingDict = GenesWithNIsoforms[1]
    unexpr_dict = GenesWithNIsoforms[2]
    MaxNumIsos = GenesWithNIsoforms[3]

    #Get ranking and dropout info for each gene
    rankingDict = Dict()
    pObsDict = Dict()
    pDropoutDict = Dict()
    genesList = []
    spearmans_arr = []
    pDropout_Rank_cor = []
    std_arr = []

    for (key, value) in workingDict
        Arr1 = FindRank(value)
        Arr2 = FindRankbyNoCells(value)
        Arr3 = ordinalrank(FindPDropout(value, Km))
        push!(std_arr, FindStd(value))

        push!(spearmans_arr,
        corspearman(convert(Array{Float64, 1}, Arr1),
        convert(Array{Float64, 1}, Arr2)))

        push!(pDropout_Rank_cor,
        corspearman(convert(Array{Float64, 1}, Arr1),
        convert(Array{Float64, 1}, Arr3)))

    end

    writedlm(path * "_spearmans_array.txt", spearmans_arr)
    writedlm(path * "_pdropout_rank_cor_array.txt", pDropout_Rank_cor)
    writedlm(path * "_std.txt", std_arr)

end

export globalPredict
export oneGenePrediction
export NegativeControl
export NegativeControlDrop
export rankingPredict
export rankingOriginal

end #end of module

#GPU structure of kernel
#find p(dropout)
#for i in 1:10000
#simulation process
#save average
#return an array of averages

max_stars_repo_name: jenni-westoby/Obstacles

max_issues_repo_name: jenni-westoby/Obstacles

max_forks_repo_name: jenni-westoby/Obstacles


text:
import NearlyNewton: minimize!
function theta(x)
    if x[1] > 0
        return atan(x[2] / x[1]) / (2.0 * pi)
    else
        return (pi + atan(x[2] / x[1])) / (2.0 * pi)
    end
end

function fletcher_powell_fg!(∇f, x)
    theta_x = theta(x)

    if !(∇f==nothing)
        if ( x[1]^2 + x[2]^2 == 0 )
            dtdx1 = 0;
            dtdx2 = 0;
        else
            dtdx1 = - x[2] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
            dtdx2 =   x[1] / ( 2 * pi * ( x[1]^2 + x[2]^2 ) );
        end
        ∇f[1] = -2000.0*(x[3]-10.0*theta_x)*dtdx1 +
            200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[1]/sqrt( x[1]^2+x[2]^2 );
        ∇f[2] = -2000.0*(x[3]-10.0*theta_x)*dtdx2 +
            200.0*(sqrt(x[1]^2+x[2]^2)-1)*x[2]/sqrt( x[1]^2+x[2]^2 );
        ∇f[3] =  200.0*(x[3]-10.0*theta_x) + 2.0*x[3];
    end

    fx = 100.0 * ((x[3] - 10.0 * theta_x)^2 + (sqrt(x[1]^2 + x[2]^2) - 1.0)^2) + x[3]^2

    return fx
end

tp_fletch_powell_fg! = TestProblem(fletcher_powell_fg!, [-1.0, 0.0, 0.0], I, NearlyNewton.OptOptions())
tp_fletch_powell_fg!_alt = TestProblem(fletcher_powell_fg!, [-0.5, 0.0, 0.0], I, NearlyNewton.OptOptions())

minimize!(tp_fletch_powell_fg!, BFGS(InverseApprox()))
minimize!(tp_fletch_powell_fg!, BFGS(DirectApprox()))

minimize!(tp_fletch_powell_fg!, SR1(InverseApprox()))
minimize!(tp_fletch_powell_fg!, SR1(DirectApprox()))

minimize!(tp_fletch_powell_fg!, DFP(InverseApprox()))
minimize!(tp_fletch_powell_fg!, DFP(DirectApprox()))

minimize!(tp_fletch_powell_fg!_alt, BFGS(InverseApprox()))
minimize!(tp_fletch_powell_fg!_alt, BFGS(DirectApprox()))

minimize!(tp_fletch_powell_fg!_alt, SR1(InverseApprox()))
minimize!(tp_fletch_powell_fg!_alt, SR1(DirectApprox()))

minimize!(tp_fletch_powell_fg!_alt, DFP(InverseApprox()))
minimize!(tp_fletch_powell_fg!_alt, DFP(DirectApprox()))

ext: jl


text:
"""
    NonstationaryArmBase

`NonstationaryArmBase` is the concrete type for all nonstationary arm models. This
includes rested, restless and evolving arm models. Any subtype should implement
`pull!()`, `tick!()` and `reset!()` methods
"""
abstract type NonstationaryArmBase <: ArmBase end

"""
    Pulse( period, changePoint, holdDuration [,isRestless] )

`Pulse` creates a pulse reward arm which switches between 0 and 1. `period`(Int) is
the period of the pulse wave, `changePoint` (Int) is the point at which the pulse
should switch from 0 to 1 from the start of period. `holdDuration` (Int) is the time
for which the pulse should stay 1. `isRestless` (bool) sets whether the arm is
restless or not.
"""
mutable struct Pulse <: NonstationaryArmBase
    step::Int             # For tracking internally
    period::Int           # Period of the pulse
    changePoint::Int      # The point at which the signal goes high
    highDuration::Int     # Time to stay HIGH
    isRestless::Bool        # Is this a restless arm

    function Pulse( _period, _cp, _hd, _ir = true )
        new(
            0,
            _period,
            _cp,
            _hd,
            _ir
        )
    end
end

function pull!( arm::Pulse )
    arm.step    = arm.step + 1
    if arm.step > arm.period
        arm.step = 1
    end

    return ((arm.step>=arm.changePoint)&&(arm.step<=arm.changePoint+arm.highDuration))  ? 1 : 0;
end

function tick!( arm::Pulse )
    if arm.isRestless
        arm.step = arm.step + 1
        if arm.step > arm.period
            arm.step = 1
        end
    end
    nothing
end

function reset!( arm::Pulse )
    arm.step    = 0
    nothing
end

"""
    Sinusoidal( period[, isRestless = true [, offset = 2*π*rand()]] )

`Sinusoidal` creates an arm which sinusoidally changes between 0 and 1 with
`period`. `isRestless` (bool) sets whether the arm is restless or not. `offset`
sets the offset for the arm (set to a random offset by default).
"""
mutable struct Sinusoidal <: NonstationaryArmBase
    step::Int         # Time step
    period::Int       # Period of the sinusoidal wave
    isRestless::Bool    # Is this arm restless - True by default
    offset::Float64     # Initial pahse, random by default

    function Sinusoidal( p::Int; isRestless::Bool=true, offset::Float64=2*π*rand() )
        new(
            0,
            p,
            isRestless,
            offset
        )
    end
end

function pull!( arm::Sinusoidal )
    arm.step = arm.step + 1
    return 1/2 + 1/2 * sin( 2 * π * (arm.step-1)/arm.period + arm.offset )  # To limit rewards between 0 and 1
end

function tick!( arm::Sinusoidal )
    if arm.isRestless
        arm.step = arm.step + 1
    end
    nothing
end

function reset!( arm::Sinusoidal )
    arm.step = 0
    nothing
end

"""
    Square( period::Int, chnagePoints::Dict{Int,Float64}[,isRestless=true] )

`Square` creates an arm which switches value at specified instances. Starts from
a default value of `0` and changes to values as per the dictionary of changePoints.
```julia-repl
arm1 = Arms.Square( 50, Dict(10=>0.5,25=>0.20,40=>1.0) )
```
"""
mutable struct Square <: NonstationaryArmBase
    step::Int                         # For internal tracking
    period::Int                       # Period at which wave repeats
    changePoints::Dict{Int,Float64}   # Dictionary of change points
    isRestless::Bool                    # Is this a restless arm

    function Square( _period::Int, _cps::Dict{Int,Float64}, _ir::Bool = true )
        new(
            0,
            _period,
            _cps,
            _ir
        )
    end
end

function pull!( arm::Square )
    arm.step    = arm.step + 1
    if arm.step > arm.period
        arm.step = 1
    end
    lastCpIndex = 0.0;
    for timeIdx in keys(arm.changePoints)
        if timeIdx <= arm.step && timeIdx > lastCpIndex
            lastCpIndex = timeIdx
        end
    end
    if lastCpIndex == 0 # No change occured
        return 0.0
    else
        return arm.changePoints[lastCpIndex]
    end
end

function tick!( arm::Square )
    if arm.isRestless
        arm.step = arm.step + 1
        if arm.step > arm.period
            arm.step = 1
        end
    end
    nothing
end

function reset!( arm::Square )
    arm.step    = 0
    nothing
end

"""
    Variational( variation, period[, isRestless=true[,offset=2π*rand()]] )

`Variational` creates a nonstationary arm as mentioned in
    Besbes, O., Gur, Y., & Zeevi, A. (2014). Optimal Exploration-Exploitation in a Multi-Armed-Bandit Problem with Non-stationary Rewards, 1–20.
"""
mutable struct Variational <: NonstationaryArmBase
    step::Int64         # Time step
    variation::Float64  # Variation of the arms
    period::Float64     # Period of the sinusoidal wave
    isRestless::Bool    # Is this arm restless - True by default
    offset::Float64     # Initial phase, random by default

    function Variational( variation::Float64, period::Float64; isRestless::Bool=true, offset::Float64=2*π*rand() )
        new( 0, variation, period, isRestless, offset )
    end
end

function pull!( arm::Variational )
    arm.step = arm.step + 1
    return 1/2 + 1/2 * sin( arm.variation * π * (arm.step-1)/arm.period + arm.offset )  # To limit rewards between 0 and 1
end

function tick!( arm::Variational )
    if arm.isRestless
        arm.step = arm.step + 1
    end
    nothing
end

function reset!( arm::Variational )
    arm.step = 0
    nothing
end


text:
@testset "mean_functions" begin
    @testset "ZeroMean" begin
        P = 3
        Q = 2
        D = 4
        # X = ColVecs(randn(rng, D, P))
        x = randn(P)
        x̄ = randn(P)
        f = ZeroMean{Float64}()

        for x in [x]
            @test AbstractGPs._map(f, x) == zeros(size(x))
            # differentiable_mean_function_tests(f, randn(rng, P), x)
        end

        # Manually verify the ChainRule. Really, this should employ FiniteDifferences, but
        # currently ChainRulesTestUtils isn't up to handling this, so this will have to do
        # for now.
        y, pb = rrule(AbstractGPs._map, f, x)
        @test y == AbstractGPs._map(f, x)
        Δmap, Δf, Δx = pb(randn(P))
        @test iszero(Δmap)
        @test iszero(Δf)
        @test iszero(Δx)
    end
    @testset "ConstMean" begin
        rng, D, N = MersenneTwister(123456), 5, 3
        # X = ColVecs(randn(rng, D, N))
        x = randn(rng, N)
        c = randn(rng)
        m = ConstMean(c)

        for x in [x]
            @test AbstractGPs._map(m, x) == fill(c, N)
            # differentiable_mean_function_tests(m, randn(rng, N), x)
        end
    end
    @testset "CustomMean" begin
        rng, N, D = MersenneTwister(123456), 11, 2
        x = randn(rng, N)
        foo_mean = x -> sum(abs2, x)
        f = CustomMean(foo_mean)

        @test AbstractGPs._map(f, x) == map(foo_mean, x)
        # differentiable_mean_function_tests(f, randn(rng, N), x)
    end
end


text:
# Transpose is treated in a special way. There are two transpose operations, namely x' and trans(x)
# A' creates a special ADtrans node. The purpose of this is that, for example A'*B can then be recompiled later in to a GEMM call GEMM('T','N',A,B) -- this is efficient since the transpose of A is not actually computed or stored in the graph. These ADtrans nodes are deleted later after optimisation.
#
# On the other hand, trans(A) actually computes and stores the transpose of A on the graph. This should be used sparingly since this will, in general, be very wasteful of storage and computation.

# f(X)=transpose(X)
import Base.transpose
transpose(A::ADnode)=ADtrans(A)
export transpose

ftranspose(A::ADnode)=ADnode(Ftranspose,A)

export ftranpose

Ftranspose(X)=(X',nothing)

function Ftranspose_inplace(value::Array,auxvalue,X::Array)
    copy!(value,X')
end

function Dtranspose(derivativeIDX,f_c,faux_c,grad_c,grad_n,X::Array)
    axpy!(1.0,grad_c',grad_n)
end

if PROC=="GPU"
    function Ftranspose_inplace(value::CudaArray,auxvalue,X::CudaArray)
        tmp=CudaArray(Float64,size(X))
        CUBLAS.geam!('T','N',1.0,X,0.0,tmp,value) # value=X'
        free(tmp)
    end

    function Dtranspose(derivativeIDX,f_c,faux_c,grad_c,grad_n,X::CudaArray)
        tmp=CudaArray(Float64,size(grad_n))
        CUBLAS.geam!('T','N',1.0,grad_c,0.0,tmp,tmp) # tmp=grad_c'
        axpy!(1.0,tmp,grad_n); free(tmp)
    end

end

Derivative[Ftranspose]=Dtranspose
Inplace[Ftranspose]=Ftranspose_inplace

export Ftranspose

import Base.*

function *(A::ADtrans,B::ADnode)
#    splice!(node,A.index)
    node[A.index]=nothing
    return ADnode(FAtransposeX,[node[A.parent] B])
end

function *(A::ADtrans,B::Real)
    node[A.index]=nothing
    return ADnode(FAtransposeX,[node[A.parent] ADconst(B)])
end

function *(A::ADnode,B::ADtrans)
    node[B.index]=nothing
    return ADnode(FAXtranspose,[A node[B.parent]])
end

function *(A::Real,B::ADtrans)
    node[B.index]=nothing
    return ADnode(FAXtranspose,[ADconst(A) node[B.parent]])
end

function *(A::ADtrans,B::ADtrans)
    node[A.index]=nothing
    node[B.index]=nothing
    return ADnode(FAtransposeXtranspose,[node[A.parent] node[B.parent]])
end

export *

import Base.+

function +(A::ADtrans,B::ADnode)
    node[A.index]=nothing
    return ADnode(FXtransposePY,[node[A.parent] B])
end

function +(A::ADnode,B::ADtrans)
    node[B.index]=nothing
    return ADnode(FXPYtranspose,[A node[B.parent]])
end

function +(A::ADtrans,B::ADtrans)
    node[A.index]=nothing
    node[B.index]=nothing
    return ADnode(FXtransposePYtranspose,[node[A.parent] node[B.parent]])
end

export +

import Base.-

function -(A::ADtrans,B::ADnode)
    node[A.index]=nothing
    return ADnode(FXtransposeMY,[node[A.parent] B])
end

function -(A::ADnode,B::ADtrans)
    node[B.index]=nothing
    return ADnode(FXMYtranspose,[A node[B.parent]])
end

function -(A::ADtrans,B::ADtrans)
    node[A.index]=nothing
    node[B.index]=nothing
    return ADnode(FXtransposeMYtranspose,[node[A.parent] node[B.parent]])
end

trans=ftranspose
export trans

export -

#TODO A'.*B, A.*B', A'.*B'

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
# This file simulates a closed system

using Causal

# Construct the model
@defmodel model begin
    @nodes begin
        gen = FunctionGenerator(readout=sin)
        adder = Adder(signs=(+,-))
        ds = ContinuousLinearSystem()
        writer = Writer(input=Inport(2))
    end
    @branches begin
        gen[1]      =>  adder[1]
        adder[1]    =>  ds[1]
        ds[1]       =>  adder[2]
        gen[1]      =>  writer[1]
        ds[1]       =>  writer[2]
    end
end

# Simualate the model
sim = simulate!(model, 0., 0.01, 10.)

# Read and plot data
t, x = read(getnode(model, :writer).component)
using Plots
plot(t, x[:, 1], label="r(t)", xlabel="t", lw=3)
plot!(t, x[:, 2], label="y(t)", xlabel="t", lw=3)
plot!(t, 6 / 5 * exp.(-2t) + 1 / 5 * (2 * sin.(t) - cos.(t)), label="Analytical Solution", lw=3)
# fileanme = "readme_example.svg"
# path = joinpath(@__DIR__, "../docs/src/assets/ReadMePlot/")
# savefig(joinpath(path, fileanme))


text:
insts = OrderedDict()
insts["minimal"] = [
    ((1, 1, 1, false, false, true),),
    ((1, 1, 1, false, true, false),),
    ((1, 1, 1, true, false, false),),
    ]
insts["fast"] = [
    ((3, 1, 3, false, false, true),),
    ((3, 1, 3, false, true, false),),
    ((3, 1, 3, true, false, false),),
    ((1, 2, 3, false, false, true),),
    ((1, 2, 3, false, true, false),),
    ((1, 2, 3, true, false, false),),
    ]
insts["various"] = [
    ((4, 2, 5, false, false, true),),
    ((4, 2, 5, false, true, false),),
    ((4, 2, 5, true, false, false),),
    ((4, 2, 10, false, true, false),),
    ((4, 2, 10, true, false, false),),
    ((4, 4, 5, false, true, false),),
    ((4, 3, 6, true, false, false),),
    ((4, 3, 6, false, true, false),),
    ]
return (NearestPolyMatJuMP, insts)


text:
pkgs = Pkg.installed()
!haskey(pkgs, "RLESUtils") && Pkg.clone("https://github.com/rcnlee/RLESUtils.jl.git", "RLESUtils")
!haskey(pkgs, "NECPP") && Pkg.clone("https://github.com/rcnlee/NECPP.jl.git", "NECPP")


text:
"""
```julia
LaplaceLikelihood(β::T=1.0)  #  Laplace likelihood with scale β
```

Laplace likelihood for regression:
```math
1/(2β) exp(-|y-f|/β)
```
see [wiki page](https://en.wikipedia.org/wiki/Laplace_distribution)
---
For the analytical solution, it is augmented via:
```math
p(y|f,ω) = N(y|f,ω⁻¹)
```
where ``ω ~ Exp(ω | 1/(2 β^2))``, and `Exp` is the [Exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution)
We use the variational distribution ``q(ω) = GIG(ω | a,b,p)``
"""
mutable struct LaplaceLikelihood{T<:Real} <: RegressionLikelihood{T}
    β::T
    a::T
    p::T
    b::Vector{T} #Variational parameter b of GIG
    θ::Vector{T} #Expected value of ω
    function LaplaceLikelihood{T}(β::T) where {T<:Real}
        new{T}(β,β^-2,0.5)
    end
    function LaplaceLikelihood{T}(β::T,b::AbstractVector{T},θ::AbstractVector{T}) where {T<:Real}
        new{T}(β,β^(-2),0.5,b,θ)
    end
end

function LaplaceLikelihood(β::T=1.0) where {T<:Real}
    LaplaceLikelihood{T}(β)
end

function init_likelihood(likelihood::LaplaceLikelihood{T},inference::Inference{T},nLatent::Int,nSamplesUsed::Int,nFeatures::Int) where T
    if inference isa AnalyticVI || inference isa GibbsSampling
        LaplaceLikelihood{T}(
        likelihood.β,
        rand(T,nSamplesUsed),
        zeros(T,nSamplesUsed))
    else
        LaplaceLikelihood{T}(likelihood.β)
    end
end

function pdf(l::LaplaceLikelihood,y::Real,f::Real)
    Distributions.pdf(Laplace(f,l.β),y) #WARNING multioutput invalid
end

function logpdf(l::LaplaceLikelihood,y::Real,f::Real)
    Distributions.logpdf(Laplace(f,l.β),y) #WARNING multioutput invalid
end

function Base.show(io::IO,model::LaplaceLikelihood{T}) where T
    print(io,"Laplace likelihood")
end

function compute_proba(l::LaplaceLikelihood{T},μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
    return μ,max.(σ²,0.0).+ 2*l.β^2
end

## Local Updates ##

function local_updates!(l::LaplaceLikelihood{T},y::AbstractVector,μ::AbstractVector,Σ::AbstractVector) where {T}
    l.b .= Σ+abs2.(μ-y)
    l.θ .= sqrt(l.a)./sqrt.(l.b)
end

function sample_local!(l::LaplaceLikelihood,y::AbstractVector,f::AbstractVector) where {T}
    l.b .= rand.(GeneralizedInverseGaussian.(1/l.β^2,abs2.(f-y),0.5))
    set_ω!(l,inv.(l.b))
    return nothing
end

@inline ∇E_μ(l::LaplaceLikelihood{T},::AOptimizer,y::AbstractVector) where {T} =  (l.θ.*y,)
@inline ∇E_Σ(l::LaplaceLikelihood{T},::AOptimizer,y::AbstractVector) where {T} = (0.5*l.θ,)

## ELBO ##
function expec_log_likelihood(l::LaplaceLikelihood{T},i::AnalyticVI,y::AbstractVector,μ::AbstractVector,diag_cov::AbstractVector) where {T}
    tot = -0.5*length(y)*log(twoπ)
    tot += 0.5*sum(log,l.θ)
    tot += -0.5*(dot(l.θ,diag_cov)+dot(l.θ,abs2.(μ))-2.0*dot(l.θ,μ.*y)+dot(l.θ,abs2.(y)))
    return tot
end

AugmentedKL(l::LaplaceLikelihood,::AbstractVector) = GIGEntropy(l)-expecExponentialGIG(l)

GIGEntropy(l::LaplaceLikelihood{T}) where {T} = GIGEntropy(l.a,l.b,l.p)

function expecExponentialGIG(l::LaplaceLikelihood{T}) where {T}
    sum(-log(2*l.β^2).-0.5*(l.a.*sqrt.(l.b)+l.b.*sqrt(l.a))./(l.a.*l.b*l.β^2))
end

## PDF and Log PDF Gradients ##

function grad_quad(likelihood::LaplaceLikelihood{T},y::Real,μ::Real,σ²::Real,inference::Inference) where {T<:Real}
    nodes = inference.nodes*sqrt(σ²) .+ μ
    Edlogpdf = dot(inference.weights,grad_log_pdf.(likelihood,y,nodes))
    Ed²logpdf =  (1/sqrt(twoπ*σ²))/(likelihood.β^2)
    return -Edlogpdf::T, Ed²logpdf::T
end

@inline grad_log_pdf(l::LaplaceLikelihood{T},y::Real,f::Real) where {T<:Real} = sign(y-f)./l.β

@inline hessian_log_pdf(l::LaplaceLikelihood{T},y::Real,f::Real) where {T<:Real} = zero(T)


text:
using Hamburg
using DataFrames
using Test

include("Covid19.jl")

@testset "using `dataset()` to load datasets" begin
    @testset "loads existing datasets into a DataFrame" begin
        df = dataset("covid-19", "infected")
        @test isa(df, DataFrame)
        @test size(df)[1] > 0
    end

    @testset "throws an error when trying to load a non-existent dataset" begin
        @test_throws ArgumentError dataset("moo", "foo")
    end
end


text:
"""
    randomamplitudes(ts::AbstractArray{T, 1} where T)

Create a random amplitude surrogate for `ts`.

A modification of the random phases surrogates from
[Theiler et al., 1992](https://www.sciencedirect.com/science/article/pii/016727899290102S)
, where amplitudes are adjusted instead of the phases after taking the Fourier
transform.

## Arguments
- ts`**: the time series for which to generate the surrogate realization.

## References

J. Theiler et al., Physica D *58* (1992) 77-94 (1992).
[https://www.sciencedirect.com/science/article/pii/016727899290102S](https://www.sciencedirect.com/science/article/pii/016727899290102S)

"""
function randomamplitudes(ts::AbstractArray{T, 1} where T)
    n = length(ts)

    # Fourier transform
    ft = fft(ts)

    # Polar coordinate representation of the Fourier transform
    r = abs.(ft)    # amplitudes
    ϕ = angle.(ft)  # phase angles

    # To randomise phases, we multiply each complex amplitude by e^(iϕ), with ϕ picked
    # on the interval [0, 2π]. Here, we pick according to a uniform distribution.
    randomised_amplitudes = r .* rand(Uniform(0, 2*pi), n)
    real.(ifft(randomised_amplitudes .* exp.(ϕ .* 1im)))
end


text:
md"""
# What is JSServe?

* a connection between Julia & Javascript
* a reactive DOM
* a server with routing
* can be used for dashboards, wrapping JS libraries & writing webpages
* Similar to React.js in some ways
"""

using JSServe, Observables
using JSServe: @js_str, Session, App, onjs, onload, Button
using JSServe: TextField, Slider, linkjs, get_server
using JSServe.DOM
using JSServe: JSON3
JSServe.browser_display()

app = App(DOM.h1("Hello World"))
display(app)

md"""
# Interaction with observables
"""

color = Observable("red")
color_css = map(x-> "color: $(x)", color)

app = App() do
    return DOM.h1("Hello World", style=map(x-> "color: $(x)", color))
end
display(app)

color[] = "green"

app = App() do
    button = DOM.div("click me", onclick=js"JSServe.update_obs($(color), 'blue')")
    return DOM.div(
        button, DOM.h1("Hello World", style=map(x-> "color: $(x)", color))
    )
end
display(app)

md"""
# Other ways to execute Javascript
"""

app = App() do session::Session
    button = DOM.div("click me", onclick=js"JSServe.update_obs($(color), 'blue')")
    onload(session, button, js"""function load(button){
        window.alert('Hi from JavaScript');
    }""")

    onjs(session, color, js"""function update(value){
        window.alert(value);
        // throw "heey!"
    }""")

    return DOM.div(
        button, DOM.h1("Hello World", style=map(x-> "color: $(x)", color))
    )
end
display(app)

md"""
# Including assets & Widgets
"""

MUI = JSServe.Asset("//cdn.muicss.com/mui-0.10.1/css/mui.min.css")
sliderstyle = JSServe.Asset(joinpath(@__DIR__, "sliderstyle.css"))
image = JSServe.Asset(joinpath(@__DIR__, "assets", "julia.png"))

app = App() do
    button = JSServe.Button("hi", class="mui-btn mui-btn--primary")
    slider = JSServe.Slider(1:10, class="slider")

    on(button) do click
        @show click
    end

    on(slider) do slidervalue
        @show slidervalue
    end
    link = DOM.a(href="/example1", "GO TO ANOTHER WORLD")
    return DOM.div(MUI, sliderstyle, link, button, slider, DOM.img(src=image))
end
display(app)

md"""
# Can also use regex here!
# ctx contains route, match,
"""

app = App() do
    button = DOM.div("click me", onclick=js"JSServe.update_obs($(color), 'blue')")
    return DOM.div(
        button, DOM.h1("Hello World", style=map(x-> "color: $(x)", color))
    )
end
display(app)

JSServe.route!(get_server(), "/example1" => app)

app = App() do session::Session
    slider = JSServe.Slider(1:10, class="slider m-4")
    squared = map(slider) do slidervalue
        return slidervalue^2
    end
    class = "p-2 rounded border-2 border-gray-600 m-4"
    v1 = DOM.div(slider.value, class=class)
    v2 = DOM.div(squared, class=class)
    dom = DOM.div(JSServe.TailwindCSS, slider, sliderstyle, v1, v2)
    # statemap for static serving
    return JSServe.record_state_map(session, (s, r)-> dom)
end

export_path = joinpath(@__DIR__, "demo")
mkdir(export_path)
JSServe.export_standalone(app, export_path, clear_folder=true)

using LiveServer
cd(export_path)
LiveServer.serve()


text:
#!/bin/dev/env julia
module MazeSolver
#Pkg.add("DataStructures")
#Pkg.add("Images")
using Images
using DataStructures
using BenchmarkTools

function getstart(img, h, w)
	for j = 1:w
		if img[1, j].r > 0.5
			return (1, j)
		end
	end
	throw(DomainError(img, "maze must have a white pixel in first row"))
end

function getend(img, h, w)
	for j = 1:w
		if img[h, j].r > 0.5
			return (h, j)
		end
	end
	throw(DomainError(img, "maze must have a white pixel in last row"))
end

function checkNeighbors(maze_map, index, h, w)
	neighbors = Char[]
	if index[1] > 1 && maze_map[index[1]-1, index[2]] != false
		push!(neighbors, 'N')
	end
	if index[2] < w && maze_map[index[1], index[2]+1] != false
		push!(neighbors, 'E')
	end
	if index[1] < h && maze_map[index[1]+1, index[2]] != false
		push!(neighbors, 'S')
	end
	if index[2] > 1 && maze_map[index[1], index[2]-1] != false
		push!(neighbors, 'W')
	end
	num = length(neighbors)
	return num, neighbors
end

function checkChar(char, array)
	length(filter(x -> x == char, array)) >= 1
end

function getNeighbors(links, maze_map, index, h, w)
	neighbors = Array{Tuple{Int, Int}, 1}()
	if checkChar('N', links)
		offset = 1
		while index[1] - offset > 0
			if maze_map[index[1] - offset, index[2]] isa Bool &&
					maze_map[index[1] - offset, index[2]]
				offset += 1
			else
				push!(neighbors, (index[1] - offset, index[2]))
				break
			end
		end
	end
	if checkChar('E', links)
		offset = 1
		while index[2] + offset <= w
			if maze_map[index[1], index[2] + offset] isa Bool &&
					maze_map[index[1], index[2] + offset]
				offset += 1
			else
				push!(neighbors, (index[1], index[2] + offset))
				break
			end
		end
	end
	if checkChar('S', links)
		offset = 1
		while index[1] + offset <= h
			if maze_map[index[1] + offset, index[2]] isa Bool &&
					maze_map[index[1] + offset, index[2]]
				offset += 1
			else
				push!(neighbors, (index[1] + offset, index[2]))
				break
			end
		end
	end
	if checkChar('W', links)
		offset = 1
		while index[2] - offset > 0
			if maze_map[index[1], index[2] - offset] isa Bool &&
					maze_map[index[1], index[2] - offset]
				offset += 1
			else
				push!(neighbors, (index[1], index[2] - offset))
				break
			end
		end
	end
	return neighbors
end

function analyse_maze(img)
	h, w = size(img)
	maze_map = Matrix{Union{UInt8, Bool, Array{Tuple{Int, Int}, 1}}}(undef, h, w)
	Startnode = getstart(img, h, w)
	Endnode = getend(img, h, w)
	counter = 0
	for j in CartesianIndices(img)
		if img[j].r <= 0.5
			maze_map[j] = false
		elseif img[j].r > 0.5
			maze_map[j] = true
			counter += 1
		end
	end
	println("	number of path tiles: " * string(counter))
	maze_map = reshape(maze_map, (h, w))
	maze_map[Startnode[1], Startnode[2]] = UInt8(5)
	maze_map[Endnode[1], Endnode[2]] = UInt8(6)
	nodes = [Startnode, Endnode]
	for i in CartesianIndices((2:h-1, 2:w-1))
		if maze_map[i[1], i[2]] != false
			num, neighbors = checkNeighbors(maze_map, (i[1], i[2]), h, w)
			if length(neighbors) == 2 && (('N' in neighbors && 'S' in neighbors)
				|| ('E' in neighbors && 'W' in neighbors))
				maze_map[i[1], i[2]] = true
			elseif num >= 1
				maze_map[i[1], i[2]] = UInt8(num + 1)
				push!(nodes, ((i[1], i[2])))
			elseif num == 0
				maze_map[i[1], i[2]] = false
			end
		end
	end
	println("	number of nodes: " * string(length(nodes)))
	for a in nodes
		neighbors = getNeighbors(checkNeighbors(maze_map, a, h, w)[2],
			maze_map, a, h, w)
		maze_map[a[1], a[2]] = neighbors
	end
	print("	took ")
	return maze_map, nodes, Startnode, Endnode
end

function a_star(maze_map, startnode::Tuple{Int, Int}, endnode::Tuple{Int, Int})
	costs = Dict{Tuple{Int, Int}, Int}()
	dict = Dict{Tuple{Int, Int},Tuple{Int, Int}}()
	openlist = PriorityQueue{Tuple{Int, Int}, Int}()
	closedlist = Set{Tuple{Int, Int}}()

	costs[startnode] = 0
	enqueue!(openlist, startnode, 0)

	while !(isempty(openlist))
		current = dequeue!(openlist)
		if current == endnode
			return unfoldpath(endnode, startnode, dict)
		end
		push!(closedlist, current)
		for edge in maze_map[current[1], current[2]]
			if !(edge in closedlist)
				tentative_g = costs[current] + cost(edge, current)
				if edge in keys(openlist) && tentative_g >= costs[edge]
					continue
				end
				dict[edge] = current
				costs[edge] = tentative_g

				f = tentative_g + heuristic(edge, endnode)
				if edge in keys(openlist)
					openlist[edge] = f
				else
					enqueue!(openlist, edge, f)
				end
			end
		end
	end
end

function cost(a, b)
	return abs(b[1] - a[1]) + abs(b[2] - a[2])
end

function heuristic(a, b)
	return floor(hypot((b[2] - a[2]), (b[1] - a[1])))
end

function bfs(maze_map, startnode::Tuple{Int, Int},
			endnode::Tuple{Int, Int})
	dict = Dict{Tuple{Int, Int},Tuple{Int, Int}}()
	visited = Set{Tuple{Int, Int}}()
	push!(visited, startnode)
	q = Queue{Tuple{Int, Int}}()
	dict[startnode] = startnode
	enqueue!(q, startnode)
	while !isempty(q)
		vertex = dequeue!(q)
		if vertex == endnode
			return unfoldpath(endnode, startnode, dict)
		end
		for edge in maze_map[vertex[1], vertex[2]]
			if !(edge in visited)
				dict[edge] = vertex
				push!(visited, edge)
				enqueue!(q, edge)
			end
		end
	end
end

function unfoldpath(lastnode, startnode, maze)
	path = Array{Tuple{Int, Int}, 1}()
	current = lastnode
	while current != startnode
		if ((current[1] - maze[current][1]) == 0)
			if current[2] - maze[current][2] > 0
				# due west
				offset = (0, -1)
			else
				# due east
				offset = (0, 1)
			end
		else
			if current[1] - maze[current][1] > 0
				# due north
				offset = (-1, 0)
			else
				# due south
				offset = (1, 0)
			end
		end
		factor = 0
		while (current[1] + (offset[1] * factor),
			current[2] + (offset[2] * factor)) != maze[current]
			push!(path, (current[1] + (offset[1] * factor),
				current[2] + (offset[2] * factor)))
			factor += 1
		end
		current = maze[current]
	end
	push!(path, current)
	return path
end

function colourpath(path, painted)
	number = 1
	for l in result
		painted[l[1], l[2]] = RGB(1.0 * (number / length(result)),
			1.0 - 1.0 * (number / length(result)),
			0.75 - 0.5 * (number / length(result)))
		number += 1
	end
	return painted
end

if (isfile(ARGS[1]))
	filepath = ARGS[1]
	(length(ARGS) == 2 && ARGS[2] == "true") ? weighted = true : weighted = false
else
	println("couldn't load " * string(ARGS[1]))
	exit()
end

img = RGB.(load(filepath))

println("analysing maze ... ")
@time maze_map, nodes, startnode, endnode = analyse_maze(img)
reshape(maze_map, size(img))
print(weighted ? "solved maze for shortest path in " : "solved maze in ")
@time result = weighted ? a_star(maze_map, startnode, endnode) : bfs(maze_map, startnode, endnode)
println("length of the solution: " * string(length(result)))

save(string(split(filepath, ".")[1] * "_solved.png"), colourpath(result, copy(img)))
println("stored maze at: " * string(split(filepath, ".")[1] * "_solved.png"))
end


text:
using Revise
using AdFem
using PyCall
using LinearAlgebra
np = pyimport("numpy")

# Domain information
NT = 50
Δt = 1/NT
n = 50
m = n
h = 1.0/n
bdnode = Int64[]
for i = 1:m+1
    for j = 1:n+1
        if i==1 || i==m+1 || j==1|| j==n+1
            push!(bdnode, (j-1)*(m+1)+i)
        end
    end
end

# Physical parameters
b = 1.0
H = [1.0 0.0 0.0
    0.0 1.0 0.0
    0.0 0.0 0.5]
Q = compute_fvm_tpfa_matrix(m, n, h)
K = compute_fem_stiffness_matrix(H, m, n, h)
L = compute_interaction_matrix(m, n, h)
M = compute_fvm_mass_matrix(m, n, h)
A = [K -b*L'
b*L/Δt 1/Δt*M-Q]
A, Abd = fem_impose_coupled_Dirichlet_boundary_condition(A, bdnode, m, n, h)

pexact = (x,y,t)->x^2*(1-x)^2*y^2*(1-y)^2*exp(-t)
uexact1 = (x,y,t)->(x^2+y^2)*t
uexact2 = (x,y,t)->(x^2-y^2)*t
g1 = (x,y,t)-> 3t + (-x^2*y^2*(1 - y)^2*(2*x - 2)*exp(-t) - 2*x*y^2*(1 - x)^2*(1 - y)^2*exp(-t))
g2 = (x,y,t)-> -t + (-x^2*y^2*(1 - x)^2*(2*y - 2)*exp(-t) - 2*x^2*y*(1 - x)^2*(1 - y)^2*exp(-t))
f = (x,y,t)->(-x^2*y^2*(x - 1)^2*(y - 1)^2 - 2*x^2*y^2*(x - 1)^2 - 2*x^2*y^2*(y - 1)^2 - 8*x^2*y*(x - 1)^2*(y - 1) - 2*x^2*(x - 1)^2*(y - 1)^2 - 8*x*y^2*(x - 1)*(y - 1)^2 - 2*y^2*(x - 1)^2*(y - 1)^2 + 2*(x - y)*exp(t))*exp(-t)

U = zeros(m*n+2(m+1)*(n+1), NT+1)
for i = 1:m
    for j = 1:n
        x = (i-1/2)*h; y = (j-1/2)*h
        U[2(m+1)*(n+1)+(j-1)*m+i, 1] = pexact(x,y,0.0)
    end
end

x = Float64[]; y = Float64[]
for j = 1:n+1
    for i = 1:m+1
        push!(x, (i-1)*h)
        push!(y, (j-1)*h)
    end
end

for i = 1:NT
    t = i*Δt
    @info i
    g1_ = eval_f_on_gauss_pts((x,y)->g1(x,y,t), m, n, h)
    g2_ = eval_f_on_gauss_pts((x,y)->g2(x,y,t), m, n, h)
    f_ = eval_f_on_gauss_pts((x,y)->f(x,y,t), m, n, h)
    rhs1 = -compute_fem_source_term(g1_, g2_, m, n, h)
    rhs2 = -compute_fvm_source_term(f_, m, n, h) +
                b*L*U[1:2(m+1)*(n+1), i]/Δt +
                M * U[2(m+1)*(n+1)+1:end,i]/Δt

    rhs = [rhs1;rhs2]

    bdval = [uexact1.(x[bdnode],y[bdnode],t); uexact2.(x[bdnode],y[bdnode],t)]
    rhs[[bdnode; bdnode.+ (m+1)*(n+1)]] = bdval
    rhs -= Abd * bdval

    U[:,i+1] = A\rhs
end

U0 = zeros(2*(m+1)*(n+1)+m*n, NT+1)
for tk = 1:NT+1
    t = (tk-1)*Δt
    for j = 1:n+1
        for i = 1:m+1
            x = (i-1)*h
            y = (j-1)*h
            U0[(j-1)*(m+1)+i, tk] = uexact1(x, y, t)
            U0[(j-1)*(m+1)+i+(m+1)*(n+1), tk] = uexact2(x, y, t)
        end
    end
    for i = 1:m
        for j = 1:n
            x = (i-1/2)*h
            y = (j-1/2)*h
            U0[2*(m+1)*(n+1)+(j-1)*m+i, tk] = pexact(x, y, t)
        end
    end
end

visualize_displacement(U, m, n, h, name="_out")
visualize_pressure(U, m, n, h, name="_out")

visualize_displacement(abs.(U-U0), m, n, h; name="_diff")
visualize_pressure(abs.(U-U0), m, n, h; name="_diff")

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
"""
```
pseudo_measurement{T<:AbstractFloat}(m::AnSchorfheide{T},
    TTT::Matrix{T}, RRR::Matrix{T}, CCC::Vector{T})
```

Assign pseudo-measurement equation (a linear combination of states):

```
x_t = ZZ_pseudo*s_t + DD_pseudo
```
"""
function pseudo_measurement{T<:AbstractFloat}(m::AnSchorfheide{T},
                                              TTT::Matrix{T},
                                              RRR::Matrix{T},
                                              CCC::Vector{T})
    endo   = m.endogenous_states
    pseudo = m.pseudo_observables

    _n_states = n_states_augmented(m)
    _n_pseudo = n_pseudo_observables(m)

    # Initialize pseudo ZZ and DD matrices
    ZZ_pseudo = zeros(_n_pseudo, _n_states)
    DD_pseudo = zeros(_n_pseudo)

    ##########################################################
    ## PSEUDO-OBSERVABLE EQUATIONS
    ##########################################################

    ## Output
    ZZ_pseudo[pseudo[:y_t],endo[:y_t]] = 1.

    ## π_t
    ZZ_pseudo[pseudo[:π_t],endo[:π_t]] = 1.
    DD_pseudo[pseudo[:π_t]]            = 100*(m[:π_star]-1);

    ## z_t
    ZZ_pseudo[pseudo[:z_t], endo[:z_t]] = 1.

    ## Nominal FFR
    ZZ_pseudo[pseudo[:NominalFFR], endo[:R_t]] = 1.
    DD_pseudo[pseudo[:NominalFFR]] = m[:π_star] + m[:rA] + 4.0*m[:γ_Q]

    ## Real FFR
    ZZ_pseudo[pseudo[:RealFFR], endo[:R_t]] = 1.
    ZZ_pseudo[pseudo[:RealFFR], endo[:π_t]] = -4.
    DD_pseudo[pseudo[:RealFFR]] = m[:π_star] + m[:rA] + 4.0*m[:γ_Q] - 4.0*(100. * (m[:π_star] - 1.))

    return PseudoMeasurement(ZZ_pseudo, DD_pseudo)
end


text: # This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
module StorageColorTypes

using FixedPointNumbers
using ColorTypes

import Base: reinterpret
import FixedPointNumbers: floattype
import ColorTypes: red, green, blue, alpha, chroma, hue, comp1, comp2, comp3

export
    RGB565LE, RGB565BE, RGB565,
    RGB555LE, RGB555BE, RGB555,
    XRGB444LE, XRGB444BE, XRGB444,
    RGBX444LE, RGBX444BE, RGBX444,
    ARGB4444LE, ARGB4444BE, ARGB4444,
    ARGB1555LE, ARGB1555BE, ARGB1555,
    RGB332, RGB242,
    HSV24, HSL24, HLS24, Lab24, Luv24, HLS

const BE = ENDIAN_BOM === 0x01020304

include("types.jl")
include("conversions.jl")
include("traits.jl")

end # module


text:
# # @__NAME__

# PREAMBLE

# PKG_SETUP

# This example is a (slightly modified) contribution by [Aykut Satici](https://github.com/symplectomorphism).

# ## Setup

using Pkg # hide
Pkg.activate(@__DIR__) # hide
Pkg.instantiate() # hide
using LinearAlgebra
using RigidBodyDynamics
using StaticArrays

# ## Model definition

# We're going to create a [four-bar linkage](https://en.wikipedia.org/wiki/Four-bar_linkage) that looks like this:
# ![fourbar](fourbar.jpg)
#
# We'll 'cut' the mechanism at joint 4: joints 1, 2, and 3 will be part of the spanning tree of the mechanism, but joint 4 will be a 'loop joint' (see e.g. Featherstone's 'Rigid Body Dynamics Algorithms'), for which the dynamics will be enforced using Lagrange multipliers.
#
# First, we'll define some relevant constants:

## gravitational acceleration
g = -9.81

## link lengths
l_0 = 1.10
l_1 = 0.5
l_2 = 1.20
l_3 = 0.75

## link masses
m_1 = 0.5
m_2 = 1.0
m_3 = 0.75

## link center of mass offsets from the preceding joint axes
c_1 = 0.25
c_2 = 0.60
c_3 = 0.375

## moments of inertia about the center of mass of each link
I_1 = 0.333
I_2 = 0.537
I_3 = 0.4

## Rotation axis: negative y-axis
axis = SVector(0., -1., 0.);

# Construct the world rigid body and create a new mechanism:

world = RigidBody{Float64}("world")
fourbar = Mechanism(world; gravity = SVector(0., 0., g))

# Next, we'll construct the spanning tree of the mechanism,
# consisting of bodies 1, 2, and 3 connected by joints 1, 2, and 3.
# Note the use of the `moment_about_com` keyword (as opposed to `moment`):

# Link 1 and joint 1:

joint1 = Joint("joint1", Revolute(axis))
inertia1 = SpatialInertia(frame_after(joint1),
    com=SVector(c_1, 0, 0),
    moment_about_com=I_1*axis*transpose(axis),
    mass=m_1)
link1 = RigidBody(inertia1)
before_joint1_to_world = one(Transform3D,
    frame_before(joint1), default_frame(world))
attach!(fourbar, world, link1, joint1,
    joint_pose = before_joint1_to_world)

# Link 2 and joint 2:

joint2 = Joint("joint2", Revolute(axis))
inertia2 = SpatialInertia(frame_after(joint2),
    com=SVector(c_2, 0, 0),
    moment_about_com=I_2*axis*transpose(axis),
    mass=m_2)
link2 = RigidBody(inertia2)
before_joint2_to_after_joint1 = Transform3D(
    frame_before(joint2), frame_after(joint1), SVector(l_1, 0., 0.))
attach!(fourbar, link1, link2, joint2,
    joint_pose = before_joint2_to_after_joint1)

# Link 3 and joint 3:

joint3 = Joint("joint3", Revolute(axis))
inertia3 = SpatialInertia(frame_after(joint3),
    com=SVector(l_0, 0., 0.),
    moment_about_com=I_3*axis*transpose(axis),
    mass=m_3)
link3 = RigidBody(inertia3)
before_joint3_to_world = Transform3D(
    frame_before(joint3), default_frame(world), SVector(l_0, 0., 0.))
attach!(fourbar, world, link3, joint3, joint_pose = before_joint3_to_world)

# Finally, we'll add joint 4 in almost the same way we did the other joints,
# with the following exceptions:
# 1. both `link2` and `link3` are already part of the `Mechanism`, so the `attach!`
#    function will figure out that `joint4` will be a loop joint.
# 2. instead of using the default (identity) for the argument that specifies the
#    transform from the successor of joint 4 (i.e., link 3) to the frame directly after
# joint 4, we'll specify a transform that incorporates the $l_3$ offset.

## joint between link2 and link3
joint4 = Joint("joint4", Revolute(axis))
before_joint4_to_joint2 = Transform3D(
    frame_before(joint4), frame_after(joint2), SVector(l_2, 0., 0.))
joint3_to_after_joint4 = Transform3D(
    frame_after(joint3), frame_after(joint4), SVector(-l_3, 0., 0.))
attach!(fourbar, link2, link3, joint4,
    joint_pose = before_joint4_to_joint2, successor_pose = joint3_to_after_joint4)

# Note the additional non-tree joint in the printed `Mechanism` summary.

# ## Simulation

# As usual, we'll now construct a `MechanismState` and `DynamicsResult` for the
# four-bar `Mechanism`. We'll set some initial conditions for a simulation, which
# were solved for a priori using a nonlinear program (not shown here).

state = MechanismState(fourbar)
result = DynamicsResult(fourbar);

set_configuration!(state, joint1, 1.6707963267948966) # θ
set_configuration!(state, joint2, -1.4591054166649482) # γ
set_configuration!(state, joint3, 1.5397303602625536) # ϕ

set_velocity!(state, joint1, 0.5)
set_velocity!(state, joint2, -0.47295)
set_velocity!(state, joint3, 0.341)

## Invalidate the cache variables
setdirty!(state)

# Next, we'll do a 3-second simulation:

ts, qs, vs = simulate(state, 3., Δt = 1e-2);

# ## Visualization

# For visualization, we'll use [`MeshCatMechanisms`](https://github.com/JuliaRobotics/MeshCatMechanisms.jl),
# an external package based on RigidBodyDynamics.jl.

using MeshCatMechanisms

# Create a `MechanismVisualizer` for the four-bar linkage and open it in a new browser tab
# (see [`MeshCat.jl`](https://github.com/rdeits/MeshCat.jl) for other options):

mvis = MechanismVisualizer(fourbar, Skeleton(inertias=false));

#-

#nb ##NBSKIP
#nb open(mvis)
#md ## open(mvis)

# And animate:
MeshCatMechanisms.animate(mvis, ts, qs; realtimerate = 1.);

max_stars_repo_path: examples/3. Four-bar linkage/3. Four-bar linkage.jl


text:
edit(pkg) =
  isdir(Pkg.dir(pkg)) ?
    run(`atom $(Pkg.dir(pkg))`) :
    error("$pkg not installed")

isuntitled(p) = ismatch(r"^(\.\\|\./)?untitled-[\d\w]+(:\d+)?$", p)

function basepath(file)
  srcdir = joinpath(JULIA_HOME,"..","..","base")
  releasedir = joinpath(JULIA_HOME,"..","share","julia","base")
  normpath(joinpath(isdir(srcdir) ? srcdir : releasedir, file))
end

realpath′(p) = ispath(p) ? realpath(p) : p

fullpath(path) =
  (isuntitled(path) || isabspath(path) ? path : basepath(path)) |> realpath′

appendline(path, line) = line > 0 ? "$path:$line" : path

baselink(path, line) =
  path == "no file" ? span(".fade", path) :
    isuntitled(path) ? link(path, line, Text(appendline("untitled", line))) :
    isabspath(path)  ?
      link(path, line, Text(pkgpath(appendline(path, line)))) :
      link(basepath(path), line, Text(normpath(joinpath("base", appendline(path, line)))))

function pkgpath(path)
  m = match(r"([^/\\]+[/\\]src[/\\].*)$", path)
  m == nothing ? basename(path) : m.captures[1]
end


text:
"""
```
struct MMAApproxModel{TApprox <: AbstractMMAApprox} <: AbstractModel
    parent::AbstractModel
    objective_ineq_constraints::AbstractFunction
    approx_objective_ineq_constraints::TApprox
    box_min::AbstractVector
    box_max::AbstractVector
end
```

An approximate restricted model that uses the `MMAApprox` or `XMMAApprox` approximation of the objective and constraint functions. The lower and upper bounds of the approximate model are more restricted than the original model. See [`updatebounds!`](@ref) for more details on the bounds.
- `parent`: the original model.
- `objective_ineq_constraints`: a function that returns the exact objective value and constraint violations as a single vector.
- `approx_objective_ineq_constraints`: a function that returns the approximate objective value and constraint violations as a single vector.
- `box_min`: the restricted lower bounds on the decision variables.
- `box_max`: the restricted upper bounds on the decision variables.
"""
@params struct MMAApproxModel{TApprox <: AbstractMMAApprox} <: AbstractModel
    parent::AbstractModel
    objective_ineq_constraints::AbstractFunction
    approx_objective_ineq_constraints::TApprox
    box_min::AbstractVector
    box_max::AbstractVector
end

"""
    MMAApproxModel(parent::AbstractModel, x::AbstractVector; extended = false, kwargs...)

Constructs an MMA approximation of the model `parent` around the point `x`. If `extended` is true, an [`XMMAApprox`](@ref) approximation is used. Otherwise, the default [`MMAApprox`](@ref) approximation is used. `kwargs` can be used to pass additional options to `XMMAApprox`, e.g. setting the coefficients.
"""
function MMAApproxModel(
    parent::AbstractModel,
    x::AbstractVector;
    extended = false,
    kwargs...,
)
    T = eltype(x)
    σ = map(1:length(x)) do j
        diff = getmax(parent, j) - getmin(parent, j)
        if !isfinite(diff)
            diff = 1000 * one(T)
        end
        return 0.5 * diff
    end
    ρ = zeros(length(getineqconstraints(parent)) + 1)
    obj_constr = getobjectiveconstraints(parent)
    approx_obj_constr = MMAApprox(
        obj_constr,
        x;
        σ = σ,
        ρ = ρ,
    )
    model = MMAApproxModel(
        parent,
        obj_constr,
        extended ? XMMAApprox(approx_obj_constr; kwargs...) : approx_obj_constr,
        copy(getmin(parent)),
        copy(getmax(parent)),
    )
    updateapprox!(model)
    return model
end

getmin(m::MMAApproxModel)= m.box_min
getmax(m::MMAApproxModel) = m.box_max

@doc docρ
getρ(model::MMAApproxModel) = getρ(model.approx_objective_ineq_constraints)

@doc docρ
function setρ!(model::MMAApproxModel, ρ)
    setρ!(model.approx_objective_ineq_constraints, ρ)
    return model
end

@doc docσ
getσ(model::MMAApproxModel) = getσ(model.approx_objective_ineq_constraints)

@doc docσ
function setσ!(model::MMAApproxModel, σ)
    setσ!(model.approx_objective_ineq_constraints, σ)
    return model
end

@doc docxk
getxk(model::MMAApproxModel) = getxk(model.approx_objective_ineq_constraints)

@doc docxk
function setxk!(model::MMAApproxModel, x::AbstractVector)
    setxk!(model.approx_objective_ineq_constraints, x)
    return model
end

@doc docfk
getfk(model::MMAApproxModel) = getfk(model.approx_objective_ineq_constraints)

@doc docfk
function setfk!(model::MMAApproxModel, f)
    setfk!(model.approx_objective_ineq_constraints, f)
    return model
end

@doc doc∇fk
get∇fk(model::MMAApproxModel) = get∇fk(model.approx_objective_ineq_constraints)

@doc doc∇fk
function set∇fk!(model::MMAApproxModel, ∇f::AbstractVecOrMat)
    set∇fk!(model.approx_objective_ineq_constraints, ∇f)
    return model
end

getparent(model::MMAApproxModel) = model.parent

"""
    getapproxfg(model::MMAApproxModel)

Returns the last cached value of the approximate objective and constraint values.
"""
getapproxfg(model::MMAApproxModel) = getapproxfg(model.approx_objective_ineq_constraints)

"""
    updatebounds!(model::MMAApproxModel, x)

Updates the bounds on the decision variables in `model` using the asymptotes of the approximation and the bounds on the original variables. The update formula is described in the original [1987 MMA paper](https://onlinelibrary.wiley.com/doi/abs/10.1002/nme.1620240207) and the [2002 paper](https://epubs.siam.org/doi/abs/10.1137/S1052623499362822).
"""
function updatebounds!(model::MMAApproxModel, x)
    L, U = getasymptotes(model.approx_objective_ineq_constraints)
    getmin(model) .= max.(getmin(getparent(model)), L .+ 0.1 .* (x .- L))
    getmax(model) .= min.(getmax(getparent(model)), U .- 0.1 .* (U .- x))
    return model
end

# Must be called every time σ or ρ are updated
# Doesn't do function evaluations
"""
    updateapprox!(approx::MMAApproxModel)

Updates the MMA approximation in `model`. This is called when `σ`, `ρ` or `xk` change.
"""
function updateapprox!(model::MMAApproxModel)
    updateapprox!(model.approx_objective_ineq_constraints)
    updatebounds!(model, getxk(model))
    return model
end

# Performs a function evaluation at x
"""
    updateapprox!(model::MMAApproxModel, x)

Updates the MMA approximation in `model` making it around the point `x`.
"""
function updateapprox!(model::MMAApproxModel, x::AbstractVector)
    updateapprox!(model.approx_objective_ineq_constraints, x)
    updatebounds!(model, x)
    return model
end

# Doesn't do function evaluations
"""
    updateapprox!(model::MMAApproxModel, x, f, ∇f)

Updates the MMA approximation in `model` making it around the point `x` where `f` is the value of the approximated function at `x` and `∇f` is its gradient if `f` is a number and the Jacobian if `f` is a vector.
"""
function updateapprox!(model::MMAApproxModel, x::AbstractVector, f, ∇f)
    updateapprox!(model.approx_objective_ineq_constraints, x, f, ∇f)
    updatebounds!(model, x)
    return model
end


text:
@testset "Check identifiability of `ODESystem` object" begin
    @parameters a01 a21 a12
    @variables t x0(t) x1(t) y1(t) [output = true]
    D = Differential(t)

    eqs = [D(x0) ~ -(a01 + a21) * x0 + a12 * x1, D(x1) ~ a21 * x0 - a12 * x1, y1 ~ x0]

    de = ODESystem(eqs, t, name=:Test)
    correct = Dict{Any,Any}(
        "nonidentifiable" => [a12, a21, substitute(x1, t=>0), a01],
        "locally_not_globally" => Num[],
        "globally" => [substitute(x0, t=>0)]
    )
    output = identifiability_ode(de)
    # have to use string here, otherwise fails on Ubutnu
    # TODO: figure out test on Ubuntu
    @test isequal(Set(correct["nonidentifiable"]), Set(output["nonidentifiable"]))
    @test isequal(Set(correct["locally_not_globally"]), Set(output["locally_not_globally"]))
    @test isequal(Set(correct["globally"]), Set(output["globally"]))

    # --------------------------------------------------------------------------
        # # --------------------------------------------------------------------------
        # # ----------

    @parameters a01 a21 a12
    @variables t x0(t) x1(t) y1(t) [output = true]
    D = Differential(t)
    using SpecialFunctions

    eqs = [
        D(x0) ~ -(a01 + a21) * SpecialFunctions.erfc(x0) + a12 * x1,
        D(x1) ~ a21 * x0 - a12 * x1,
        y1 ~ x0
    ]

    de = ODESystem(eqs, t, name=:Test)

    @test_throws ArgumentError identifiability_ode(de)
    # ----------
    # @parameters a b c
    # @variables t x1(t) x2(t) y(t) [output = true]
    # D = Differential(t)

    # eqs = [D(x1) ~ -a * x1 + x2 * b / (x1 + b / (c^2 - x2)), D(x2) ~ x2 * c^2 + x1, y ~ x2]
    # de = ODESystem(eqs, t, name=:Test)
    # correct = Dict(a => :globally, b => :globally, c => :locally)
    # to_check = [a, b, c]
    # @test isequal(correct, assess_identifiability(de, to_check))
end

max_stars_repo_name: iliailmer/SIAN-Julia

max_issues_repo_name: iliailmer/SIAN-Julia

max_forks_repo_name: iliailmer/SIAN-Julia


text:
using Distributions
using DelimitedFiles
using Distances
using DataFrames
using Statistics
# using StatsBase

# Compute DIC for ABC rejection on simulated data directly based on simulated summary statistics under different models

# For a given set of parameters, simulate 10 target datasets
ntest = parse(UInt32, ARGS[1])
# For a simulated target data, compute DIC for 100 times
ndic = parse(UInt32, ARGS[2])
num_sample = parse(UInt32, ARGS[3])
ncell = parse(UInt32, ARGS[4])
s3_weight = parse(Float64, ARGS[5])
target_fitness = 0
if(size(ARGS,1) > 5)
    target_fitness = parse(Float64, ARGS[6])
end

# Simulate target data
# ntest = 1
# ndic = 1
# num_sample = 10
# ncell = 20
# target_fitness = -0.2

println(ncell)
println(target_fitness)
println(s3_weight)

# agg_func = Statistics.median
# agg_func = median
agg_func = mean

# true values of the parameters (varing fitness values)
mu_sim1 = 0.1
mu_sim2 = 0.1
brate_sim = 0.4

# 200 used in the DIC paper
nsample = num_sample
msample = num_sample

nsims=100000
acc=0.001

max_mu = 0.5

min_brate = 0.3
max_brate = 0.5

# very unlikely values
LNL_MIN = -1e20

SMALL_VAL = 1e-9

# change directory accordingly
dir = "../../"
progSim = dir * "/code/simpdo"
simdir = dir * "/data/power_analysis/"

rsuffix = "_logs3_weight" * string(s3_weight) * ".tsv"

skip = 3

############################## Simulate target data ##############################
function getTargetdata_neutral(seed, mu_sim1, mu_sim2, brate_sim)
    mu1 = mu_sim1
    mu2 = mu_sim2
    birth_rate = brate_sim
    death_rate = 0

    suffix=string(seed)

    cmdSim=`$progSim -o $dir -e $ncell --chr_prob $mu1 --arm_prob $mu2 --birth_rate $birth_rate --death_rate $death_rate --skip $skip --seed $seed --suffix $suffix`

    res = chomp(read(cmdSim, String))
    arr = split(res,"\t")
    targetdata = map(x->parse(Float64,x), arr)

    targetdata = targetdata[1:3]

    return targetdata
end

function getTargetdata_withsel(seed, mu_sim1, mu_sim2, brate_sim, fitness)
    mu1 = mu_sim1
    mu2 = mu_sim2
    birth_rate = brate_sim
    death_rate = 0
    suffix=string(seed)

    fitval = "1 " * string(fitness)
    model = 1

    cmdSim=`$progSim -o $dir -e $ncell --model $model --fitness "$fitval" --chr_prob $mu1 --arm_prob $mu2 --birth_rate $birth_rate --death_rate $death_rate --skip $skip --seed $seed --suffix $suffix`

    res = chomp(read(cmdSim, String))
    arr = split(res,"\t")
    targetdata = map(x->parse(Float64,x), arr)

    targetdata = targetdata[1:3]

    return targetdata
end

############################## standard routine to compute DIC ##############################
# Compute surrogate likelihood
function compute_surrogate_likelihood(simdata, targetdata, scales)
    # likelihood = exp.(-0.5 * sqeuclidean(log.(targetdata), log.(simdata))) / sqrt(2 * pi)
    likelihood = exp.(-0.5 * sqeuclidean((targetdata), (simdata))) / sqrt(2 * pi)
    return likelihood
end

function compute_lnl_nsample(params, nsample, targetdata, scales)
  #println(params)
  mu1 = params[1]
  mu2 = params[2]
  birth_rate = params[3]
  death_rate = 0

  if size(params, 1) > 3
      fitval = "1 " * string(params[4])
      model = 1
  else
      fitval = "1 0"
      model = 0
  end
  # println(model)
  # println(fitval)

  sum_likelihood = 0.0
  for i in 1:nsample
      seed=rand(UInt64, 1)[1]
      suffix=string(seed)

      # call simulation program to generate data
      cmdSim=`$progSim -o $dir -e $ncell --model $model --fitness "$fitval" --chr_prob $mu1 --arm_prob $mu2 --birth_rate $birth_rate --death_rate $death_rate --skip $skip --seed $seed --suffix $suffix`

      res = chomp(read(cmdSim, String))
      arr = split(res,"\t")
      simdata = map(x->parse(Float64,x), arr)
      simdata = simdata[1:3]
      # println(simdata)
      # simdata = vcat(simdata[1:2], log(simdata[3]/ncell))

      # println("original simdata ", simdata)
      # simdata[simdata.==0] .= SMALL_VAL
      # simdata = vcat(log(simdata[1] + 1)/scales[1], log(simdata[2] + 1)/scales[2], log(simdata[3])/scales[3])
      # simdata = vcat(log(simdata[1])/scales[1], log(simdata[2])/scales[2], log(simdata[3])/scales[3])
      simdata = vcat((simdata[1]), (simdata[2]), s3_weight * log(simdata[3]))
      # println("log simdata ", simdata)

      ll = compute_surrogate_likelihood(simdata, targetdata, scales)
      #println(string(simdata) * " with likelihood " * string(ll))

      sum_likelihood += ll
  end

  #println("sum_likelihood of nsample: " * string(sum_likelihood))
  if sum_likelihood > 0
      lnl = log(sum_likelihood / nsample)
  else
      lnl = LNL_MIN
  end

  return lnl
end

function compute_dev_param(params, nsample, targetdata, scales)
    dev = -2 * compute_lnl_nsample(params, nsample, targetdata, scales)
    return dev
end

function compute_devbar(abcres, type, msample, nsample, targetdata, scales)
    sum_likelihood = 0.0
    resampled = sample(1:size(abcres,1), msample)

    for i in 1:msample
        sel = abcres[resampled[i],:]

        if type == "neutral"
            params = [sel[1], sel[2], sel[3]]
        else
            params = [sel[1], sel[2], sel[3], sel[4]]
        end
        lnl = compute_lnl_nsample(params, nsample, targetdata, scales)
        sum_likelihood += lnl
        # println("sum_likelihood for msample " * string(i) * " is " * string(sum_likelihood))
    end

    devbar = -2 * sum_likelihood / msample

    return devbar
end

function compute_dic(devbar, dev_param)
    dic = 2 * devbar - dev_param
    return dic
end

# type: the type under which ABC is run
# agg_func: can be mean or meadian
function get_dic_by_type(type, abcres, msample, nsample, targetdata, agg_func, scales)
    # Compute posterior mean
    pst_mu1 = agg_func(abcres[:,1])
    pst_mu2 = agg_func(abcres[:,2])
    pst_brate = agg_func(abcres[:,3])

    if type == "neutral"
        pparams = [pst_mu1, pst_mu2, pst_brate]
    else
        pst_fitness = agg_func(abcres[:,4])
        pparams = [pst_mu1, pst_mu2, pst_brate, pst_fitness]
    end
    # println(pparams)
    # println("original targetdata ", targetdata)
    # targetdata[targetdata.==0] .= SMALL_VAL
    # targetdata = vcat(log(targetdata[1] + 1)/scales[1], log(targetdata[2] + 1)/scales[2], log(targetdata[3])/scales[3])
    # targetdata = vcat(log(targetdata[1])/scales[1], log(targetdata[2])/scales[2], log(targetdata[3])/scales[3])
    targetdata = vcat((targetdata[1]), (targetdata[2]), s3_weight * log(targetdata[3]))
    # println("log targetdata ", targetdata)

    dev_param = compute_dev_param(pparams, nsample, targetdata, scales)
    devbar = compute_devbar(abcres, type, msample, nsample, targetdata, scales)
    dic = compute_dic(devbar, dev_param)
    # println(devbar)
    # println(dev_param)
    # println(dic)

    return dic
end

############################## Get DIC with ABC REJ ##############################

# Get variance of summary statistics from prior distributions (simulations)
# use log scale due to skewness and huge range of s3
function get_sstat_var(simdata_sstat,)
    simdata_all[simdata_sstat.==0] .= SMALL_VAL
    # sstart = 6
    # var1 = std(log.(simdata_all[:, sstart] .+ 1))
    # var2 = std(log.(simdata_all[:, sstart+1] .+ 1))
    var1 = std(log.(simdata_sstat[:, 1]))
    var2 = std(log.(simdata_sstat[:, 2]))
    # var1 = 1
    # var2 = 1
    var3 = std(log.(simdata_sstat[:, 3]))

    return  [var1, var2, var3]
end

# acceptance rate: acc
# should generate simulations for different number of cells
# separate abc from dic computation to check posterior distribution
function get_abc_res(type, simdata_all, sstart, scales, ncell, nsims, msample, nsample, targetdata, agg_func, acc)
    # println("running abc rej")
    # targetdata = vcat(log(targetdata[1] + 1)/scales[1], log(targetdata[2] + 1)/scales[2], log(targetdata[3])/scales[3])
    targetdata = vcat((targetdata[1]), (targetdata[2]), s3_weight * log(targetdata[3]))

    dim_stat = size(simdata_all, 2)  # dimension of simulation results
    # Compute the distance between simulated data and targetdata
    dists = zeros(Float64, nsims)
    for i in 1:nsims
        simdata1 = simdata_all[i, :]

        simdata = simdata1[sstart:sstart+2]
        # println("original simdata ", simdata)
        # simdata = vcat(log(simdata[1] + 1)/scales[1], log(simdata[2] + 1)/scales[2], log(simdata[3])/scales[3])
        simdata = vcat((simdata[1]), (simdata[2]), s3_weight * log(simdata[3]))

        # just used for sorting, no need to use log distance. Use original distance to be consistent with previous results. Biased with longer branches
        # println("log targetdata ", targetdata)
        # println("log simdata ", simdata)
        dist = evaluate(Euclidean(1e-16), (targetdata), (simdata))
        # dist = evaluate(Euclidean(1e-16), log.(targetdata), log.(simdata))
        dists[i] = dist
    end

    simdata_all = hcat(simdata_all, dists)

    # Sort data by distance and take top acceptance rate as rejection results
    ndim = dim_stat + 1
    naccept = Int(size(simdata_all, 1) * acc)
    abcres = sortslices(simdata_all, dims = 1, by = x -> x[ndim])[1:naccept,:]

    return abcres
end

# write results line by line to reduce memory usage and get results in time
fout = simdir * "/dic_ncell" * string(ncell) * "_fitness" * string(target_fitness) * "_ntest" * string(ntest) * "_ndic" * string(ndic) * "_nsample" * string(num_sample) * "_f" * string(agg_func) * "_acc" * string(acc) *  rsuffix

fsim_suffix = "_ncell" * string(ncell)  * "_n" * string(nsims) * "_maxmu" * string(max_mu) * "_minb" * string(min_brate) * "_maxb" * string(max_brate) * ".tsv"

type = "neutral"
fsim = simdir * "/sim_data_" * type * fsim_suffix
simdata_all_neutral = readdlm(fsim)
sstart_neutral = 5  # start index of summary statistics, from 5 when estimating birth rate

type = "selection"
fsim = simdir * "/sim_data_" * type * fsim_suffix
simdata_all_selection = readdlm(fsim)
sstart_selection = 6

# no need to do scaling when assuming neutral model
scales_neutral = [1, 1, 1]
scales_selection = [1, 1, 1]

if target_fitness == 0
    println("target data under neutral evolution")
    # res_all = zeros(Float64, ntest * ndic, 10)
    open(fout, "a") do f
        for i in 1:ntest
            seed = rand(UInt64, 1)[1]
            targetdata = getTargetdata_neutral(seed, mu_sim1, mu_sim2, brate_sim)
            # println(targetdata)

            for j in 1:ndic
                type = "neutral"
                abcres_neutral = get_abc_res(type, simdata_all_neutral, sstart_neutral, scales_neutral, ncell, nsims, msample, nsample, targetdata, agg_func, acc)
                dic_neutral = get_dic_by_type(type, abcres_neutral, msample, nsample, targetdata, agg_func, scales_neutral)

                type = "selection"
                abcres_withsel = get_abc_res(type, simdata_all_selection, sstart_selection, scales_selection, ncell, nsims, msample, nsample, targetdata, agg_func, acc)
                dic_withsel = get_dic_by_type(type, abcres_withsel, msample, nsample, targetdata, agg_func, scales_selection)

                # println(string(dic_neutral) * " " * string(dic_withsel))
                res = vcat(ncell, mu_sim1, mu_sim2, brate_sim, target_fitness, seed, targetdata, j, dic_neutral, dic_withsel)
                # println(res)

                k = (i - 1) * ndic + j
                println(k)
                # res_all[k,:] = res
                writedlm(f, transpose(res))
            end
        end
    end
else
    println("target data with selection")
    # res_all = zeros(Float64, ntest * ndic, 11)
    open(fout, "a") do f
        for i in 1:ntest
            seed = rand(UInt64, 1)[1]
            targetdata = getTargetdata_withsel(seed, mu_sim1, mu_sim2, brate_sim, target_fitness)
            # println(targetdata)
            for j in 1:ndic
                type = "neutral"
                abcres_neutral = get_abc_res(type, simdata_all_neutral, sstart_neutral, scales_neutral, ncell, nsims, msample, nsample, targetdata, agg_func, acc)
                dic_neutral = get_dic_by_type(type, abcres_neutral, msample, nsample, targetdata, agg_func, scales_neutral)

                type = "selection"
                abcres_withsel = get_abc_res(type, simdata_all_selection, sstart_selection, scales_selection, ncell, nsims, msample, nsample, targetdata, agg_func, acc)
                dic_withsel = get_dic_by_type(type, abcres_withsel, msample, nsample, targetdata, agg_func, scales_selection)

                # println(string(dic_neutral) * " " * string(dic_withsel))
                res = vcat(ncell, mu_sim1, mu_sim2, brate_sim, target_fitness, seed, targetdata, j, dic_neutral, dic_withsel)
                # println(res)

                k = (i - 1) * ndic + j
                println(k)
                # res_all[k,:] = res
                writedlm(f, transpose(res))
            end
        end
    end
end


text:
using ConcurrentPlayground
using Test

@testset "ConcurrentPlayground.jl" begin
    # Write your tests here.
end


text:
include("plotter_setup.jl")
using LinearAlgebra

# visualization cl
cla()
fig, axes = subplots(1, 3, figsize=[6.,1.8])

for (ax, a, b) = zip(axes, [.01, .01, .01], [.2, .5, .99])
    @assert a < 3*b
    @assert a + b <= 1

    l1 = ax.plot([0, 1], [(1+b)/2, 1/2], "-k", label="norm. Laplacian")
    l2 = ax.plot([0, 1], [1/2+b, 1/2], "--r", label="adjacency matrix")
    l3 = ax.plot([0, 1], [(1+b), 1], "b", label="classical", linestyle="dashdot")
    setp(ax, xlabel=L"degree exponent $a+ \frac{i}{n}b$")
    ax.hlines([0.5, 1.], -0.05, 1.05, linestyle=":", color="k", linewidth=1)
    ax.text(.55, 1.75, latexstring("b = $b"))
end

setp(axes, yticks=0:0.5:2, yticklabels=fill("", 5), xticks=[0,1], xticklabels=[L"a", L"a+b"])
setp(axes[1], ylabel="time exponent", yticklabels=0:0.5:2)
setp(axes, ylim=[-0.0,2.0], xlim=[-0.0,1.0])

handles, labels = axes[1].get_legend_handles_labels()
leg = fig.legend(handles, labels, loc="upper center", ncol=3,framealpha=0)
leg.get_frame().set_linewidth(0.0)
savefig("plots/complex_search/cl_visualization.pdf", bbox_inches="tight", metadata = Dict("CreationDate" => nothing))

## BA

m0 = 3
nodes = [1, :n]
repeating = 20 #may be smaller than maximal
data_n = 100:100:5_000
matrices = ["norm. laplacian", "laplacian", "adjacency matrix"]

for node=nodes
    cla()
    fig, axs = subplots(2, 3, figsize=(6,3.))

    for (i, matrix) = enumerate(matrices), m = 1:repeating
        color = 3/4* (1-m/repeating) .* (1,1,1)
        data_raw = npzread("complex_search/BA-$m0-$node-$(matrix[1:3])-$m.npz")
        axs[1,i].plot(data_n, data_raw[:,end] , "-k", linewidth=.1, color=color)
        axs[2,i].loglog(data_n, data_raw[:,end-1] ./ data_raw[:,end], "-", linewidth=.1, color=color)
        axs[1,i].set_title(replace(matrix, "laplacian" => "Laplacian"))
    end

    axs[1, 1].set(ylabel=L"succ. prob. $p(T)$")
    axs[2, 1].set(ylabel=L"exp. time $T/p(T)$")

    axs[1,2].tick_params(labelleft=false)
    axs[1,3].tick_params(labelleft=false)
    axs[2,2].tick_params(labelleft=false)
    axs[2,3].tick_params(labelleft=false)

    for i=1:3
        axs[1,i].set_xscale("log")
        axs[1,i].tick_params(labelbottom=false)
        axs[1,i].set_xlim(95,6_000)
        axs[2,i].set_xlim(95,6_000)
        axs[2,i].set(xlabel=L"graph order $n$")
        axs[2,i].tick_params(which="minor", left=true, labelleft=false)
        axs[2,i].set_yticks(vcat([ (2*k):(2*k):(9*k) for k=10 .^ (0:4)]...), minor=true)
        axs[2,i].set_yticks(10 .^ (0:4))
        axs[1,i].set_yticks(0:0.2:1.)
        axs[1,i].set_ylim(-0.05,1.05)
        axs[2,i].set_ylim(2, 15_000)
    end

    savefig("plots/complex_search/ba_$m0-$node.pdf", bbox_inches="tight", metadata = Dict("CreationDate" => nothing))
end

repeating = 200 # consider all!
for node=nodes
    cla()
    fig, axs = subplots(2, 3, figsize=(6,3.))
    for (i,matrix) = enumerate(matrices)
        data_alpha = zeros(repeating)
        data_reg_qual = zeros(repeating)

        for m = 1:repeating
            data_raw = npzread("complex_search/BA-$m0-$node-$(matrix[1:3])-$m.npz")
            data_y = log.(data_raw[:,end-1] ./ data_raw[:,end])
            beta, alpha = linreg(log.(data_n), data_y)
            data_alpha[m] = alpha
            data_reg_qual[m] = norm(data_y .- alpha*log.(data_n) .- beta)
        end
        @show maximum(data_reg_qual)
        axs[1,i].hist(data_alpha, bins=-0.0:.05:2.0, color="k", linewidth=.1)
        axs[2,i].plot(data_alpha, data_reg_qual, "kx")
        axs[1,i].set_title(replace(matrix, "laplacian" => "Laplacian"))
    end

    axs[1,2].tick_params(labelleft=false)
    axs[1,3].tick_params(labelleft=false)
    axs[2,2].tick_params(labelleft=false)
    axs[2,3].tick_params(labelleft=false)

    axs[1,1].set_ylabel("frequency")
    axs[2,1].set_ylabel("regression quality")
    for i=1:3
        axs[1,i].tick_params(labelbottom=false)
        axs[2,i].set(xlabel=L"exponent $\alpha$")
        axs[1,i].tick_params(which="minor", left=true, labelleft=false)
        axs[1,i].set_xticks(.25:.5:1.75, minor=true)
        axs[1,i].set_xticks(0.:.5:2.)
        axs[2,i].tick_params(which="minor", left=true, labelleft=false)
        axs[2,i].set_xticks(.25:.5:1.75, minor=true)
        axs[2,i].set_xticks(0.:.5:2.)

        axs[1,i].set_xlim(-0.1,2.1)
        axs[2,i].set_xlim(-0.1,2.1)
        axs[1,i].set_ylim(0, 130) #TODO fill
        axs[2,i].set_ylim(-0.2, 10.2) #TODO fill
    end

    savefig("plots/complex_search/ba_$m0-$node-exponent.pdf", bbox_inches="tight", metadata = Dict("CreationDate" => nothing))
end