text:
## NIScope built-in high level functions
# Configure vertical
"""
niScope_ConfigureVertical(obj::T;ch::String="1",vrang=5.0,voffs=0.5,cpl="dc",att="10x",st="on") where {T<:NIScope}
"""
function niScope_ConfigureVertical(obj::T;ch::String="1",vrang=5.0,voffs=0.5,
	cpl="dc",att="10x",st="on") where {T<:NIScope}
	chlist = [UInt8.(collect("0")); UInt8(0)] # terminate with NULL char
	cpl = NIScope_Vertical_Dict[cpl]
	att = NIScope_Vertical_Dict[att]
	st = NIScope_Vertical_Dict[st]
	status = niScope_ConfigureVertical(obj.handle, chlist, vrang, voffs, cpl, att, st)
end

# Configure horizontal
"""
status = niScope_ConfigureHorizontalTiming(obj::T;minSampleRate=2e5,minNumPts=1000,refPosition=0.1,numRecords=1,enforceRealtime=Int16(1)) where {T<:NISCOPE}
"""
function niScope_ConfigureHorizontalTiming(obj::T;minSampleRate=2e5,
	minNumPts=1000,refPosition=0.1,numRecords=1,enforceRealtime=Int16(1)) where {T<:NIScope}
	status = niScope_ConfigureHorizontalTiming(obj.handle, minSampleRate, minNumPts, refPosition, numRecords, enforceRealtime)
end


text:
using package_example

package_example.TestC.print_vars()

ext: jl


text:
### WARNING This is an old ODE solver which is not based on an iterator implementation. We keep it until LSODA has an iterator implementation

"""
Same as the `solve` for `CHV(::DiffEqBase.DEAlgorithm)` but for `CHV(::Symbol)`. This is an old implementation of the CHV algorithm which can be used with `:lsoda`. For all other solvers, use the the new solver.
"""
function solve(problem::PDMPProblem, algo::CHV{Tode}; verbose::Bool = false, ind_save_d = -1:1, ind_save_c = -1:1, n_jumps = Inf64, reltol = 1e-7, abstol = 1e-9, save_positions = (false, true), save_rate = false, finalizer = finalize_dummy) where {Tode <: Symbol}
	verbose && println("#"^30)
	ode = algo.ode
	@assert ode in [:cvode, :lsoda, :adams, :BDF]
	verbose && printstyled(color=:red,"--> Start CHV method (algo::Symbol)\n")

	# table to use DiffEqBase
	odeTable = Dict(:lsoda => lsoda(), :BDF => CVODE_BDF(), :adams => CVODE_Adams(), :cvode => CVODE_BDF())

	# initialise the problem. If I call twice this solve function, it should give the same result...
	init!(problem)

	# we declare the characteristics for convenience
	caract = problem.caract
	ratecache = caract.ratecache

	ti, tf = problem.tspan
	n_jumps += 1 # to hold initial vector
	nsteps   = 1 # index for the current jump number

	xc0 = caract.xc0
	xd0 = caract.xd0

	# Set up initial simulation time
	t = ti

	X_extended = similar(xc0, length(xc0) + 1)
	for ii in eachindex(xc0)
		X_extended[ii] = xc0[ii]
	end
	X_extended[end] = ti

	#useful to use the same array, as it can be used in CHV(ode)
	Xd = caract.xd
	if ind_save_c[1] == -1
		ind_save_c = 1:length(xc0)
	end

	if ind_save_d[1] == -1
		ind_save_d = 1:length(xd0)
	end
	xc_hist = VectorOfArray([copy(xc0)[ind_save_c]])
	xd_hist = VectorOfArray([copy(xd0)[ind_save_d]])

	res_ode = zeros(length(X_extended))

	nsteps += 1

	probExtLsoda = ODEProblem((du, u, p, _t) -> algo(du, u, caract, _t), copy(X_extended), (ti, tf))

	function Flow(_X0, _Xd, Δt, _r; _alg = odeTable[ode])
		prob = DiffEqBase.remake(probExtLsoda; tspan = (0, Δt))
		prob.u0 .= _X0
		sol = solve(prob, _alg; abstol = abstol, reltol = reltol, save_everystep = false)
		return sol.u[end]
	end

	# we use the first time interval from the one generated by the constructor PDMPProblem
	δt = problem.simjptimes.tstop_extended

	# Main loop
	while (t < tf) && (nsteps < n_jumps)

		verbose && println("--> t = $t, δt = $δt, nstep =  $nsteps")

		res_ode .= Flow(X_extended, Xd, δt, ratecache.rate)

		verbose && println("--> ode solve has been performed!")

		if (res_ode[end] < tf) && nsteps < n_jumps
			# this is the next jump time
			t = res_ode[end]

			# this holds the new state of the continuous component
			@inbounds for ii in eachindex(X_extended)
				X_extended[ii] = res_ode[ii]
			end

			caract.R(ratecache.rate, X_extended, Xd, caract.parms, t, false)

			# Update event
			ev = pfsample(ratecache.rate)

			# we perform the jump, it changes Xd and (possibly) X_extended
			affect!(caract.pdmpjump, ev, X_extended, Xd, caract.parms, t)

			verbose && println("--> Which reaction? => ", ev)
			verbose && println("--> xd = ", Xd)

			# save state, post-jump
			if save_positions[2] || (nsteps == n_jumps - 1)
				pushTime!(problem, t)
				push!(xc_hist, copy(X_extended[ind_save_c]))
				push!(xd_hist, copy(Xd[ind_save_d]))
			end

			save_rate && push!(problem.rate_hist, caract.R(ratecache.rate, X_extended, Xd, caract.parms, t, true)[1])

			finalizer(ratecache.rate, caract.xc, caract.xd, caract.parms, t)

			δt = - log(rand())

		else
			probLast = ODEProblem((du, u, p, _t) -> caract.F(du, u, Xd, caract.parms, _t), X_extended[1:end-1], (t, tf))
			res_ode_last = solve(probLast, odeTable[ode]; abstol = 1e-9, reltol = 1e-7, save_everystep = false)

			t = tf

			# save state
			pushTime!(problem, tf)
			push!(xc_hist, copy(res_ode_last[end][ind_save_c]))
			push!(xd_hist, copy(Xd[ind_save_d]))
		end
		nsteps += 1
	end
	verbose && println("--> Done")
	if verbose && save_positions[2]
		println("--> xc = ", xd_hist[:, 1:nsteps-1])
	end
	return PDMPResult(copy(problem.time), xc_hist, xd_hist, problem.rate_hist, save_positions, length(problem.time), 0)
end


text:
# External modules
using SparseArrays
using LinearAlgebra
using Random

######################
# NOTE: This is not a long term solution.
# Once the package is set up, we can put this dependency in the .toml
try
    using SymEngine
catch e
    println("Julia SymEngine is not yet installed. Installing now.");
    using Pkg
    Pkg.add("SymEngine")
    using SymEngine
end
try
    using WriteVTK
catch e
    println("WriteVTK is not yet installed. Installing now.");
    using Pkg
    Pkg.add("WriteVTK")
    using WriteVTK
end
# try
#     using Latexify
# catch e
#     println("Latexify is not yet installed. Installing now.");
#     using Pkg
#     Pkg.add("Latexify")
#     using Latexify
# end
######################

# include these first
include("femshop_constants.jl");
include("femshop_config.jl");
include("femshop_prob.jl");

# include these next
include("logging.jl");

include("refel.jl");
include("mesh_data.jl");
include("mesh_read.jl");
include("mesh_write.jl")
include("grid.jl");
include("simple_mesh.jl");
include("recursive_ordering.jl");
include("tiled_ordering.jl");
include("ef_ordering.jl");
include("bad_ordering.jl");

include("general_utils.jl");
include("function_utils.jl");
include("symexpression.jl");
include("variables.jl");
include("coefficient.jl");
include("parameter.jl");
include("indexer.jl");
include("geometric_factors.jl");
include("dg_utils.jl");
include("time_steppers.jl");
include("fv_utils.jl");
include("grid_parent_child.jl");
include("polyharmonic_interp.jl");

include("output_data.jl");
include("cachsim_solve.jl");

# Femshop submodules
include("cachesim/CachesimOut.jl");
using .CachesimOut
include("SymbolicParser.jl")
using .SymbolicParser
include("CodeGenerator.jl");
using .CodeGenerator
include("DGSolver.jl");
using .DGSolver
include("CGSolver.jl");
using .CGSolver
include("FVSolver.jl");
using .FVSolver


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

"""
    umatrix(weights::EsomWeights{Float64})

Generate a U-matrix for the given ESOM weights.
"""
function umatrix(weights::EsomWeights{Float64}, settings::Settings = defaultSettings)
    (_, k, m) = size(weights)
    result = Array{Float64, 2}(undef, k, m)
    for index in CartesianIndices((k,m))
        weight = weights[:,index]
        neighbours = directNeighbours(index, settings)
        dist(w) = settings.distance(w, weight)
        result[index] = mean(map(dist, Slices(weights[:,neighbours], 1)))
    end
    result
end

"""
    pmatrix(data::AbstractMatrix{Float64}, weights::EsomWeights{Float64}; radius = nothing)

Generate a P-matrix for the given data and ESOM weights.

If no pareto radius is given, a suitable one is estimated.
"""
function pmatrix(data::AbstractMatrix{Float64}, weights::EsomWeights{Float64},
                 settings::Settings = defaultSettings;
                 radius::Union{Real, Nothing} = nothing)
    @assert size(data, 2) == size(weights, 1)

    if radius == nothing
        # described in 'Visualization and 3D Printing of Multivariate Data of
        # Biomarkers', Thrun et al., 2016a
        # produces rather flat maps, may chose a different algo?
        distances = begin
            d = pairwise(settings.distance, data, dims=1)
            (x,y) = size(d)
            [d[i,j] for j in 1:x for i in (j+1):y]
        end
        p20 = quantile(distances, 20/100)
        abc = ABCanalysis(distances)
        v = maximum(distances[abc.c_indices]) / minimum(distances[abc.a_indices])
        radius = v * p20
    end

    (_, k, m) = size(weights)
    result = fill(0, (k, m))
    for index in CartesianIndices((k,m))
        weight = weights[:,index]
        for row in eachrow(data)
            if settings.distance(row, weight) <= radius
                result[index] += 1
            end
        end
    end
    return result
end

"""
    ustarmatrix(um::AbstractMatrix{Float64}, pm::AbstractMatrix{Int})

Generate a U*-matrix from the given U- and P-matrices.
"""
function ustarmatrix(um::AbstractMatrix{Float64}, pm::AbstractMatrix{Int},
                     settings::Settings = defaultSettings)
    @assert size(um) == size(pm)
    # Using the calculation form 'U*-matrix: a tool to visualize clusters in
    # high dimensional data' (Ultsch 2003)
    scaleFactor(p::Int) = (p - mean(pm)) / (mean(pm) - maximum(pm)) + 1
    return um .* scaleFactor.(pm)
end

"""
    ustarmatrix(data::AbstractMatrix{Float64}, weights::EsomWeights{Float64})

Generate a U*-matrix for the given data and ESOM weights.
"""
function ustarmatrix(data::AbstractMatrix{Float64}, weights::EsomWeights{Float64},
                     settings::Settings = defaultSettings)
    um = umatrix(weights, settings)
    pm = pmatrix(data, weights, settings)
    return ustarmatrix(um, pm, settings)
end

# compatability with DataIO.LRNData
for f in (:pmatrix, :ustarmatrix)
    @eval @inline ($f)(data::LRNData, args...; kwargs...) = ($f)(data.data, args...; kwargs...)
end


text:
include("model.jl")
include("sampler.jl")
include("simulations.jl")
include("ABC_Gillespie.jl")
include("distances.jl")

using DataFrames
using Gadfly
using Compose
using Colors
using JLD
using MultivariateStats
using RCall
using EntropyEstimators

""" save results of ABC rejection or ABC SMC (each population is saved in a separated directory)

"""
function save_results(results::Array{ABC_population,1},model::ABC_Model,path::ASCIIString)

	#save the priors into a txt and jld file
    writedlm(string(path,"/priors.txt"),model.priors)
    jldopen(string(path,"/priors.jld"),"w") do file
        write(file,"priors",model.priors)
    end

    #save the particles from each population into a file
    for i in 1:length(results)

        #save_pop(i,results[i],model,path)
        compare_prior_posterior(path,pop_num=i)

    end

end

"""
save the results of one population
"""
function save_pop(pop_num::Int64,pop::ABC_population,model::ABC_Model,path::ASCIIString)

    #create population directory (only if non existing)
    pop_dir = string(path,"population",pop_num)
    run(`mkdir -p $pop_dir`)

    #save the variables of the current population into a .jld file that can be later loaded
    jldopen(string(pop_dir,"/pop$pop_num.jld"),"w") do file
       write(file,"distances",pop.distances)
       write(file,"weights",pop.weights)
       write(file,"particles",pop.particles)
       write(file,"epsilon",pop.epsilon)
       write(file,"num_steps",pop.num_steps)
    end

    #open new file for writing
    pop_file = open(string(pop_dir,"/particles.txt"), "w")

    particles = pop.particles

    for j in 1:length(particles)
        write(pop_file,string(particles[j]), "\n")
    end

    close(pop_file)

    #save weights and particles
    weights_file = string(pop_dir,"/weights.txt")
    dis_file = string(pop_dir,"/distances.txt")
    writedlm(weights_file,pop.weights)
    writedlm(dis_file,pop.distances)

    df = DataFrame()
    all_plots = Array{Gadfly.Plot,1}(model.num_params*model.num_params)
    #for models with many parameters, create a plot with the histograms only
    hist_plots = Array{Gadfly.Plot,1}(model.num_params)
    df_cols = collect(1:model.num_params)

    hist_plot_num = 1
    scatter_plot_num = 1

    cols = round(Int,model.num_params/2)
    if (model.num_params%2)==1
        cols+=1
    end

    rows = 2

    if model.num_params==3
    	rows=1
    end

    for j in 1:model.num_params

        #save all the plots of the current parameter into one plot
        par_plots = Array{Gadfly.Plot,1}(model.num_params)

        #get all the current parameter values
        p_vals = [particles[k][j] for k in 1:length(particles)]

        #add them to the dataframe
        insert!(df,j,p_vals,symbol("parameter$j"))

        #generate a histogram plot
        plot_j = plot(x=p_vals,Geom.histogram(bincount=10),
                        Guide.xlabel(string("Parameter ",j)),
                        Guide.ylabel("Frequency"),
                        Guide.title("Parameters Histogram")
                        )
        all_plots[hist_plot_num] = plot_j
        hist_plots[j] = plot_j
        par_plots[j] = plot_j
        hist_plot_num += (model.num_params+1)

        #generate scatter plots with all other parameters
        for n in 1:model.num_params
            if n!=j

                #only create scatter plots with other parameters
                n_vals = [particles[k][n] for k in 1:length(particles)]
                scatter_j_n = plot(x=p_vals,y=n_vals,Geom.point,
                                   Guide.xlabel(string("Parameter ",j)),
                                   Guide.ylabel(string("Parameter ",n)))
                all_plots[scatter_plot_num] = scatter_j_n
                par_plots[n] = scatter_j_n
                scatter_plot_num +=1
            else
                scatter_plot_num +=1
            end
        end

        #save the plots of the current parameter into one plot
        cs = transpose(reshape([Context[render(par_plots[i]) for i in 1:length(par_plots)]], cols,rows))
        p = gridstack(cs)
        draw(PDF(string(pop_dir,"/parameter_",j,"_plots.pdf"),18inch,10inch),p)

    end

    #save the plots into one plot
    cs = transpose(reshape([Context[render(all_plots[i]) for i in 1:length(all_plots)]], model.num_params,model.num_params))
    p = gridstack(cs)
    draw(PDF(string(pop_dir,"/parameters_histograms.pdf"),12inch,12inch),p)

    #save only histogram plots
    cs = transpose(reshape([Context[render(hist_plots[i]) for i in 1:length(hist_plots)]], cols,rows))
    h = gridstack(cs)
    draw(PDF(string(pop_dir,"/parameters_histograms_only.pdf"),18inch,10inch),h)

    #save the dataframe into a file
    writetable(string(pop_dir,"/particles_df.csv"),df)

    #plot the weighted posteriror distributions
    plot_weighted_posterior(pop_dir)
    #TODO: change so that this is done during the loop, create a function for making the plots and call it twice

    #PCA analysis
    #robustness_analysis(pop_dir)
end

"""
plot the data such that each variable is shown over time
"""
function plot_data(data::Array{Array{Float64,1},1},time_points::FloatRange{Float64},num_params::Int64,path::ASCIIString)

    df = DataFrame()
    df_cols = collect(1:num_params)
    p_layers = Layer[]
    p_colors = distinguishable_colors(num_params)
    times = collect(time_points)

    """p = plot(Guide.XLabel("Time"),
                 Guide.YLabel("Value"),
                 Guide.Title("Data"))
    """
    for j in 1:num_params
        #get all the current parameter values
        vals = [data[k][j] for k in 1:length(data)]

        #add them to the dataframe
        insert!(df,j,vals,symbol("var$j"))

        #create the layer
        l = layer(x=times,y=vals,Geom.point,Theme(default_point_size=8px,default_color=color(p_colors[j])))
        push!(p_layers,l[1])
        #Gadfly.add_plot_element!(p,l)
    end

    p = plot(p_layers,Guide.XLabel("Time"),
                     Guide.YLabel("Value"),
                     Guide.Title("Data"),
                     Guide.manual_color_key("Variable", [string("var",i) for i in 1:num_params], p_colors))

    draw(PDF(string(path,"data.pdf"),12inch,12inch),p)

end

"""
plot random walks
"""
function plot_random_walks(walks::Array{Random_Walk,1},path::ASCIIString;walks_title="Random Walks")

    p_layers = Array{Layer,1}(length(walks))
    walk_color = distinguishable_colors(length(walks))

    for i in 1:length(walks)

        w = walks[i]
        points= w.steps_points

        point_x = [points[i][1] for i in 1:length(points)]
        point_y = [points[i][2] for i in 1:length(points)]

        #l = layer(x=point_x,y=point_y,Geom.line)
        l = layer(x=point_x,y=point_y,Geom.path(),Theme(default_color=walk_color[i]))
        p_layers[i] = l[1]
    end

    p = plot(p_layers,Guide.XLabel("X"),
                     Guide.YLabel("Y"),
                     Guide.Title(walks_title),
                     Coord.Cartesian(xmin=-5, xmax=5, ymin=-10, ymax=5)
                     )

                    #Coord.Cartesian(xmin=-5, xmax=5, ymin=-20, ymax=0)

    draw(PDF(string(path,"random_walks.pdf"),12inch,12inch),p)

    return p

end

"""plotting function for transition matrix"""
function my_spy(M::AbstractMatrix, elements::Gadfly.ElementOrFunction...; mapping...)
    is, js, values = Gadfly._findnz(x->!isnan(x), M)
    n,m = size(M)
    df = DataFrames.DataFrame(i=is, j=js, value=values)
    plot(df, x=:j, y=:i, color=:value,
        Coord.cartesian(fixed=true, xmin=0.5, xmax=m+.5, ymin=0.5, ymax=n+.5),
        #Scale.color_continuous(minvalue=0.0, maxvalue=0.015),
        Scale.color_continuous,
        Geom.rectbin,
        Scale.x_continuous,
        Scale.y_continuous,
        elements...; mapping...)
end

"""
plot transition matrixes on the same figure
"""
function plot_TMs(TM_arr::Array{Array{Float64,2},1},path::ASCIIString,particle_num::Int64,pop_num::Int64;num_groups::Int64=4)

    #create population directory (only if non existing)
    pop_dir = string(path,"population",pop_num)
    run(`mkdir -p $pop_dir`)

    TM_plots = Array{Gadfly.Plot,1}(length(TM_arr))

    for i in 1:length(TM_arr)

        #plot the transition matrix
        T_plot = my_spy(TM_arr[i],Guide.xlabel(""),Guide.ylabel(""))
        TM_plots[i] = T_plot

    end

    #plot all the TMs in one plot
    cs = transpose(reshape([Context[render(TM_plots[i]) for i in 1:length(TM_plots)]], num_groups,num_groups))
    p = gridstack(cs)
    title = compose(context(0, 0, 1w, 0.25inch),text(0.5, 1.0, "Time after Injury / Distance from Wound", hcenter, vbottom))
    draw(PDF(string(pop_dir,"/TMs_particle",particle_num,".pdf"),12inch,12inch),vstack(title, p))

end

"""
robustness analysis - PCA on correlation matrix
"""
function robustness_analysis(path::ASCIIString)

    # particles = string(path,"particles_df.csv")

    # weights = string(path,"weights.txt")

    # #open the DF of the particles
    # post_dist_df = readtable(particles)

    # #open the DF of the weights
    # weights_df = readtable(weights,header=false)
    # weights_vec = Array(weights_df[:,1])

    # #get the wieghted posterior
    # post_dist_weighted = weighted_posterior(post_dist_df,weights_vec)

    # """
    # #get the correlation matrix from R - using the weights
    # cor_mat = rcopy(R"cov.wt($post_dist_df,wt=$weights_vec,cor=TRUE)$cor")

    # #run PCA using R
    # #cor_mat_pca_sdev = rcopy(R"princomp($cor_mat,cor=T)$sdev")
    # cor_mat_pca_sdev = rcopy(R"prcomp($cor_mat,cor=T)$sdev")

    # #cor_mat_pca_loadings = rcopy(R"princomp($cor_mat,cor=T)$loadings")
    # cor_mat_pca_loadings = rcopy(R"prcomp($cor_mat,cor=T)$rotation")

    # PCA_vars = cor_mat_pca_sdev

    # PCA_comps = cor_mat_pca_loadings
    # """

    #get the weighted posterior
    post_dist_weighted = readtable(string(path,"/weighted_particles_df.csv"))

    #get the correlation matrix
    post_dist_cor = cor(DataMatrix(post_dist_weighted))

    #run PCA on the correlation matrix
    post_dist_PCA = fit(PCA,Array(post_dist_cor))

    #get the variances of the components
    PCA_vars = principalvars(post_dist_PCA)

    #get the components
    PCA_comps = projection(post_dist_PCA)

    #plot the variances for the first 5 components
    comps = ["PC1","PC2","PC3","PC4","PC5"]
    p = plot(x=comps,y=PCA_vars[1:5],Geom.bar,
                     Guide.XLabel("Principal Components"),
                     Guide.YLabel("Variance"),
                     Guide.Title("Correlation Matrix PCA")
                     )

    draw(PDF(string(path,"/PCA_variance.pdf"),12inch,12inch),p)

    #parameters
    pars = ["w","mb","b0","mp","p0","A","D","tau0","Kd","R0"]

    #plot the first component for all parameters
    pc1=PCA_comps[:,1]
    p = plot(x=pars,y=pc1,Geom.bar,
                     Guide.XLabel(""),
                     Guide.YLabel("PC1"),
                     Guide.Title("Sloppy")
                     )

    draw(PDF(string(path,"/PCA_PC1.pdf"),12inch,12inch),p)

    #plot the fifth component for all parameters
    pc5=PCA_comps[:,5]
    p = plot(x=pars,y=pc5,Geom.bar,
                     Guide.XLabel(""),
                     Guide.YLabel("PC5"),
                     Guide.Title("Stiff")
                     )

    draw(PDF(string(path,"/PCA_PC5.pdf"),12inch,12inch),p)

    """
    #get the projection
    amount = zeros(size(PCA_comps))
    for i in 1:size(PCA_comps)[1]
        #get the sum of the ith row
        #PC_sum = sum(PCA_comps[i,:].^2)
        for j in 1:size(PCA_comps)[2]
            #get the sum of the jth col
            PC_sum = sum(PCA_comps[:,j].^2)
            amount[i,j] = (PCA_comps[i,j]^2) / PC_sum
        end
    end

    #plot the first and fifth projections

    #plot the first component for all parameters
    pc1=amount[:,1]
    p = plot(x=pars,y=pc1,Geom.bar,
                     Guide.XLabel(""),
                     Guide.YLabel("PC1"),
                     Guide.Title("Sloppy")
                     )

    draw(PDF(string(path,"PCA_PC1_projection.pdf"),12inch,12inch),p)

    #plot the fifth component for all parameters
    pc5=amount[:,5]
    p = plot(x=pars,y=pc5,Geom.bar,
                     Guide.XLabel(""),
                     Guide.YLabel("PC5"),
                     Guide.Title("Stiff")
                     )

    draw(PDF(string(path,"PCA_PC5_projection.pdf"),12inch,12inch),p)
    """

end

"""
sample from posterior population accodring to weights
"""
function weighted_posterior(particles::DataFrames.DataFrame,weights::Array{Float64,1})

        sampled_particles = DataFrame()

        for j in 1:size(particles)[2]
            #sample each of the parameters vector
            par_vec = particles[:,j]
            sampled = StatsBase.sample(par_vec,WeightVec(weights),size(particles)[1])
            insert!(sampled_particles,j,sampled,symbol("parameter$j"))
        end

        return sampled_particles
end

"""
plot weighted posterior distributions
"""
function plot_weighted_posterior(path::ASCIIString)

    #get the particles and weights
    particles = readtable(string(path,"/particles_df.csv"))
    weights_df = readtable(string(path,"/weights.txt"),header=false)
    weights= Array(weights_df[:,1])

    #get the weighted posterior
    weighted_particles = weighted_posterior(particles,weights)

    #save the weighted particles
    writetable(string(path,"/weighted_particles_df.csv"),weighted_particles)

    num_params = size(particles)[2]

    all_plots = Array{Gadfly.Plot,1}(num_params*num_params)
    #for models with many parameters, create a plot with the histograms only
    hist_plots = Array{Gadfly.Plot,1}(num_params)

    hist_plot_num = 1
    scatter_plot_num = 1

    cols = round(Int,num_params/2)
    if (num_params%2)==1
        cols+=1
    end

    rows=2
    if num_params==3
    	rows=1
    end

    for j in 1:num_params

        #save all the plots of the current parameter into one plot
        par_plots = Array{Gadfly.Plot,1}(num_params)

        #get all the current parameter values
        p_vals = weighted_particles[:,j]

        #generate a histogram plot
        plot_j = plot(x=p_vals,Geom.histogram(bincount=10),
                        Guide.xlabel(string("Parameter ",j)),
                        Guide.ylabel("Frequency"),
                        Guide.title("Weighted Parameters Histogram")
                        )
        all_plots[hist_plot_num] = plot_j
        hist_plots[j] = plot_j
        par_plots[j] = plot_j
        hist_plot_num += (num_params+1)

        #generate scatter plots with all other parameters
        for n in 1:num_params
            if n!=j

                #only create scatter plots with other parameters
                n_vals = weighted_particles[:,n]
                scatter_j_n = plot(x=p_vals,y=n_vals,Geom.point,
                                   Guide.xlabel(string("Parameter ",j)),
                                   Guide.ylabel(string("Parameter ",n)))
                all_plots[scatter_plot_num] = scatter_j_n
                par_plots[n] = scatter_j_n
                scatter_plot_num +=1
            else
                scatter_plot_num +=1
            end
        end

        #save the plots of the current parameter into one plot
        cs = transpose(reshape([Context[render(par_plots[i]) for i in 1:length(par_plots)]], cols,rows))
        p = gridstack(cs)
        draw(PDF(string(path,"/weighted_parameter_",j,"_plots.pdf"),18inch,10inch),p)

    end

    #save the plots into one plot
    cs = transpose(reshape([Context[render(all_plots[i]) for i in 1:length(all_plots)]], num_params,num_params))
    p = gridstack(cs)
    draw(PDF(string(path,"/weighted_parameters_histograms.pdf"),12inch,12inch),p)

    #save only histogram plots
    cs = transpose(reshape([Context[render(hist_plots[i]) for i in 1:length(hist_plots)]], cols,rows))
    h = gridstack(cs)
    draw(PDF(string(path,"/weighted_parameters_histograms_only.pdf"),18inch,10inch),h)

end

"""
compare prior and posterior distribution using Entropy
"""
function compare_prior_posterior(path::ASCIIString;pop_num=5,estimator="shrinkage")

    #get the prior distributions
    priors = jldopen(string(path,"/priors.jld"),"r") do file
        read(file, "priors")
    end

    pop_dir = string(path,"/population",pop_num,"/")

    #get the particles
    particles = readtable(string(pop_dir,"/particles_df.csv"))

    #get the weighted particles
    weighted_particles = readtable(string(pop_dir,"/weighted_particles_df.csv"))

    #TODO: check what happends with weighted particles

    entropy_arr = zeros(size(particles)[2])

    #go over each parameter
    for i in 1:size(particles)[2]
        #create the prior distribution
        prior_dist = rand(priors[i],size(particles)[1])

        #get the posterior distribution
        post_dist = particles[:,i]

        #get the entropy using the requsted estimate
        #prior_post_entropy = get_entropy(prior_dist, post_dist)
        prior_ent = get_entropy(prior_dist)
        post_ent = get_entropy(post_dist)
        entropy_arr[i] = prior_ent - post_ent

    end

    #plot the entropy values
    pars = ["w","mb","b0","mp","p0","A","D","tau0","Kd","R0"]
    p = plot(x=pars,y=entropy_arr,Geom.bar,
                     Guide.XLabel(""),
                     Guide.YLabel("Entropy"),
                     Guide.Title("Prior Posterior Entropy")
                     )

    draw(PDF(string(pop_dir,"/entropy.pdf"),12inch,12inch),p)

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct TargetResource <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id
    resourceName::Any # spec type: Union{ Nothing, String } # spec name: resourceName
    resourceType::Any # spec type: Union{ Nothing, String } # spec name: resourceType

    function TargetResource(;id=nothing, resourceName=nothing, resourceType=nothing)
        o = new()
        validate_property(TargetResource, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        validate_property(TargetResource, Symbol("resourceName"), resourceName)
        setfield!(o, Symbol("resourceName"), resourceName)
        validate_property(TargetResource, Symbol("resourceType"), resourceType)
        setfield!(o, Symbol("resourceType"), resourceType)
        o
    end
end # type TargetResource

const _property_map_TargetResource = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"), Symbol("resourceName")=>Symbol("resourceName"), Symbol("resourceType")=>Symbol("resourceType"))
const _property_types_TargetResource = Dict{Symbol,String}(Symbol("id")=>"String", Symbol("resourceName")=>"String", Symbol("resourceType")=>"String")
Base.propertynames(::Type{ TargetResource }) = collect(keys(_property_map_TargetResource))
Swagger.property_type(::Type{ TargetResource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_TargetResource[name]))}
Swagger.field_name(::Type{ TargetResource }, property_name::Symbol) =  _property_map_TargetResource[property_name]

function check_required(o::TargetResource)
    true
end

function validate_property(::Type{ TargetResource }, name::Symbol, val)
end


text:
import Base.ndims, Base.*, Base.==, Base.convert, Base.getindex

export
	AbstractInterpolation,
	SimpleInterpolation,
	CompoundInterpolation,
	InterpolationSupport

abstract AbstractIndexManipulator

abstract AbstractInterpolation{N} <: AbstractIndexManipulator
AbstractInterpolation(N=1) = AbstractInterpolation{N}
ndims{N}(::Type{AbstractInterpolation{N}}) = N
ndims{I<:AbstractInterpolation}(::Type{I}) = ndims(I.super)

abstract SimpleInterpolation <: AbstractInterpolation{1}

abstract CompoundInterpolation{N, INS<:Tuple{Vararg{AbstractInterpolation}}} <: AbstractInterpolation{N}
CompoundInterpolation{INS<:Tuple{Vararg{AbstractInterpolation}}}(::Type{INS}) = CompoundInterpolation{sum(map(ndims, INS.parameters)), INS}
CompoundInterpolation(ins...) = CompoundInterpolation(Tuple{ins...})
ndims{INC<:CompoundInterpolation}(::Type{INC}) = INC.parameters[1]

using Base.Cartesian
@inline @generated function getindex{N, INS}(interps::Type{CompoundInterpolation{N, INS}}, x::Vararg{TypeVar(:T), N})
	bsym = [Symbol("b_", i) for i in 1:N]
	ibsym = [Symbol("ib_", i) for i in 1:N]

	bsymdefs = Expr(:block, [Expr(:(=), bsym[i], :(getindex($(INS.parameters[i]), x[$i]))) for i in 1:N]...)

	flesh = Expr(:tuple, Expr(:call, :*, [:($i[1]) for i in ibsym]...), [:($i[2]) for i in ibsym]...)
	flesh = Expr(:generator, flesh, [Expr(:(=), ibsym[i], bsym[i]) for i in 1:N]...)

	quote
		$bsymdefs
		return vec($(Expr(:comprehension, flesh)))
	end
end

@inline @generated function getindex{AT, N, INS, IT<:Number, IM<:Tuple}(arr::AbstractArray{AT,N}, ::Type{IM}, ::Type{CompoundInterpolation{N, INS}}, x::Vararg{IT, N})
	symbolic = [getindex_symbolic(INS.parameters[i], :(x[$i])) for i in 1:N]

	ibsym = [[Symbol("ib_", i, "_", j) for j in 1:length(symbolic[i][2])] for i in 1:N]
	cbsym = [[Symbol("cb_", i, "_", j) for j in 1:length(symbolic[i][2])] for i in 1:N]

	setups = Expr(:block, [symbolic[i][1] for i in 1:N]...)
	bsymdefs = Expr(:block, [Expr(:block, [Expr(:block, Expr(:(=), cbsym[i][j], symbolic[i][2][j][1]), Expr(:(=), ibsym[i][j], symbolic[i][2][j][2])) for j in 1:length(symbolic[i][2])]...) for i in 1:N]...)

	indices_prod = Base.product(ibsym...)
	indices = map(a -> Expr(:ref, :arr, IM.parameters..., a...), indices_prod)

	for i in 1:N
		indices = mapslices(a->Expr(:call, :+, [Expr(:call, :*, ic, ia) for (ic,ia) in zip(cbsym[i], a)]...), indices, i)
	end

	quote
		$(Expr(:meta, :inline))
		$setups
		$bsymdefs
		return $(indices[1])
	end
end

@inline getindex{AT, N, IT<:Number, IN<:SimpleInterpolation, IM<:Tuple}(arr::AbstractArray{AT, N}, ::Type{IM}, ::Type{IN}, x::Vararg{IT, N}) = getindex(arr, IM, CompoundInterpolation{N, NTuple{N, IN}}, x...)
@inline getindex{AT, N, IT<:Number, IN<:AbstractInterpolation}(arr::AbstractArray{AT,N}, ::Type{IN}, x::Vararg{IT, N}) = getindex(arr, Tuple{}, IN, x...)

@inline @generated function inplaceadd!{AT, N, INS, VT, IT<:Number, IM<:Tuple}(arr::AbstractArray{AT,N}, ::Type{IM}, ::Type{CompoundInterpolation{N, INS}},  v::VT, x::Vararg{IT, N})
	symbolic = [getindex_symbolic(INS.parameters[i], :(x[$i])) for i in 1:N]

	ibsym = [[Symbol("ib_", i, "_", j) for j in 1:length(symbolic[i][2])] for i in 1:N]
	cbsym = [[Symbol("cb_", i, "_", j) for j in 1:length(symbolic[i][2])] for i in 1:N]

	setups = Expr(:block, [symbolic[i][1] for i in 1:N]...)
	bsymdefs = Expr(:block, [Expr(:block, [Expr(:block, Expr(:(=), cbsym[i][j], symbolic[i][2][j][1]), Expr(:(=), ibsym[i][j], symbolic[i][2][j][2])) for j in 1:length(symbolic[i][2])]...) for i in 1:N]...)

	coeff_prod = Base.product(cbsym...)
	indices_prod = Base.product(ibsym...)

	flesh = Expr(:block, map((c,i) -> Expr(:call, :inplaceadd!, :arr, IM.parameters..., Expr(:call, :*, c..., :v), i...), coeff_prod, indices_prod)...)

	quote
		$(Expr(:meta, :inline))
		$setups
		$bsymdefs
		$flesh
		return
	end
end

@inline inplaceadd!{AT, N, IN<:SimpleInterpolation, VT, IM<:Tuple}(arr::AbstractArray{AT,N}, ::Type{IM}, ::Type{IN},  v::VT, x::Vararg{TypeVar(:T), N})  = inplaceadd!(arr, IM, CompoundInterpolation{N, NTuple{N, IN}}, v, x...)
@inline inplaceadd!{AT, N, IN<:AbstractInterpolation, VT}(arr::AbstractArray{AT,N}, ::Type{IN},  v::VT, x::Vararg{TypeVar(:T), N})  = inplaceadd!(arr, Tuple{}, IN, v, x...)

max_stars_repo_name: JKrehl/Interpol

max_issues_repo_name: JKrehl/Interpol

max_forks_repo_name: JKrehl/Interpol


text:
# mri/z-test.jl

using Test: @test

@test exp_mult(:test)
@test exp_xform(:test)
@test image_geom_mri(:test)
@test ir_mri_coil_compress(:test)
@test ir_mri_kspace_ga_radial(:test)
@test ir_mri_sensemap_sim(:test)
@test mri_objects(:test)


text:
using Random

@testset "Get clusters" begin
    cluster_name = ENV["BANYAN_CLUSTER_NAME"]

    clusters = get_clusters()
    get_cluster_s3_bucket_name(cluster_name)
    running_clusters = get_running_clusters()

    @test haskey(clusters, cluster_name)
    @test all(c -> c[2].status == :running, running_clusters)

end

@testset "Update clusters" begin
    cluster_name = ENV["BANYAN_CLUSTER_NAME"]

    update_cluster(cluster_name)
    cluster_status = get_cluster_status(cluster_name)

    @test cluster_status == :updating

    while cluster_status == :updating
        sleep(5)
        cluster_status = get_cluster_status(cluster_name)
    end
    @test cluster_status == :running
end

function bucket_exists(s3_bucket_name)
    ispath(S3Path("s3://$(s3_bucket_name)", config=Banyan.get_aws_config()))
end

@testset "Create clusters" begin
    Random.seed!()
    cluster_name = "cluster-$(Random.randstring(['a':'z'; '0':'9'], 6))"
    create_cluster(
        name=cluster_name,
        instance_type="t3.xlarge",  # 4 vCPUs
        max_num_workers=32,
        initial_num_workers=8,
        min_num_workers=2,
        scaledown_time=30,
    )
    @test c.status == :running
    delete_cluster(cluster_name)
end

@testset "Destroy and delete clusters with $s3_bucket S3 bucket" for s3_bucket in [
        "default", "user-provided"
    ]
    Random.seed!()
    cluster_name = "cluster-$(Random.randstring(['a':'z'; '0':'9'], 6))"
    if s3_bucket == "default"
        s3_bucket = nothing
    elseif s3_bucket == "user-provided"
        s3_bucket = Random.randstring(['a':'z'; '0':'9'], 6)
        s3_create_bucket(Banyan.get_aws_config(), s3_bucket)
    end

    # Create a cluster (at least initiate) and check that S3 bucket exists
    c = create_cluster(
        name=cluster_name,
        instance_type="t3.large",
        s3_bucket_name=s3_bucket,
        nowait=true
    )
    sleep(30) # Just to ensure that cluster creation has initiated
    s3_bucket_name = get_cluster_s3_bucket_name(cluster_name)
    s3_bucket_exists = bucket_exists(s3_bucket_name)
    if !isnothing(s3_bucket)
        @test s3_bucket == s3_bucket_name
    end
    @test s3_bucket_exists

    # Destroy cluster and check that S3 bucket still exists
    destroy_cluster(cluster_name)
    s3_bucket_exists = bucket_exists(s3_bucket_name)
    @test s3_bucket_exists
    sleep(30) # Just to ensure that cluster destruction is complete

    # Re-create cluster and check that S3 bucket exists and is same as before
    while get_cluster_status(cluster_name) != :terminated
        sleep(15)
    end
    c_r = create_cluster(
        name=cluster_name,
        nowait=true
    )
    s3_bucket_name_r = get_cluster_s3_bucket_name(cluster_name)
    s3_bucket_exists = bucket_exists(s3_bucket_name_r)
    @test s3_bucket_exists
    @test s3_bucket_name == s3_bucket_name_r

    # Delete cluster
    delete_cluster(cluster_name)
    sleep(30)  # Just to ensure that bucket has been deleted
    s3_bucket_exists = bucket_exists(s3_bucket_name_r)
    @test !s3_bucket_exists

    # Check that the cluster cannot be created again
    @test_throws ErrorException create_cluster(name=cluster_name, nowait=true)
end

@testset "Benchmark create_cluster with $instance_type instance type" for instance_type in [
    "t3.xlarge", "t3.2xlarge", "c5.2xlarge", "m4.4xlarge", "m4.10xlarge"
]
    Random.seed!()
    cluster_name = "cluster-$(Random.randstring(['a':'z'; '0':'9'], 6))"
    t = @elapsed begin
        c = create_cluster(
            name=cluster_name,
            instance_type=instance_type,
            max_num_workers=16,
            initial_num_workers=1
        )
    end
    delete_cluster(cluster_name)

    # Save results to file
    open("create_cluster_times.txt", "a") do f
        write(f, "$(instance_type)\t$(string(t/60))\n")
    end

    # Verify that cluster was spun up
    @test c.status == :running
end

@testset "Upload file to S3 from $src_type" for src_type in [
    "local",
    "http",
    "s3"
]
    if src_type == "local"
        src_path = "file://data/iris.csv"
        dst_name = "iris_from_local.csv"
    elseif src_type == "http"
        src_path = "https://raw.githubusercontent.com/banyan-team/banyan-julia/v0.1.3/BanyanDataFrames/test/res/iris.csv"
        dst_name = "iris_from_http.csv"
    elseif src_type == "s3"
        s3_bucket = Random.randstring(['a':'z'; '0':'9'], 6)
        dst_name = "data_from_s3"
        src_path = "s3://$s3_bucket/$dst_name"
        # Create a bucket and upload data
        s3_create_bucket(Banyan.get_aws_config(), s3_bucket)
        s3_put(Banyan.get_aws_config(), s3_bucket, dst_name, "some file contents")
    end

    cluster_name = ENV["BANYAN_CLUSTER_NAME"]
    upload_to_s3(src_path; dst_name=dst_name, cluster_name=cluster_name)

    cluster_s3_bucket = get_cluster_s3_bucket_name(cluster_name)
    @test ispath(S3Path("s3://$cluster_s3_bucket/$dst_name"))

    # Cleanup
    s3_delete(Banyan.get_aws_config(), cluster_s3_bucket, dst_name)
    if src_type == "s3"
        s3_delete(Banyan.get_aws_config(), s3_bucket, dst_name)
        s3_delete_bucket(Banyan.get_aws_config(), s3_bucket)
    end
end

@testset "Upload dir to S3 from $src_type" for src_type in [
    "local",
]
    if src_type == "local"
        src_path = "data"
    end

    cluster_name = ENV["BANYAN_CLUSTER_NAME"]
    dst_name = upload_to_s3(src_path; cluster_name=cluster_name)

    cluster_s3_bucket = get_cluster_s3_bucket_name(cluster_name)
    @test ispath(S3Path("s3://$cluster_s3_bucket/$dst_name"))
    for f_name in readdir(src_path)
        @test ispath(S3Path("s3://$cluster_s3_bucket/$dst_name/$f_name"))
    end

    # Cleanup
    for f_name in readdir(src_path)
        s3_delete(Banyan.get_aws_config(), cluster_s3_bucket, "$dst_name/$f_name")
    end
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
module NeuroImages

export NeuroImage, phys_ndims, phys_size, location, transform_to, push_transformation!, mni, talairach, voxel, scanner, coordsys, read_nifti

include("reference.jl")
include("transformation.jl")
include("image.jl")
include("io.jl")

end # module


text:
abstract type AbstractAssimilation end

abstract type AbstractAssimilationC <: AbstractAssimilation end
abstract type AbstractAssimilationCN <: AbstractAssimilation end
abstract type AbstractAssimilationCNP <: AbstractAssimilation end

@columns struct AssimilationC{NSB,KD,YDE,NC,NX} <: AbstractAssimilationC
    N_SB::NSB   | mol/mol   | "Mol transporters per mol of biomass carbon"
    K_D::KD     | mol/m^3   | "Reference half-saturation constant"
    y_DE::YDE   | mol/mol   | "Yield on assimilation"
    N_C::NC     | _         | "Number of C atoms per monomer"
    N_X::NX     | _         | "Monomer id for return-on-investment"
end

for fn in fieldnames(AssimilationC)
    @eval $fn(p::AssimilationC) = p.$fn
end

assimilation!(J_DE::Vector{Float64}, p::AssimilationC, D::Vector{Float64}, V::Vector{Float64}) = begin
    n_substrates = size(D,1)
    n_consumers  = size(V,1)
    k2p = 180.0*60*60*ones(n_substrates)
    ECA  = ECA_kinetics!(zeros(n_substrates,n_consumers), D, V, K_D(p), k2p, N_SB(p))
    J_DE = vcat(sum((1.0 .- y_DE(p)).*N_C(p).*ECA[:, 1:n_consumers], dims=1)...)
end

uptake!(J_D::Vector{Float64}, p::AbstractAssimilation, D::Vector{Float64}, V::Vector{Float64}) = begin
    n_substrates = size(D,1)
    n_consumers  = size(V,1)
    k2p = 180.0*60*60*ones(n_substrates)
    ECA  = ECA_kinetics!(zeros(n_substrates,n_consumers), D, V, K_D(p), k2p, N_SB(p))
    J_D = vcat(sum(ECA[:, 1:n_consumers], dims=2)...).*N_C(p)
end

assimilation_production!(J_DE_CO2::Vector{Float64}, p::AbstractAssimilation, D::Vector{Float64}, V::Vector{Float64}) = begin
    n_substrates = size(D,1)
    n_consumers  = size(V,1)
    k2p = 180.0*60*60*ones(n_substrates)
    ECA  = ECA_kinetics!(zeros(n_substrates,n_consumers), D, V, K_D(p), k2p, N_SB(p))
    J_DE_CO2 = vcat(sum(y_DE(p).*N_C(p).*ECA[:, 1:n_consumers], dims=1)...)
end

max_stars_repo_name: giannamars/DEBMicroTrait

max_issues_repo_name: giannamars/DEBMicroTrait

max_forks_repo_name: giannamars/DEBMicroTrait


text:
using Documenter, DemoCards
using Augmentor
using Random
using MLDatasets

try
    using ISICArchive
catch
    using Pkg
    Pkg.add(url="https://github.com/Evizero/ISICArchive.jl.git", rev="master")
    using ISICArchive
end

ENV["DATADEPS_ALWAYS_ACCEPT"] = true # MLDatasets

op_templates, op_theme = cardtheme("grid")
operations, operations_cb = makedemos("operations", op_templates)

format = Documenter.HTML(edit_link = "master",
                         prettyurls = get(ENV, "CI", nothing) == "true",
                         assets = [
                             joinpath("assets", "favicon.ico"),
                             joinpath("assets", "style.css"),
                             op_theme
                        ]
)

About = "Introduction" => "index.md"

GettingStarted = "gettingstarted.md"

UserGuide = "User's guide" => [
        "interface.md",
        operations
    ]

DevGuide = "Developer's guide" => [
        "wrappers.md"
    ]

Examples = "Examples" => [
        "examples/flux.md"
    ]

License = "License" => "license.md"

PAGES = [
    About,
    GettingStarted,
    UserGuide,
    DevGuide,
    Examples,
    License
    ]

makedocs(
    modules = [Augmentor],
    sitename = "Augmentor.jl",
    authors = "Christof Stocker",
    format = format,
    checkdocs = :exports,
    pages = PAGES
)

operations_cb()

deploydocs(repo = "github.com/Evizero/Augmentor.jl.git")


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
function add_outage_results(m, p, r::Dict)
	# TODO with many outages the dispatch arrays are so large that it can take hours to create them
	# (eg. 8760 * 12 hour outages with PV, storage and diesel makes 7*12*8760 = 735,840 values)
	# For now the outage dispatch outputs are not created (commented out below). Perhaps make a new
	# function to optionally get the outage dispatch values so that we don't slow down returning the
	# other results.
	r["expected_outage_cost"] = value(m[:ExpectedOutageCost])
	r["max_outage_cost_per_outage_duration"] = value.(m[:dvMaxOutageCost]).data
	r["dvUnservedLoad"] = value.(m[:dvUnservedLoad]).data
	S = length(p.elecutil.scenarios)
	T = length(p.elecutil.outage_start_timesteps)
	unserved_load_per_outage = Array{Float64}(undef, S, T)
	for s in 1:S, t in 1:T
		unserved_load_per_outage[s, t] = sum(r["dvUnservedLoad"][s, t, ts] for
											 ts in 1:p.elecutil.outage_durations[s])
	# need the ts in 1:p.elecutil.outage_durations[s] b/c dvUnservedLoad has unused values in third dimension
	end
	r["unserved_load_per_outage"] = round.(unserved_load_per_outage, digits=2)
	r["mg_storage_upgrade_cost"] = value(m[:dvMGStorageUpgradeCost])
	r["dvMGDischargeFromStorage"] = value.(m[:dvMGDischargeFromStorage]).data

	for t in p.techs
		r[t * "_upgraded"] = value(m[:binMGTechUsed][t])
	end
	r["storage_upgraded"] = value(m[:binMGStorageUsed])

	if !isempty(p.pvtechs)
		for t in p.pvtechs

			# need the following logic b/c can have non-zero mg capacity when not using the capacity
			# due to the constraint for setting the mg capacities equal to the grid connected capacities
			if Bool(round(r[t * "_upgraded"], digits=1))
				r[string(t, "mg_kw")] = round(value(m[:dvMGsize][t]), digits=4)
			else
				r[string(t, "mg_kw")] = 0
			end
			r[string("mg_", t, "_upgrade_cost")] = round(value(m[:dvMGTechUpgradeCost][t]), digits=2)

			if !isempty(p.storage.types)
				PVtoBatt = (m[:dvMGProductionToStorage][t, s, tz, ts] for
					s in p.elecutil.scenarios,
					tz in p.elecutil.outage_start_timesteps,
					ts in p.elecutil.outage_timesteps)
			else
				PVtoBatt = []
			end
			r[string("mg", t, "toBatt")] = round.(value.(PVtoBatt), digits=3)

			PVtoCUR = (m[:dvMGCurtail][t, s, tz, ts] for
				s in p.elecutil.scenarios,
				tz in p.elecutil.outage_start_timesteps,
				ts in p.elecutil.outage_timesteps)
			r[string("mg", t, "toCurtail")] = round.(value.(PVtoCUR), digits=3)

			PVtoLoad = (
				m[:dvMGRatedProduction][t, s, tz, ts] * p.production_factor[t, tz+ts]
						* p.levelization_factor[t]
				- m[:dvMGCurtail][t, s, tz, ts]
				- m[:dvMGProductionToStorage][t, s, tz, ts] for
					s in p.elecutil.scenarios,
					tz in p.elecutil.outage_start_timesteps,
					ts in p.elecutil.outage_timesteps
			)
			r[string("mg", t, "toLoad")] = round.(value.(PVtoLoad), digits=3)
		end
	end

	if !isempty(p.gentechs)
		for t in p.gentechs

			# need the following logic b/c can have non-zero mg capacity when not using the capacity
			# due to the constraint for setting the mg capacities equal to the grid connected capacities
			if Bool(round(r[t * "_upgraded"], digits=1))
				r[string(t, "_mg_kw")] = round(value(m[:dvMGsize][t]), digits=4)
			else
				r[string(t, "mg_kw")] = 0
			end

			r[string("mg_", t, "_fuel_used")] = value.(m[:dvMGFuelUsed][t, :, :]).data
			r[string("mg_", t, "_upgrade_cost")] = round(value(m[:dvMGTechUpgradeCost][t]), digits=2)

			if !isempty(p.storage.types)
				GenToBatt = (m[:dvMGProductionToStorage][t, s, tz, ts] for
					s in p.elecutil.scenarios,
					tz in p.elecutil.outage_start_timesteps,
					ts in p.elecutil.outage_timesteps)
			else
				GenToBatt = []
			end
			r[string("mg", t, "toBatt")] = round.(value.(GenToBatt), digits=3)

			GENtoCUR = (m[:dvMGCurtail][t, s, tz, ts] for
				s in p.elecutil.scenarios,
				tz in p.elecutil.outage_start_timesteps,
				ts in p.elecutil.outage_timesteps)
			r[string("mg", t, "toCurtail")] = round.(value.(GENtoCUR), digits=3)

			GENtoLoad = (
				m[:dvMGRatedProduction][t, s, tz, ts] * p.production_factor[t, tz+ts]
						* p.levelization_factor[t]
				- m[:dvMGCurtail][t, s, tz, ts]
				- m[:dvMGProductionToStorage][t, s, tz, ts] for
					s in p.elecutil.scenarios,
					tz in p.elecutil.outage_start_timesteps,
					ts in p.elecutil.outage_timesteps
			)
			r[string("mg", t, "toLoad")] = round.(value.(GENtoLoad), digits=3)
		end
	end
end

max_stars_repo_name: brtietz/REoptLite

max_issues_repo_name: brtietz/REoptLite

max_forks_repo_name: brtietz/REoptLite