text:
# Types
struct ALPHA end
struct BETA end
struct MULTIALPHA end
struct MULTIBETA end
struct NTDALPHA end
struct NTDBETA end

# # Check NaN value (only GD)
# function checkNaN(W::AbstractArray, cca::GD)
#     l = length(W)
#     for i=1:l
#         if any(isnan, W[i])
#             error("NaN values are generated. Select other stepsize")
#         end
#     end
# end

# # Check NaN value (other CCA)
# function checkNaN(N::AbstractArray, s::Number, n::Number, W::AbstractArray, evalfreq::Number, cca::Union{SGD,RSGD,HORST,ORTHITER})
#     l = length(W)
#     for i=1:l
#         if mod((N[i]*(s-1)+n), evalfreq) == 0
#             if any(isnan, W[i])
#                 error("NaN values are generated. Select other stepsize")
#             end
#         end
#     end
# end

# # Output the result of CCA
# function output(outdir::AbstractString, out::Tuple, cca::NIPALS)
#     writecsv(joinpath(outdir, "T.csv"), out[1])
#     writecsv(joinpath(outdir, "P.csv"), out[2])
#     writecsv(joinpath(outdir, "U.csv"), out[3])
#     writecsv(joinpath(outdir, "Q.csv"), out[4])
#     writecsv(joinpath(outdir, "B.csv"), out[5])
# end

# # Output the result of CCA
# function output(outdir::AbstractString, out::Tuple, cca::Union{GD,SGD,RSGD,HORST,ORTHITER})
#     l = length(out[1])
#     for i = 1:l
#         writecsv(joinpath(outdir, "Eigen_vectors"*string(i)*".csv"), out[1][i])
#         writecsv(joinpath(outdir, "Eigen_values"*string(i)*".csv"), out[2][i])
#         writecsv(joinpath(outdir, "LatentVariables"*string(i)*".csv"), out[3][i])
#         writecsv(joinpath(outdir, "ExpCorVar"*string(i)*".csv"), out[4][i])
#         writecsv(joinpath(outdir, "TotalCorVar"*string(i)*".csv"), out[5][i])
#         writecsv(joinpath(outdir, "Objective"*string(i)*".csv"), out[6][i])
#         writecsv(joinpath(outdir, "AvgLatentCor"*string(i)*".csv"), out[7][i])
#     end
#     if out[8] == 1
#         touch(joinpath(outdir, "Converged"))
#     end
# end

# # Output the result of CCA
# function output(outdir::AbstractString, out::Tuple, cca::OOCMCCA)
#     l = length(out[1])
#     for i = 1:l
#         writecsv(joinpath(outdir, "Eigen_vectors"*string(i)*".csv"), out[1][i])
#         writecsv(joinpath(outdir, "Eigen_values"*string(i)*".csv"), out[2][i])
#         writecsv(joinpath(outdir, "LatentVariables"*string(i)*".csv"), out[3][i])
#         writecsv(joinpath(outdir, "ExpCorVar"*string(i)*".csv"), out[4][i])
#         writecsv(joinpath(outdir, "TotalCorVar"*string(i)*".csv"), out[5][i])
#         writecsv(joinpath(outdir, "Objective"*string(i)*".csv"), out[6][i])
#         writecsv(joinpath(outdir, "AvgLatentCor"*string(i)*".csv"), out[7][i])
#     end
# end

writecsv(filename::AbstractString, data) = writedlm(filename, data, ',')
readcsv(filename::AbstractString) = readdlm(filename, ',')
readcsv(filename::AbstractString, ::Type{T}) where {T} = readdlm(filename, ',', T)

# # Parse command line options
# function parse_commandline(cca::OOCMCCA)
#     s = ArgParseSettings()

#     @add_arg_table s begin
#         "--input", "-i"
#             help = "Julia Binary file generated by `OnlineCCA.csv2bin` function."
#             arg_type = AbstractString
#             required = true
#         "--outdir", "-o"
#             help = "The directory specified the directory you want to save the result."
#             arg_type = Union{Nothing,AbstractString}
#             default = "."
#             required = false
#         "--scale"
#             help = "{log,ftt,raw}-scaling of the value."
#             arg_type = AbstractString
#         "--pseudocount", "-p"
#             help = "The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated."
#             arg_type = Union{Number,AbstractString}
#         "--colmeanlist", "-m"
#             help = "The mean of each row of matrix. The CSV file is generated by `OnlineCCA.sumr` functions."
#             arg_type = AbstractString
#             default = ""
#             required = false
#         "--colvarlist", "-v"
#             help = "The variance of each row of matrix. The CSV file is generated by `OnlineCCA.sumr` functions."
#             arg_type = AbstractString
#             default = ""
#             required = false
#         "--dim", "-d"
#             help = "The number of dimension of CCA."
#             arg_type = Union{Number,AbstractString}
#             default = 3
#         "--numepoch", "-e"
#             help = "The number of epoch."
#             arg_type = Union{Number,AbstractString}
#             default = 5
#         "--lower"
#             help = "Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)"
#             arg_type = Union{Number,AbstractString}
#             default = 0
#         "--upper"
#             help = "Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)"
#             arg_type = Union{Number,AbstractString}
#             default = 1.0f+38
#         "--offsetStoch"
#             help = "Off set value for avoding overflow when calculating stochastic gradient"
#             arg_type = Union{Number,AbstractString}
#             default = 1.0f-20
#         "--initW"
#             help = "The CSV file saving the initial values of eigenvectors."
#             arg_type = Union{Nothing,AbstractString,}
#             default = nothing
#         "--perm"
#             help = "Whether the data matrix is shuffled at random"
#             arg_type = Union{Bool,AbstractString}
#             default = false
#     end

#     return parse_args(s)
# end

# # Parse command line options
# function parse_commandline(cca::Union{HORST,ORTHITER})
#     s = ArgParseSettings()

#     @add_arg_table s begin
#         "--input", "-i"
#             help = "Julia Binary file generated by `OnlineCCA.csv2bin` function."
#             arg_type = AbstractString
#             required = true
#         "--outdir", "-o"
#             help = "The directory specified the directory you want to save the result."
#             arg_type = Union{Nothing,AbstractString}
#             default = "."
#             required = false
#         "--scale"
#             help = "{log,ftt,raw}-scaling of the value."
#             arg_type = AbstractString
#         "--pseudocount", "-p"
#             help = "The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated."
#             arg_type = Union{Number,AbstractString}
#         "--colmeanlist", "-m"
#             help = "The mean of each row of matrix. The CSV file is generated by `OnlineCCA.sumr` functions."
#             arg_type = AbstractString
#             default = ""
#             required = false
#         "--colvarlist", "-v"
#             help = "The variance of each row of matrix. The CSV file is generated by `OnlineCCA.sumr` functions."
#             arg_type = AbstractString
#             default = ""
#             required = false
#         "--dim", "-d"
#             help = "The number of dimension of CCA."
#             arg_type = Union{Number,AbstractString}
#             default = 3
#         "--numepoch", "-e"
#             help = "The number of epoch."
#             arg_type = Union{Number,AbstractString}
#             default = 5
#         "--lower"
#             help = "Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)"
#             arg_type = Union{Number,AbstractString}
#             default = 0
#         "--upper"
#             help = "Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)"
#             arg_type = Union{Number,AbstractString}
#             default = 1.0f+38
#         "--evalfreq"
#             help = "Evaluation Frequency of Reconstruction Error"
#             arg_type = Union{Number,AbstractString}
#             default = 5000
#         "--offsetStoch"
#             help = "Off set value for avoding overflow when calculating stochastic gradient"
#             arg_type = Union{Number,AbstractString}
#             default = 1.0f-20
#         "--initW"
#             help = "The CSV file saving the initial values of eigenvectors."
#             arg_type = Union{Nothing,AbstractString,}
#             default = nothing
#         "--logdir", "-l"
#             help = "The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration."
#             arg_type = Union{Nothing,AbstractString}
#             default = nothing
#         "--perm"
#             help = "Whether the data matrix is shuffled at random"
#             arg_type = Union{Bool,AbstractString}
#             default = false
#     end

#     return parse_args(s)
# end

# # Parse command line options
# function parse_commandline(cca::NIPALS)
#     s = ArgParseSettings()

#     @add_arg_table s begin
#         "--input", "-i"
#             help = "Julia Binary file generated by `OnlineCCA.csv2bin` function."
#             arg_type = AbstractString
#             required = true
#         "--outdir", "-o"
#             help = "The directory specified the directory you want to save the result."
#             arg_type = Union{Nothing,AbstractString}
#             default = "."
#             required = false
#         "--scale"
#             help = "{log,ftt,raw}-scaling of the value."
#             arg_type = AbstractString
#         "--pseudocount", "-p"
#             help = "The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated."
#             arg_type = Union{Number,AbstractString}
#         "--colmeanlist", "-m"
#             help = "The mean of each row of matrix. The CSV file is generated by `OnlineCCA.sumr` functions."
#             arg_type = AbstractString
#             default = ""
#             required = false
#         "--colvarlist", "-v"
#             help = "The variance of each row of matrix. The CSV file is generated by `OnlineCCA.sumr` functions."
#             arg_type = AbstractString
#             default = ""
#             required = false
#         "--dim", "-d"
#             help = "The number of dimension of CCA."
#             arg_type = Union{Number,AbstractString}
#             default = 3
#         "--numepoch", "-e"
#             help = "The number of epoch."
#             arg_type = Union{Number,AbstractString}
#             default = 5
#         "--lower"
#             help = "Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)"
#             arg_type = Union{Number,AbstractString}
#             default = 0
#         "--upper"
#             help = "Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)"
#             arg_type = Union{Number,AbstractString}
#             default = 1.0f+38
#         "--offsetStoch"
#             help = "Off set value for avoding overflow when calculating stochastic gradient"
#             arg_type = Union{Number,AbstractString}
#             default = 1.0f-20
#     end

#     return parse_args(s)
# end

# # Parse command line options
# function parse_commandline(cca::Union{GD,SGD,RSGD})
#     s = ArgParseSettings()

#     @add_arg_table s begin
#         "--input", "-i"
#             help = "Julia Binary file generated by `OnlineCCA.csv2bin` function."
#             arg_type = AbstractString
#             required = true
#         "--outdir", "-o"
#             help = "The directory specified the directory you want to save the result."
#             arg_type = Union{Nothing,AbstractString}
#             default = "."
#             required = false
#         "--scale"
#             help = "{log,ftt,raw}-scaling of the value."
#             arg_type = AbstractString
#         "--pseudocount", "-p"
#             help = "The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated."
#             arg_type = Union{Number,AbstractString}
#         "--colmeanlist", "-m"
#             help = "The mean of each row of matrix. The CSV file is generated by `OnlineCCA.sumr` functions."
#             arg_type = AbstractString
#             default = ""
#             required = false
#         "--colvarlist", "-v"
#             help = "The variance of each row of matrix. The CSV file is generated by `OnlineCCA.sumr` functions."
#             arg_type = AbstractString
#             default = ""
#             required = false
#         "--dim", "-d"
#             help = "The number of dimension of CCA."
#             arg_type = Union{Number,AbstractString}
#             default = 3
#         "--stepsize", "-s"
#             help = "The parameter used in every iteration."
#             arg_type = Union{Number,AbstractString}
#             default = 1.0f3
#         "--numepoch", "-e"
#             help = "The number of epoch."
#             arg_type = Union{Number,AbstractString}
#             default = 5
#         "--scheduling"
#             help = "Learning parameter scheduling. `robbins-monro`, `momentum`, `nag`, and `adagrad` are available."
#             arg_type = AbstractString
#             default = "robbins-monro"
#         "-g"
#             help = "The parameter that is used when scheduling is specified as nag."
#             arg_type = Union{Number,AbstractString}
#             default = 0.9f0
#         "--epsilon"
#             help = "The parameter that is used when scheduling is specified as adagrad."
#             arg_type = Union{Number,AbstractString}
#             default = 1.0f-8
#         "--lower"
#             help = "Stopping Criteria (When the relative change of error is below this value, the calculation is terminated)"
#             arg_type = Union{Number,AbstractString}
#             default = 0
#         "--upper"
#             help = "Stopping Criteria (When the relative change of error is above this value, the calculation is terminated)"
#             arg_type = Union{Number,AbstractString}
#             default = 1.0f+38
#         "--evalfreq"
#             help = "Evaluation Frequency of Reconstruction Error"
#             arg_type = Union{Number,AbstractString}
#             default = 5000
#         "--offsetFull"
#             help = "Off set value for avoding overflow when calculating full gradient"
#             arg_type = Union{Number,AbstractString}
#             default = 1.0f-20
#         "--offsetStoch"
#             help = "Off set value for avoding overflow when calculating stochastic gradient"
#             arg_type = Union{Number,AbstractString}
#             default = 1.0f-20
#         "--initW"
#             help = "The CSV file saving the initial values of eigenvectors."
#             arg_type = Union{Nothing,AbstractString,}
#             default = nothing
#         "--logdir", "-l"
#             help = "The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration."
#             arg_type = Union{Nothing,AbstractString}
#             default = nothing
#         "--perm"
#             help = "Whether the data matrix is shuffled at random"
#             arg_type = Union{Bool,AbstractString}
#             default = false
#     end

#     return parse_args(s)
# end

# # Return N, M
# function nm(input::AbstractArray)
#     l = length(input)
#     N = zeros(UInt32, l)
#     M = zeros(UInt32, l)
#     tmpN = zeros(UInt32, 1)
#     tmpM = zeros(UInt32, 1)
#     for i=1:l
#         open(input[i]) do file
#             stream = ZstdDecompressorStream(file)
#             read!(stream, tmpN)
#             N[i] = tmpN[1]
#             read!(stream, tmpM)
#             M[i] = tmpM[1]
#             close(stream)
#         end
#     end
#     return N, M
# end

# # Initialization
# function init(input::AbstractArray, pseudocount::AbstractArray, dim::Number, colmeanlist::AbstractArray, colvarlist::Union{Nothing,AbstractString,AbstractArray}, initW::Union{Nothing,AbstractArray}, cca::Union{OOCMCCA}, offsetStoch::Number, scale::AbstractArray)
#     if !(length(input) >= 2)
#         error("Please specify the input as length-2 or longer array (e.g. [\"X.zst\", \"Y.zst\", \"Z.zst\"]).")
#     end
#     check1 = length(input) == length(pseudocount)
#     check2 = length(input) == length(colmeanlist)
#     check3 = true
#     if !(colvarlist == "")
#         check3 = length(input) == length(colvarlist)
#     end
#     check4 = length(input) == length(scale)
#     if !(check1 && check2 && check3 && check4)
#         error("The length of vectors input, pseudocount, colmeanlist, colvarlist, and scale have to be same.")
#     end
#     N, M = nm(input)
#     tmpN = zeros(UInt32, 1)
#     tmpM = zeros(UInt32, 1)
#     l = length(input)
#     tmppseudocount = zeros(Float32, l)
#     for i=1:l
#         tmppseudocount[i] = Float32(pseudocount[i])
#     end
#     pseudocount = tmppseudocount
#     offsetStoch = Float32(offsetStoch)
#     W = Any[]
#     if initW != nothing
#         for i=1:l
#             tmpW = readcsv(initW[i], Float32)
#             push!(W, tmpW) # Eigen vectors
#         end
#     end
#     v = Any[]
#     D = Any[]
#     colmeanvec = Any[]
#     colvarvec = Any[]
#     for i=1:l
#         push!(v, zeros(Float32, M[i], dim)) # Temporal Vector (Same length
#         if initW == nothing
#             tmpW = zeros(Float32, M[i], dim)
#             for j=1:dim
#                 tmpW[(i-1)*dim+j, j] = 1
#             end
#             push!(W, tmpW) # Eigen vectors
#         end
#         push!(D, Diagonal(reverse(1:dim))) # Diagonaml Matrix
#         push!(colmeanvec, readcsv(colmeanlist[i], Float32))
#         if length(colvarlist) != 0
#             push!(colvarvec, readcsv(colvarlist[i], Float32))
#         else
#             push!(colvarvec, zeros(Float32, N[i], 1))
#         end
#     end
#     if colvarlist == ""
#         colvarlist=Any[]
#         for i=1:l
#             push!(colvarlist, "")
#         end
#     end
#     TotalCorVar = zeros(Float32, l)
#     for i=1:l
#         x = zeros(UInt32, M[i])
#         normx = zeros(Float32, M[i])
#         open(input[i]) do file
#             stream = ZstdDecompressorStream(file)
#             read!(stream, tmpN)
#             read!(stream, tmpM)
#             for n = 1:N[i]
#                 # Data Import
#                 read!(stream, x)
#                 normx = normalizex(x, n, stream, scale[i], pseudocount[i], colmeanlist[i], colmeanvec[i], colvarlist[i], colvarvec[i])
#                 TotalCorVar[i] = TotalCorVar[i] .+ (normx'normx)[1]
#             end
#             close(stream)
#         end
#         TotalCorVar[i] = TotalCorVar[i] / M[i]
#     end
#     return pseudocount, W, v, D, colmeanlist, colmeanvec, colvarlist, colvarvec, N, M, TotalCorVar, offsetStoch
# end

# # Initialization
# function init(input::AbstractArray, pseudocount::AbstractArray, dim::Number, colmeanlist::AbstractArray, colvarlist::Union{Nothing,AbstractString,AbstractArray}, initW::Union{Nothing,AbstractArray}, logdir::Union{Nothing,AbstractString}, cca::Union{HORST,ORTHITER}, lower::Number, upper::Number, evalfreq::Number, offsetStoch::Number, scale::AbstractArray)
#     if !(length(input) >= 2)
#         error("Please specify the input as length-2 or longer array (e.g. [\"X.zst\", \"Y.zst\", \"Z.zst\"]).")
#     end
#     check1 = length(input) == length(pseudocount)
#     check2 = length(input) == length(colmeanlist)
#     check3 = true
#     if !(colvarlist == "")
#         check3 = length(input) == length(colvarlist)
#     end
#     check4 = length(input) == length(scale)
#     if !(check1 && check2 && check3 && check4)
#         error("The length of vectors input, pseudocount, colmeanlist, colvarlist, and scale have to be same.")
#     end
#     N, M = nm(input)
#     tmpN = zeros(UInt32, 1)
#     tmpM = zeros(UInt32, 1)
#     l = length(input)
#     tmppseudocount = zeros(Float32, l)
#     for i=1:l
#         tmppseudocount[i] = Float32(pseudocount[i])
#     end
#     pseudocount = tmppseudocount
#     lower = Float32(lower)
#     upper = Float32(upper)
#     evalfreq = Int64(evalfreq)
#     offsetStoch = Float32(offsetStoch)
#     W = Any[]
#     if initW != nothing
#         for i=1:l
#             tmpW = readcsv(initW[i], Float32)
#             push!(W, tmpW) # Eigen vectors
#         end
#     end
#     v = Any[]
#     D = Any[]
#     colmeanvec = Any[]
#     colvarvec = Any[]
#     for i=1:l
#         push!(v, zeros(Float32, M[i], dim)) # Temporal Vector (Same length
#         if initW == nothing
#             tmpW = zeros(Float32, M[i], dim)
#             for j=1:dim
#                 tmpW[(i-1)*dim+j, j] = 1
#             end
#             push!(W, tmpW) # Eigen vectors
#         end
#         push!(D, Diagonal(reverse(1:dim))) # Diagonaml Matrix
#         push!(colmeanvec, readcsv(colmeanlist[i], Float32))
#         if length(colvarlist) != 0
#             push!(colvarvec, readcsv(colvarlist[i], Float32))
#         else
#             push!(colvarvec, zeros(Float32, N[i], 1))
#         end
#     end
#     if colvarlist == ""
#         colvarlist=Any[]
#         for i=1:l
#             push!(colvarlist, "")
#         end
#     end
#     TotalCorVar = zeros(Float32, l)
#     for i=1:l
#         x = zeros(UInt32, M[i])
#         normx = zeros(Float32, M[i])
#         open(input[i]) do file
#             stream = ZstdDecompressorStream(file)
#             read!(stream, tmpN)
#             read!(stream, tmpM)
#             for n = 1:N[i]
#                 # Data Import
#                 read!(stream, x)
#                 normx = normalizex(x, n, stream, scale[i], pseudocount[i], colmeanlist[i], colmeanvec[i], colvarlist[i], colvarvec[i])
#                 TotalCorVar[i] = TotalCorVar[i] .+ (normx'normx)[1]
#             end
#             close(stream)
#         end
#         TotalCorVar[i] = TotalCorVar[i] / M[i]
#     end
#     # directory for log file
#     if logdir isa String
#         if !isdir(logdir)
#             mkdir(logdir)
#         end
#     end
#     return pseudocount, W, v, D, colmeanlist, colmeanvec, colvarlist, colvarvec, N, M, TotalCorVar, lower, upper, evalfreq, offsetStoch
# end

# # Initialization
# function init(input::AbstractArray, pseudocount::AbstractArray, dim::Number, colmeanlist::AbstractArray, colvarlist::Union{Nothing,AbstractString,AbstractArray}, cca::NIPALS, lower::Number, upper::Number, offsetStoch::Number, scale::AbstractArray)
#     if !(length(input) == 2)
#         error("Please specify the input as length-2 array (e.g. [\"X.zst\", \"Y.zst\"]).")
#     end
#     check1 = length(input) == length(pseudocount)
#     check2 = length(input) == length(colmeanlist)
#     check3 = true
#     if !(colvarlist == "")
#         check3 = length(input) == length(colvarlist)
#     end
#     check4 = length(input) == length(scale)
#     if !(check1 && check2 && check3 && check4)
#         error("The length of vectors input, pseudocount, colmeanlist, colvarlist, and scale have to be same.")
#     end
#     N, M = nm(input)
#     tmpN = zeros(UInt32, 1)
#     tmpM = zeros(UInt32, 1)
#     l = length(input)
#     tmppseudocount = zeros(Float32, l)
#     for i=1:l
#         tmppseudocount[i] = Float32(pseudocount[i])
#     end
#     pseudocount = tmppseudocount
#     lower = Float32(lower)
#     upper = Float32(upper)
#     offsetStoch = Float32(offsetStoch)
#     W = Any[]
#     v = Any[]
#     D = Any[]
#     colmeanvec = Any[]
#     colvarvec = Any[]
#     for i=1:l
#         push!(v, zeros(Float32, M[i], dim)) # Temporal Vector (Same length
#         tmpW = zeros(Float32, M[i], dim)
#         for j=1:dim
#             tmpW[(i-1)*dim+j, j] = 1
#         end
#         push!(W, tmpW) # Eigen vectors
#         push!(D, Diagonal(reverse(1:dim))) # Diagonaml Matrix
#         push!(colmeanvec, readcsv(colmeanlist[i], Float32))
#         if length(colvarlist) != 0
#             push!(colvarvec, readcsv(colvarlist[i], Float32))
#         else
#             push!(colvarvec, zeros(Float32, N[i], 1))
#         end
#     end
#     if colvarlist == ""
#         colvarlist=Any[]
#         for i=1:l
#             push!(colvarlist, "")
#         end
#     end
#     TotalCorVar = zeros(Float32, l)
#     for i=1:l
#         x = zeros(UInt32, M[i])
#         normx = zeros(Float32, M[i])
#         open(input[i]) do file
#             stream = ZstdDecompressorStream(file)
#             read!(stream, tmpN)
#             read!(stream, tmpM)
#             for n = 1:N[i]
#                 # Data Import
#                 read!(stream, x)
#                 normx = normalizex(x, n, stream, scale[i], pseudocount[i], colmeanlist[i], colmeanvec[i], colvarlist[i], colvarvec[i])
#                 TotalCorVar[i] = TotalCorVar[i] .+ (normx'normx)[1]
#             end
#             close(stream)
#         end
#         TotalCorVar[i] = TotalCorVar[i] / M[i]
#     end
#     return pseudocount, W, v, D, colmeanlist, colmeanvec, colvarlist, colvarvec, N, M, TotalCorVar, lower, upper, offsetStoch
# end

# # Initialization
# function init(input::AbstractArray, pseudocount::AbstractArray, stepsize::Number, g::Number, epsilon::Number, dim::Number, colmeanlist::AbstractArray, colvarlist::Union{AbstractString,AbstractArray}, initW::Union{Nothing,AbstractArray}, logdir::Union{Nothing,AbstractString}, cca::Union{GD,SGD,RSGD}, lower::Number, upper::Number, evalfreq::Number, offsetFull::Number, offsetStoch::Number, scale::AbstractArray)
# 	if !(length(input) >= 2)
#         error("Please specify the input as length-2 or longer array (e.g. [\"X.zst\",\"Y.zst\",\"Z.zst\"]).")
# 	end
#     check1 = length(input) == length(pseudocount)
#     check2 = length(input) == length(colmeanlist)
#     check3 = true
#     if !(colvarlist == "")
#         check3 = length(input) == length(colvarlist)
#     end
#     check4 = length(input) == length(scale)
#     if !(check1 && check2 && check3 && check4)
#         error("The length of vectors input, pseudocount, colmeanlist, colvarlist, and scale have to be same.")
#     end
#     N, M = nm(input)
#     tmpN = zeros(UInt32, 1)
#     tmpM = zeros(UInt32, 1)
#     l = length(input)
#     tmppseudocount = zeros(Float32, l)
#     for i=1:l
#         tmppseudocount[i] = Float32(pseudocount[i])
#     end
#     pseudocount = tmppseudocount
#     stepsize = Float32(stepsize)
#     g = Float32(g)
#     epsilon = Float32(epsilon)
#     lower = Float32(lower)
#     upper = Float32(upper)
#     evalfreq = Int64(evalfreq)
#     offsetFull = Float32(offsetFull)
#     offsetStoch = Float32(offsetStoch)
#     W = Any[]
#     if initW != nothing
#         for i=1:l
#             tmpW = readcsv(initW[i], Float32)
#             push!(W, tmpW) # Eigen vectors
#         end
#     end
#     v = Any[]
#     D = Any[]
#     colmeanvec = Any[]
#     colvarvec = Any[]
#     for i=1:l
#         push!(v, zeros(Float32, M[i], dim)) # Temporal Vector (Same length
#         tmp = zeros(Float32, M[i], dim)
#         for j=1:dim
#             tmp[j, j] = 1
#         end
#         if initW == nothing
#             push!(W, tmp) # Eigen vectors
#         end
#         push!(D, Diagonal(reverse(1:dim))) # Diagonaml Matrix
#         push!(colmeanvec, readcsv(colmeanlist[i], Float32))
#         if length(colvarlist) != 0
#             push!(colvarvec, readcsv(colvarlist[i], Float32))
#         else
#             push!(colvarvec, zeros(Float32, N[i], 1))
#         end
#     end
#     if colvarlist == ""
#         colvarlist=Any[]
#         for i=1:l
#             push!(colvarlist, "")
#         end
#     end
#     TotalCorVar = zeros(Float32, l)
#     for i=1:l
#         x = zeros(UInt32, M[i])
#         normx = zeros(Float32, M[i])
#         open(input[i]) do file
#             stream = ZstdDecompressorStream(file)
#             read!(stream, tmpN)
#             read!(stream, tmpM)
#             for n = 1:N[i]
#                 # Data Import
#                 read!(stream, x)
#                 normx = normalizex(x, n, stream, scale[i], pseudocount[i], colmeanlist[i], colmeanvec[i], colvarlist[i], colvarvec[i])
#                 TotalCorVar[i] = TotalCorVar[i] .+ (normx'normx)[1]
#             end
#             close(stream)
#         end
#         TotalCorVar[i] = TotalCorVar[i] / M[i]
#     end
#     # directory for log file
#     if logdir isa String
#         if !isdir(logdir)
#             mkdir(logdir)
#         end
#     end
#     return pseudocount, stepsize, g, epsilon, W, v, D, colmeanlist, colmeanvec, colvarlist, colvarvec, N, M, TotalCorVar, lower, upper, evalfreq, offsetFull, offsetStoch
# end

# # Eigen value, Loading
# function WλV(W::AbstractArray, input::AbstractArray, dim::Number, scale::AbstractArray, pseudocount::AbstractArray, colmeanlist::AbstractArray, colmeanvec::AbstractArray, colvarlist::AbstractArray, colvarvec::AbstractArray, TotalCorVar::AbstractArray)
#     N, M = nm(input)
#     l = length(W)
#     λ = Any[]
#     V = Any[]
#     x = Any[]
#     normx = Any[]
#     for i=1:l
#         push!(V, zeros(N[i], dim))
#         push!(x, zeros(UInt32, M[i]))
#         push!(normx, zeros(UInt32, M[i]))
#     end
#     ExpCorCov = zeros(Float32, l)
#     tmpN = zeros(UInt32, 1)
#     tmpM = zeros(UInt32, 1)
#     Objective = Any[]
#     AvgLatentCor = Any[]
#     for i=1:l
#         open(input[i]) do file
#             stream = ZstdDecompressorStream(file)
#             read!(stream, tmpN)
#             read!(stream, tmpM)
#             for n=1:N[i]
#                 # Data Import
#                 read!(stream, x[i])
#                 normx[i] = normalizex(x[i], n, stream, scale[i], pseudocount[i], colmeanlist[i], colmeanvec[i], colvarlist[i], colvarvec[i])
#                 V[i][n, :] = normx[i]'W[i]
#             end
#             close(stream)
#         end
#         # Eigen value
#         σ = Float32[norm(V[i][:, x]) for x=1:dim]
#         for n = 1:dim
#             V[i][:, n] ./= σ[n]
#         end
#         push!(λ, σ .* σ ./ M[i])
#         # Sort by Eigen value
#         idx = sortperm(λ[i], rev=true)
#         W[i] .= W[i][:, idx]
#         λ[i] .= λ[i][idx]
#         V[i] .= V[i][:, idx]
#         ExpCorCov[i] = sum(λ[i]) / TotalCorVar[i]
#     end
#     # The value of objective function/The correletion between latent variables
#     for i=1:l
#         for j=1:l
#             if i != j
#                 push!(Objective, sum(diag(V[i]'V[j])))
#                 tmpAvgLatentCor = 0.0
#                 for k=1:dim
#                     tmpAvgLatentCor += cor(V[i][:,k].+1e-10, V[j][:,k].+1e-10)
#                 end
#                 push!(AvgLatentCor, tmpAvgLatentCor/(dim*(l-1)))
#             end
#         end
#     end
#     # Return
#     return W, λ, V, ExpCorCov, TotalCorVar, Objective, AvgLatentCor
# end

# # Output log file （only GD）
# function outputlog(s::Number, input::AbstractArray, dim::Number, logdir::AbstractString, W::AbstractArray, cca::GD, TotalCorVar::AbstractArray, scale::AbstractArray, pseudocount::AbstractArray, colmeanlist::AbstractArray, colmeanvec::AbstractArray, colvarlist::AbstractArray, colvarvec::AbstractArray, lower::Number, upper::Number, stop::Number)
#     REs = RecError(W, input, TotalCorVar, scale, pseudocount, colmeanlist, colmeanvec, colvarlist, colvarvec, dim)
#     l = length(W)
#     RelChange = Any[]
#     for i = 1:l
#         outREs = ["E"=>REs[1][2][i], "AE"=>REs[2][2][i],
#         "RMSE"=>REs[3][2][i], "ARE"=>REs[4][2][i],
#         "Explained Correlation/Covariance"=>REs[5][2][i],
#         "Total Correlation/Variance"=>REs[6][2][i],
#         "Objective"=> REs[7][2][i], "AvgLatentCor"=> REs[8][2][i]]
#         writecsv(joinpath(logdir, "RecError_Epoch"*string(s)*"_"*string(i)*".csv"), outREs)
#         writecsv(joinpath(logdir, "W_Epoch"*string(s)*"_"*string(i)*".csv"), W[i])
#     end
#     if s != 1
#         for i = 1:l
#             old_Objective = readcsv(joinpath(logdir, "RecError_Epoch"*string(s-1)*"_"*string(i)*".csv"))
#             push!(RelChange, abs(REs[7][2][i] - old_Objective[7,2]) / abs(REs[7][2][i]))
#         end
#         if sum(RelChange) < lower
#             println("Relative change of reconstruction error is below the lower value (no change)")
#             stop = 1
#         end
#         if sum(RelChange) > upper
#             println("Relative change of reconstruction error is above the upper value (unstable)")
#             stop = 2
#         end
#     end
#     return stop
# end

# # Output log file (other CCA)
# function outputlog(N::AbstractArray, s::Number, n::Number, input::AbstractArray, dim::Number, logdir::AbstractString, W::AbstractArray, cca::Union{SGD,RSGD}, TotalCorVar::AbstractArray, scale::AbstractArray, pseudocount::AbstractArray, colmeanlist::AbstractArray, colmeanvec::AbstractArray, colvarlist::AbstractArray, colvarvec::AbstractArray, lower::Number, upper::Number, stop::Number, evalfreq::Number)
#     if(mod((N[1]*(s-1)+n), evalfreq) == 0)
#         REs = RecError(W, input, TotalCorVar, scale, pseudocount, colmeanlist, colmeanvec, colvarlist, colvarvec, dim)
#         l = length(W)
#         RelChange = Any[]
#         for i = 1:l
#             outREs = ["E"=>REs[1][2][i], "AE"=>REs[2][2][i],
#             "RMSE"=>REs[3][2][i], "ARE"=>REs[4][2][i],
#             "Explained Correlation/Covariance"=>REs[5][2][i],
#             "Total Correlation/Variance"=>REs[6][2][i],
#             "Objective"=> REs[7][2][i], "AvgLatentCor"=> REs[8][2][i]]
#             writecsv(joinpath(logdir, "RecError_"*string((N[1]*(s-1)+n))*"_"*string(i)*".csv"), outREs)
#             writecsv(joinpath(logdir, "W_"*string((N[1]*(s-1)+n))*"_"*string(i)*".csv"), W[i])
#         end

#         if n != evalfreq && (N[1]*(s-1)+(n-evalfreq)) != 0
#             for i = 1:l
#                 old_Objective = readcsv(joinpath(logdir, "RecError_"*string((N[1]*(s-1)+(n-evalfreq)))*"_"*string(i)*".csv"))
#                 push!(RelChange, abs(REs[7][2][i] - old_Objective[7,2]) / abs(REs[7][2][i]))
#             end
#             if sum(RelChange) < lower
#                 println("Relative change of reconstruction error is below the lower value (no change)")
#                 stop = 1
#             end
#             if sum(RelChange) > upper
#                 println("Relative change of reconstruction error is above the upper value (unstable)")
#                 stop = 2
#             end
#         end
#     end
#     return stop
# end

# # Reconstuction Error
# function RecError(W::AbstractArray, input::AbstractArray, TotalCorVar::AbstractArray, scale::AbstractArray, pseudocount::AbstractArray, colmeanlist::AbstractArray, colmeanvec::AbstractArray, colvarlist::AbstractArray, colvarvec::AbstractArray, dim::Number)
#     N, M = nm(input)
#     l = length(W)
#     λ = Any[]
#     V = Any[]
#     x = Any[]
#     normx = Any[]
#     for i=1:l
#         push!(V, zeros(N[i], dim))
#         push!(x, zeros(UInt32, M[i]))
#         push!(normx, zeros(UInt32, M[i]))
#     end
#     ExpCorCov = zeros(Float32, l)
#     tmpN = zeros(UInt32, 1)
#     tmpM = zeros(UInt32, 1)
#     Objective = Any[]
#     AvgLatentCor = Any[]
#     E = zeros(Float32, l)
#     AE = zeros(Float32, l)
#     RMSE = zeros(Float32, l)
#     ARE = zeros(Float32, l)
#     for i=1:l
#         open(input[i]) do file
#             stream = ZstdDecompressorStream(file)
#             read!(stream, tmpN)
#             read!(stream, tmpM)
#             for n = 1:N[i]
#                 # Data Import
#                 read!(stream, x[i])
#                 normx[i] = normalizex(x[i], n, stream, scale[i], pseudocount[i], colmeanlist[i], colmeanvec[i], colvarlist[i], colvarvec[i])
#                 pc = W[i]'normx[i]
#                 E[i] = E[i] + dot(normx[i], normx[i]) - dot(pc, pc)
#                 V[i][n, :] = pc
#             end
#             close(stream)
#         end
#         # Eigen value
#         σ = Float32[norm(V[i][:, x]) for x=1:dim]
#         for n = 1:dim
#             V[i][:, n] ./= σ[n]
#         end
#         push!(λ, σ .* σ ./ M[i])
#         # Sort by Eigen value
#         idx = sortperm(λ[i], rev=true)
#         W[i] .= W[i][:, idx]
#         λ[i] .= λ[i][idx]
#         V[i] .= V[i][:, idx]
#         ExpCorCov[i] = sum(λ[i]) / TotalCorVar[i]
#         ExpCorCov[i] = sum(λ[i]) / TotalCorVar[i]
#         AE[i] = E[i] / M[i]
#         RMSE[i] = sqrt(E[i] / (N[i] * M[i]))
#         ARE[i] = sqrt(E[i] / TotalCorVar[i])
#         @assert E[i] isa Float32
#         @assert AE[i] isa Float32
#         @assert RMSE[i] isa Float32
#         @assert ARE[i] isa Float32
#     end
#     # The value of objective function/The correletion between latent variables
#     for i=1:l
#         tmpObjective = 0.0
#         tmpAvgLatentCor = 0.0
#         for j=1:l
#             if i != j
#                 tmpObjective += sum(diag(V[i]'V[j]))
#                 for k=1:dim
#                     tmpAvgLatentCor += cor(V[i][:,k].+1e-10, V[j][:,k].+1e-10)
#                 end
#             end
#         end
#         push!(Objective, tmpObjective)
#         push!(AvgLatentCor, tmpAvgLatentCor/(dim*(l-1)))
#     end
#     # Return
#     return ["E"=>E, "AE"=>AE, "RMSE"=>RMSE, "ARE"=>ARE, "Explained Correlation/Covariance"=>ExpCorCov, "Total Correlation/Variance"=>TotalCorVar, "Objective"=> Objective, "AvgLatentCor"=> AvgLatentCor]
# end

# # Row vector
# function normalizex(x::Array{UInt32,1}, n::Number, stream, scale::AbstractString, pseudocount::Number, colmeanlist::AbstractString, colmeanvec::AbstractArray, colvarlist::AbstractString, colvarvec::AbstractArray)
#     # Input
#     if !(scale in ["log", "ftt", "raw"])
#         error("scale must be specified as log, ftt, or raw")
#     end

#     # Logscale, FTTscale, Raw
#     if scale == "log"
#         pc = UInt32(pseudocount)
#         xx = Vector{Float32}(undef, length(x))
#         @inbounds for i in 1:length(x)
#             xx[i] = log10.(x[i] .+ pc)
#         end
#     end
#     if scale == "ftt"
#         xx = Vector{Float32}(undef, length(x))
#         @inbounds for i in 1:length(x)
#             xx[i] = sqrt.(x[i]) .+ sqrt.(x[i] .+ 1.0f0)
#         end
#     end
#     if scale == "raw"
#         xx = convert(Vector{Float32}, x)
#     end

#     # Centering, Normalization
#     if (colmeanlist != "") && (colvarlist == "")
#         xx = xx .- colmeanvec
#     end
#     if (colmeanlist != "") && (colvarlist != "")
#         xx = xx .- colmeanvec ./ colvarvec
#     end
#     # Return
#     return xx
# end

# # Full Gradient
# function ∇f(W::AbstractArray, input::AbstractArray, D::AbstractArray, scale::AbstractArray, pseudocount::AbstractArray, colmeanlist::AbstractArray, colmeanvec::AbstractArray, colvarlist::AbstractArray, colvarvec::AbstractArray, stepsize::Number, offsetFull::Number, offsetStoch::Number, perm::Bool)
#     N, M = nm(input)
#     tmpN = zeros(UInt32, 1)
#     tmpM = zeros(UInt32, 1)
#     l = length(input)
#     tmpW = Any[]
#     x = Any[]
#     normx = Any[]
#     for i=1:l
#         push!(tmpW, zeros(Float32, size(W[i])))
#         push!(x, zeros(UInt32, M[i]))
#         push!(normx, zeros(Float32, M[i]))
#     end
#     stream = Any[]
#     # Stream
#     for i=1:l
#         push!(stream, ZstdDecompressorStream(open(input[i])))
#     end
#     # tmpN, tmpN
#     for i=1:l
#         read!(stream[i], tmpN)
#         read!(stream[i], tmpM)
#     end
#     # x, norm
#     for n=1:N[1]
#         for i=1:l
#             read!(stream[i], x[i])
#             normx[i] = normalizex(x[i], n, stream[i], scale[i], pseudocount[i], colmeanlist[i], colmeanvec[i], colvarlist[i], colvarvec[i])
#             if perm
#                 normx[i] .= normx[randperm(length(normx[i]))]
#             end
#         end
#         for i=1:l
#             # Full Gradient
#             tmpW[i] .+= offsetFull * ∇fn(i, W, normx, D, M, stepsize, offsetStoch)
#         end
#     end
#     # GD
#     for i=1:l
#         W[i] .= tmpW[i] / offsetFull
#     end
#     # close
#     for i=1:l
#         close(stream[i])
#     end
#     return W
# end

# # Stochastic Gradient
# function ∇fn(i::Number, W::AbstractArray, x::AbstractArray, D::AbstractArray, M::AbstractArray, stepsize::Number, offsetStoch::Number)
#     l = length(W)
#     tmp∇fn = zeros(Float32, size(W[i]))
#     for j=1:l
#         if i != j
#             tmp∇fn .= x[i] * (offsetStoch * x[j]'W[j] * D[j])
#         end
#     end
#     return 1/offsetStoch * stepsize * tmp∇fn * Float32(2 / M[i])
# end

# # sym
# function sym(Y::AbstractArray)
#     return (Y + Y') / 2
# end

# # Riemannian Gradient
# function Pw(Z::AbstractArray, W::AbstractArray)
#     l = length(W)
#     out = Any[]
#     for i=1:l
#         push!(out, Z[i] - W[i] * sym(W[i]'Z[i]))
#     end
#     return out
# end

# # Retraction
# function Retraction(W::AbstractArray)
#     l = length(W)
#     for i=1:l
#         W[i] = Array(qr!(W[i]).Q)
#     end
#     return W
# end


text:
# Julia wrapper for header: /usr/local/include/libavutil/cast5.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_cast5_alloc,
    av_cast5_init,
    av_cast5_crypt,
    av_cast5_crypt2

function av_cast5_alloc()
    ccall((:av_cast5_alloc,libavutil),Ptr{AVCAST5},())
end

function av_cast5_init(ctx,key,key_bits::Integer)
    ccall((:av_cast5_init,libavutil),Cint,(Ptr{AVCAST5},Ptr{UInt8},Cint),ctx,key,key_bits)
end

function av_cast5_crypt(ctx,dst,src,count::Integer,decrypt::Integer)
    ccall((:av_cast5_crypt,libavutil),Void,(Ptr{AVCAST5},Ptr{UInt8},Ptr{UInt8},Cint,Cint),ctx,dst,src,count,decrypt)
end

function av_cast5_crypt2(ctx,dst,src,count::Integer,iv,decrypt::Integer)
    ccall((:av_cast5_crypt2,libavutil),Void,(Ptr{AVCAST5},Ptr{UInt8},Ptr{UInt8},Cint,Ptr{UInt8},Cint),ctx,dst,src,count,iv,decrypt)
end

ext: jl

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
module ReliablePowerSimulations

#################################################################################
# Exports
export ThermalStandardUCOutages
export ThermalDispatchOutages
export HydroDispatchRunOfRiverOutage
export HydroDispatchReservoirStorageOutage
export HydroDispatchPumpedStoragewReservationOutage
export BookKeepingwReservationOutage
export EndOfPeriodEnergyTargetOutage
export RenewableOutageDispatch

#################################################################################
# Imports
using PowerSystems
import InfrastructureSystems
import Dates
import PowerSimulations
import PowerModels
import JuMP
import ParameterJuMP

const PSY = PowerSystems
const IS = InfrastructureSystems
const PM = PowerModels
const PSI = PowerSimulations
const PJ = ParameterJuMP

#################################################################################
# Includes
include("./core/variables.jl")
include("./core/initial_condition_types.jl")
include("./core/update_initial_conditions.jl")

include("device_models/devices/common/constraints_structs.jl")
include("device_models/devices/common/duration_constraint.jl")
include("device_models/devices/common/outage_constraint.jl")

include("./device_models/devices/thermal_generation.jl")
include("./device_models/devices/renewable_generation.jl")
include("./device_models/devices/hydro_generation.jl")
include("./device_models/devices/storage.jl")

include("./device_models/device_constructors/thermalgeneration_constructor.jl")
include("./device_models/device_constructors/hydrogeneration_constructor.jl")
include("./device_models/device_constructors/renewablegeneration_constructor.jl")
include("./device_models/device_constructors/storage_constructor.jl")

end # module


text:
#=
XYZ reader
This module contains functions that can read an XYZ file and return a
Python dictionary with its contents.

=#

"""
read_xyz(datapath_or_datastring, is_datafile=true)

  Read data in .xyz format, from either a file or a raw string.
    `datapath_or_datastring`: Either the path to the XYZ file (can be relative
                                            or absolute), or a string corresponding to the content
                                            of an XYZ file (including newline characters).
    `is_datafile`: Either True (default) if passing the filepath to the data,
                                         or False if passing a string of raw data.
    `return`: A list of the atoms in the order that
                   they appear on the file, stored in
                   objects with keys "symbol", "x", "y",
                   and "z".
"""
function read_xyz(datapath_or_datastring; is_datafile=true)
  if is_datafile
    records = decode(read(datapath_or_datastring),"UTF-8")
  else
    req = Base.download(datapath_or_datastring)
    records = decode(read(req),"UTF-8")
  end
  lines = split(records, "\n")

  atoms =Vector{Dict}(undef,0)

  for line in lines
    rs = match(
      r"^\s*([\w]+)\s+([\w\.\+\-]+)\s+([\w\.\+\-]+)\s+([\w\.\+\-]+)\s*", string(line)
    )
    if !(rs isa Nothing) && (length(rs.captures) == 4)
      atom = Dict(
        "symbol" => rs.captures[1],
        "x" => parse(Float64,rs.captures[2]),
        "y" => parse(Float64,rs.captures[3]),
        "z" => parse(Float64,rs.captures[4])
      )
      push!(atoms, atom)
    end
  end

  return atoms
end


text:
using Plots
using DifferentialEquations
using Catalyst
using ModelingToolkit
using PWS.Trajectories

sn = @reaction_network begin
    0.005, S --> ∅
    0.25, ∅ --> S
end

u0 = [50]
tspan = (0., 10.)
discrete_prob = DiscreteProblem(sn, u0, tspan)
jump_prob = JumpProblem(sn, discrete_prob, Direct())
sol = solve(jump_prob, SSAStepper())
signal = Trajectories.trajectory(sol)

plot(signal)
signal.u[3:end] .= [SVector{1}(0)]
signal.u

rn = @reaction_network begin
    0.01 * θ, S --> X + S
    0.01 * 50.0 * (1-θ), ∅ --> X
    0.01, X --> ∅
end θ

mutable struct SignalStepper
    signal::Trajectory
    index::Int
end

function (signal_stepper::SignalStepper)(integrator)
    if integrator.t == zero(integrator.t)
        signal_stepper.index = 1
    end
    integrator.u[1] = signal_stepper.signal.u[signal_stepper.index][1]
    signal_stepper.index += 1
end
stepper = SignalStepper(signal, 1)

signal_cb = PresetTimeCallback(sol.t, stepper)

u0r = [50, 50]
p = [1.0]
response_prob = DiscreteProblem(rn, u0r, tspan, p)
jump_prob_r = JumpProblem(rn, response_prob, Direct())
sol_r = solve(jump_prob_r, SSAStepper(), callback=signal_cb, tstops=sol.t)
plot(sol_r)

Trajectories.trajectory(sol_r)

react = reactions(rn)[1]
ratelaw = jumpratelaw(react)
s_rate = build_function(ratelaw, species(rn), params(rn))
s_ratefun = eval(s_rate)
s_ratefun(0.5, 1.0)

sol_r.t[end]
sol_r[end]

function build_rate_functions(rn::ReactionSystem)
    totalrate = sum((jumpratelaw(react, rxvars=[]) for react in reactions(rn)))
    print(totalrate)
    totalrate_fun = eval(build_function(totalrate, species(rn), params(rn), expression=Val{true}))
    rate_funs = [eval(build_function(jumpratelaw(react, rxvars=[]), species(rn), params(rn), expression=Val{true})) for react in reactions(rn)]
    (totalrate_fun, rate_funs)
end

function log_likelihood(solution::ODESolution, rn::ReactionSystem, totalrate_fun, rate_funs; θ=1.0)
    result = 0.0
    for i in Base.Iterators.drop(eachindex(solution), 1)
        dt = solution.t[i] - solution.t[i - 1]
        u = solution[i - 1]
        result += - dt * totalrate_fun(u, θ)
    end

    symbol_map = speciesmap(rn)

    reaction_props = zeros(length(solution) - 1)
    for (rate_fun, react) in zip(rate_funs, reactions(rn))
        net_change = zeros(Int64, numspecies(rn))
        for (species, change) in react.netstoich
            net_change[symbol_map[species]] = change
        end

        du = solution[2] - solution[1]
        for i in Base.Iterators.drop(eachindex(solution), 1)
            du .= solution.u[i] .- solution.u[i - 1]
            u = solution.u[i - 1]

            if du == net_change
                reaction_props[i-1] += rate_fun(u, θ) / totalrate_fun(u, θ)
            end
        end
    end

    for r in reaction_props
        if r > 0.0
            result += log(r)
        end
    end

    result
end

totalrate, rates = build_rate_functions(rn)
log_likelihood(sol_r, rn, totalrate, rates; θ=0.0)

totalrate([54, 50], 0.99)

function signal_response_pair(signal_problem, response_problem; θ=1.0)
    signal = solve(signal_problem, SSAStepper())
    signal_cb = PresetTimeCallback(signal.t, (integrator) -> integrator.u[1] = signal(integrator.t)[1])

    problem = remake(response_problem; p=[θ])
    response = solve(response_problem, SSAStepper(), callback=signal_cb, tstops=signal.t)
    response
end

result = signal_response_pair(jump_prob, jump_prob_r, θ=0.0)
log_likelihood(result, rn, totalrate, rates, 0.0)
plot(result)

ll_vals = Vector{Float64}()
θ_vals = rand(100000)
for θ in θ_vals
    result = signal_response_pair(jump_prob, jump_prob_r, θ=θ)
    ll = log_likelihood(result, rn, totalrate, rates; θ=θ)
    push!(ll_vals, ll)
end

scatter(θ_vals, ll_vals)

ll_vals_cond = Vector{Float64}()
for θ in ones(100000)
    result = signal_response_pair(jump_prob, jump_prob_r, θ=θ)
    ll = log_likelihood(result, rn, totalrate, rates, θ)
    push!(ll_vals_cond, ll)
end

histogram([ll_vals_cond, ll_vals], alpha=0.5)

mean(ll_vals) - mean(ll_vals_cond)S


text:
using StrideArrays
using ArrayInterface
using LoopVectorization
using YaoLocations
using BenchmarkTools
using BQCESubroutine
using StrideArrays
using BQCESubroutine.Utilities
using LoopVectorization
using ThreadingUtilities
using BQCESubroutine.Utilities: BitSubspace

function subspace_mul_generic!(S::AbstractStrideArray{<:Any, <:Any, T, 3}, indices, C, U::AbstractStrideArray{<:Any, <:Any, T, 3}, subspace, offset=0) where {T <: Base.HWReal}
    Bmax = size(S, 2)
    for s in 1:length(subspace)
        k = subspace[s]
        # idx = k .+ indices
        # _k = k - 1
        for _b ∈ 0:(Bmax-1) >>> 3
            b =    _b << 3;
            bmax = b + 8
            if bmax ≤ Bmax # full static block
                @avx for n ∈ 1:8, m ∈ axes(U, 2)
                    C_re_m_n = zero(T)
                    C_im_m_n = zero(T)
                    for i ∈ axes(U, 3)
                        j = k + indices[i] + offset
                        C_re_m_n += U[1,m,i] * S[1,n+b,j] - U[2,m,i] * S[2,n+b,j]
                        C_im_m_n += U[1,m,i] * S[2,n+b,j] + U[2,m,i] * S[1,n+b,j]
                    end
                    C[1,m,n] = C_re_m_n
                    C[2,m,n] = C_im_m_n
                end
                @avx for n ∈ 1:8, m ∈ axes(U, 2)
                    j = k + indices[m] + offset
                    S[1,n+b,j] = C[1,m,n]
                    S[2,n+b,j] = C[2,m,n]
                end
                # AmulB!(C_re, C_im, U_re, U_im,
            else # dynamic block
                Nmax = 8 + Bmax - bmax
                @avx for n ∈ 1:Nmax, m ∈ axes(U, 2)
                    C_re_m_n = zero(T)
                    C_im_m_n = zero(T)
                    for i ∈ axes(U, 3)
                        j = k + indices[i] + offset
                        C_re_m_n += U[1,m,i] * S[1,n+b,j] - U[2,m,i] * S[2,n+b,j]
                        C_im_m_n += U[1,m,i] * S[2,n+b,j] + U[2,m,i] * S[1,n+b,j]
                    end
                    C[1,m,n] = C_re_m_n
                    C[2,m,n] = C_im_m_n
                end
                @avx for n ∈ 1:Nmax, m ∈ axes(U, 2)
                    j = k + indices[m] + offset
                    S[1,n+b,j] = C[1,m,n]
                    S[2,n+b,j] = C[2,m,n]
                end
            end
        end
    end
    return S
end

S = @StrideArray rand(Float64, 2, 100, 1<<20)
U = @StrideArray rand(Float64, 2, 8, 8)
C = StrideArray{Float64}(undef, (StaticInt{2}(), ArrayInterface.size(U, 2), StaticInt{8}()))
locs = Locations((1, 3, 5))
subspace = bsubspace(20, locs)
comspace = bcomspace(20, locs)
indices = StrideArray{Int}(undef, (StaticInt{length(comspace)}(), ))
@simd ivdep for i in eachindex(indices)
    indices[i] = comspace[i] + 1
end

@benchmark subspace_mul_generic!($(copy(S)), indices, C, U, subspace)

Sr = subspace_mul_generic!(copy(S), indices, U, subspace)
So = reinterpret(reshape, ComplexF64, S);
Uo = reinterpret(reshape, ComplexF64, U);
S2 = BQCESubroutine.subspace_mul_generic!(copy(So), indices, Uo, subspace)
reinterpret(reshape, ComplexF64, Sr) ≈ S2

@benchmark BQCESubroutine.subspace_mul_generic!($(rand(ComplexF64, 100, 1<<20)), $indices, $(rand(ComplexF64, 1<<3, 1<<3)), $subspace)
@benchmark subspace_mul_generic!($(copy(S)), indices, C, U, subspace)

S = @StrideArray rand(Float64, 2, 100, 1<<20)
U = @StrideArray rand(Float64, 2, 8, 8)

ArrayInterface.size(S)
ArrayInterface.size(U)
Base.@gc_preserve begin
@allocated C = StrideArray{Float64}(undef, (StaticInt{2}(), ArrayInterface.size(U, 2), StaticInt{8}()))
end

StrideArray{Float64}(undef, (StaticInt{2}(), ArrayInterface.size(U, 2), StaticInt{8}()))


text:
"""
```julia
using DifferentialEquations, Plots
function lorenz!(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)
sol = solve(prob)
plot(sol)
```
"""
module DifferentialEquations

  using Reexport

  @reexport using DiffEqBase
  @reexport using DiffEqNoiseProcess
  @reexport using RecursiveArrayTools

  @reexport using SteadyStateDiffEq
  @reexport using StochasticDiffEq
  @reexport using OrdinaryDiffEq
  @reexport using BoundaryValueDiffEq
  using Sundials
  @reexport using DelayDiffEq

  @reexport using DiffEqCallbacks
  @reexport using DiffEqJump

  @reexport using DiffEqFinancial
  @reexport using MultiScaleArrays
  @reexport using DiffEqPhysics

  @reexport using DimensionalPlotRecipes
  @reexport using ParameterizedFunctions

  using LinearAlgebra

  import DiffEqBase: solve

  include("default_solve.jl")
  include("default_arg_parsing.jl")
  include("ode_default_alg.jl")
  include("sde_default_alg.jl")
  include("dae_default_alg.jl")
  include("dde_default_alg.jl")
  include("discrete_default_alg.jl")
  include("rode_default_alg.jl")
  include("steady_state_default_alg.jl")
  include("bvp_default_alg.jl")

end # module


text:
using Documenter, eQuilibrator

makedocs(
    modules = [eQuilibrator],
    clean = false,
    sitename = "eQuilibrator.jl",
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        assets = ["assets/favicon.ico"],
    ),
    authors = "The developers of eQuilibrator.jl",
    linkcheck = !("skiplinks" in ARGS),
    pages = [
        "Home" => "index.md",
        "Installation" => "installation.md",
        "Example usage" => "examples.md",
        "Function reference" => "functions.md",
    ],
)

deploydocs(
    repo = "github.com/stelmo/eQuilibrator.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "main",
    push_preview = true,
)


text:
# This is a minimal working example for running praff on several
# workers in Julia and saving them into a shared matrix

using Random
@everywhere using RAFF

# Set Debug for Logging
@everywhere using Logging
@everywhere using Base.CoreLogging

@everywhere global_logger(ConsoleLogger(stdout, Logging.Error))

@everywhere gmodel!(x, t_, g) = begin
    g[1] = exp(t_ * x[2])
    g[2] = t_ * x[1] * exp(t_ * x[2]);
end

@everywhere model(x, t) = x[1] * exp(t * x[2])

n = 2

np = 100

p = 70

# Set the seed for generating the same data
Random.seed!(123456789)

data, xSol, = RAFF.generateNoisyData(model, n, np, p)

praff(model, gmodel!, data, n, MAXMS=2)

println("True sol: $(xSol).")


text:
import BS
using Test

s = 42
k = 40
r = 0.1
sigma = 0.2
t = 0.5

@test BS.price(BS.EuropeanCall(), s, k, t, r, sigma) ≈ 4.759422392871542
@test BS.price(BS.EuropeanPut(), s, k, t, r, sigma) ≈ 0.8085993729000975
@test BS.delta(BS.EuropeanCall(), s, k, t, r, sigma) ≈ 0.7791312909426691
@test BS.delta(BS.EuropeanPut(), s, k, t, r, sigma) ≈ -0.22086870905733091
@test BS.theta(BS.EuropeanCall(), s, k, t, r, sigma) ≈ -4.559092194592626
@test BS.theta(BS.EuropeanPut(), s, k, t, r, sigma) ≈ -0.7541744965897705
@test BS.gamma(BS.EuropeanCall(), s, k, t, r, sigma) ≈ 0.04996267040591185
@test BS.gamma(BS.EuropeanPut(), s, k, t, r, sigma) ≈ 0.04996267040591185
@test BS.vega(BS.EuropeanCall(), s, k, t, r, sigma) ≈ 8.813415059602853
@test BS.vega(BS.EuropeanPut(), s, k, t, r, sigma) ≈ 8.813415059602853
@test BS.rho(BS.EuropeanCall(), s, k, t, r, sigma) ≈ 13.982045913360281
@test BS.rho(BS.EuropeanPut(), s, k, t, r, sigma) ≈ -5.042542576653999
@test BS.impvol(BS.EuropeanCall(), 4.12, s, k, t, r) ≈ 0.1135753892186858
@test BS.impvol(BS.EuropeanPut(), 4.12, s, k, t, r) ≈ 0.5246966268060681


text:
# sample from the Muller potential using ABOBA
using Plots
using Printf
using Random
using LinearAlgebra
using ForwardDiff

push!(LOAD_PATH,"../src/")

using JuBasicMD

include("potentials.jl")

β = 0.1;
γ = 0.1;
M = [1.5, 0.75];
q₀ = [-0.75, 1.5];
p₀ = [0.0, 0.0];
Δt = 1e-2
n_iters = 10^6;

V = x->Muller(x);
cfg = ForwardDiff.GradientConfig(V, q₀);
gradV! = (gradV, x)-> ForwardDiff.gradient!(gradV, V, x, cfg);

sampler = ABOBA(gradV!, β, γ, M, Δt);

Random.seed!(100);
Q₀ = copy(q₀);
P₀ = copy(p₀);
sample_trajectory!([Q₀, P₀], sampler, options=Options(n_iters=n_iters));
@printf("In Place Q after %d iterations: (%g, %g)\n",n_iters, Q₀[1], Q₀[2])
@printf("In Place P after %d iterations: (%g, %g)\n",n_iters, P₀[1], P₀[2])

Random.seed!(100);
Xvals = sample_trajectory([q₀, p₀], sampler, options=Options(n_iters=n_iters,n_save_iters=n_iters));
Q = Xvals[end][1]
P = Xvals[end][2]
@printf("Q after %d iterations: (%g, %g)\n",n_iters, Q[1], Q[2])
@printf("P after %d iterations: (%g, %g)\n",n_iters, P[1], P[2])

Random.seed!(100);
Xvals = sample_trajectory([q₀, p₀], sampler, options=Options(n_iters=n_iters));
Qvals = [X[1] for X in Xvals];
histogram2d([Q[1] for Q in Qvals], [Q[2] for Q in Qvals],normalize=true,color=:viridis)
xlims!(-1.5,1.5)
ylims!(-0.5, 2.0)
xlabel!("x")
ylabel!("y")

max_stars_repo_name: cortner/JuBasicMD

max_issues_repo_name: cortner/JuBasicMD

max_forks_repo_name: cortner/JuBasicMD


text:
using DataStructures
using Statistics

"""
A vantage-point tree. Implemented from Samet, H. (2006). *Foundations of
multidimensional and metric data structures*. San Francisco, California: Morgan
Kaufmann. This data structure allows for searching a metric space in a more
efficient way than comparing all points to each other. The "text" portion of the
name is to imply that this tree will work with text-based data, and, in fact,
is what it was tested to work with. In theory, however, it should work for
numerical data as well.

# Parameters
* **pivot** The pivot point chosen for the tree for a particular node; the
    remaining data points will be assined to either the left or right tree based
    on the median distance between them and this pivot point
* **d** A distance function that will produce a similarity score between two
    data points
* **left** The left subtree containing all the points that are less than the
    median distance away from ``pivot``
* **right** The right subtree containing all the points that are greater than or
    equal to the median distance away from ``pivot``
"""
mutable struct TextVPTree{T}
    pivot::T
    d
    r::Number
    left::TextVPTree{T}
    right::TextVPTree{T}

    """
        TextVPTree{T}(pivot::T, d) where T

    Inner constructor for a `TextVPTree`. There are two methods. The first
    leaves `r` undefined, as would be the case for a leaf node. The other takes
    a value for `r`, which would be used for nodes that aren't leaf nodes.
    """
    function TextVPTree{T}(pivot::T, d) where T
        new(pivot, d)
    end

    function TextVPTree{T}(pivot::T, d, r::Number) where T
        new(pivot, d, r)
    end
end

"""
    TextVPTree(items::Array, d)

Outer constructor for a `TextVPTree`. Takes in an array of items `items` and a
distance function `d` and proceeds to build a vantage-point tree from them.
"""
function TextVPTree(items::Array, d)

    if length(items) == 1 # leaf node
        @inbounds return TextVPTree{typeof(items[1])}(items[1], d)
    end

    # select a pivot
    pivotI = rand(1:length(items))
    @inbounds pivot = items[pivotI]
    @inbounds items = items[1:end .!= pivotI]

    # determine value of r
    distances = [d(pivot, item) for item in items]
    r = median(distances)

    # Assign items to left and right subtrees
    S1 = Vector()
    S2 = Vector()
    for (item, dist) in zip(items, distances)
        if dist < r
            push!(S1, item)
        else
            push!(S2, item)
        end
    end

    # create tree
    tree = TextVPTree{typeof(pivot)}(pivot, d, r)

    # recursively define left and right subtrees
    if length(S1) > 0
        tree.left = TextVPTree(S1, d)
    end

    if length(S2) > 0
        tree.right = TextVPTree(S2, d)
    end

    return tree
end

"""
    radiusSearch(tree::TextVPTree, query, epsilon)

Performs a search for all items in a VP tree `tree` that are within a radius
`epsilon` from a query `query`.

# Returns

A `Vector` of items that are within the given radius `epsilon`
"""
function radiusSearch(tree::TextVPTree, query, epsilon)

    results = Vector()
    d = tree.d
    r = isdefined(tree, :r) ? tree.r : -1
    p = tree.pivot
    q = query

    # add root node if possible
    if q != p && d(q, p) <= epsilon
        push!(results, p)
    end

    # at leaf node, so can't recurse any further
    if r == -1
        return results
    end

    # check if we need to recurse into left and right subtrees using
    # inequalities from the implementation text
    if isdefined(tree, :left) && max(d(q, p) - r, 0) <= epsilon
        results = vcat(results, radiusSearch(tree.left, query, epsilon))
    end

    if isdefined(tree, :right) && max(r - d(q, p), 0) <= epsilon
        results = vcat(results, radiusSearch(tree.right, query, epsilon))
    end

    return results
end

"""
    nneighbors(tree::TextVPTree, query, n)

Find the `n` nearest neighbors in a VP tree `tree` to a given query `query`.

# Returns

* A `PriorityQueue` of items where the keys are the items themselves and the values are the distances from the items to `query`; the `PriorityQueue` is defined such that small values have higher priorities than large ones
"""
function nneighbors(tree::TextVPTree, query, n)

    # Helper function that adds items to the `PriorityQueue` that will
    # eventually be returned
    function nneighbors_(tree::TextVPTree, query, n, pq)

        d = tree.d
        r = isdefined(tree, :r) ? tree.r : -1
        p = tree.pivot
        q = query

        # if the `PriorityQueue` is not yet full, add the item
        # otherwise, add the item if it's closer to `query` than the furthest
        # item in `pq`
        if length(pq) < n && ! haskey(pq, p) && p != q
            enqueue!(pq, p=>d(p, q))
        elseif d(p, q) < peek(pq)[2] && ! haskey(pq, p) && p != q
            enqueue!(pq, p=>d(p, q))
            dequeue!(pq)
        end

        # check if the left subtree may contain items closer to `query` than
        # are in the queue; if so, check the left subtree for neighbors
        if isdefined(tree, :left)
            rlo = 0
            rhi = r
            d1 = max(d(q, p) - rhi, rlo - d(q, p), 0)

            if length(pq) < n || d1 <= peek(pq)[2]
                pq = nneighbors_(tree.left, q, n, pq)
            end
        end

        # check if the right subtree may contain items closer to `query` than
        # are in the queue; if so, check the right subtree for neighbors
        if isdefined(tree, :right)
            rlo = r
            rhi = Inf
            d1 = max(d(q, p) - rhi, rlo - d(q, p), 0)

            if length(pq) < n || d1 <= peek(pq)[2]
                pq = nneighbors_(tree.right, q, n, pq)
            end
        end

        return pq
    end

    d = tree.d
    r = isdefined(tree, :r) ? tree.r : -1
    p = tree.pivot
    q = query

    # if this first element is not the same as `query`, add it to `pq`
    # otherwise, instantiate `pq` as an empty `PriorityQueue`
    local pq
    if q != p
        pq = PriorityQueue(Base.Order.Reverse, p=>d(p, q))
    else
        pq = PriorityQueue(Base.Order.Reverse)
    end

    # check if the left subtree may contain items closer to `query` than
    # are in the queue; if so, check the left subtree for neighbors
    if isdefined(tree, :left)
        rlo = 0
        rhi = r
        d1 = max(d(q, p) - rhi, rlo - d(q, p), 0)

        if length(pq) < n || d1 <= peek(pq)[2]
            pq = nneighbors_(tree.left, q, n, pq)
        end
    end

    # check if the right subtree may contain items closer to `query` than
    # are in the queue; if so, check the right subtree for neighbors
    if isdefined(tree, :right)
        rlo = r
        rhi = Inf
        d1 = max(d(q, p) - rhi, rlo - d(q, p), 0)

        if length(pq) < n || d1 <= peek(pq)[2]
            pq = nneighbors_(tree.right, q, n, pq)
        end
    end

    return pq
end


text:
import IncrementalInference: getSample
import Base: convert, show
# import DistributedFactorGraphs: getManifold

using .Images

export ScatterAlignPose2, PackedScatterAlignPose2

export overlayScatter, overlayScatterMutate

"""
$TYPEDEF

This is but one incarnation for how radar alignment factor could work, treat it as a starting point.

Notes
- Stanard `cvt` argument is lambda function to convert incoming images to user convention of image axes,
  - default `cvt` flips image rows so that Pose2 xy-axes corresponds to img[x,y] -- i.e. rows down and across from top left corner.
- Use rescale to resize the incoming images for lower resolution (faster) correlations
- Both images passed to the construct must have the same type some matrix of type `T`.

Example
-------
```julia
arp2 = ScatterAlignPose2(img1, img2, 2) # e.g. 2 meters/pixel
```

See also: [`overlayScanMatcher`](@ref)
"""
struct ScatterAlignPose2{H1<:HeatmapGridDensity,H2<:HeatmapGridDensity} <: IIF.AbstractManifoldMinimize
  """ reference image for scan matching. """
  hgd1::H1
  """ test image to scan match against the reference image. """
  hgd2::H2
  """ Common grid scale for both images -- i.e. units/pixel.
  Constructor uses two arguments `gridlength`*`rescale=1`=`gridscale`.
  Arg 0 < `rescale` ≤ 1 is also used to rescale the images to lower resolution for speed. """
  gridscale::Float64
  """ how many heatmap sampled particles to use for mmd ailgnment """
  sample_count::Int
  """ bandwidth to use for mmd """
  bw::Float64
end

  # replace inner constructor with transform on image
  ScatterAlignPose2(im1::AbstractMatrix{T},
                    im2::AbstractMatrix{T},
                    domain::Tuple{<:AbstractVector{<:Real},<:AbstractVector{<:Real}};
                    sample_count::Integer=75,
                    bw::Real=5e-5, # from a sensitivity analysis with marine radar data (50 or 100 samples)
                    rescale::Real=1,
                    N::Integer=1000,
                    cvt = (im)->reverse(Images.imresize(im,trunc.(Int, rescale.*size(im))),dims=1)
                  ) where {T} = ScatterAlignPose2(HeatmapGridDensity(cvt(im1),domain,N=N),
                                                  HeatmapGridDensity(cvt(im2),domain,N=N),
                                                  float(rescale),
                                                  sample_count, bw  )
#

getManifold(::IIF.InstanceType{<:ScatterAlignPose2}) = getManifold(Pose2Pose2)

function getSample( cf::CalcFactor{<:ScatterAlignPose2} )

  pts1, = sample(cf.factor.hgd1.densityFnc, cf.factor.sample_count)
  pts2, = sample(cf.factor.hgd1.densityFnc, cf.factor.sample_count)

  M = getManifold(Pose2)
  e0 = ProductRepr(SA[0 0.0], SMatrix{2,2}(1.0, 0, 0, 1)) # identity_element(M)

  # precalc SE2 points
  R0 = e0.parts[2]

  ## FIXME REMOVE DEBUG + [20.0;-10]
  pts2_ = map(pt->ProductRepr(pt, R0), pts2)

  return pts1, pts2_, M #, e0
end

function (cf::CalcFactor{<:ScatterAlignPose2})(Xtup, p, q)
  #

  # will move to CalcFactor, follow IIF #1415
  M = Xtup[3]

  pts1, pts2_ = Xtup[1], Xtup[2]

  # get the current relative transform estimate
  pTq = Manifolds.compose(M, inv(M,p), q)

  # move other points with relative transform
  pts2__ = map(i->Manifolds.compose(M, pTq, pts2_[i]).parts[1], 1:length(pts1))

  # cf._residual[1] =
  return Float64[mmd(M.manifold[1], pts1, pts2__, bw=SA[cf.factor.bw;]);]
end

function Base.show(io::IO, sap::ScatterAlignPose2{H1,H2}) where {H1,H2}
  printstyled(io, "ScatterAlignPose2{", bold=true, color=:blue)
  println(io)
  printstyled(io, "    H1 = ", color=:magenta)
  println(io, H1)
  printstyled(io, "    H2 = ", color=:magenta)
  println(io, H2)
  printstyled(io, " }", color=:blue, bold=true)
  println(io)
  println(io, "  size(.data):     ", size(sap.hgd1.data))
  println(io, "    min/max:       ", round(minimum(sap.hgd1.data),digits=3), "/", round(maximum(sap.hgd1.data),digits=3))

  println(io, "  size(.data):     ", size(sap.hgd2.data))
  println(io, "    min/max:       ", round(minimum(sap.hgd2.data),digits=3), "/", round(maximum(sap.hgd2.data),digits=3))

  println(io, "  gridscale:       ", sap.gridscale)
  println(io, "  sample_count:    ", sap.sample_count)
  println(io, "  bw:              ", sap.bw)
  nothing
end

Base.show(io::IO, ::MIME"text/plain", sap::ScatterAlignPose2) = show(io, sap)
Base.show(io::IO, ::MIME"application/juno.inline", sap::ScatterAlignPose2) = show(io, sap)

"""
    $SIGNATURES

Overlay the two images from `AlignRadarPose2` with the first (red) fixed and transform the second image (blue) according to `tf`.

Notes:
- `tf` is a Manifolds.jl type `::ProductRepr` (or newer `::ArrayPartition`) to represent a `SpecialEuclidean(2)` manifold point.
"""
function overlayScatter(sap::ScatterAlignPose2,
                        trans::AbstractVector{<:Real}=[0;0.0],
                        rot::Real=0.0;
                        user_coords = [trans; rot],
                        offsetTrans::AbstractVector{<:Real}=[0;0.0],
                        user_offset = [offsetTrans;0.0],
                        score=Ref(0.0),
                        sample_count::Integer=sap.sample_count,
                        showscore::Bool=true,
                        findBest::Bool=true  )
  #

  pts1_, = sample(sap.hgd1.densityFnc, sample_count)
  pts2_, = sample(sap.hgd2.densityFnc, sample_count)

  # artificial offset for testing
  # TODO must still add the rotation user_offset component (not critical)
  for pt in pts2_
    pt .+= user_offset[1:2]
  end

  M = getManifold(Pose2Pose2)
  e0 = identity_element(M)
  # not efficient, but okay for here
  pTq(xyr=user_coords) = exp(M, e0, hat(M, e0, xyr))

  R0 = e0.parts[2]
  _pts2_nt = map(pt->ProductRepr(pt, R0), pts2_)

  # add user_offset
  uo = exp(M, e0, hat(M, e0, user_offset))
  _pts2_ = map(pt->Manifolds.compose(M, uo, pt), _pts2_nt)
  for (i,pt) in enumerate(_pts2_)
    pt.parts[2] .= R0
    # overwrite original for updated return values
    pts2_[i] .= pt.parts[1]
  end

  # take p as reference identity
  pts2T_ = map(pt->Manifolds.compose(M, pTq(), pt).parts[1], _pts2_)
  best_coords = zeros(3)

  if findBest
    # keep tfg separately so that optim can be more efficient
    tfg = initfg()
    ev_(xyr) = calcFactorResidualTemporary( sap, (Pose2,Pose2),
                                            (pts1_,_pts2_,M),
                                            (e0,pTq(xyr));
                                            tfg )[1]
    #
    score[] = ev_(user_coords)
    best = Optim.optimize(ev_, user_coords)
    best_coords .= round.(best.minimizer,digits=3)

    if showscore
      @info "overlayScatter score" score[] string(best_coords) best.minimum user_offset
    end
  end

  pts2T_best = map(pt->Manifolds.compose(M, pTq(best_coords), pt).parts[1], _pts2_)

  @cast pts1[i,j] := pts1_[j][i]
  @cast pts2[i,j] := pts2_[j][i]
  @cast pts2T[i,j] := pts2T_[j][i]
  @cast pts2Tb[i,j] := pts2T_best[j][i]

  (;pPp=pts1, qPq=pts2, pPq=pts2T, pPq_best=pts2Tb, best_coords, user_coords, user_offset)
end

function overlayScatterMutate(sap_::ScatterAlignPose2;
                              sample_count::Integer = sap_.sample_count,
                              bw::Real = sap_.bw,
                              kwargs... )
  #
  sap = ScatterAlignPose2(sap_.hgd1, sap_.hgd2, sap_.gridscale, sample_count, float(bw))
  Caesar.overlayScatter(sap; kwargs...);
end

## =========================================================================================
## Factor serialization below
## =========================================================================================

Base.@kwdef struct PackedScatterAlignPose2 <: PackedInferenceType
  hgd1::PackedHeatmapGridDensity
  hgd2::PackedHeatmapGridDensity
  gridscale::Float64 = 1.0
  sample_count::Int = 50
  bw::Float64 = 0.01
end

function convert(::Type{<:PackedScatterAlignPose2}, arp::ScatterAlignPose2)
  PackedScatterAlignPose2(;
    hgd1 = convert(PackedHeatmapGridDensity,arp.hgd1),
    hgd2 = convert(PackedHeatmapGridDensity,arp.hgd2),
    gridscale = arp.gridscale,
    sample_count = arp.sample_count,
    bw = arp.bw )
end

function convert(::Type{<:ScatterAlignPose2}, parp::PackedScatterAlignPose2)
  ScatterAlignPose2(
    covert(HeatmapGridDensity,parp.hgd1),
    covert(HeatmapGridDensity,parp.hgd2),
    parp.gridscale,
    parp.sample_count,
    parp.bw )
end

#


text:
using TransformVariables
using TransformVariables: AbstractTransform, CallableTransform, CallableInverse

export Pushforward
export Pullback

struct Pushforward{F,M} <: AbstractMeasure
    f::F
    μ::M
    logjac::Bool
end

Pushforward(f,μ) = Pushforward(f, μ, true)

struct Pullback{F,M} <: AbstractMeasure
    f::F
    ν::M
    logjac::Bool
end

Pullback(f,ν) = Pullback(f, ν, true)

function logdensity(pb::Pullback{F}, x) where {F <: CallableTransform}
    f = pb.f
    ν = pb.ν
    if pb.logjac
        y, logJ = transform_and_logjac(f.t, x)
        return logdensity(ν, y) + logJ
    else
        y = f(x)
        return logdensity(ν, y)
    end
end

function logdensity(pf::Pushforward{F}, y) where {F <: CallableTransform}
    f = pf.f
    μ = pf.μ
    x = inverse(f.t)(y)
    if pf.logjac
        _, logJ = transform_and_logjac(f.t, x)
        return logdensity(μ, x) - logJ
    else
        return logdensity(μ, x)
    end
end

Pullback(f::AbstractTransform, ν, logjac::Bool=true) = Pullback(transform(f), ν, logjac)
Pushforward(f::AbstractTransform, ν, logjac::Bool=true) = Pushforward(transform(f), ν, logjac)

Pullback(f::CallableInverse, ν, logjac::Bool=true) = Pushforward(transform(f.t), ν, logjac)

Pushforward(f::CallableInverse, ν, logjac::Bool=true) = Pullback(transform(f.t), ν, logjac)

Base.rand(rng::AbstractRNG, T::Type, ν::Pushforward) = ν.f(rand(rng, T, ν.μ))

Base.rand(rng::AbstractRNG, T::Type, μ::Pullback) = μ.f(rand(rng, T, μ.ν))

testvalue(ν::Pushforward) = transform(ν.f, testvalue(ν.μ))

testvalue(μ::Pullback) = transform(inverse(μ.f), testvalue(μ.ν))

basemeasure(μ::Pullback) = Pullback(μ.f, basemeasure(μ.ν), false)

basemeasure(ν::Pushforward) = Pushforward(ν.f, basemeasure(ν.μ), false)

TV.as(ν::Pushforward) = ν.f ∘ as(ν.μ)

TV.as(μ::Pullback) = inverse(μ.f) ∘ μ.ν

TV.as(::Lebesgue) = asℝ

# TODO: Make this work for affine embeddings
TV.as(d::Affine) = _as_affine(_firstval(d))

_firstval(d::Affine) = first(values(getfield(getfield(d, :f), :par)))
_as_affine(x::Real) = asℝ
_as_affine(x::AbstractArray) = as(Vector, size(x,1))

basemeasure(::Pushforward{TV.CallableTransform{T}, Lebesgue{ℝ}}) where {T <: TV.ScalarTransform} = Lebesgue(ℝ)
basemeasure(::Pullback{TV.CallableTransform{T}, Lebesgue{ℝ}}) where {T <: TV.ScalarTransform} = Lebesgue(ℝ)
# t = as𝕀
# μ = Normal()
# ν = Pushforward(t, μ)
# x = rand(μ)
# julia> logdensity(μ, x) ≈ logdensity(Pushforward(inverse(t), ν), x)
# true


text:
using .KernelFunctions

function DeterminantalPointProcess(kernel::Kernel, X::AbstractVector; parallelize=false)
    return DeterminantalPointProcess(kernelmatrix(kernel, X); parallelize=parallelize)
end

function DeterminantalPointProcess(kernel::Kernel, X::AbstractMatrix; obsdim=1, parallelize=false)
    return DeterminantalPointProcess(kernel, KernelFunctions.vec_of_vecs(X; obsdim=obsdim); parallelize=parallelize)
end


text:
#
function psi_compute(T::Int64, q::Int64; itrend::Int64=1)
  fvec = π*(1:1:q);
  tvec = (0.5:1:T)/T;
  psi = (sqrt(2))*(cos(tvec*fvec'))./(sqrt(T));
  fvec_2t = fvec./(2*T)
  iota = sin(fvec_2t)./fvec_2t
  if itrend == 1
    psi = psi.*iota'
  end
return psi
end
#


text:
#=
DPM.jl

10/27/2015
Adham Beyki, odinay@gmail.com
=#

#############################################
###### Dirichlet Process Mixture Model ######
#############################################
type DPM{T} <: MixtureModel
    component::T
    KK::Int
    aa::Float64
    aa1::Float64
    aa2::Float64

    DPM{T}(c::T, KK::Int, aa::Float64, aa1::Float64, aa2::Float64) = new(c, KK, aa, aa1, aa2)
end
DPM{T}(c::T, KK::Int, aa::Real, aa1::Real, aa2::Real) = DPM{typeof(c)}(c, KK, convert(Float64, aa),
    convert(Float64, aa1), convert(Float64, aa2))

function Base.show(io::IO, dpm::DPM)
    println(io, "Dirichlet Process Mixture Model with $(dpm.KK) $(typeof(dpm.component)) components")
end

"""
resampling concentration parameter based on Escobar and West 1995
"""
function sample_hyperparam!(dpm::DPM, NN::Int, iters::Int)

    for n = 1:iters
        eta_pp = rand(Distributions.Beta(dpm.aa+1, NN))
        rr = (dpm.aa1 + dpm.KK - 1) / (n*(dpm.aa2 - log(eta_pp)))
        pi_eta = rr / (1.0 + rr)

        if rand() < pi_eta
            dpm.aa = rand(Distributions.Gamma(dpm.aa1+dpm.KK)) / (dpm.aa2-log(eta_pp))
        else
            dpm.aa = rand(Distributions.Gamma(dpm.aa1+dpm.KK-1)) / (dpm.aa2-log(eta_pp))
        end
    end
end

count_active_clusters(nn::Vector{Int}) = length(find(x -> x>0, nn))

function storesample{T}(
    dpm::DPM{T},
    KK_list::Vector{Int},
    KK_dict::Dict{Int, Vector{Int}},
    sample_n::Int,
    filename::ASCIIString)

    if endswith(filename, "_")
        dummy_filename = string(filename, sample_n, ".h5")
    else
        dummy_filename = string(filename, "_", sample_n, ".h5")
    end

    KK_dict_keys = collect(keys(KK_dict))

    println("storing on disk...")
    HDF5.h5open(dummy_filename, "w") do file
        for kk in KK_dict_keys
            HDF5.write(file, "zz/$(kk)", KK_dict[kk])
        end
        HDF5.write(file, "sample_n", sample_n)
    end
end

function collapsed_gibbs_sampler{T1, T2}(
    dpm::DPM{T1},
    xx::Vector{T2},
    zz::Vector{Int},
    n_burnins::Int, n_lags::Int, n_samples::Int,
    sample_hyperparam::Bool=true, n_internals::Int=10,
    store_every::Int=250, results_path="", filename::ASCIIString="DPM_results")

    n_iterations = n_burnins + (n_samples)*(n_lags+1)
    NN = length(xx)
    nn = zeros(Int, dpm.KK)
    log_likelihood = 0.0

    components = Array(typeof(dpm.component), dpm.KK)
    for kk = 1:dpm.KK
        components[kk] = deepcopy(dpm.component)
    end

    KK_list = Int[]
    KK_dict = Dict{Int, Vector{Int}}()

    # initializing the model
    tic()
    for ii = 1:NN
        kk = zz[ii]
        additem!(components[kk], xx[ii])
        nn[kk] += 1
        log_likelihood += loglikelihood(components[kk], xx[ii])
    end
    push!(KK_list, dpm.KK)
    elapsed_time = toq()

    #=
    starting the MCMC chain
    1. remote the data point
    2. sample for zz
    3. add the datapoint back
    =#

    for iteration = 1:n_iterations

        # verbose
        if iteration < n_burnins
            print_with_color(:blue, "Burning... ")
        end
        println(@sprintf("iteration: %d, KK=%d, KK mode=%d, aa=%.2f, time=%.2f, likelihood=%.2f",
            iteration, dpm.KK, indmax(hist(KK_list, .5:maximum(KK_list)+0.5)[2]),
            dpm.aa, elapsed_time, log_likelihood))
        log_likelihood = 0.0

        tic()
        @inbounds for ii = randperm(NN)

            # 1
            kk = zz[ii]
            delitem!(components[kk], xx[ii])
            nn[kk] -= 1

            # remove the cluster if it is empty
            if nn[kk] == 0
                println("\tcomponent $kk has become inactive")
                splice!(nn, kk)
                splice!(components, kk)
                dpm.KK -= 1

                # shifting the labels one cluster back
                idx = find(x -> x>kk, zz)
                zz[idx] -= 1
            end

            # 2
            pp = zeros(Float64, dpm.KK+1)
            @inbounds for kk = 1:dpm.KK
                pp[kk] = log(nn[kk]) + logpredictive(components[kk], xx[ii])
            end
            pp[dpm.KK+1] = log(dpm.aa) + logpredictive(dpm.component, xx[ii])
            lognormalize!(pp)
            kk = sample(pp)

            # instanciate a new component if needed
            if kk == dpm.KK+1
                println("\tcomponent $(kk) activated.")
                push!(components, deepcopy(dpm.component))
                push!(nn, 0)
                dpm.KK += 1
            end

            # 3
            zz[ii] = kk
            nn[kk] += 1
            additem!(components[kk], xx[ii])
            log_likelihood += loglikelihood(components[kk], xx[ii])
        end # ii

        # sample the hyperparameter
        if sample_hyperparam
            sample_hyperparam!(dpm, NN, n_internals)
        end
        elapsed_time = toq()

        # save the sample
        if (iteration-n_burnins) % (n_lags+1) == 0 &&  iteration > n_burnins
            sample_n = convert(Int, (iteration-n_burnins)/(n_lags+1))
            push!(KK_list, dpm.KK)
            KK_dict[dpm.KK] = deepcopy(zz)

            if sample_n % store_every == 0
                normalized_filename = normpath(joinpath(results_path, filename))
                storesample(dpm, KK_list, KK_dict, sample_n, normalized_filename)
            end
        end
    end # iteration

    sample_n = convert(Int, (n_iterations-n_burnins)/(n_lags+1))
    normalized_filename = normpath(joinpath(results_path, filename))
    storesample(dpm, KK_list, KK_dict, sample_n, normalized_filename)

    KK_list, KK_dict
end

function truncated_gibbs_sampler{T1, T2}(
    dpm::DPM{T1},
    KK_truncation::Int,
    xx::Vector{T2},
    zz::Vector{Int},
    n_burnins::Int, n_lags::Int, n_samples::Int,
    sample_hyperparam::Bool=true, n_internals::Int=10,
    store_every::Int=100, results_path="", filename::ASCIIString="DPM_results_")

    n_iterations = n_burnins + (n_samples)*(n_lags+1)
    NN = length(xx)
    nn = zeros(Int, KK_truncation)
    pp = zeros(Float64, KK_truncation)
    log_likelihood = 0.0

    components = Array(typeof(dpm.component), KK_truncation)
    for kk = 1:KK_truncation
        components[kk] = deepcopy(dpm.component)
    end

    KK_list = Int[]
    KK_dict = Dict{Int, Vector{Int}}()

    tic()
    for ii = 1:NN
        kk = zz[ii]
        additem!(components[kk], xx[ii])
        nn[kk] += 1
        log_likelihood += loglikelihood(components[kk], xx[ii])
    end
    push!(KK_list, count_active_clusters(nn))

    # beta variables for mixing proprtions, Ishwaran and James 2001
    pi_tilde = ones(Float64, KK_truncation)
    for kk = 1:KK_truncation-1
        gamma1 = 1 + nn[kk]
        gamma2 = dpm.aa + sum(nn[kk+1:KK_truncation])
        pi_tilde[kk] = gamma1 / (gamma1 + gamma2)
    end
    pi_tilde[KK_truncation] = 1.0

    elapsed_time = toq()

    #########################################################
    ################ starting the MCMC chain ################
    #########################################################

    for iteration = 1:n_iterations

        # Verbose
        if iteration < n_burnins
            print_with_color(:blue, "Burning... ")
        end
        println(@sprintf("iteration: %d, KK=%d, KK mode=%d, aa=%.2f, time=%.2f, likelihood=%.2f",
            iteration, dpm.KK, indmax(hist(KK_list, .5:maximum(KK_list)+0.5)[2]), dpm.aa, elapsed_time, log_likelihood))
        log_likelihood = 0.0

        tic()
        for ii = randperm(NN)

            # 1
            # remove the point
            kk = zz[ii]
            delitem!(components[kk], xx[ii])
            nn[kk] -= 1

            # 2
            # sample zz
            pp = zeros(Float64, KK_truncation)
            for kk = 1:KK_truncation
                pp[kk] = log(pi_tilde[kk]) + sum(log(1-pi_tilde[1:kk-1])) + logpredictive(components[kk], xx[ii])
            end
            lognormalize!(pp)
            kk = sample(pp)

            # 3
            # add the point to the newly resampled cluster
            zz[ii] = kk
            nn[kk] += 1
            additem!(components[kk], xx[ii])
            log_likelihood += loglikelihood(components[kk], xx[ii])
        end # ii

        # 4
        # sample pi_tilde
        for kk = 1:KK_truncation-1
            gamma1 = 1 + nn[kk]
            gamma2 = dpm.aa + sum(nn[kk+1:KK_truncation])
            pi_tilde[kk] = rand(Distributions.Beta(gamma1, gamma2))
        end
        pi_tilde[KK_truncation] = 1.0

        dpm.KK = count_active_clusters(nn)

        # 5
        # sample the hyperparameter
        if sample_hyperparam
            sample_hyperparam!(dpm, NN, n_internals)
        end
        elapsed_time = toq()

        # save the sample
        if (iteration-n_burnins) % (n_lags+1) == 0 &&  iteration > n_burnins
            sample_n = convert(Int, (iteration-n_burnins)/(n_lags+1))

            push!(KK_list, dpm.KK)
            KK_dict[dpm.KK] = deepcopy(zz)

            if sample_n % store_every == 0
                normalized_filename = normpath(joinpath(results_path, filename))
                storesample(dpm, KK_list, KK_dict, sample_n, normalized_filename)
            end
        end

    end # iteration

    sample_n = convert(Int, (n_iterations-n_burnins)/(n_lags+1))
    normalized_filename = normpath(joinpath(results_path, filename))
    storesample(dpm, KK_list, KK_dict, sample_n, normalized_filename)

    KK_list, KK_dict
end

#######################
###    Posterior    ###
#######################

# posterior for DPM with collapsed Gibbs sampling
function posterior{T1, T2}(
    dpm::DPM{T1},
    xx::Vector{T2},
    zz::Vector{Int},
    KK::Int)

    components = Array(typeof(dpm.component), KK)
    for kk = 1:KK
        components[kk] = deepcopy(dpm.component)
    end

    NN = length(xx)
    nn = zeros(Int, KK)

    # since in truncated Gibbs sampling we might have a case where the sampler jumps over
    # a cluster and doesn't allocate any point to it, we have to map them back
    zz_unq = sort(unique(zz))
    for ii = 1:NN
        kk = find(x -> x==zz[ii], zz_unq)[1]
        additem!(components[kk], xx[ii])
        nn[kk] += 1
    end

    return([posterior(components[kk]) for kk=1:KK], nn)
end


text:
struct Householder <: SRSolver end

function impliedVolatilitySRHalley(
    isCall::Bool,
    price::T,
    f::T,
    strike::T,
    tte::T,
    df::T,
    ftolrel::T,
    maxEval::Int,
    solver::Householder,
)::T where {T}
    c, ex = normalizePrice(isCall, price, f, strike, df)
    if c >= 1 / ex || c <= 0
        throw(DomainError(c, string("Price out of range, must be < ", 1 / ex, " and > 0")))
    end
    x = log(ex)
    guess = T(impliedVolatilitySRGuessUndiscountedCall(c, ex,x))
    b = guess
    xtolrel = 32 * eps(T)
    xtol = T(0)
    ftolrel = max(ftolrel, eps(T))
    ftol = ftolrel * max(1, c)
    logc = log(c)
    fb, fbOverfpb, fp2bOverfpb, fp3bOverfpb = objectiveHouseholderLog(x, ex, b, logc)
    if abs(fb) < ftol
        return b / sqrt(tte)
    end
    for iteration = 1:maxEval
        x0 = b
        lf = (fbOverfpb) * (fp2bOverfpb)
        hn = -fbOverfpb
        num = (1 + fp2bOverfpb * hn / 2)
        denom = (1 + fp2bOverfpb * hn + fp3bOverfpb / 6 * hn * hn)
        x1 = x0 + hn * num / denom
        a = x0
        fb, fbOverfpb, fp2bOverfpb, fp3bOverfpb = objectiveHouseholderLog(x, ex, x1, logc)
        b = x1
        xtol_ = xtol + max(1, abs(b)) * xtolrel
        if abs(b - a) <= xtol_ || abs(fb) < ftol
            break
        end
    end
    return b / sqrt(tte)
end

function objectiveHouseholderLog(x::T, ex::T, v::T, logc::T) where {T}
    v = abs(v)
    h = x / v
    t = v / 2
    h2 = h^2
    t2 = t^2
    sqrt2 = sqrt(T(2))
    Np = erfcx(-(h + t) / sqrt2)
    Nm = erfcx(-(h - t) / sqrt2)
    eh2t2 = exp(-(h2 + t2) / 2)
    norm = 1 / (2 * sqrt(ex)) * eh2t2
    cEstimate = norm * (Np - Nm)
    logcEstimate = log(cEstimate)
    twopi = 2 * T(pi)
    logvega = (2 /sqrt(twopi)) / (Np - Nm)  #u'/u
    volgaOverVega = (h + t) * (h - t) / v
    logvolgaOverVega = volgaOverVega - logvega  #u'' u - u'2 /u^2
    h2mt2 = h2 - t2
    c3OverVega = (-3 * h2 - t2 + h2mt2 * h2mt2) / (v^2)   #
    logc3overVega = c3OverVega - 3 * logvega * volgaOverVega + 2 * logvega * logvega
    return logcEstimate - logc, (logcEstimate - logc) / logvega, logvolgaOverVega, logc3overVega
end

function objectiveHouseholder(x::T, ex::T, v::T, c::T) where {T}
    v = abs(v)
    h = x / v
    t = v / 2
    h2 = h^2
    t2 = t^2
    sqrt2 = sqrt(T(2))
    Np = erfcx(-(h + t) / sqrt2)
    Nm = erfcx(-(h - t) / sqrt2)
    eh2t2 = exp(-(h2 + t2) / 2)
    norm = 1 / (2 * sqrt(ex)) * eh2t2
    cEstimate = norm * (Np - Nm)
    eht = exp(-h * t)
    twopi = 2 * T(pi)
    vega = eh2t2 * eht / sqrt(twopi)    #normpdf(h + t)  #this is exp(-0.5* (h^2+t^2 +2*h*t)) = eh2t2*exp(-h*t)
    cEstimateOverVega = (Np - Nm) / (2 * eht * sqrt(ex / (twopi))) - c / vega
    volgaOverVega = (h + t) * (h - t) / v
    h2mt2 = h2 - t2

    c3OverVega = (-3 * h2 - t2 + h2mt2 * h2mt2) / (v^2)
    return cEstimate - c, cEstimateOverVega, volgaOverVega, c3OverVega
end


text:
# power functions
function load_pwr_data(set)
    data_net = PowerModels.parse_file("data/$(set[:pwr_case])")
    gen_data = CSV.read("data/gen_data.csv", DataFrame; header=1)
    load_dyn = CSV.read("data/load_factor.csv", DataFrame; header=1)[:,1]
    # Network size
    G = size(gen_data,1)
    N = length(data_net["bus"])
    E = length(data_net["branch"])
    D = length(data_net["load"])

    # order bus indexing
    bus_keys=collect(keys(data_net["bus"]))
    bus_key_dict = Dict()
    for i in 1:N
        push!(bus_key_dict, i => bus_keys[i])
    end
    node(key) = [k for (k,v) in bus_key_dict if v == key][1]

    # Load generation data
    p̅ = zeros(G); p̲ = zeros(G); c1 = zeros(G); c2 = zeros(G); M_p = zeros(N,G)
    for g in 1:G
        p̅[g] = gen_data[g,:P_max]
        p̲[g] = gen_data[g,:P_min] * 0
        c1[g] = gen_data[g,:C]
        c2[g] = gen_data[g,:C_2]
        M_p[gen_data[g,:node],g] = 1
    end

    # Load demand data
    load_key=collect(keys(data_net["load"]))
    d = zeros(D); M_d = zeros(N,D)
    for h in load_key
        d[parse(Int64,h)] = data_net["load"][h]["pd"]*data_net["baseMVA"] + 1e-3
        M_d[node(string(data_net["load"][h]["load_bus"])),parse(Int64,h)] = 1
    end
    load_dyn = ones(set[:T])

    # Load transmission data
    line_key=collect(keys(data_net["branch"]))
    β = zeros(E); f̅ = zeros(E); n_s = trunc.(Int64,zeros(E)); n_r = trunc.(Int64,zeros(E))
    for l in line_key
        β[data_net["branch"][l]["index"]] = -imag(1/(data_net["branch"][l]["br_r"] + data_net["branch"][l]["br_x"]im))
        n_s[data_net["branch"][l]["index"]] = data_net["branch"][l]["f_bus"]
        n_r[data_net["branch"][l]["index"]] = data_net["branch"][l]["t_bus"]
        f̅[data_net["branch"][l]["index"]] = data_net["branch"][l]["rate_a"]*data_net["baseMVA"]
    end
    # merge parallel lines
    ff = zeros(N,N); ββ = zeros(N,N)
    for l in line_key
        ff[node(string(n_s[data_net["branch"][l]["index"]])),node(string(n_r[data_net["branch"][l]["index"]]))] += f̅[data_net["branch"][l]["index"]]
        ff[node(string(n_r[data_net["branch"][l]["index"]])),node(string(n_s[data_net["branch"][l]["index"]]))] += f̅[data_net["branch"][l]["index"]]
        ββ[node(string(n_s[data_net["branch"][l]["index"]])),node(string(n_r[data_net["branch"][l]["index"]]))]  = β[data_net["branch"][l]["index"]]
        ββ[node(string(n_r[data_net["branch"][l]["index"]])),node(string(n_s[data_net["branch"][l]["index"]]))]  = β[data_net["branch"][l]["index"]]
    end
    # find all parallel lines
    parallel_lines = []
    for l in line_key, e in line_key
        if l != e && node(string(n_s[data_net["branch"][l]["index"]])) == node(string(n_s[data_net["branch"][e]["index"]])) && node(string(n_r[data_net["branch"][l]["index"]])) == node(string(n_r[data_net["branch"][e]["index"]]))
            push!(parallel_lines,l)
        end
    end
    # update number of edges
    E = E - Int(length(parallel_lines)/2)
    # get s and r ends of all edge
    n_s = trunc.(Int64,zeros(E)); n_r = trunc.(Int64,zeros(E))
    ff = LowerTriangular(ff)
    for l in 1:E
        n_s[l] = findall(!iszero, ff)[l][1]
        n_r[l] = findall(!iszero, ff)[l][2]
    end
    β = zeros(E); f̅ = zeros(E);
    likely_congested_lines = [2 25 34 76 77 112 116 174]
    for l in 1:E
        β[l] = ββ[n_s[l],n_r[l]]
        l in likely_congested_lines ? f̅[l] = ff[n_s[l],n_r[l]]*3 : f̅[l] = ff[n_s[l],n_r[l]]*3
    end

    # Find reference node
    ref = 0
    for n in 1:N
        if sum(M_p[n,:]) == 0 &&  sum(M_d[n,:]) == 0 == 0
            ref = n
        end
    end

    # Compute PTDF matrix
    B_line = zeros(E,N); B̃_bus = zeros(N,N); B = zeros(N,N)
    for n in 1:N
        for l in 1:E
            if n_s[l] == n
                B[n,n] += β[l]
                B_line[l,n] = β[l]
            end
            if n_r[l] == n
                B[n,n] += β[l]
                B_line[l,n] = -β[l]
            end
        end
    end
    for l in 1:E
        B[Int(n_s[l]),Int(n_r[l])] = - β[l]
        B[Int(n_r[l]),Int(n_s[l])] = - β[l]
    end
    B̃_bus = remove_col_and_row(B,ref)
    B̃_bus = inv(B̃_bus)
    B̃_bus = build_B̆(B̃_bus,ref)
    PTDF = B_line*B̃_bus

    # power-gas data
    gas_bus = [10 12 15 25 26 61 65 66 69]
    power_node = [45 33 36 31 30 40 39 38 25]
    Λ = zeros(48,N)
    Λ = zeros(48,N)
    for i in 1:length(power_node)
        Λ[power_node[i],gas_bus[i]] = 0.25
    end

    @info("done loading power system data")

    # safe network data
    net = Dict(
    # transmission data
    :f̅ => f̅, :n_s => n_s, :n_r => n_r, :T => round.(PTDF,digits=8),
    # load data
    :d => round.(d,digits=5), :M_d => M_d, :k_l => load_dyn,
    # generation data
    :p̅ => round.(p̅,digits=5), :p̲ => round.(p̲,digits=5), :M_p => M_p,
    :c1 => round.(c1,digits=5), :c2 => round.(c2,digits=5), :c̀2 => round.(sqrt.(c2),digits=5),
    # graph data
    :N => N, :E => E, :G => G, :D => D, :ref => ref,
    :gas_bus => [10 12 15 25 26 61 65 66 69],
    :Λ => Λ
    )
    return net
end
function solve_chance_constrained_OPF(net,set,ξ)
    # DC-OPF definition
    model = Model(optimizer_with_attributes(Mosek.Optimizer, "LOG" => 0))
    # model variables
    P =  [@variable(model, [1:net[:G], 1:set[:k_t][t]]) for t = 1:set[:T]]
    v =  [@variable(model, [1:net[:N]]) for t = 1:set[:T]]
    # auxiliary variables
    A =  [@variable(model, [1:net[:N], 1:set[:k_t][t]]) for t = 1:set[:T]]
    B =  [@variable(model, [1:net[:E], 1:set[:k_t][t]]) for t = 1:set[:T]]
    y_f =  [@variable(model, [1:net[:E]]) for t = 1:set[:T]]
    x_f =  [@variable(model, [1:net[:E]]) for t = 1:set[:T]]
    y_p =  [@variable(model, [1:net[:G]]) for t = 1:set[:T]]
    x_p =  [@variable(model, [1:net[:G]]) for t = 1:set[:T]]
    # model objective (exp cost + var of non-CCGT units)
    @objective(model, Min, sum(net[:c1]'P[t]*S(t)*ξ[:μ̂] + 10000 * sum(v[t][i] for i in setdiff(1:net[:N],net[:gas_bus])) for t in 1:set[:T]))
    # prioritize CCGT plants to balance renewables
    @constraint(model, ϕ[t=1:set[:T],i=setdiff(1:net[:N],net[:gas_bus])],
        [v[t][i];ξ[:Σ̅]*(net[:M_p]*P[t]*S(t))[i,:]] in SecondOrderCone())
    # OPF equations
    @constraint(model, α[t=1:set[:T]], A[t] .== (net[:M_p]*P[t] .+ ξ[:M_w]*W(t) .- net[:M_d]*D(t)*net[:k_l][t]))
    @constraint(model, β[t=1:set[:T]], B[t] .== net[:T] * A[t])
    @constraint(model, λ[t=1:set[:T]], ones(net[:N])' * A[t] .== 0)
    # double-sided reformulation of power flow chance constraints
    @constraint(model, μ_1[t=1:set[:T], l=1:net[:E]], [sqrt(set[:ε̅_p]) * (net[:f̅] .- x_f[t])[l]; vcat(ξ[:Σ̅]*(B[t]*S(t))[l,:],y_f[t][l])] in SecondOrderCone())
    @constraint(model, μ_2[t=1:set[:T], l=1:net[:E]], (B[t] * S(t) * ξ[:μ̂])[l] <= y_f[t][l] + x_f[t][l])
    @constraint(model, μ_3[t=1:set[:T], l=1:net[:E]], (B[t] * S(t) * ξ[:μ̂])[l] >= - y_f[t][l] - x_f[t][l])
    @constraint(model, μ_4[t=1:set[:T], l=1:net[:E]], net[:f̅][l] >= x_f[t][l] >= 0)
    @constraint(model, μ_5[t=1:set[:T], l=1:net[:E]], y_f[t][l] >= 0)
    # double-sided generation (reformulated) chance constraints
    @constraint(model, κ_1[t=1:set[:T], i=1:net[:G]], [sqrt(set[:ε̅_p]) * (1/2*(net[:p̅] - net[:p̲]) - x_p[t])[i]; vcat(ξ[:Σ̅]*(P[t]*S(t))[i,:],y_p[t][i])] in SecondOrderCone())
    @constraint(model, κ_2[t=1:set[:T], i=1:net[:G]], (P[t]*S(t)*ξ[:μ̂])[i] - 1/2*(net[:p̅] - net[:p̲])[i] <= y_p[t][i] + x_p[t][i])
    @constraint(model, κ_3[t=1:set[:T], i=1:net[:G]], (P[t]*S(t)*ξ[:μ̂])[i] - 1/2*(net[:p̅] - net[:p̲])[i] >= - y_p[t][i] - x_p[t][i])
    @constraint(model, κ_4[t=1:set[:T], i=1:net[:G]], 1/2*(net[:p̅] - net[:p̲])[i] >= x_p[t][i] >= 0)
    @constraint(model, κ_5[t=1:set[:T], i=1:net[:G]], y_p[t][i] >= 0)
    # solve model
    optimize!(model)
    @info("done solving chance-constrained OPF in $(solve_time(model)) sec. with status $(termination_status(model))")
    sol = Dict(:status => termination_status(model),
                :obj => JuMP.objective_value(model),
                :P => [JuMP.value.(P[t]) for t = 1:set[:T]],
                :CPUtime => solve_time(model))
    return sol
end
function wind_forecast(set)
    # power system data
    w_f = CSV.read("data/wind.csv", DataFrame; header=0)[:,1]
    bus = [3 8 11 20 24 26 31 38 43 49 53]
    W_ = length(bus)
    M_w = zeros(net[:N],W_)
    for i in 1:W_
        M_w[bus[i],i] = 1
    end
    # uncertainty data
    k = set[:k_t][set[:T]]
    # mean and covariance
    μ = ones(k-1)
    Σ = zeros(k-1,k-1)
    for i in 1:k-1
        Σ[i,i] = set[:σ]
    end
    # first two moments of altered dimentions
    μ̂ = vcat(1,μ)
    # E[ξξᵀ] = Cov + μμᵀ
    Σ̂ = zeros(k,k)
    Σ̂[1,1] = 1; Σ̂[2:end,1] .= μ; Σ̂[1,2:end] .= μ;
    Σ̂[2:end,2:end] .=  Σ .+ μ*μ'
    # factorization of E[ξξᵀ]
    Σ̅ = cholesky(Σ̂).L .+ zeros(k,k)
    Σ̅[:,1] .= 0
    ξ̂ = ones(k,set[:Ns])
    ξ̂[2:end,:] = rand(MvNormal(μ,Σ),set[:Ns])
    return Dict(:w_f => w_f, :bus => bus, :W => W_, :M_w => M_w, :μ => μ, :Σ => Σ, :μ̂ => μ̂, :Σ̂ => Σ̂, :Σ̅ => Σ̅, :ξ̂ => ξ̂)
end