text:
using JuMP
using StatsBase
using DataFrames
using CSV
using Plots
include("SupportFunction.jl")
println("Leeeeroooy Jenkins")

industry = "NoDur"
folder = "12"
folderArr = [12, 24, 36, 48, 120, 240]
industryArr = ["Durbl", "Enrgy", "HiTec", "Hlth", "Manuf", "NoDur", "Other", "Shops", "Telcm", "Utils"]
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"*industry*"/"*folder*"-1/"

# Common parameters within 1 industry over time but from different datasets. Loop over datasets
# Hypothesis is that some paraemeters are always important no matter what the expansions of the dataset we are dealing with.
j=2
trainingSize=240
path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"
#Fetch summary file
println("Predicting industry ",industryArr[j])
summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)

testMat = Array{Array{Float64}}(24)

for i = 1:24
    #Fetch dataset from summary
    fileName = summary[i,1]

    #Fetch best gamma from summary, based on R^2
    bestGamma = summary[i,7]

    nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

    #Fetch bSolved matrix
    testMat[i] = bSolveMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
        nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))
end

#=Filenames:
1,     2     3      4       5       6       7        8
Clean, TA,   Exp,   ExpTA,  Time,   TimeTA, TimeExp, TimeExpTA
=#
### Raw datasets ###
# Check for all datasets if the raw input are used
# Are any of the raw inputs always used?
# How many percent of the 8 datasets uses the prediction term in period t
getRaw()
function getRaw()
    avgUsage = zeros(10,10)
    signifMat = Array{String}(10,10)
    for j = 1:10
        println("Predicting industry ", industryArr[j])
        summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)
        testMat = zeros(845,10)
        for i = 1:1

            #Fetch dataset from summary
            fileName = summary[i,1]

            #Fetch best gamma from summary, based on R^2
            bestGamma = summary[i,7]

            nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

            #Fetch bSolved matrix
            tempMat = bSolveMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
                nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))

            testMat += Array(tempMat)[:,1:10]
            # 10 industries predicting 10 industries
        end
        testMat = testMat/1
        testMat[find(testMat)] = 1

        for k = 1:10
            avgUsage[j,k] = countnz(testMat[:,k])/845
        end
    end

    realmean = mean(avgUsage)
    realStd = std(avgUsage)
    for j = 1:10
        for k = 1:10
            z_score = (avgUsage[j,k]-realmean)/(realStd)
            if z_score > 2.326347874
                signifMat[j,k] = "***"
            elseif z_score > 1.644853627
                signifMat[j,k] = "**"
            elseif z_score > 1.281551566
                signifMat[j,k] = "*"
            else
                signifMat[j,k] = " "
            end
        end
    end

    println(signifMat)
    println(avgUsage)

    # 1.959964 --> 0.95
    # 2.575829 --> 0.99
    # 3.290527 --> 0.999
    writedlm(path*"ParametersAnalysis/industryRawAverageSignif.csv",signifMat,",")
    writedlm(path*"ParametersAnalysis/industryRawAverageUsage.csv",avgUsage,",")
end

getRawCorrectPredictions()
function getRawCorrectPredictions()
    avgUsage = zeros(10,10)
    avgCorrectUsage = zeros(10,10)
    signifMat = Array{String}(10,10)
    for j = 1:10
        println("Predicting industry ", industryArr[j])
        summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)
        testMat = zeros(845,10)
        for i = 1:1
            #Fetch dataset from summary
            fileName = summary[i,1]

            #Fetch best gamma from summary, based on R^2
            bestGamma = summary[i,7]

            nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

            #Fetch bSolved matrix
            tempMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
                nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))

            realValues = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_real.csv",
                nullable=false, types=[Int64, Int64, Int64, Float64])

            predictedValues = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_predicted.csv",
                nullable=false, types=vcat([Int64, Int64, Int64],fill(Float64,10)))

            for t = 1:845
                if sign(predictedValues[t,3+bestGamma]) == sign(realValues[t,4])
                    for pt = 1:10
                        if tempMat[t,pt] != 0
                            avgCorrectUsage[j,pt] += 1
                        end
                    end
                end
            end

            testMat += Array(tempMat)[:,1:10]
        end

        testMat = testMat/1
        testMat[find(testMat)] = 1

        for k = 1:10
            avgUsage[j,k] = countnz(testMat[:,k])/845
        end
    end

    avgCorrectUsage = avgCorrectUsage/845

    avgCorrectUsagePercent = zeros(10,10)
    for i=1:10
        for j=1:10
            avgCorrectUsagePercent[i,j] = avgCorrectUsage[i,j]/avgUsage[i,j]
        end
    end
    # 1.959964 --> 0.95
    # 2.575829 --> 0.99
    # 3.290527 --> 0.999
    writedlm(path*"ParametersAnalysis/RawIndustryCorrectUsage.csv",avgCorrectUsagePercent,",")
end

getMacroTimeTA()
function getMacroTimeTA()
    path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"
    avgUsage = zeros(10,35)
    signifMat = Array{String}(10,35)
    for j = 1:10
        println("Predicting industry ", industryArr[j])
        summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)
        testMat = zeros(845,359)
        for i = 14
            #Fetch dataset from summary
            fileName = summary[i,1]

            #Fetch best gamma from summary, based on R^2
            bestGamma = summary[i,7]

            nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

            #Fetch bSolved matrix
            tempMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
                nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))

            testMat += Array(tempMat)
            # 10 industries predicting 10 industries
        end
        testMat[find(testMat)] = 1

        for k = 1:27
            for t = 0:12
                avgUsage[j,k] += countnz(testMat[:,k+27*t])/845
            end
            avgUsage[j,k] = avgUsage[j,k]/13
        end

        for k = 28:35
            avgUsage[j,k] = countnz(testMat[:, end-35+k])/845
        end
    end

    realmean = mean(avgUsage)
    realStd = std(avgUsage)
    for j = 1:10
        for k = 1:35
            z_score = (avgUsage[j,k]-realmean)/(realStd)
            if z_score > 2.326347874
                signifMat[j,k] = "***"
            elseif z_score > 1.644853627
                signifMat[j,k] = "**"
            elseif z_score > 1.281551566
                signifMat[j,k] = "*"
            else
                signifMat[j,k] = " "
            end
        end
    end

    println(signifMat)
    println(avgUsage)

    writedlm(path*"ParametersAnalysis/MacroTimeTAAverageSignif.csv",signifMat,",")
    writedlm(path*"ParametersAnalysis/MacroTimeTAAverageUsage.csv",avgUsage,",")
end

getMacroTimeTAMostUsedVariables()
function getMacroTimeTAMostUsedVariables()
    path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"
    avgUsage = zeros(13,27)
    signifMat = Array{String}(13,27)

    allMat = zeros(845,359)
    for j = 1:10
        println("Predicting industry ", industryArr[j])
        summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)
        testMat = zeros(845,359)
        for i = 14
            #Fetch dataset from summary
            fileName = summary[i,1]

            #Fetch best gamma from summary, based on R^2
            bestGamma = summary[i,7]

            nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

            #Fetch bSolved matrix
            tempMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
                nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))

            testMat += Array(tempMat)
            # 10 industries predicting 10 industries
        end
        testMat[find(testMat)] = 1

        allMat += testMat

        #find(x -> x>0.9, mean(testMat,1))

        for k = 1:27
            for t = 0:12
                avgUsage[t+1,k] += countnz(testMat[:,k+27*t])/845
            end
        end
    end

    allMat = allMat/10
    allMatMean = mean(allMat,1)
    mean(allMatMean)
    std(allMatMean)

    n = 500

    stepVector = zeros(n+1,3)
    for i = 0:n
        stepVector[i+1,1] = i/n
        x = length(find(x-> x <= i/n, mean(allMat,1)))
        stepVector[i+1,3] = x
        y = length(find(x-> x >= i/n && x < (i+1)/n, mean(allMat,1)))
        stepVector[i+1,2] = y
    end
    writedlm(path*"ParametersAnalysis/MacroTimeTAParameterDistribution.CSV",stepVector,",")

    x = find(x-> x< 0.05 && x>0.01, mean(allMat,1))
    x = find(x-> x> 0.5, mean(allMat,1))

    for i=x
        if floor(i/27) < 13
            println("T",Int64(floor(i/27)),", ",i%27)
        else
            println("TA term ",i%27)
        end
    end
    mean(allMat,1)[x]
end

getVixMacroTimeTA()
function getVixMacroTimeTA()
    path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"
    avgUsage = zeros(10,36)
    signifMat = Array{String}(10,36)
    for j = 1:10
        println("Predicting industry ", industryArr[j])
        summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)
        testMat = zeros(845,372)
        for i = 22
            #Fetch dataset from summary
            fileName = summary[i,1]

            #Fetch best gamma from summary, based on R^2
            bestGamma = summary[i,7]

            nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

            #Fetch bSolved matrix
            tempMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
                nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))

            testMat += Array(tempMat)
            # 10 industries predicting 10 industries
        end
        testMat[find(testMat)] = 1

        for k = 1:28
            for t = 0:12
                avgUsage[j,k] += countnz(testMat[:,k+27*t])/845
            end
        end

        for k = 29:36
            avgUsage[j,k] = countnz(testMat[:, end-36+k])/845
        end
    end

    realmean = mean(avgUsage)
    realStd = std(avgUsage)
    for j = 1:10
        for k = 1:10
            z_score = (avgUsage[j,k]-realmean)/(realStd)
            if z_score > 2.326347874
                signifMat[j,k] = "***"
            elseif z_score > 1.644853627
                signifMat[j,k] = "**"
            elseif z_score > 1.281551566
                signifMat[j,k] = "*"
            else
                signifMat[j,k] = " "
            end
        end
    end

    println(signifMat)
    println(avgUsage)

    writedlm(path*"ParametersAnalysis/VIXMacroTimeTAAverageSignif.csv",signifMat,",")
    writedlm(path*"ParametersAnalysis/VIXMacroTimeTAAverageUsage.csv",avgUsage,",")
end

getVixPeriodMacroTimeTA()
function getVixPeriodMacroTimeTA()
    path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"
    avgUsage = zeros(10,36)
    signifMat = Array{String}(10,36)
    for j = 1:10
        println("Predicting industry ", industryArr[j])
        summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)
        testMat = zeros(323,372)
        for i = 22
            #Fetch dataset from summary
            fileName = summary[i,1]

            #Fetch best gamma from summary, based on R^2
            bestGamma = summary[i,7]

            nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

            #Fetch bSolved matrix
            tempMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
                nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))

            testMat += Array(tempMat)[end-322:end,:]
            # 10 industries predicting 10 industries
        end
        testMat[find(testMat)] = 1

        for k = 1:28
            for t = 0:12
                avgUsage[j,k] += countnz(testMat[:,k+27*t])/323
            end
        end

        for k = 29:36
            avgUsage[j,k] = countnz(testMat[:, end-36+k])/323
        end
    end

    realmean = mean(avgUsage)
    realStd = std(avgUsage)
    for j = 1:10
        for k = 1:36
            z_score = (avgUsage[j,k]-realmean)/(realStd)
            if z_score > 3.090232306
                signifMat[j,k] = "***"
            elseif z_score > 2.326347874
                signifMat[j,k] = "**"
            elseif z_score > 1.644853627
                signifMat[j,k] = "*"
            else
                signifMat[j,k] = " "
            end
        end
    end

    println(signifMat)
    println(avgUsage)

    writedlm(path*"ParametersAnalysis/VIXPeriodMacroTimeTAAverageSignif.csv",signifMat,",")
    writedlm(path*"ParametersAnalysis/VIXPeriodMacroTimeTAAverageUsage.csv",avgUsage,",")
end

getVixMacroTimeTAMostUsedVariables()
function getVixMacroTimeTAMostUsedVariables()
    path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"
    avgUsage = zeros(13,28)
    signifMat = Array{String}(13,28)

    allMat = zeros(845,372)
    for j = 1:10
        println("Predicting industry ", industryArr[j])
        summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)
        testMat = zeros(845,372)
        for i = 22
            #Fetch dataset from summary
            fileName = summary[i,1]

            #Fetch best gamma from summary, based on R^2
            bestGamma = summary[i,7]

            nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

            #Fetch bSolved matrix
            tempMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
                nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))

            testMat += Array(tempMat)
            # 10 industries predicting 10 industries
        end
        testMat[find(testMat)] = 1

        allMat += testMat

        #find(x -> x>0.9, mean(testMat,1))

        for k = 1:28
            for t = 0:12
                avgUsage[t+1,k] += countnz(testMat[:,k+27*t])/845
            end
        end
    end

    allMat = allMat/10
    allMatMean = mean(allMat,1)
    mean(allMatMean)
    std(allMatMean)

    stepVector = zeros(10001,3)
    for i = 0:10000
        stepVector[i+1,1] = i/10000
        x = length(find(x-> x <= i/10000, mean(allMat,1)))
        stepVector[i+1,3] = x
        y = length(find(x-> x >= i/10000 && x < (i+1)/10000, mean(allMat,1)))
        stepVector[i+1,2] = y
    end
    writedlm(path*"ParametersAnalysis/VIXMacroTimeTAParameterDistribution.CSV",stepVector,",")

    x = find(x-> x< 0.05 && x>0.01, mean(allMat,1))
    x = find(x-> x> 0.25, mean(allMat,1))

    for i=x
        if floor(i/28) < 13
            println("T",Int64(floor(i/28)),", ",i%28)
        else
            println("TA term ",i%28)
        end
    end
    mean(allMat,1)[x]
end

function getSuperSummary()
    tempData = zeros(25,8)
    path = "/Users/SkipperAfRosenborg/Google Drive/DTU/10. Semester/Thesis/GitHubCode/Results/IndexData/LassoTest/"
    i="NoDur"
    for i in industryArr
        tempData += Array{Float64}(CSV.read(path*i*"/Summary New240.csv",nullable=false)[:,2:end])
    end
    tempData = tempData/10
    writedlm(path*"SuperSummaryNew.csv",tempData,",")
end

rawInputUsage()
#NoDur is clearly the strongest predictor and also shows the strongest prediction power in recent terms
function rawInputUsage()
    println("Predicting dataset ",industryArr[j])
    industryCounter = zeros(845,10)
    for t = 1:845
        for curIndustry = 1:10
            for dataset = 1:10
                if testMat[dataset][t,curIndustry] != 0
                    industryCounter[t,curIndustry] += 1
                end
            end
        end
    end
    industryCounter = industryCounter/10

    avgUsage = mean(industryCounter,1)
    for i=1:10
        println(industryArr[i], " is used \t", round(avgUsage[i],3),"% of the time across all datasets")
    end
    println()

    avgUsageIfActive = zeros(1,10)
    for curIndustry = 1:10
        iter = 0
        for t = 1:845
            if industryCounter[t,curIndustry] != 0
                iter += 1
                avgUsageIfActive[1,curIndustry] += industryCounter[t,curIndustry]
            end
        end
        avgUsageIfActive[1,curIndustry] = avgUsageIfActive[1,curIndustry]/iter
    end

    for i=1:10
        println(industryArr[i], " is used \t", round(avgUsageIfActive[i],3),"% of the time across all datasets when it's active")
    end

    #Create label rowvector
    industryArr2 = Array{String}(1,10)
    for i = 1:10
        industryArr2[1,i] = industryArr[i]
    end

    # Graph it
    plotlyjs()
    plot(industryCounter, title="Industry Overview", linewidth=2, label=industryArr2, yaxis=("% used in period t",(0,1),0:0.1:1),xaxis=("Months",0:100:845))
end

getRawTime()
function getRawTime()
    avgUsage = zeros(10,10)
    signifMat = Array{String}(10,10)
    for j = 1:10
        println("Predicting industry ", industryArr[j])
        summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)
        testMat = zeros(845,10)
        for i = 5:5
            #Fetch dataset from summary
            fileName = summary[i,1]

            #Fetch best gamma from summary, based on R^2
            bestGamma = summary[i,7]

            nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

            #Fetch bSolved matrix
            tempMat = bSolveMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
                nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))

            testMat = Array(tempMat)
            # 10 industries predicting 10 industries
        end
        countnz(testMat[8,:])
        testMat[find(testMat)] = 1

        for k = 1:10
            avgUsage[j,k] = countnz(testMat[:,[k,k+10,k+20,k+30,k+40,k+50,k+60,k+70,k+80,k+90,k+100,k+110,k+120]])/845
        end
    end

    realmean = mean(avgUsage)
    realStd = std(avgUsage)
    for j = 1:10
        for k = 1:10
            z_score = (avgUsage[j,k]-realmean)/(realStd)
            if z_score > 3.090232306
                signifMat[j,k] = "***"
            elseif z_score > 2.326347874
                signifMat[j,k] = "**"
            elseif z_score > 1.644853627
                signifMat[j,k] = "*"
            else
                signifMat[j,k] = " "
            end
        end
    end

    println(signifMat)
    println(avgUsage)

    # 1.959964 --> 0.95
    # 2.575829 --> 0.99
    # 3.290527 --> 0.999
    writedlm(path*"ParametersAnalysis/RawTimeIndustryAverageSignif.csv",signifMat,",")
    writedlm(path*"ParametersAnalysis/RawTimeIndustryAverageUsage.csv",avgUsage,",")
end

getRawExp()
function getRawExp()
    avgUsage = zeros(10,10)
    signifMat = Array{String}(10,10)
    for j = 1:10
        println("Predicting industry ", industryArr[j])
        summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)
        testMat = zeros(845,10)
        for i = 3
            #Fetch dataset from summary
            fileName = summary[i,1]

            #Fetch best gamma from summary, based on R^2
            bestGamma = summary[i,7]

            nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

            #Fetch bSolved matrix
            tempMat = bSolveMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
                nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))

            testMat = Array(tempMat)
            # 10 industries predicting 10 industries
        end
        countnz(testMat[8,:])
        testMat[find(testMat)] = 1

        for k = 1:10
            avgUsage[j,k] = countnz(testMat[:,[k,k+10,k+20,k+30]])/845
        end
    end

    realmean = mean(avgUsage)
    realStd = std(avgUsage)
    for j = 1:10
        for k = 1:10
            z_score = (avgUsage[j,k]-realmean)/(realStd)
            if z_score > 3.090232306
                signifMat[j,k] = "***"
            elseif z_score > 2.326347874
                signifMat[j,k] = "**"
            elseif z_score > 1.644853627
                signifMat[j,k] = "*"
            else
                signifMat[j,k] = " "
            end
        end
    end

    println(signifMat)
    println(avgUsage)

    # 1.959964 --> 0.95
    # 2.575829 --> 0.99
    # 3.290527 --> 0.999
    writedlm(path*"ParametersAnalysis/RawExpIndustryAverageSignif.csv",signifMat,",")
    writedlm(path*"ParametersAnalysis/RawExpIndustryAverageUsage.csv",avgUsage,",")
end

getRawTimeExp()
function getRawExp()
    avgUsage = zeros(10,10)
    signifMat = Array{String}(10,10)
    for j = 1:10
        println("Predicting industry ", industryArr[j])
        summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)
        testMat = zeros(845,10)
        for i = 7
            #Fetch dataset from summary
            fileName = summary[i,1]

            #Fetch best gamma from summary, based on R^2
            bestGamma = summary[i,7]

            nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

            #Fetch bSolved matrix
            tempMat = bSolveMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
                nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))

            testMat = Array(tempMat)[:,:]
            # 10 industries predicting 10 industries
        end
        countnz(testMat[8,:])
        testMat[find(testMat)] = 1

        lRange = vcat(Array{Int64}(15:15+11),Array{Int64}(28:28+11),Array{Int64}(41:41+11))

        for k = 1:10
            for l = lRange
                avgUsage[j,k] += countnz(testMat[:,[k+10*(l-1)]])/845
            end
        end
    end

    realmean = mean(avgUsage)
    realStd = std(avgUsage)
    for j = 1:10
        for k = 1:10
            z_score = (avgUsage[j,k]-realmean)/(realStd)
            if z_score > 3.090232306
                signifMat[j,k] = "***"
            elseif z_score > 2.326347874
                signifMat[j,k] = "**"
            elseif z_score > 1.644853627
                signifMat[j,k] = "*"
            else
                signifMat[j,k] = " "
            end
        end
    end

    println(signifMat)
    println(avgUsage)

    # 1.959964 --> 0.95
    # 2.575829 --> 0.99
    # 3.290527 --> 0.999
    writedlm(path*"ParametersAnalysis/RawTimeExpRawTimeExpIndustryAverageSignif.csv",signifMat,",")
    writedlm(path*"ParametersAnalysis/RawTimeExpRawTimeExpIndustryAverageUsage.csv",avgUsage,",")
end

getMacro()
function getMacro()
    avgUsage = zeros(10,17)
    signifMat = Array{String}(10,17)
    for j = 1:10
        println("Predicting industry ", industryArr[j])
        summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)
        testMat = zeros(845,17)
        for i = 17-8:24-8

            #Fetch dataset from summary
            fileName = summary[i,1]

            #Fetch best gamma from summary, based on R^2
            bestGamma = summary[i,7]

            nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

            #Fetch bSolved matrix
            tempMat = bSolveMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
                nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))

            testMat += Array(tempMat)[:,11:27]
            # 10 industries predicting 10 industries
        end
        testMat = testMat/8
        testMat[find(testMat)] = 1

        for k = 1:17
            avgUsage[j,k] = countnz(testMat[:,k])/845
        end
    end

    realmean = mean(avgUsage)
    realStd = std(avgUsage)
    for j = 1:10
        for k = 1:17
            z_score = (avgUsage[j,k]-realmean)/(realStd)
            if z_score > 3.090232306
                signifMat[j,k] = "***"
            elseif z_score > 2.326347874
                signifMat[j,k] = "**"
            elseif z_score > 1.644853627
                signifMat[j,k] = "*"
            else
                signifMat[j,k] = " "
            end
        end
    end

    println(signifMat)
    println(avgUsage)

    writedlm(path*"ParametersAnalysis/MacroAverageSignif.csv",signifMat,",")
    writedlm(path*"ParametersAnalysis/MacroAverageUsage.csv",avgUsage,",")
end

getTA()
function getTA()
    avgUsage = zeros(10,8)
    signifMat = Array{String}(10,8)
    for j = 1:10
        println("Predicting industry ", industryArr[j])
        summary = CSV.read(path*industryArr[j]*"/"*"Summary "*string(trainingSize)*".csv", nullable=false)
        testMat = zeros(845,8)
        for i = [2,6,8,2+8,6+8,8+8,2+16,6+16,8+16]

            #Fetch dataset from summary
            fileName = summary[i,1]

            #Fetch best gamma from summary, based on R^2
            bestGamma = summary[i,7]

            nPredictionTerms = [10, 18, 40, 48, 130, 138, 520, 528, 27, 35, 108, 116, 351, 359, 1404, 1412, 28, 36, 112, 120, 364, 372, 1456, 1464]

            #Fetch bSolved matrix
            tempMat = bSolveMat = CSV.read(path*industryArr[j]*"/"*string(trainingSize)*"-1/240_"*fileName*"_bMatrix240_"*string(bestGamma)*".csv",
                nullable=false, header=false, datarow=1, types=fill(Float64,nPredictionTerms[i]))

            testMat += Array(tempMat)[:,end-7:end]
            # 10 industries predicting 10 industries
        end
        testMat = testMat/8
        testMat[find(testMat)] = 1

        for k = 1:8
            avgUsage[j,k] = countnz(testMat[:,k])/845
        end
    end

    realmean = mean(avgUsage)
    realStd = std(avgUsage)
    for j = 1:10
        for k = 1:8
            z_score = (avgUsage[j,k]-realmean)/(realStd)
            if z_score > 3.090232306
                signifMat[j,k] = "***"
            elseif z_score > 2.326347874
                signifMat[j,k] = "**"
            elseif z_score > 1.644853627
                signifMat[j,k] = "*"
            else
                signifMat[j,k] = " "
            end
        end
    end

    println(signifMat)
    println(avgUsage)

    writedlm(path*"ParametersAnalysis/TAAverageSignif.csv",signifMat,",")
    writedlm(path*"ParametersAnalysis/TAAverageUsage.csv",avgUsage,",")
end

rawNonLinearUsage()
# Are non linear transformation used?
# Check how many and which non linear transformation is used.
# Only x^2, log(x), sqrt(x)?
function rawNonLinearUsage()
    avgUsage = zeros(1,40)

    #Load Exp dataset
    for i=1:40
        avgUsage[1,i] = countnz(testMat[3][:,i])/845
    end

    rawAvgUse = mean(avgUsage[1:10])
    SquaredAvgUse = mean(avgUsage[11:20])
    LogAvgUse = mean(avgUsage[21:30])
    SqrtAvgUse = mean(avgUsage[31:40])

    println("Raw usage average usage ",round(rawAvgUse,3))
    for i=1:10
        println(industryArr[i],"is used ", round(avgUsage[i],3),"% on average")
    end

    println("\nSquared usage average usage ", round(SquaredAvgUse,3))
    for i=11:20
        println(industryArr[i-10],"is used ", round(avgUsage[i],3),"% on average")
    end

    println("\nLog usage average usage ", round(LogAvgUse,3))
    for i=21:30
        println(industryArr[i-20],"is used ", round(avgUsage[i],3),"% on average")
    end

    println("\nSqrt usage average usage ", round(SqrtAvgUse,3))
    for i=31:40
        println(industryArr[i-30],"is used ", round(avgUsage[i],3),"% on average")
    end
end

# Are time series used?
    # Are the time series usage based on a single period (t-2, t-3... ) across all industries
    # or a certain industry with multiple timeseries.
function rawTimeSeriesUsage()
    avgUsage = zeros(10,13)

    #Load Exp dataset
    for industryIndex=1:10
        for i=1:13
            avgUsage[industryIndex,i] = countnz(testMat[5][:,industryIndex+(i-1)*10])/845
        end
    end

    for i=0:12
        println("Time series ", i, " is on average used ", round(mean(avgUsage[:,i+1]),3))
    end

    for industryIndex=1:10
        for i=0:12
            println("Time series ", i, " is on average used ", round(avgUsage[industryIndex, i+1],3), " in industry ",industryArr[industryIndex])
        end
        println()
    end
end

# Are TA used?
    # Are the TA input used?
function rawTAUsage()
    avgUsage = zeros(8)

    for dataset=[2,6,8]
        #Load TA dataset
        for i=1:8
            avgUsage[i] = countnz(testMat[dataset][:,end-8+i:end-8+i])/845
        end

        for i=1:8
            println("TA  ", i, " is on average used ", round(avgUsage[i],3), " in dataset ",dataset)
        end
        println("")
    end
end

# Are non linear transformation of times series used?
    # Important!
function rawTimeExpUsage()
    avgUsage = zeros(10)

    for i=1:10
        avgUsage[i] = countnz(testMat[1][:,i])
    end

    realmean = mean(avgUsage)
    degOfFree = 519

    stdev = std(avgUsage)

    for i=1:10
        zscore = (avgUsage[i]-realmean)/stdev
        if zscore > 1.965 || zscore < -1.965
            println("Paremeter ", i, " is significant, and had an average of ", avgUsage[i])
        end
    end
end

#Are the combination of parameters used within the past 10 years seen earlier?

for i = 9:9+8
    macroUsage(i)
end
# Macro dataset
function macroUsage(datIndx)
    dataset = testMat[datIndx]

    nCols = size(dataset)[2]
    nCols = 27

    avgUsage = zeros(nCols)

    for i=1:nCols
        avgUsage[i] = countnz(dataset[:,i])
    end

    realmean = mean(avgUsage)
    degOfFree = nCols-1

    stdev = std(avgUsage)

    for i=1:nCols
        zscore = (avgUsage[i]-realmean)/stdev
        if zscore > 1.965 || zscore < -1.965
            #println("Paremeter ", i, " is significant, and had an average of ", avgUsage[i])
        end
    end

    println(find(x -> x < 845*0.05, avgUsage))
    #println(avgUsage[11:27])
end

function macroUsage()
    nCols = 27

    avgUsage = zeros(8,nCols)

    for j=9:9+7
        for i=1:nCols
            avgUsage[j-8,i] = countnz(testMat[j][:,i])
        end
    end

    println(find(x -> x < 845*0.05, mean(avgUsage,1)))
    #println(avgUsage[11:27])
end
    # Are non linear transformation used
    # Are time series used?
    # Are TA used?
    # Are non linear transformation of times series used?
    # Are any of the raw inputs always used?
    # Is VIX actually used?

max_stars_repo_name: skipperafrosenborg/Thesis

max_issues_repo_name: skipperafrosenborg/Thesis

max_forks_repo_name: skipperafrosenborg/Thesis


text:
function dft(atoms  ::Atoms,                # atoms
             basis  ::Basis,                # molecular basis
             g      ::MultiGrid,            # integration grid
             N      ::Float64,              # number of particles
             vembed ::AbstractArray,        # embedding potential
             xc     ::ldaxc,                # exchange-correlation functional
             scf    ::SCF,
             mixer  ::AbstractMixer,        # density mixing
             lmax   ::Int,                  # max l-quantum number in hartree multipole expansion
             nh     ::Int,                  # number of hartree radial points
             nspl   ::Int,                  # number of spline points
             T      ::Float64,              # temperature
             logger ::Callable)
    local E,Etot,Enuc,vs,vh,δrho_mp,δvh_mp,evs,efs,f,mu,H,S,rho_out

    stats = Stats()
    init_mixer(mixer,g)
    init_scf(scf,g)

    @stats "* total" begin

        # initialize
        @stats "* evaluate vext"                vext     = sum([at[:v](g) for at in atoms])                 # external potential
                                                vext    += vembed                                           # embedding potential
        @stats "* evaluate rho_free"            rho_free = sum([at[:rho](g) for at in atoms])               # free atom density
        @stats "* evaluate vh_free"             vh_free  = sum([at[:vh](g) for at in atoms])                # free atom hartree potential
                                                rho_in   = isempty(scf.rho) ? rho_free : scf.rho[end]       # initial guess for density

                                                sort!(basis, by=φ->φ.R)                                     # sort basis by center
        @stats "* evaluate basis"               φ,∆φ = eval_basis_lapl(basis,g)                             # evaluate basis on grid
        @stats "* overlap"                      wφT  = (g.w.*φ)'                                            # grid weights times basis
        @stats "* overlap"                      S    = Symmetric(wφT*φ)                                     # overlap matrix
                                                Enuc = sum([μ != ν ? μ.Z*ν.Z/dist(μ.R,ν.R)/2 : 0.           # nuclear attraction
                                                    for μ in atoms, ν in atoms])

                                                lmax = max(lmax,maximum(Int[f.l for f in basis]))
        @stats "* sphharm basis"                realsphharmbasis!(atoms,lmax)                               # atomic spherical harmonics basis
        @stats "* sphharm tables"               sphharmtab,disttab = compute_sphharm_dist_tab(basis,g,lmax) # spherical harmonics, distances

        # storage
        vx = similar(rho_in); vc = similar(rho_in)
        ex = similar(rho_in); ec = similar(rho_in)
        vh = similar(rho_in); vs = similar(rho_in)
        hφ = similar(φ)

        while !scf.converged
                                                δrho = rho_in - rho_free
            @stats "poisson multipole"          δrho_mp,δvh_mp = poissonmultipole(g,δrho,atoms,         # solve hartree potential
                                                    nh=nh,nspl=nspl)
            @stats "eval hartree multipole"     vh[:] = isempty(δvh_mp) ? vh_free : vh_free + eval_basis_sum(δvh_mp,g,sphharmtab,disttab)

            @stats "exchange-correlation"       xc(rho_in,ex,vx,ec,vc)                                  # exchange-correlation
                                                vs[:] = vext + vx + vc + vh                             # Kohn-Sham potential

            @stats "hamiltonian"                h!(hφ,φ,∆φ,vs)                                          # apply Hamiltonian to basis
            @stats "hamiltonian"                H = Symmetric(wφT*hφ)                                   # Hamiltonian matrix
            @stats "KS equations"               evs,efs = LinAlg.eig(H,S)                               # compute eigenspectrum

                                                f,mu = compute_occupations(evs,2,N,T)
                                                ζ = f .> 1e-12                                          # nonzero occupations
                                                efs = efs[:,ζ]; evs = evs[ζ]; f = f[ζ]                  # discard unused eigenspectrum

            # @stats "normalize efs"              efs ./= sqrt(integrate(g,abs2(φ*efs)))'                 # normalize eigenfunctions

            @stats "compute density"            rho_out = abs2(φ*efs)*f                                 # compute density
                                                E = dot(f,evs)                                          # sum of eigenvalues
            @stats "Harris functional"          Etot = E + integrate(g, rho_in.*(ex+ec - (vx+vc) - vh/2)) + Enuc    # Harris functional

            @stats "check convergence"          scf(Etot,rho_out)                                       # check convergence
            @stats "density mixing"             rho_in = mixer(rho_in,rho_out)                          # mix densities

            if !is(logger,Void) logger(scf) end

        end  # scf cycle

        Ex   = integrate(g, rho_in.*ex)                         # exchange energy
        Ec   = integrate(g, rho_in.*ec)                         # correlation energy
        Exc  = Ex + Ec                                          # exchange-correlation energy
        Eh   = integrate(g, rho_in.*vh/2)                       # Hartree energy
        Epot = integrate(g, rho_in.*vext)                       # potential energy (external)
        Ekin = Etot - Epot - Eh - Ex - Ec                       # kinetic energy

        # functions for visualization
        rho_free_ = r -> sum([at[:rho](r) for at in atoms])     # free atom density
        rho_      = r -> (abs2(basis(r)'*efs)*f)[]              # total density
        vh_free_  = r -> sum([at[:vh](r) for at in atoms])      # free atom hartree potential
        vh_       = r -> vh_free_(r) + sum(δvh_mp(r))           # hartree potential
        vx_       = r -> xc.x.v([rho_(r)])[]                    # exchange potential
        vc_       = r -> xc.c.v([rho_(r)])[]                    # correlation potential
        vxc_      = r -> vx_(r) + vc_(r)                        # XC potential
        vhxc_     = r -> vh_(r) + vxc_(r)                       # hartree + XC potential
        ex_       = r -> xc.x.e([rho(r)])[]                     # exchange energy density
        ec_       = r -> xc.c.e([rho(r)])[]                     # correlation energy density
        exc_      = r -> ex_(r) + ec_(r)                        # XC energy density

    end  # total

    struct(Etot=Etot,Ekin=Ekin,Epot=Epot,Eh=Eh,Ex=Ex,Ec=Ec,Exc=Exc,Enuc=Enuc,rho=rho_in,
        rho_=rho_,rho_free_=rho_free_,vs=vs,vh=vh,vh_=vh_,vh_free_=vh_free_,vx_=vx_,vc_=vc_,
        vxc_=vxc_,vhxc_=vhxc_,ex_=ex_,ec_=ec_,exc_=exc_,drho_mp=δrho_mp,dvh_mp=δvh_mp,
        evs=evs,efs=efs,f=f,mu=mu,stats=stats)
end

function dft!(m::Molecule; params...)
    merge!(m,Dict{Symbol,Any}(params))

    res = dft(
        m.atoms,                                                # atoms
        m.basis,                                                # molecular basis
        m.g,                                                    # integration grid
        get!(m, :N,      sum([at[:N] for at in m.atoms])),      # particle number
        get!(m, :vembed, zeros(length(m.g))),                   # embedding potential
        get!(m, :xc,     ldaxcunpol()),                         # xc functional
        get!(m, :scf,    SCF()),
        get!(m, :mixer,  PulayMixer()),                         # density mixer
        get!(m, :lmax,   4),                                    # max l-quantum number in hartree multipole expansion
        get!(m, :nh,     300),                                  # number of hartree radial points
        get!(m, :nspl,   100),                                  # number of spline points
        get!(m, :T,      1e-6),                                 # temperature
        get!(m, :logger, Void)
    )

    merge!(m,res)
    return
end

function dft(m::Molecule; params...)
    let m = deepcopy(m)
        dft!(m,params...)
        return m
    end
end

function dft(; name::AbstractString="", atoms::Atoms=Atom[], params...)
    m = Molecule(;name=name,atoms=atoms,params...)
    dft!(m)
    m
end

# sparse version
# function dft!(sys::Molecule; params...)
#     local E,Etot,Enuc,vs,vh,δρ_mp,δvh_mp,evs,efs,f,mu

#     merge!(sys, Dict{Symbol,Any}(params))

#     atoms = sys[:atoms]     ::Atoms
#     basis = sys[:basis]     ::Basis             # molecular basis
#     g     = sys[:g]         ::MultiGrid         # integration grid
#     stats = Stats()

#     @defaults sys begin
#         N::Float64           = :N in keys(sys) ? at[:N] : sum([at[:N] for at in atoms])  # number of particles
#         xc::ldaxc            = ldaxcunpol()                     # exchange-correlation functional
#         scf::SCF             = SCF()
#         mixer::AbstractMixer = PulayMixer()                     # density mixing
#         lmax::Int            = 4                                # max l-quantum number in hartree multipole expansion
#         num::Int             = 300                              # num of radial grid points to solve hartree potential
#         T::Float64           = 1e-6                             # temperature
#         logger               = nothing
#     end

#     if :grid in fieldnames(mixer) mixer.grid = g end
#     scf.grid = g
#     scf.converged = false

#     @stats "* total" begin

#         # initialize
#         @stats "* evaluate vext"                vext     = sum([at[:v](g) for at in atoms])                 # external potential
#         @stats "* evaluate rho_free"            rho_free = sum([at[:rho](g) for at in atoms])               # free atom density
#         @stats "* evaluate vh_free"             vh_free  = sum([at[:vh](g) for at in atoms])                # free atom hartree potential
#                                                 rho_in   = isempty(scf.rho) ? rho_free : scf.rho[end]       # initial guess for density

#         @stats "* evaluate basis"               φ,∆φ = eval_basis_v1(basis,g)                               # evaluate basis on grid
#         @stats "* overlap"                      ζ     = abs(φ) .< 1e-12
#         @stats "* overlap"                      φ[ζ]  = 0.
#         @stats "* overlap"                      ∆φ[ζ] = 0.
#         @stats "* overlap"                      φ     = sparse(φ)
#         @stats "* overlap"                      ∆φ    = sparse(∆φ,φ)
#         @stats "* overlap"                      wφT   = sparse((g.w.*φ)')
#         @stats "* overlap"                      S     = Symmetric(full(wφT*φ))

#         # storage
#         vx = similar(rho_in); vc = similar(rho_in)
#         ex = similar(rho_in); ec = similar(rho_in)
#         vh = similar(rho_in); vs = similar(rho_in)

#         while !scf.converged
#                                                 vh[:] = vh_free
#                                                 # δrho = rho_in - rho_free
#             # @stats "hartree"                    δρ_mp,δvh_mp = poissonmultipole!(vh,g,δrho,atoms;lmax=lmax,num=num)  # Hartree potential

#             @stats "exchange-correlation"       xc(rho_in,ex,vx,ec,vc)                              # exchange-correlation
#                                                 vs[:] = vext + vx + vc + vh                         # Kohn-Sham potential

#             @stats "hamiltonian"                hφ = h(φ,∆φ,vs)                                     # apply Hamiltonian to basis
#             @stats "hamiltonian"                H = Symmetric(full(wφT*hφ))                         # Hamiltonian matrix
#             @stats "KS equations"               evs,efs = LinAlg.eig(H,S)                           # compute eigenspectrum

#             @stats "KS equations"               efs ./= sqrt(integrate(g,abs2(φ*efs)))              # normalize eigenfunctions

#                                                 f,mu = compute_occupations(evs,2,N,T)

#             @stats "compute density"            rho_out = abs2(φ*efs)*f                             # output density
#                                                 E = dot(f,evs)                                      # sum of eigenvalues
#                                                 Enuc = sum([μ != ν ? μ[:Z]*ν[:Z]/dist(μ[:R],ν[:R])/2 : 0.           # nuclear attraction
#                                                     for μ in atoms, ν in atoms])
#             @stats "harris energy"              Etot = E + integrate(g, rho_in.*(ex+ec - (vx+vc) - vh/2)) + Enuc    # Harris functional

#             @stats "check convergence"          scf(Etot, rho_out)                                  # check convergence
#             @stats "density mixing"             rho_in = mixer(rho_in, rho_out)                     # mix densities

#             if !is(logger, nothing) logger(scf) end

#         end  # scf cycle

#         Ex = integrate(g, rho_in.*ex)                           # exchange energy
#         Ec = integrate(g, rho_in.*ec)                           # correlation energy
#         Exc = Ex + Ec                                           # exchange-correlation energy
#         Eh = integrate(g, rho_in.*vh/2)                         # Hartree energy
#         Epot = integrate(g, rho_in.*vext)                       # potential energy (external)
#         Ekin = Etot - Epot - Eh - Ex - Ec                       # kinetic energy

#         # functions for visualization
#         rho_      = r -> (abs2(basis(r)'*efs)*f)[]              # total density
#         rho_free_ = r -> sum([at[:rho](r) for at in atoms])     # free atom density

#     end  # total

#     res = struct(Etot=Etot,Ekin=Ekin,Epot=Epot,Eh=Eh,Ex=Ex,Ec=Ec,Exc=Exc,Enuc=Enuc,rho=rho_in,
#         rho_=rho_,rho_free_=rho_free_,vs=vs,vh=vh,evs=evs,efs=efs,f=f,mu=mu,stats=stats)

#     merge!(sys, res)
#     return
# end

ext: jl

max_stars_repo_name: johncsnyder/Element


text:
# defining constants
amu = 931.4940954E6   # eV / c^2
ħc = 0.19732697E-6    # eV m
boltz = 8.617333262E-5 # eV/K

export Potential, pot_from_dict, filter_sample_points!, fit_pot!, solve_pot!, find_crossing, pot_from_file, cleave_pot
export Plotter
# export solve1D_ev_amu
# export sqwell, harmonic, double_well, polyfunc, morse
# export get_bspline, get_spline

"""
Stores a potential in one-dimensional space Q, with discreet points (E0, Q0) and fitting function func.

## Fields

- `name` -- the name of Potential.
- `QE_data`   -- the (n X 2) DataFrame of data points (Q vs Energy).
- `E0`, `Q0`  -- the minimum point of the Potential [`Q0`, `E0`].
- `func_type`     -- the type of fitting function ("bspline", "spline", "harmonic", "polyfunc", "morse_poly", "morse").
- `params`    -- the list of hyper parameters for the fitting function.
- `Q`, `E`  -- `Q` and `E`=`func(Q, p_opt; params)`. They are not the same as QE_data, they are the data points of the fitted function.
- `nev`  -- the number of eigenvalues to be evaluated.
- `ϵ` -- the list of eigenvalues
- `T` -- temperature (only used for filtering sample points)
- `T_weight` -- turn on temperature-dependent Boltzman weighting

## Constructor
    Potential()

"""
mutable struct Potential
    name::String
    QE_data::DataFrame
    E0::Float64; Q0::Float64
    func_type::String
    func # can be either Function or Interpolations.ScaledInterpolation
    params::Dict{String, Any}
    Q::Array{Float64,1}; E::Array{Float64,1}
    nev::Int64
    # ϵ includes E0
    ϵ::Array{Float64,1}; χ::Array{Float64,2}
    # options for thermally accessible fit
    T::Float64
    T_weight::Bool
    Potential() = new("", DataFrame([0 0], [:Q, :E]), 0, 0,
                      "harmonic_fittable", x->0, Dict(),
                      [], [],
                      0, [], Array{Float64}(undef, 0, 2), 293.15,false)
end

# extend copy to work with potentials, a cleaner implementation would be welcome
function Base.copy(pot::Potential)
    copy_pot = Potential()
    copy_pot.name = pot.name
    copy_pot.QE_data = pot.QE_data
    copy_pot.E0 = pot.E0
    copy_pot.Q0 = pot.Q0
    copy_pot.func_type = pot.func_type
    copy_pot.func = pot.func
    copy_pot.params = pot.params
    copy_pot.Q = pot.Q
    copy_pot.E = pot.E
    copy_pot.nev = pot.nev
    copy_pot.ϵ = pot.ϵ
    copy_pot.χ = pot.χ
    copy_pot.T = pot.T
    return copy_pot
end

"""
    pot_from_file(filename::String, resolution::Int64 = 3001)

Parse a two column file with data for reaction coordinate and energy.
Lines beginning with a hash are ignored. The optional argument `resolution`
determines how many points are included in the interpolation between the
read-in data points, pot.Q.

Additionally, any pot.QE_data.Q not included in pot.Q are added.

"""
function pot_from_file(filename::String, resolution::Int64 = 3001)
    Q_dat = Float64[]
    E_dat = Float64[]
    # read each line
    for line in eachline(filename)
        # skip comments
        if line[begin] != '#'
            line_split=split(line)
            append!(Q_dat,parse(Float64,line_split[1]))
            append!(E_dat,parse(Float64,line_split[2]))
        end
    end
    pot = Potential()
    pot.QE_data = DataFrame(Q = Q_dat[:], E = E_dat[:])
    # range of points
    interpol_points = range(Q_dat[1], stop=Q_dat[end], length=resolution)
    # add sample points from QE_data, sort and remove duplicates
    pot.Q = unique(sort(vcat(interpol_points, Q_dat[:])))

    # default harmonic fittable
    pot.func_type = "spline"

    return pot
end

"""
    filter_sample_points!(pot::Potential)

Remove all fit points after a thermally unsurmountable point in increasing order.
This function is only robust if

"""
function filter_sample_points!(pot::Potential, thermal_energy::Float64)
    # find minima
    min_ind_array = findall(pot.QE_data.E .== minimum(pot.QE_data.E))
    len_min_ind = length(min_ind_array)
    # check if there are multiple maxima
    if len_min_ind > 1
        throw(ArgumentError("The fitted potential has several minima."))
    end

    min_ind = len_min_ind[1]

    # select the interconnected data points below the energy threshold and
    # containing the minimum energy point we call these points an island
    below_thresh = pot.QE_data.E .<= thermal_energy + pot.E0

    island = []

    for (ind, bool) in enumerate(below_thresh)
        # if the data point is not part of an island
        if bool == 0
            # if there is a non-0 island in memory, add it and clear memory
            if length(island) > 0
                if min_ind in island
                    break
                else
                island = []
                end
            end
        # if the point is part of an island
        else
            # add point to island
            append!(island, ind)
        end
    end

    pot.QE_data = DataFrame(Q = pot.QE_data.Q[island], E = pot.QE_data.E[island])

end

"""
    function filter_sample_points!(pot::Potential, thermal_energy::Float64)

Same as above but use the temperature associated with the potential to filter points

"""
function filter_sample_points!(pot::Potential)
    filter_sample_points!(pot, pot.T*boltz)
end

"""
    fit_pot!(pot::Potential)

Fit a function `pot.func_type` to `pot.QE_data` on the domain `pot.Q`.

## parameters

- `pot`: `Potential`
    - `pot.func_type`: the fitting function; `{"spline" (preferred), "bspline", "harmonic", "polyfunc", "morse_poly", "morse"}`.
    - `pot.params`: the hyperparameters.

### Hyperparameters `params`

- spline (preferred)

    Spline interpolation. See more detail in [Dierckx.jl](https://github.com/kbarbary/Dierckx.jl).
    - `order`: spline order (between 1 and 5; default 2).
    - `smoothness`: the amount of smoothness is determined by the condition that `sum((w[i]*(y[i]-spline(x[i])))**2) <= s`
    - `weight`: the weight applied to each `QE_data` point (length m 1-d array).

- bspline

    Basic spline interpolation. See more detail in [Interpolations.jl](https://github.com/JuliaMath/Interpolations.jl).

- harmonic

    Harmonic function whose minimum is at [`pot.Q0`, `pot.E0`].
    - `hw`: the energy quanta of the harmonic oscillator.

- polyfunc

    Polynomial function;
        `y = E₀ + Σ coeffs[i].* (x .- Q₀) .^(i-1)`.
    - `poly_order`: the maximum order of polynomials.
    - `p0`: the initial parameters for the fitting function.

## Example
- Spline fit
```julia
nev = 60
name = "DX-+h"

Q1 = [30.66721918 29.860133 29.05306268 28.24612992 27.43911044 26.63197583 25.82507425 25.01797262 24.21096115 23.40387798 22.59690043 21.78991117 20.98281579 20.17570172 19.36884219 18.56169249 17.75463179 16.94772679 16.14061031 15.33347439 14.52663309 13.71945696 12.91240658 12.10544441 11.29841968 10.4913338 9.684370388 8.877289725 8.070184138]
E1 = [8.0902 7.5970 7.0749 6.5242 5.9461 5.3451 4.7300 4.1147 3.5182 2.9637 2.4769 2.0819 1.7972 1.6315 1.5800 1.6237 1.7301 1.8586 1.9687 2.0283 2.0190 2.0673 1.9910 2.0528 1.9730 2.0857 2.4550 3.1653 4.3448]

pot = Potential(); pot.name = name
pot.nev = nev
pot.Q0 = Q1[findmin(E1)[2]]; pot.E0 = 1.69834
pot.QE_data = DataFrame(Q = Q1[:], E = E1[:])
pot.QE_data.E .+= - minimum(pot.QE_data.E) + pot.E0
pot.Q = Q

pot.func_type = "spline"
# spline fitting parameters
params = Dict()
params["order"] = 4
params["smoothness"] = 0.001
params["weight"] = [1 1 1 1 1 0.5 0.4 0.4 0.5 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1]
fit_pot!(pot)
```

- Harmonic
```julia
nev = 40
name = "SnIII+h"

pot = Potential(); pot.name = name
pot.Q0 = 1.90291674728; pot.E0 = 0.585005
pot.nev = nev
pot.func_type = "harmonic"
params = Dict()
params["hw"] = 0.0281812646475
fit_pot!(pot)
```

"""
function fit_pot!(pot::Potential)
    # pot.params["E0"] = pot.E0
    # pot.params["Q0"] = pot.QE_data.Q[findmin(pot.QE_data.E)[2]]

    E_CUT = 2 # defines upper energy limit on data used for fitting
    e_cut_ind = pot.QE_data.E .< E_CUT + pot.E0 # boolean

    println("Potential fitting: $(pot.name)")

    if pot.func_type == "bspline"
        println("=========bspline=========\n")

        bspline = get_bspline(pot.QE_data.Q, pot.QE_data.E)
        pot.E = bspline(pot.Q)
        pot.func = bspline

    elseif pot.func_type == "spline"
        println("=========spline==========\n")

        weight = get(pot.params, "weight", nothing)
        if isa(weight, String)
            weight = parse.(Float64, split(weight))
        end
        if weight ≠ nothing
            weight = vec(weight)
        end
        smoothness = get(pot.params, "smoothness", 0)
        order = get(pot.params, "order", 2)

        spline = get_spline(pot.QE_data.Q, pot.QE_data.E; weight = weight, smoothness = smoothness, order = order)
        pot.E = spline(pot.Q)
        pot.func = spline

    elseif pot.func_type == "harmonic"
        println("========harmonic=========\n")
        # func = x -> harmonic(x, params["hw"]; param = params)
        pot.E = harmonic.(pot.Q, pot.params["hw"]; E₀ = pot.E0, Q₀ = pot.Q0)
        pot.func = x -> harmonic(x, pot.params["hw"]; E₀ = pot.E0, Q₀ = pot.Q0)

    # potentials that need fitting
    else
        if pot.func_type == "polyfunc"
            println("========polynomial========\n")
            poly_order = get(pot.params, "poly_order", 4)
            poly_params = get(pot.params, "p0", zeros(poly_order))
            func = (x, p) -> polyfunc(x, p; E₀ = pot.E0, Q₀ = pot.Q0, poly_order = poly_order)
        elseif pot.func_type == "morse_poly"
            println("========morse polynomial========\n")
            poly_params = get(pot.params, "p0", zeros(poly_order))
            poly_order = get(pot.params, "poly_order", 4)
            func = (x, p) -> morse_poly(x, p; E₀ = pot.E0, Q₀ = pot.Q0, poly_order = poly_order)
        elseif pot.func_type == "morse"
            poly_params = zeros(2)
            func = (x, p) -> morse(x, p; E₀ = pot.E0, Q₀ = pot.Q0)
        elseif pot.func_type == "harmonic_fittable"
            println("========harmonic fit========\n")
            poly_params = zeros(1)
            func = (x, p) -> harmonic_fittable(x, p; E₀ = pot.E0, Q₀ = pot.Q0)
        end
        if pot.T_weight
            # Boltzman weighting
            weights = exp.((-(pot.QE_data.E[e_cut_ind].-minimum(pot.QE_data.E[e_cut_ind]))./(pot.T*boltz)))
            print("Weights: ", weights)
            fit = curve_fit(func, pot.QE_data.Q[e_cut_ind], pot.QE_data.E[e_cut_ind], weights, Float64.(poly_params)) # curve_fit : model, x data, y data, weights, p0
        else
            fit = curve_fit(func, pot.QE_data.Q[e_cut_ind], pot.QE_data.E[e_cut_ind], Float64.(poly_params)) # curve_fit : model, x data, y data, p0
        end

        pot.E = func.(pot.Q, Ref(fit.param)) # when calling function, need to Ref() so that it can deal with array mismatch
        pot.func = x -> func(x, fit.param)

        println("===========Fit===========")
        println("Function: $(pot.func_type)")
        println("Best fit: $(fit.param)")
        println("=========================\n")
    end
end

"""
    solve_pot!(pot::Potential)

Solve 1D Shrödinger equation for `Potential`. `pot.ϵ` and `pot.χ` store the eigenvalues and eigenvectors, respectively.
"""
function solve_pot!(pot::Potential)
    # pot.ϵ, pot.χ = solve1D_ev_amu(pot.func, pot.Q; nev=pot.nev)
    # Q = pot.Q
    # NQ=length(Q); Qi=minimum(Q); Qf=maximum(Q)
    pot.ϵ, pot.χ = solve1D_ev_amu(pot.func, pot.Q; nev=pot.nev)
end

"""
Solve 1D Shrödinger equation. The Brooglie wrapper in the unit of `eV` and `amu`.
"""
function solve1D_ev_amu(func, Q::AbstractArray{<:Real,1}; nev=30, maxiter=nev*length(Q))
    NQ=length(Q); Qi=minimum(Q); Qf=maximum(Q)
    return solve1D_ev_amu(func, NQ, Qi, Qf; nev=nev, maxiter=maxiter)
end

function solve1D_ev_amu(func, NQ::Int, Qi::Real, Qf::Real; nev=30, maxiter=nev*length(Q))
    factor = (1/amu) * (ħc*1E10)^2

    ϵ1, χ1 = Brooglie.solve(x -> func.(x*factor^0.5);
                  N=NQ, a=Qi/factor^0.5, b=Qf/factor^0.5, m=1, nev=nev, maxiter=maxiter)
    return ϵ1, vcat(χ1'...)/factor^0.25
end

"""
    find_crossing(pot_1::Potential, pot_2::Potential)

Find the crossing point between two potential energy surfaces `pot1` and `pot2`.
`Qx, Ex = find_crossing(pot1, pot2)`.
"""
function find_crossing(pot_1::Potential, pot_2::Potential)
    # find root of pot_1 - pot_2 = 0
    diff_func = x-> pot_1.func(x) - pot_2.func(x)
    Q = pot_1.Q
    rts = find_zero(diff_func, Q[length(Q)÷2]) # start search at the midpoint of the Q vector
    return rts, pot_1.func.(rts)
end

function sqwell(x, width, depth; x0=0.)
    x = x .- x0
    pot_well = (x .< -width/2.) .| (x .> width/2.)
    return pot_well .* depth
end

function harmonic(x, ħω::Real; E₀, Q₀)
    a = amu/2 * (ħω/ħc/1E10)^2
    return a*(x - Q₀)^2 + E₀
end

function harmonic_fittable(x, coeff; E₀, Q₀)
    # an alternate syntax for harmonic() in order to make it fittable
    y = (0 .* x) .+ E₀ .+ coeff[1] * (x .- Q₀) .^2
    return y
end

function double_well(x, ħω1::Real, ħω2::Real; param)
    a1 = amu / 2. * (ħω1/ħc/1E10)^2
    a2 = amu / 2. * (ħω2/ħc/1E10)^2
    return - a1*x.*x + a2*x.*x.*x.*x
end

# Set of potentials for fitting
function polyfunc(x, coeffs; E₀::Real, Q₀::Real, poly_order::Int)
    # E0 = param["E0"]
    # Q₀ = param["Q0"]
    # poly_order = Int(param["poly_order"])
    y = 0 .* x # array of size x
    y = y .+ E₀

    for i = 2:poly_order + 1
        y = y .+ coeffs[i].* (x .- Q₀) .^(i-1)
    end

    return y
end

function morse(x, coeffs; E₀, Q₀)
    y = 0 .* x # array of size x
    y = y .+ E₀

    y = y .+ coeffs[1] .* ((1 .- exp.(-coeffs[2].*(x.-Q₀))).^2)

    return y
end

function morse_poly(x, coeffs; E₀, Q₀, poly_order)
    # TODO: CHECK # of coeffs and param["poly_order"], and give a WARNING
    # E₀ = param["E0"]
    # Q₀ = param["Q0"]
    # po = param["poly_order"]

    orders = if (typeof(po) == Int64) Array([po]) else parse.(Int, split(po)) end

    A = abs(coeffs[1])
    a = coeffs[2]

    r₀ = coeffs[3]

    morse = x -> A*(exp(-2*a*(x-Q₀-r₀)) - 2*exp(-a*(x-Q₀-r₀)))

    coeff = zeros(maximum(orders)+1)
    for (index, poly_order) in enumerate(orders)
        coeff[poly_order+1] = coeffs[3+index]
    end

    coeff[end] = abs(coeff[end])

    poly = Poly(coeff)
    tot_derv = polyder(poly, 1) -  A*(-2a)*(exp(2a*r₀) - exp(a*r₀))

    r₁ = Polynomials.roots(tot_derv)
    r₁ = real.(r₁[imag.(r₁) .== 0])
    r₁ = r₁[argmin(abs.(r₁))]

    return morse(x) + poly(x-Q₀-r₁) + E₀ - morse(Q₀) - poly(-r₁)
end

function get_spline(Qs, Es; weight = nothing, smoothness = -1, order = 2)
    if Qs[1] > Qs[end]
        Qs = reverse(Qs)
        Es = reverse(Es)
    end

    if weight == nothing
        spl = Spline1D(Qs, Es, k = order, bc = "extrapolate", s = smoothness)
    else
        spl = Spline1D(Qs, Es, w = weight, k = order, bc = "extrapolate", s = smoothness)
    end
    return spl
end

function get_bspline(Qs, Es)
    # BSplines assume your data is uniformly spaced on the grid
    # Qs, Es have to be equally-spaced (Range)
    if Qs[1] > Qs[end]
        Qs = reverse(Qs)
        Es = reverse(Es)
    end
    Qs_range = range(Qs[1], stop = Qs[end], length = length(Qs))

    itp = interpolate(Es, BSpline(Quadratic(Line(OnCell()))))
    etpf = extrapolate(itp, Line())
    setpf = scale(etpf, Qs_range)

    return setpf
end

"""
    function cleave_pot(pot::Potential)

Return two potentials with the `QE_data` of the original potential split by the apex.
If there is more than one maximum, throw an error

"""
function cleave_pot(pot::Potential; discard_max=false)
    # find maxima
    max_ind_array = findall(pot.E .== maximum(pot.E))
    len_max_ind = length(max_ind_array)
    # check if there are multiple maxima
    if len_max_ind > 1
        throw(ArgumentError("The fitted potential has several maxima."))
    end

    # cast unique maximum to int
    max_ind = max_ind_array[1]

    # find the corresponding Q value
    Q_limit = pot.Q[max_ind]

    # define two sub-potentials
    pot_a = copy(pot)
    pot_a.E = []
    pot_a.func = x->0

    pot_b = copy(pot_a)

    # assign corresponding data points
    Q_data = pot.QE_data.Q
    E_data = pot.QE_data.E

    # discard maxima
    if discard_max
        mask = E_data.<maximum(E_data)
        E_data = E_data[mask]
        Q_data = Q_data[mask]
    end

    pot_a.QE_data = DataFrame(Q = Q_data[Q_data.<Q_limit], E = E_data[Q_data.<Q_limit])
    pot_b.QE_data = DataFrame(Q = Q_data[Q_data.>Q_limit], E = E_data[Q_data.>Q_limit])

    # assign the ground states as initial and final points
    pot_a.E0 = pot.QE_data.E[1]
    pot_a.Q0 = pot.QE_data.Q[1]

    pot_b.E0 = pot.QE_data.E[end]
    pot_b.Q0 = pot.QE_data.Q[end]

    # make cleft potentials harmonic
    pot_a.func_type = "harmonic_fittable"
    pot_b.func_type = "harmonic_fittable"

    return pot_a, pot_b
end


text:
module LTC

using Reexport

using Distributions
using DiffEqBase
using OrdinaryDiffEq
using DiffEqSensitivity
import DiffEqFlux: initial_params, paramlength, FastChain, FastDense, sciml_train
using GalacticOptim
using ModelingToolkit
using Flux
using NNlib: sigmoid

rand_uniform(TYPE, lb,ub,dims...) = TYPE.(rand(Uniform(lb,ub),dims...))
rand_uniform(TYPE, lb,ub) = rand_uniform(TYPE, lb,ub,1)[1]

#Zygote.@nograd rand_uniform, reshape

include("layers.jl")
include("mtk_recur.jl")
include("optimization.jl")
include("losses.jl")
include("variables.jl")
# include("mkt_sysstruct.jl")
# include("zygote.jl")

include("systems/systems.jl")
include("systems/ncp/ncp_sys_gen.jl")
include("systems/ncp/wiring.jl")

export MTKRecur, MTKCell, Mapper, Broadcaster, get_bounds
export initial_params, paramlength
export reset_state!
export optimize
export loss_seq

export Wiring, FWiring, Net
end


text:
"""
    qc_age_consistency(data, age_var)

Identifies cases where the patient's age is less than ages of onset
"""
function qc_age_consistency(data::DataFrame, age_var::Symbol)
    df = data[Not(ismissing.(data.age .- data[:, age_var])),:]
    df = df[df.age .< df[:, age_var], :]
    df = df[:, [:family_id, :person_id, :centre]]
    #df[:,"Problem"] .= missing
    if age_var == :age_of_onset
        df[:,"Problem"] .= "Age < Age of Onset"
    elseif age_var == :onset_d
        df[:, "Problem"] .= "Age < Onset D"
    elseif age_var == :onset_m
        df[:, "Problem"] .= "Age < Onset Hypomania"
    elseif age_var == :onset_M
        df[:, "Problem"] .= "Age < Onset Mania"
    elseif age_var == :age_first_suicide_attempt
        df[:, "Problem"] .= "Age < Age at first SA"
    end
    return df
end

"""
    filter_qc_exclusions(data, qc_exclusions)

Removes the subjects whose data were flagged in qc_exclusions.
"""
function filter_qc_exclusions(data::DataFrame, qc_exclusions::DataFrame)
    for row ∈ eachrow(qc_exclusions)
        fam = data.family_id .== row[1]
        pers = data.person_id .== row[2]
        site = data.centre .== row[3]
        data = data[Not(fam .& pers .& site),:]
    end
    return data
end


text:
"""
    calc_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{FreeSlip}, ::ModelSetup)
apply free slip boundary condition for velocity
apply zero numerical flux in the normal direction
"""
function calc_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{FreeSlip},
    ::ModelSetup,
    ::Nothing,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρu = state⁻.ρu

    return nothing
end

function calc_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{FreeSlip},
    ::ModelSetup,
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρu = state⁻.ρu
    gradflux⁺.ν∇u = n⁻ * (@SVector [-0, -0, -0])'

    return nothing
end

"""
    calc_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{NoSlip}, ::ModelSetup)
apply no slip boundary condition for velocity
sets ghost point to have no numerical flux on the boundary for U
"""
@inline function calc_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{NoSlip},
    ::ModelSetup,
    ::Nothing,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρu = -state⁻.ρu

    return nothing
end

@inline function calc_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Impenetrable{NoSlip},
    ::ModelSetup,
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρu = -state⁻.ρu
    gradflux⁺.ν∇u = gradflux⁻.ν∇u

    return nothing
end

"""
    calc_boundary_state!(::NumericalFluxSecondOrder, ::Insulating, ::HBModel)

apply insulating boundary condition for velocity
sets ghost point to have no numerical flux on the boundary for κ∇θ
"""
@inline function calc_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Insulating,
    ::ModelSetup,
    ::Nothing,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρθ = state⁻.ρθ

    return nothing
end

@inline function calc_boundary_state!(
    ::NumericalFluxSecondOrder,
    ::Insulating,
    ::ModelSetup,
    ::ConstantViscosity,
    state⁺,
    gradflux⁺,
    hyperflux⁺,
    aux⁺,
    n⁻,
    state⁻,
    gradflux⁻,
    hyperflux⁻,
    aux⁻,
    t,
    _...,
)
    state⁺.ρθ = state⁻.ρθ
    gradflux⁺.κ∇θ = n⁻ * -0

    return nothing
end

# """
#     calc_boundary_state!(::NumericalFluxSecondOrder, ::Penetrable{FreeSlip}, ::ModelSetup)
# apply free slip boundary condition for velocity
# apply zero numerical flux in the normal direction
# """
# @inline function calc_boundary_state!(
#     ::NumericalFluxSecondOrder,
#     ::Penetrable{FreeSlip},
#     ::ModelSetup,
#     ::ConstantViscosity,
#     state⁺,
#     gradflux⁺,
#     aux⁺,
#     n⁻,
#     state⁻,
#     gradflux⁻,
#     aux⁻,
#     t,
#     args...,
# )
#     state⁺.ρu = state⁻.ρu
#     gradflux⁺.ν∇u = n⁻ * (@SVector [-0, -0, -0])'

#     return nothing
# end

# """
#     calc_boundary_state!(::NumericalFluxSecondOrder, ::Impenetrable{MomentumFlux}, ::HBModel)
# apply kinematic stress boundary condition for velocity
# sets ghost point to have specified flux on the boundary for ν∇u
# """
# @inline function calc_boundary_state!(
#     ::NumericalFluxSecondOrder,
#     bc::Impenetrable{<:MomentumFlux},
#     model::ModelSetup,
#     ::ConstantViscosity,
#     state⁺,
#     gradflux⁺,
#     aux⁺,
#     n⁻,
#     state⁻,
#     gradflux⁻,
#     aux⁻,
#     t,
#     args...,
# )
#     state⁺.ρu = state⁻.ρu
#     gradflux⁺.ν∇u = n⁻ * bc.drag.flux(state⁻, aux⁻, t)'

#     return nothing
# end

# """
#     calc_boundary_state!(::NumericalFluxSecondOrder, ::Penetrable{MomentumFlux}, ::HBModel)
# apply kinematic stress boundary condition for velocity
# sets ghost point to have specified flux on the boundary for ν∇u
# """
# @inline function calc_boundary_state!(
#     ::NumericalFluxSecondOrder,
#     bc::Penetrable{<:MomentumFlux},
#     shallow::ModelSetup,
#     ::ConstantViscosity,
#     state⁺,
#     gradflux⁺,
#     aux⁺,
#     n⁻,
#     state⁻,
#     gradflux⁻,
#     aux⁻,
#     t,
#     args...,
# )
#     state⁺.ρu = state⁻.ρu
#     gradflux⁺.ν∇u = n⁻ * bc.drag.flux(state⁻, aux⁻, t)'

#     return nothing
# end

# """
#     calc_boundary_state!(::NumericalFluxSecondOrder, ::TemperatureFlux, ::HBModel)

# apply insulating boundary condition for velocity
# sets ghost point to have specified flux on the boundary for κ∇θ
# """
# @inline function calc_boundary_state!(
#     ::NumericalFluxSecondOrder,
#     bc::TemperatureFlux,
#     ::ModelSetup,
#     state⁺,
#     gradflux⁺,
#     aux⁺,
#     n⁻,
#     state⁻,
#     gradflux⁻,
#     aux⁻,
#     t,
# )
#     state⁺.ρθ = state⁻.ρθ
#     gradflux⁺.κ∇θ = n⁻ * bc.flux(state⁻, aux⁻, t)

#     return nothing
# end

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
using DynamicAxisWarping, Plots

fs = 70
t  = range(0,stop=1,step=1/fs)
y0 = sin.(2pi .*t)
y1 = sin.(3pi .*t)
y  = [y0;y1[2:end]] .+ 0.01 .* randn.()
q  = [y0;y0[2:end]] .+ 0.01 .* randn.()

y[10:15] .+= 0.5
q[13:25] .+= 0.5

##
f1 = plot([q y])
f2 = dtwplot(q,y,lc=:green, lw=1)
f3 = matchplot(q,y,ds=3,separation=1)
plot(f1,f2,f3, legend=false, layout=3, grid=false)

##

svec = exp10.(LinRange(-3, 1, 500))
res = map(svec) do s
    y  = [y0;y1[2:end]] .+ s .* randn.()
    q  = [y0;y0[2:end]] .+ s .* randn.()
    c,i1,i2 = dtw(q,y, transportcost=1.05)
    ni = length(y)
    length(i1) - findfirst(==(ni), i1) + 1
end
scatter(svec, res, xscale=:log10, xlabel="Noise std", ylabel="N points matched to last")


text:
using AtomEnergyLevels
using PyPlot
pygui(true);

function cobalt()
    E = []
    for q = 7:0.1:9
        push!(E, lda(27, conf = conf_enc("[Ar] 3d$q 4s$(9 - q)")).energy.total)
    end
    return E
end

E = cobalt()

begin
    title(L"Energy of the cobalt atom with configuration $3d^q\,4s^{9-q}$")
    xlabel(L"q")
    ylabel(L"E - E_{min}")

    ax1 = PyPlot.axes()
    ax1.set_xlim([7,9])

    plot(7:0.1:9, E .- minimum(E))
    grid("on")
end


text:
"""
Zhu, F., Ye, N., Yu, W., Xu, S. and Li, G., 2014.
Boundary detection and sample reduction for one-class support vector machines.
Neurocomputing, 123, pp.166-173.

Neighbors’ Distribution Properties and Sample Reduction algorithm (NDPSR)
Parameter guideline from the paper, given a data set of size N:
k = 20
none for eps (on benchmark data eps ∈ [12-20])
"""
mutable struct NDPSR <: Sampler
    k::Int
    eps::Real
    function NDPSR(k::Int=20, eps::Real=15)
        k > 0 || throw(ArgumentError("Invalid number of neighbors to consider $(k)."))
        eps >= 0 || throw(ArgumentError("Invalid threshold $(eps)."))
        return new(k, eps)
    end
end

function sample(sampler::NDPSR, data::DataSet, labels::Vector{Symbol})::BitArray{1}
    sampler.k = min(size(data, 2) - 1, sampler.k)
    tree = KDTree(data)
    sample_mask = falses(size(data, 2))

    for i in 1:size(data, 2)
        idx, dist = knn(tree, data[:, i], sampler.k + 1, true)
        x_k = mean(data, dims=2)[:, 1]
        d_ik = normalize(data[:, i] .- x_k)
        d_ij = mapslices(normalize, data[:, i] .- data[:, idx[2:end]], dims=1)
        c_sum0 = sum(dot(d_ik, d_ij[:, j]) for j in 1:sampler.k)
        if c_sum0 >= sampler.eps
            sample_mask[i] = true
        end
    end
    return sample_mask
end


text:
using Batched, LinearAlgebra, Test

batched_tr(rand(2, 2, 10))


text:
#=
  Changepoint detection
  From https://nishanthu.github.io/articles/ChangePointAnalysis.html

  See ~/jags/change_point.R
  The JAGS model
  model {
     for (i in 1:N) {
       mu0[i] = ifelse(i < cp, mu1, mu2)
       x[i] ~ dnorm(mu0[i], 1)
     }
     mu1 ~ dunif(0, mu)
     mu2 ~ dunif(mu, 10*mu)
     cp ~ dunif(0, N)
   }
"""
out$cp %>% as.vector %>% mean %>% round(1)
[1] 13.5
"""

   Compare with
     - changepoint_coal_miners.jl
     - changepoint_text_messages.jl

=#

# using Memoization
using Turing, StatsPlots, DataFrames
using ReverseDiff, Zygote, Tracker
using Printf
# Turing.setadbackend(:reversediff)
# Turing.setadbackend(:zygote)
# Turing.setadbackend(:tracker)

include("jl_utils.jl")

@model function changepoint(ys)
    N = length(ys)
    mu = mean(ys)

    mu1 ~ Uniform(0, mu)
    mu2 ~ Uniform(mu, 10*mu)
    cp ~ DiscreteUniform(1,N)

    for i in 1:N
        ys[i] ~ i < cp ? Normal(mu1,1) : Normal(mu2,1)
    end

end

# Generated by R:
# > c(rnorm(13, mean=0, sd=0.5), rnorm(20, mean=4, sd=0.5))
# ys = [0.90685440,0.38219559,0.13508408,-0.96054590,0.78553871,0.06791336,-0.09750806,-0.85431057,
#      0.33478539,0.76388689,0.52911522,-0.41226421,1.53302740,
#      5.01118838,4.75164473,3.79149212,3.97544567,4.33985561,4.52511631,4.21723600,
#      4.22733898,3.80925885,3.65504091,3.25267865,3.85339052,3.40136529,4.21462117,
#      3.94033196,4.00900359,4.44148942,4.23449516,4.52032000,4.41335270];

# But we use Julia here...
# Note:
#  - cp should be 14
#  - mu1 should be 0
#  - mu2 should be 4
ys = vcat(rand(Normal(0,0.5),13), rand(Normal(4,0.5),20) )

# Larger instance (130 + 200 elements)
# cp: 131
# mu1: 1
# mu2: 14
# ys = vcat(rand(Normal(1,0.5),130), rand(Normal(14,0.5),200) )

model = changepoint(ys)

num_chns = 4

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

# chns = sample(model, MH(), MCMCThreads(), 40_000, num_chns)
chns = sample(model, MH(), 100_000)

# chns = sample(model, SMC(), 10_000)
# chns = sample(model, PG(20), 1_000)

display(chns[[:cp,:mu1,:mu2]]) # We are not interested in :ys

show_var_dist_pct(chns, :cp)
show_var_dist_pct(chns, :mu1)
show_var_dist_pct(chns, :mu2)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using FromFile
using SymbolicUtils
using SymbolicUtils: Chain, If, RestartedChain, IfElse, Postwalk, Fixpoint, @ordered_acrule, isnotflat, flatten_term, needs_sorting, sort_args, is_literal_number, hasrepeats, merge_repeats, _isone, _iszero, _isinteger, istree, symtype, is_operation, has_trig, polynormalize
@from "Core.jl" import Options
@from "InterfaceSymbolicUtils.jl" import SYMBOLIC_UTILS_TYPES

function multiply_powers(eqn::T)::SYMBOLIC_UTILS_TYPES where {T<:Union{<:Number,SymbolicUtils.Sym{<:Number}}}
	return eqn
end

function multiply_powers(eqn::T, op::F)::SYMBOLIC_UTILS_TYPES where {F,T<:SymbolicUtils.Term{<:Number}}
	args = SymbolicUtils.arguments(eqn)
	nargs = length(args)
	if nargs == 1
		return op(multiply_powers(args[1]))
	elseif op == ^
		l = multiply_powers(args[1])
		n::Int = args[2]
		if n == 1
			return l
		elseif n == -1
			return 1.0 / l
		elseif n > 1
			return reduce(*, [l for i=1:n])
		elseif n < -1
			return reduce(/, vcat([1], [l for i=1:abs(n)]))
		else
			return 1.0
		end
	elseif nargs == 2
		return op(multiply_powers(args[1]), multiply_powers(args[2]))
	else
		return mapreduce(multiply_powers, op, args)
	end
end

function multiply_powers(eqn::T)::SYMBOLIC_UTILS_TYPES where {T<:SymbolicUtils.Term{<:Number}}
	op = SymbolicUtils.operation(eqn)
	return multiply_powers(eqn, op)
end

# Operators required for each rule:
function get_simplifier(binops::A, unaops::B) where {A,B}
    PLUS_RULES = [
       rule for (required_ops, rule) in [
       ((+,), @rule(~x::isnotflat(+) => flatten_term(+, ~x))),
       ((+,), @rule(~x::needs_sorting(+) => sort_args(+, ~x))),
       ((+,), @ordered_acrule(~a::is_literal_number + ~b::is_literal_number => ~a + ~b)),
       ((+,), @acrule(*(~~x) + *(~β, ~~x) => *(1 + ~β, (~~x)...))),
       ((+,), @acrule(*(~α, ~~x) + *(~β, ~~x) => *(~α + ~β, (~~x)...))),
       ((+,), @acrule(*(~~x, ~α) + *(~~x, ~β) => *(~α + ~β, (~~x)...))),
       ((+,), @acrule(~x + *(~β, ~x) => *(1 + ~β, ~x))),
       ((+,), @acrule(*(~α::is_literal_number, ~x) + ~x => *(~α + 1, ~x))),
       ((+,), @rule(+(~~x::hasrepeats) => +(merge_repeats(*, ~~x)...))),
       ((+,), @ordered_acrule((~z::_iszero + ~x) => ~x)),
       ((+,), @rule(+(~x) => ~x))]
       if all([(op in binops || op in unaops) for op in required_ops])
    ]
    TIMES_RULES = [
       rule for (required_ops, rule) in [
       ((*,), @rule(~x::isnotflat(*) => flatten_term(*, ~x))),
       ((*,), @rule(~x::needs_sorting(*) => sort_args(*, ~x))),
       ((*,), @ordered_acrule(~a::is_literal_number * ~b::is_literal_number => ~a * ~b)),
       ((*,), @rule(*(~~x::hasrepeats) => *(merge_repeats(^, ~~x)...))),
       ((*,), @acrule((~y)^(~n) * ~y => (~y)^(~n+1))),
       ((*,), @ordered_acrule((~x)^(~n) * (~x)^(~m) => (~x)^(~n + ~m))),
       ((*,), @ordered_acrule((~z::_isone  * ~x) => ~x)),
       ((*,), @ordered_acrule((~z::_iszero *  ~x) => ~z)),
       ((*,), @rule(*(~x) => ~x))]
       if all([(op in binops || op in unaops) for op in required_ops])
    ]
    POW_RULES =[
       rule for (required_ops, rule) in [
       ((*,), @rule(^(*(~~x), ~y::_isinteger) => *(map(a->SymbolicUtils.pow(a, ~y), ~~x)...))),
       ((*,), @rule((((~x)^(~p::_isinteger))^(~q::_isinteger)) => (~x)^((~p)*(~q)))),
       ((*,), @rule(^(~x, ~z::_iszero) => 1)),
       ((*,), @rule(^(~x, ~z::_isone) => ~x)),
       ((*, /,), @rule(inv(~x) => ~x ^ -1))]
       if all([(op in binops || op in unaops) for op in required_ops])
    ]
    ASSORTED_RULES =[
       rule for (required_ops, rule) in [
       ((), @rule(identity(~x) => ~x)),
       ((*,), @rule(-(~x) => -1*~x)),
       ((*, -,), @rule(-(~x, ~y) => ~x + -1(~y))),
       ((/, *,), @rule(~x / ~y => ~x * SymbolicUtils.pow(~y, -1))),
       ((), @rule(one(~x) => one(symtype(~x)))),
       ((), @rule(zero(~x) => zero(symtype(~x))))]
       if all([(op in binops || op in unaops) for op in required_ops])
    ]
    TRIG_RULES = [
       rule for (required_ops, rule) in [
       ((sin, cos, *, +,), @acrule(sin(~x)^2 + cos(~x)^2 => one(~x))),
       ((sin, cos, *, +,), @acrule(sin(~x)^2 + -1        => cos(~x)^2)),
       ((sin, cos, *, +,), @acrule(cos(~x)^2 + -1        => sin(~x)^2)),
       ((tan, sec, *, +,), @acrule(tan(~x)^2 + -1*sec(~x)^2 => one(~x))),
       ((tan, sec, *, +,), @acrule(tan(~x)^2 +  1 => sec(~x)^2)),
       ((tan, sec, *, +,), @acrule(sec(~x)^2 + -1 => tan(~x)^2)),
       ((cot, csc, *, +,), @acrule(cot(~x)^2 + -1*csc(~x)^2 => one(~x))),
       ((cot, csc, *, +,), @acrule(cot(~x)^2 +  1 => csc(~x)^2)),
       ((cot, csc, *, +,), @acrule(csc(~x)^2 + -1 => cot(~x)^2))]
       if all([(op in binops || op in unaops) for op in required_ops])
    ]
    function number_simplifier()
        rule_tree = [If(istree, Chain(ASSORTED_RULES)),
                     If(is_operation(+),
                        Chain(PLUS_RULES)),
                     If(is_operation(*),
                        Chain(TIMES_RULES)),
                     If(is_operation(^),
                        Chain(POW_RULES))] |> RestartedChain

        rule_tree
    end
    trig_simplifier(;kw...) = Chain(TRIG_RULES)
    function default_simplifier(; kw...)
        IfElse(has_trig,
               Postwalk(Chain((number_simplifier(),
                               trig_simplifier())),
                        ; kw...),
               Postwalk(number_simplifier())
                        ; kw...)
    end
    # reduce overhead of simplify by defining these as constant
    serial_simplifier = If(istree, Fixpoint(default_simplifier()))
    serial_polynormal_simplifier = If(istree,
                                      Fixpoint(Chain((polynormalize,
                                                      Fixpoint(default_simplifier())))))
    return serial_polynormal_simplifier
end

function custom_simplify(init_eqn::T, options::Options)::SYMBOLIC_UTILS_TYPES where {T<:SYMBOLIC_UTILS_TYPES}
    if !istree(init_eqn) #simplifier will return nothing if not a tree.
        return init_eqn
    end
    simplifier = get_simplifier(options.binops, options.unaops)
    eqn = simplifier(init_eqn)::SYMBOLIC_UTILS_TYPES #simplify(eqn, polynorm=true)

	# Remove power laws
    eqn = multiply_powers(eqn::SYMBOLIC_UTILS_TYPES)
	return eqn
end


text:
module Completion

include("unify.jl")
include("complete.jl")
include("orders.jl")

end # module


text:
function validate(options::Comonicon)
    if options.sysimg !== nothing
        isabspath(options.sysimg.path) && error("system image path must be project relative")
    end

    if options.application !== nothing
        isabspath(options.application.path) && error("application build path must project relative")
    end
    return
end

"""
    find_comonicon_toml(path::String, files=["Comonicon.toml", "JuliaComonicon.toml"])

Find `Comonicon.toml` or `JuliaComonicon.toml` in given path.
"""
function find_comonicon_toml(path::String, files=["Comonicon.toml", "JuliaComonicon.toml"])
    # user input file path
    basename(path) in files && return path

    # user input dir path
    for file in files
        path = joinpath(path, file)
        if ispath(path)
            return path
        end
    end
    return
end

"""
    read_toml(path::String)

Read `Comonicon.toml` or `JuliaComonicon.toml` in given path.
"""
function read_toml(path::String)
    file = find_comonicon_toml(path)
    file === nothing && return Dict{String,Any}()
    return TOML.parsefile(file)
end

"""
    read_toml(mod::Module)

Read `Comonicon.toml` or `JuliaComonicon.toml` in given module's project path.
"""
function read_toml(mod::Module)
    return read_toml(pkgdir(mod))
end

function has_comonicon_toml(m::Module)
    path = pkgdir(m)
    isnothing(path) && return false
    !isnothing(find_comonicon_toml(path))
end

"""
    read_options(comonicon; kwargs...)

Read in Comonicon build options. The argument `comonicon` can be:

- a module of a Comonicon CLI project.
- a path to a Comonicon CLI project that contains either `JuliaComonicon.toml` or `Comonicon.toml`.
- a path to a Comonicon CLI build configuration file named either `JuliaComonicon.toml` or `Comonicon.toml`.

In some cases, you might want to change the configuration written in the TOML file temporarily, e.g for writing
build tests etc. In this case, you can modify the configuration using corresponding keyword arguments.

keyword arguments of [`Application`](@ref) and [`SysImg`](@ref) are the same, thus keys like `filter_stdlibs`
are considered ambiguous in `read_options`, but you can specifiy them by specifiy the specific [`Application`](@ref)
or [`SysImg`](@ref) object, e.g

```julia
read_options(MyCLI; sysimg=SysImg(filter_stdlibs=false))
```

See also [`Comonicon`](@ref), [`Install`](@ref), [`SysImg`](@ref), [`Application`](@ref),
[`Download`](@ref), [`Precompile`](@ref).
"""
function read_options(m::Union{Module,String}; kwargs...)
    d = read_toml(m)
    if !haskey(d, "name")
        d["name"] = default_cmd_name(m)
    end

    options = from_dict(Comonicon, d; kwargs...)
    validate(options)
    return options
end

default_cmd_name(m::Module) = lowercase(string(nameof(m)))
default_cmd_name(m) = error("command name is not specified")


text:
using PastaQ
using ITensors
using Random

Random.seed!(1234)

N = 4     # Number of qubits
depth = 10 # Depth of the quantum circuit

# Generate random quantum circuit built out of
# random 2-qubit unitaries
println("Random circuit of depth $depth on $N qubits:")
circuit = randomcircuit(N, depth)

# Define custom function to measure an observable, in this
# case a Pauli operator on `site`
function measure_pauli(ψ::MPS, site::Int, pauli::String)
  ψ = orthogonalize!(copy(ψ), site)
  ϕ = ψ[site]
  obs_op = gate(pauli, firstsiteind(ψ, site))
  T = noprime(ϕ * obs_op)
  return real((dag(T) * ϕ)[])
end

# pauli X on site 2
σx2(ψ::MPS) = measure_pauli(ψ, 2, "X")
σz(ψ::MPS) = [measure_pauli(ψ, j, "Z") for j in 1:length(ψ)]

# define the Circuit observer
obs = Observer([
  "χs" => linkdims,      # bond dimension at each bond
  "χmax" => maxlinkdim,  # maximum bond dimension
  "σˣ(2)" => σx2,        # pauli X on site 2
  "σᶻ" => σz,
])

# run the circuit
ψ = runcircuit(circuit; (observer!)=obs)

# collect the measurements
println("Bond dimensions at each layer:")
display(results(obs, "χs"))
println()

println("Maximum bond dimension at each layer:")
display(results(obs, "χmax")')
println()

println("⟨ψ|σˣ(2)|ψ⟩ at each layer:")
display(results(obs, "σˣ(2)"))
println()

println("⟨ψ|σᶻ(n)|ψ⟩ at each layer:")
display(results(obs, "σᶻ"))
println()


text:
import TimeSeries.TimeArray

struct APIResponse
    data::String
    http_resp::HTTP.Messages.Response
end

"""
Description

    The yahoo() method is a wrapper for downloading historical stock prices from Yahoo.

Usage

    AAPL = yahoo("AAPL)
    SPX = yahoo()

Method Signature(s)

    yahoo(data::ASCIIString="^GSPC")

Details

    The yahoo method takes a stock name in the form of a string and returns a TimeSeries.TimeArray data structure
    corresponding to the Yahoo Finance ticker. With no argument, the default historical time series is the S&P 500.

References

    http://www.finance.yahoo.com

See Also

    fred() which accesses the St. Louis Federal Reserve financial and economic data sets.
"""
function yahoo(data::String="^GSPC")
    Base.depwarn("Yahoo Finance API changed, this function may not work anymore", :yahoo)
    url = "http://ichart.yahoo.com/table.csv?s=$data"
    http_resp = HTTP.request("GET", url)
    resp = APIResponse(data, http_resp)
    TimeArray(resp)
end

"""
Description

    The fred() method is a wrapper to download financial and economic time series data from the St. Louis Federal Reserve (FRED).

Usage

    DGS = fred("DGS10")
    CPI = fred()

Method Signature(s)

    fred(data::String="CPIAUCNS")

Details

    The fred() method takes a string argument that corresponds to a series code from the St. Louis Federal
    Reserve (FRED) database. It returns the data in the TimeSeries.TimeArray data structure.  When no argument
    is provided, the default data set is the Consumer Price Index for All Urban Consumers: All Items (CPIAUCNS).

References

    https://research.stlouisfed.org/fred2

See Also

    yahoo() which is a wrapper from downloading financial time series for stocks from Yahoo Finance.
"""
function fred(data::String="CPIAUCNS")
    url = "http://research.stlouisfed.org/fred2/series/$data/downloaddata/$data.csv"
    http_resp = HTTP.request("GET", url)
    resp = APIResponse(data, http_resp)
    TimeArray(resp)
end

function TimeArray(resp::APIResponse)
    #This function transform the Response object into a TimeArray
    # Split the data on every "\n"
    raw_data = String(resp.http_resp.body)
    data = split(raw_data, "\n")
    # Extract the head and body of the data
    head = strip(data[1])
    body = data[2:end]
    # Parse body
    body[end] == "" ? pop!(body) : nothing # remove trailing empty string if it's there
    body      = [split(line, ",") for line in body] # split on comma
    ######### Timestamp
    # take the first row (assuming it's date)
    # TODO: regex query needed to catch edge cases
    dates     = [line[1] for line in body]
    timestamp = Date[Date(d) for d in dates] # parse dates
    ######### Values
    svals = [line[2:end] for line in body] # get rows 2 to the end
    fvals = zeros(length(svals),length(svals[1]))
    for r in 1:size(fvals,1)
        for c in 1:size(fvals,2)
            # is not empty and is not equal to FRED's iconic "." sentinel for missingness
            if ~isempty(svals[r][c]) && ~isequal(svals[r][c],".\r")
            fvals[r,c] = parse(Float64, svals[r][c])
            else
            # captures FRED's "." sentinel
            fvals[r,c] = NaN
            end
        end
    end
    ######### Column names
    names = split(head, ",")[2:end] # Won't need the Date name (fist column) for TimeArray
    names = String[name for name in names]
    return TimeArray(timestamp, fvals, names, resp.data)
end


text:
#######################################
#
# Filtering
#
#######################################

"""
    filter_highpass(a::SSR; cutOff::Real=2, fs::Real=samplingrate(a), order::Int=3, tolerance::Real=0.01, kwargs...)

Applly a high pass filter.

A zero phase high pass filter is applied to the data using `filtfilt`.
A check is performed to ensure the filter does not affect the modulation rate.
The filter coefficents are stored in the processing field.

# Examples
```julia
a = read_SSR(fname)
b = filter_highpass(a)
c = filter_highpass(a, cutOff = 1)
```

"""
function filter_highpass(
    a::SSR;
    cutOff::Union{typeof(1.0u"Hz"),typeof(1u"Hz")} = 2u"Hz",
    order::Int = 3,
    tolerance::Real = 0.01,
    phase::AbstractString = "zero-double",
    kwargs...,
)
    Wn = (cutOff |> u"Hz" |> ustrip) / (samplingrate(Float64, a) / 2)
    f = digitalfilter(Highpass(Wn), Butterworth(order))
    _filter_check(f, modulationrate(a), samplingrate(Float64, a), tolerance)
    a = filter(a, f, phase)
    _append_filter(a, f)
end

"""
    filter_lowpass(a::SSR; cutOff::Real=150, fs::Real=samplingrate(a), order::Int=3, tolerance::Real=0.01, kwargs...)

Applly a low pass filter.

A zero phase high pass filter is applied to the data using `filtfilt`.
A check is performed to ensure the filter does not affect the modulation rate.
The filter coefficents are stored in the processing field.

# Examples
```julia
a = read_SSR(fname)
b = filter_lowpass(a)
c = filter_lowpass(a, cutOff = 1)
```

"""
function filter_lowpass(
    a::SSR;
    cutOff::Union{typeof(1.0u"Hz"),typeof(1u"Hz")} = 150u"Hz",
    order::Int = 3,
    tolerance::Real = 0.01,
    phase::AbstractString = "zero-double",
    kwargs...,
)
    Wn = (cutOff |> u"Hz" |> ustrip) / (samplingrate(Float64, a) / 2)
    f = digitalfilter(Lowpass(Wn), Butterworth(order))
    _filter_check(f, modulationrate(a), samplingrate(Float64, a), tolerance)
    a = filter(a, f, phase)
    _append_filter(a, f)
end

"""
    filter_bandpass(a::SSR; lower::Number=modulationrate(a) - 1, upper::Number=modulationrate(a) + 1, fs::Real=samplingrate(a), n::Int=24, rp::Number = 0.0001, tolerance::Real=0.01, kwargs...)

Applly a band pass filter.
A check is performed to ensure the filter does not affect the modulation rate.
The filter coefficents are stored in the processing field.

# Examples
```julia
a = read_SSR(fname)
a = filter_bandpass(a)
```
"""
function filter_bandpass(
    a::SSR;
    lower::Union{typeof(1.0u"Hz"),typeof(1u"Hz")} = (modulationrate(a) - 1) * 1.0u"Hz",
    upper::Union{typeof(1.0u"Hz"),typeof(1u"Hz")} = (modulationrate(a) + 1) * 1.0u"Hz",
    n::Int = 24,
    rp::Number = 0.0001,
    tolerance::Number = 0.01,
    phase::AbstractString = "zero-double",
    kwargs...,
)

    Wn_lower = (lower |> u"Hz" |> ustrip) / (samplingrate(Float64, a) / 2)
    Wn_upper = (upper |> u"Hz" |> ustrip) / (samplingrate(Float64, a) / 2)

    f = digitalfilter(Bandpass(Wn_lower, Wn_upper), Chebyshev1(n, rp))
    _filter_check(f, modulationrate(a), samplingrate(Float64, a), tolerance)
    a = filter(a, f, phase)
    _append_filter(a, f)
end

function _filter_check(
    f::FilterCoefficients,
    mod_freq::Number,
    fs::Number,
    tolerance::Number,
)
    #
    # Ensure that the filter does not alter the modulation frequency greater than a set tolerance
    #

    mod_change = abs.(freqresp(f, mod_freq * ((2pi) / fs)))

    if mod_change > 1 + tolerance || mod_change < 1 - tolerance
        @warn(
            "Filtering has modified modulation frequency greater than set tolerance: $mod_change"
        )
    end
    @debug("Filter magnitude at modulation frequency: $(mod_change)")
end

function _append_filter(a::SSR, f::FilterCoefficients; name::AbstractString = "filter")
    #
    # Put the filter information in the SSR processing structure
    #

    key_name = new_processing_key(a.processing, name)
    merge!(a.processing, Dict(key_name => f))

    return a
end

#######################################
#
# Downsample
#
#######################################

"""
    downsample(s::SSR, ratio::Rational)

Downsample signal by specified ratio.

"""
function downsample(s::SSR, ratio::Rational)

    @info("Downsampling SSR by ratio $ratio")

    dec_filter = DSP.FIRFilter([1], ratio)

    new_data =
        zeros(typeof(s.data[1, 1]), round.(Int, size(s.data, 1) * ratio), size(s.data, 2))

    for c = 1:size(s.data, 2)

        new_data[:, c] = DSP.filt(dec_filter, vec(s.data[:, c]))
    end

    s.data = new_data

    s.triggers["Index"] = round.(Int, s.triggers["Index"] .* ratio)
    if s.triggers["Index"][1] == 0
        s.triggers["Index"][1] = 1
    end

    s.samplingrate = float(samplingrate(Float64, s) * ratio) * 1.0u"Hz"

    return s
end


text:
module AdaptiveToleranceABC_MCMC
# Adaptive tolerance ABC-MCMC algorithm with post-correction

export ABCCutoff, ABCSimpleCutoff, ABCGenCutoff, ABCGaussCutoff, ABCEpaCutoff,
       abc_mcmc, abc_postprocess, ABCdefaultDist, plot_abc

using Statistics, LinearAlgebra, Random, AdaptiveMCMC

import MonteCarloMarkovKernels, Plots

# Abstract 'cutoff function' type
abstract type ABCCutoff end

# Simple cutoff at 1.0
mutable struct ABCSimpleCutoff{F,T} <: ABCCutoff where {F<:Function,T}
    tol::T
    fun::F
end

# General cutoff function
mutable struct ABCGenCutoff{F,T} <: ABCCutoff where {F<:Function,T}
    tol::T
    fun::F
end

# Outer constructors for the abstract type:
# if no function supplied, it's simple...
function ABCCutoff(tol::FT, fun=nothing) where {FT <: AbstractFloat}
    if fun isa Nothing
        return ABCSimpleCutoff(tol, x -> (x <= one(FT)))
    else
        return  ABCGenCutoff(tol, fun)
    end
end
# Gaussian
function ABCGaussCutoff(tol::FT) where {FT <: AbstractFloat}
    ABCGenCutoff(tol, x -> exp(-FT(0.5) * x^2))
end
# Epanechnikov
function ABCEpaCutoff(tol::FT) where {FT <: AbstractFloat}
    ABCGenCutoff(tol, x -> x >= one(FT) ? zero(FT) : one(FT) - x^2)
end

# Convenience functions
@inline function (c::ABCGenCutoff)(x)
    c.fun(x/c.tol)
end
@inline function (c::ABCSimpleCutoff)(x)
    c.fun(x/c.tol)
end

@inline function ABCdefaultDist(x::T, y::T) where {
    FT <: AbstractFloat, T <: AbstractVector{FT}}

    s = zero(FT)
    for i in 1:length(x)
        s += (x[i] - y[i])^2
    end
    isnan(s) ? FT(Inf) : sqrt(s)
end

include("abc_mcmc.jl")
include("post_correct.jl")
include("regress_correct.jl")
include("plot_abc.jl")

end # module


text:
module SetProg

using LinearAlgebra

include("Sets/Sets.jl")

import MutableArithmetics
const MA = MutableArithmetics

import Reexport
Reexport.@reexport using MultivariateBases
Reexport.@reexport using JuMP
const MOI = JuMP.MOI
Reexport.@reexport using Polyhedra

using SumOfSquares
using MultivariateMoments
import DynamicPolynomials
import MultivariatePolynomials
const MP = MultivariatePolynomials
const SpaceVariable = DynamicPolynomials.PolyVar{true}

export Polytope, Ellipsoid, PolySet
export nth_root, L1_heuristic

@enum Space Undecided PrimalSpace DualSpace

include("utilities.jl")
include("spaces.jl")
include("macros.jl") # need to be before `variables.jl` and `constraints.jl` as they use `⊆` in `@constraint`
include("variables.jl")
include("map.jl")
include("constraints.jl")
include("objective.jl")
include("data.jl")
include("optimize.jl")

end # module


text:
import Base.convert

convert(::Type{Int}, v::SubString{String}) = parse(Int, v)
convert(::Type{Float64}, v::SubString{String}) = parse(Float64, v)