text:
# XXX prints dataframe to csv file
"""
```julia
printObject(print_df::DataFrame, model_object::anyModel)
```

Writes a DataFrame of parameters, constraints, or variables to a `.csv` file in readable format (strings instead of ids). See [Individual elements](@ref).
"""
function printObject(print_df::DataFrame,anyM::anyModel; fileName::String = "", rtnDf::Tuple{Vararg{Symbol,N} where N} = (:csv,), filterFunc::Function = x -> true)

	sets = anyM.sets
	options = anyM.options

	colNam_arr = namesSym(print_df)
    cntCol_int = size(colNam_arr,1)

	# filters values according to filter function,
	print_df = copy(filter(filterFunc,print_df))

	# converts variable column to value of variable
	if :var in colNam_arr
		print_df[!,:var] = value.(print_df[!,:var])
	end

    for i = 1:cntCol_int
        lookUp_sym = Symbol(split(String(colNam_arr[i]),"_")[1])
		if !(lookUp_sym in keys(sets)) && lookUp_sym == :eqn
			print_df[!,i] = string.(print_df[!,i])
        elseif lookUp_sym in keys(sets)
			print_df[!,i] = map(x -> createFullString(x,sets[lookUp_sym]),print_df[!,i])
        end
    end

	# rename columns
	colName_dic = Dict(:Ts_dis => :timestep_dispatch, :Ts_exp => :timestep_expansion, :Ts_expSup => :timestep_superordinate_expansion, :Ts_disSup => :timestep_superordinate_dispatch,
															:R => :region, :R_dis => :region_dispatch, :R_exp => :region_expansion, :R_to => :region_to, :R_from => :region_from, :C => :carrier, :Te => :technology,
																:cns => :constraint, :var => :variable)

	rename!(print_df,map(x -> x in keys(colName_dic) ? colName_dic[x] : x, namesSym(print_df)) )
	if :csv in rtnDf
    	CSV.write("$(options.outDir)/$(fileName)_$(options.outStamp).csv",  print_df)
	end

	if :csvDf in rtnDf return print_df end
end

# <editor-fold desc="report results to csv files"
"""
```julia
reportResults(reportType::Symbol, model_object::anyModel; rtnOpt::Tuple = (:csv,))
```

Writes results to `.csv` file with content depending on `reportType`. Available types are `:summary`, `:exchange`, and `:costs`. See [Analysed results](@ref).
"""
reportResults(reportType::Symbol,anyM::anyModel; kwargs...) = reportResults(Val{reportType}(),anyM::anyModel; kwargs...)

# XXX summary of all capacity and dispatch results
function reportResults(objGrp::Val{:summary},anyM::anyModel; wrtSgn::Bool = true, rtnOpt::Tuple{Vararg{Symbol,N} where N} = (:csv,))

    techSym_arr = collect(keys(anyM.parts.tech))
	allData_df = DataFrame(Ts_disSup = Int[], R_dis = Int[], Te = Int[], C = Int[], variable = Symbol[], value = Float64[])

	# XXX get demand values
	dem_df = copy(anyM.parts.bal.par[:dem].data)
	if !isempty(dem_df)
		dem_df[!,:lvlR] = map(x -> anyM.cInfo[x].rDis, :C in namesSym(dem_df) ? dem_df[!,:C] : filter(x -> x != 0,getfield.(values(anyM.sets[:C].nodes),:idx)))

		# aggregates demand values

		# artificially add dispatch dimensions, if none exist
		if :Ts_dis in namesSym(dem_df)
			ts_dic = Dict(x => anyM.sets[:Ts].nodes[x].lvl == anyM.supTs.lvl ? x : getAncestors(x,anyM.sets[:Ts],:int,anyM.supTs.lvl)[end] for x in unique(dem_df[!,:Ts_dis]))
			dem_df[!,:Ts_disSup] = map(x -> ts_dic[x],dem_df[!,:Ts_dis])
		else
			dem_df[!,:Ts_disSup] .= anyM.supTs.step
			dem_df = flatten(dem_df,:Ts_disSup)
		end

		dem_df[!,:val] = dem_df[!,:val]	.* getResize(dem_df,anyM.sets[:Ts],anyM.supTs) ./ anyM.options.redStep

		allR_arr = :R_dis in namesSym(dem_df) ? unique(dem_df[!,:R_dis]) : getfield.(getNodesLvl(anyM.sets[:R],1),:idx)
		allLvlR_arr = unique(dem_df[!,:lvlR])
		r_dic = Dict((x[1], x[2]) => (anyM.sets[:R].nodes[x[1]].lvl < x[2] ? getDescendants(x[1], anyM.sets[:R],false,x[2]) : getAncestors(x[1],anyM.sets[:R],:int,x[2])[end]) for x in Iterators.product(allR_arr,allLvlR_arr))
		if :R_dis in namesSym(dem_df)
			dem_df[!,:R_dis] = map(x -> r_dic[x.R_dis,x.lvlR],eachrow(dem_df[!,[:R_dis,:lvlR]]))
		else
			dem_df[!,:R_dis] .= 0
		end

		dem_df = combine(groupby(dem_df,[:Ts_disSup,:R_dis,:C]),:val => ( x -> sum(x) / 1000) => :value)
		dem_df[!,:Te] .= 0
		dem_df[!,:variable] .= :demand
		if wrtSgn dem_df[!,:value] = dem_df[!,:value] .* -1 end
		allData_df = vcat(allData_df,dem_df)
	end

	# XXX get expansion and capacity variables

	for t in techSym_arr
		part = anyM.parts.tech[t]
		tech_df = DataFrame(Ts_disSup = Int[], R_dis = Int[], Te = Int[], C = Int[], variable = Symbol[], value = Float64[])

		# get installed capacity values
		for va in intersect(keys(part.var),(:expConv, :expStIn, :expStOut, :expStSize, :expExc, :capaConv, :capaStIn, :capaStOut,  :capaStSize, :oprCapaConv, :oprCapaStIn, :oprCapaStOut, :oprCapaStSize))
			capa_df = copy(part.var[va])
			if va in (:expConv, :expStIn, :expStOut, :expStSize)
				capa_df = flatten(capa_df,:Ts_expSup)
				select!(capa_df,Not(:Ts_disSup))
				rename!(capa_df,:Ts_expSup => :Ts_disSup)
			end
			# set carrier column to zero for conversion capacities and add a spatial dispatch column
			if va in (:expConv,:capaConv,:oprCapaConv)
				capa_df[!,:C] .= 0
				capa_df[!,:R_dis] = map(x -> getAncestors(x,anyM.sets[:R],:int,part.balLvl.ref[2])[end],capa_df[!,:R_exp])
			else
				capa_df[!,:R_dis] = map(x -> getAncestors(x.R_exp,anyM.sets[:R],:int,anyM.cInfo[x.C].rDis)[end],eachrow(capa_df))
			end

			select!(capa_df,Not(:R_exp))
			# aggregate values and add to tech data frame
			capa_df = combine(groupby(capa_df,[:Ts_disSup,:R_dis,:C,:Te]),:var => ( x -> value.(sum(x))) => :value)
			capa_df[!,:variable] .= va
			tech_df = vcat(tech_df,capa_df)
		end

		# add tech dataframe to overall data frame
		allData_df = vcat(allData_df,tech_df)
	end

	# XXX get dispatch variables
	for va in (:use, :gen, :stIn, :stOut, :stExtIn, :stExtOut, :stIntIn, :stIntOut, :emission, :crt, :lss, :trdBuy, :trdSell)
		# get all variables, group them and get respective values
		allVar_df = getAllVariables(va,anyM)
		if isempty(allVar_df) continue end

		disp_df = combine(groupby(allVar_df,intersect(intCol(allVar_df),[:Ts_disSup,:R_dis,:C,:Te])),:var => (x -> value(sum(x))) => :value)
		# scales values to twh (except for emissions of course)
		if va != :emission disp_df[!,:value] = disp_df[!,:value]  ./ 1000 end
		disp_df[!,:variable] .= va

		# add empty values for non-existing columns
		for dim in (:Te,:C)
			if !(dim in namesSym(disp_df))
				disp_df[:,dim] .= 0
			end
		end

		# adjust sign, if enabled
		if wrtSgn && va in (:use,:stIn,:stIntIn,:stExtIn,:crt,:trdSell) disp_df[!,:value] = disp_df[!,:value] .* -1 end

		allData_df = vcat(allData_df,disp_df)
	end

	# XXX get exchange variables aggregated by import and export
	allExc_df = getAllVariables(:exc,anyM)
	if !isempty(allExc_df)
	    # add losses to all exchange variables
	    allExc_df = getExcLosses(convertExcCol(allExc_df),anyM.parts.exc.par,anyM.sets)
	    # compute export and import of each region, losses are considered at import
	    excFrom_df = rename(combine(groupby(allExc_df,[:Ts_disSup,:R_a,:C]),:var => ( x -> value(sum(x))/1000) => :value),:R_a => :R_dis)
	    excFrom_df[!,:variable] .= :export; excFrom_df[!,:Te] .= 0
		if wrtSgn excFrom_df[!,:value] = excFrom_df[!,:value] .* -1 end

	    excTo_df = rename(combine(x -> (value = value(dot(x.var,(1 .- x.loss)))/1000,),groupby(allExc_df,[:Ts_disSup,:R_b,:C])),:R_b => :R_dis)
	    excTo_df[!,:variable] .= :import; excTo_df[!,:Te] .= 0

	    allData_df = vcat(allData_df,vcat(excFrom_df,excTo_df))
	end

	# XXX get full load hours for conversion, storage input and storage output
	if anyM.options.decomm == :none
		flh_dic = Dict(:capaConv => :flhConv, :capaStIn => :flhStIn, :capaStOut => :flhStOut)
	else
		flh_dic = Dict(:oprCapaConv => :flhConv, :oprCapaStIn => :flhStIn, :oprCapaStOut => :flhStOut)
	end

	flhAss_dic = Dict(:capaConv => [:use,:stIntOut,:gen,:stIntIn],:oprCapaConv => [:use,:stIntOut,:gen,:stIntIn], :capaStIn => [:stIntIn,:stExtIn],:oprCapaStIn => [:stIntIn,:stExtIn], :capaStOut => [:stIntOut,:stExtOut],:oprCapaStOut => [:stIntOut,:stExtOut])

	for flhCapa in collect(keys(flh_dic))
		# gets relevant capacity variables
		capaFlh_df = rename(filter(x -> x.variable == flhCapa && x.value > 0.0, allData_df),:variable => :varCapa, :value => :valCapa)
		if isempty(capaFlh_df) continue end

		# expand with relevant dispatch variables
		capaFlh_df[!,:varDisp] = map(x -> flhAss_dic[x.varCapa],eachrow(capaFlh_df))
		capaFlh_df = flatten(capaFlh_df,:varDisp)
		intCol(capaFlh_df,:varDisp)
		# omit carrier in case of conversion capacity
		if flhCapa in (:capaConv,:oprCapaConv) select!(capaFlh_df,Not([:C])) end

		# joins capacity and with relevant dispatch variables
		allFLH_df = innerjoin(capaFlh_df,rename(allData_df,:variable => :varDisp), on = intCol(capaFlh_df,:varDisp))

		# remove row for gen variables in case gen exists as well and compute full load hours for conversion
		allFLH_df = combine(combine(y -> filter(x -> :use in y[!,:varDisp] ? x.varDisp != :gen : true,y), groupby(allFLH_df, intCol(capaFlh_df,:varCapa)), ungroup = false), AsTable([:value,:valCapa]) => (x -> 1000*abs(sum(x.value)/x.valCapa[1])) => :value)
		allFLH_df[!,:variable] = map(x -> flh_dic[x], allFLH_df[!,:varCapa])
		# adds carrier again in case of conversion capacity
		allFLH_df[!,:C] .= 0

		allData_df = vcat(allData_df, select(allFLH_df,Not([:varCapa])))
	end

	# XXX comptue storage cycles

	if anyM.options.decomm == :none
		cyc_dic = Dict(:capaStIn => :cycStIn, :capaStOut => :cycStOut)
	else
		cyc_dic = Dict(:oprCapaStIn => :cycStIn, :oprCapaStOut => :cycStOut)
	end

	cycAss_dic = Dict(:oprCapaStIn => [:stIntIn,:stExtIn], :capaStIn => [:stIntIn,:stExtIn], :oprCapaStOut => [:stIntOut,:stExtOut], :capaStOut => [:stIntOut,:stExtOut])

	for cycCapa in collect(keys(cyc_dic))
		capaCyc_df = rename(filter(x -> x.variable == :capaStSize && x.value > 0.0, allData_df),:variable => :varCapa, :value => :valCapa)
		if isempty(capaCyc_df) continue end

		# expand with relevant dispatch variables
		capaCyc_df[!,:varDisp] = map(x -> cycAss_dic[cycCapa],eachrow(capaCyc_df))
		capaCyc_df = flatten(capaCyc_df,:varDisp)

		# joins capacity and with relevant dispatch variables
		capaCyc_df = innerjoin(capaCyc_df,rename(allData_df,:variable => :varDisp), on = intCol(capaCyc_df,:varDisp))

		# compute cycling value and add to overall
		capaCyc_df = combine(groupby(capaCyc_df, intCol(capaCyc_df,:varCapa)), AsTable([:value,:valCapa]) => (x -> 1000*abs(sum(x.value)/x.valCapa[1])) => :value)
		capaCyc_df[!,:variable] = map(x -> cyc_dic[cycCapa], capaCyc_df[!,:varCapa])

		allData_df = vcat(allData_df,select(capaCyc_df,Not([:varCapa])))
	end

	# return dataframes and write csv files based on specified inputs
	if :csv in rtnOpt || :csvDf in rtnOpt
		csvData_df = printObject(allData_df,anyM, fileName = "results_summary",rtnDf = rtnOpt)
	end

	if :raw in rtnOpt
		CSV.write("$(anyM.options.outDir)/results_summary_$(anyM.options.outStamp).csv", allData_df)
	end

	if :rawDf in rtnOpt && :csvDf in rtnOpt
		return allData_df, csvData_df
	else
		if :rawDf in rtnOpt return allData_df end
		if :csvDf in rtnOpt return csvData_df end
	end
end

# XXX results for costs
function reportResults(objGrp::Val{:costs},anyM::anyModel; rtnOpt::Tuple{Vararg{Symbol,N} where N} = (:csv,))
	# prepare empty dataframe
	allData_df = DataFrame(Ts_disSup = Int[], R = Int[], Te = Int[], C = Int[], variable = Symbol[], value = Float64[])

	# loops over all objective variables with keyword "cost" in it
	for cst in filter(x -> occursin("cost",string(x)),keys(anyM.parts.obj.var))
		cost_df = copy(anyM.parts.obj.var[cst])
		# rename all dispatch and expansion regions simply to region
		if !isempty(intersect([:R_dis,:R_exp],namesSym(cost_df)))
			rename!(cost_df,:R_dis in namesSym(cost_df) ? :R_dis : :R_exp => :R)
		end
		# add empty column for non-existing dimensions
		for dim in (:Te,:C,:R)
			if !(dim in namesSym(cost_df))
				cost_df[:,dim] .= 0
			end
		end
		# obtain values and write to dataframe
		cost_df[:,:variable] .= string(cst)
		cost_df[:,:value] = value.(cost_df[:,:var])
        if :Ts_exp in namesSym(cost_df) cost_df = rename(cost_df,:Ts_exp => :Ts_disSup) end
		allData_df = vcat(allData_df,cost_df[:,Not(:var)])
	end

	# return dataframes and write csv files based on specified inputs
	if :csv in rtnOpt || :csvDf in rtnOpt
		csvData_df = printObject(allData_df,anyM, fileName = "results_costs", rtnDf = rtnOpt)
	end

	if :raw in rtnOpt
		CSV.write("$(anyM.options.outDir)/results_costs_$(anyM.options.outStamp).csv", allData_df)
	end

	if :rawDf in rtnOpt && :csvDf in rtnOpt
		return allData_df, csvData_df
	else
		if :rawDf in rtnOpt return allData_df end
		if :csvDf in rtnOpt return csvData_df end
	end
end

# XXX results for exchange
function reportResults(objGrp::Val{:exchange},anyM::anyModel; rtnOpt::Tuple{Vararg{Symbol,N} where N} = (:csv,))
	allData_df = DataFrame(Ts_disSup = Int[], R_from = Int[], R_to = Int[], C = Int[], variable = Symbol[], value = Float64[])
	if isempty(anyM.parts.exc.var) error("No exchange data found") end

	# XXX expansion variables
	if :expExc in keys(anyM.parts.exc.var)
		exp_df = copy(anyM.parts.exc.var[:expExc]) |> (x -> vcat(x,rename(x,:R_from => :R_to, :R_to => :R_from)))
		exp_df = flatten(exp_df,:Ts_expSup)
		select!(exp_df,Not(:Ts_disSup))
		rename!(exp_df,:Ts_expSup => :Ts_disSup)

		exp_df = combine(groupby(exp_df,[:Ts_disSup,:R_from,:R_to,:C]), :var => (x -> value.(sum(x))) => :value)
		exp_df[!,:variable] .= :expExc
	else
		exp_df = DataFrame(Ts_disSup = Int[], R_from = Int[], R_to = Int[], C = Int[], variable = Symbol[], value = Float64[])
	end

	# XXX capacity variables
	capa_df = copy(anyM.parts.exc.var[:capaExc])
	capa_df = vcat(capa_df,rename(filter(x -> x.dir == 0, capa_df),:R_from => :R_to, :R_to => :R_from))
	capa_df = combine(groupby(capa_df,[:Ts_disSup,:R_from,:R_to,:C]), :var => (x -> value.(sum(x))) => :value)
	capa_df[!,:variable] .= :capaExc

	if anyM.options.decomm != :none
		oprCapa_df = copy(anyM.parts.exc.var[:oprCapaExc])
		oprCapa_df = vcat(oprCapa_df,rename(filter(x -> x.dir == 0, oprCapa_df),:R_from => :R_to, :R_to => :R_from))
		oprCapa_df = combine(groupby(oprCapa_df,[:Ts_disSup,:R_from,:R_to,:C]), :var => (x -> value.(sum(x))) => :value)
		oprCapa_df[!,:variable] .= :oprCapaExc
		capa_df = vcat(capa_df,oprCapa_df)
	end

	# XXX dispatch variables
	disp_df = getAllVariables(:exc,anyM)
	disp_df = combine(groupby(disp_df,[:Ts_disSup,:R_from,:R_to,:C]), :var => (x -> value.(sum(x)) ./ 1000) => :value)
	disp_df[!,:variable] .= :exc

	# XXX get full load hours
	capaExt_df = replCarLeafs(copy(capa_df),anyM.sets[:C])
	flh_df = innerjoin(rename(select(capaExt_df,Not(:variable)),:value => :capa),rename(select(disp_df,Not(:variable)),:value => :disp),on = [:Ts_disSup,:R_from,:R_to,:C])
	flh_df[!,:value] = flh_df[!,:disp] ./ flh_df[!,:capa] .* 1000
	flh_df[!,:variable] .= :flhExc

	# XXX merge and print all data
	allData_df = vcat(exp_df,capa_df,disp_df,select(flh_df,Not([:capa,:disp])))

	# return dataframes and write csv files based on specified inputs
	if :csv in rtnOpt || :csvDf in rtnOpt
		csvData_df = printObject(allData_df,anyM, fileName = "results_exchange", rtnDf = rtnOpt)
	end

	if :raw in rtnOpt
		CSV.write("$(anyM.options.outDir)/results_exchange_$(anyM.options.outStamp).csv", allData_df)
	end

	if :rawDf in rtnOpt && :csvDf in rtnOpt
		return allData_df, csvData_df
	else
		if :rawDf in rtnOpt return allData_df end
		if :csvDf in rtnOpt return csvData_df end
	end
end

# XXX print time series for in and out into separate tables
"""
```julia
reportTimeSeries(car_sym::Symbol, model_object::anyModel)
```

Writes elements of energy balance for carrier specified by `car_sym` to `.csv` file. See [Time-series](@ref).
"""
function reportTimeSeries(car_sym::Symbol, anyM::anyModel; filterFunc::Function = x -> true, unstck::Bool = true, signVar::Tuple = (:in,:out), minVal::Number = 1e-3, mergeVar::Bool = true, rtnOpt::Tuple{Vararg{Symbol,N} where N} = (:csv,))

	# XXX converts carrier named provided to index
	node_arr = filter(x -> x.val == string(car_sym),collect(values(anyM.sets[:C].nodes)))
	if length(node_arr) != 1
		error("no carrier named $car_sym defined")
		return
	end
	c_int = node_arr[1].idx

	# XXX initialize dictionary to save data
	allData_dic = Dict{Symbol,DataFrame}()
	for signItr in signVar
		allData_dic[signItr] = DataFrame(Ts_disSup = Int[], Ts_dis = Int[], R_dis = Int[], variable = String[], value = Float64[])
	end

	# XXX initialize relevant dimensions and carriers
	relDim_df = filter(filterFunc,createPotDisp([c_int],anyM))
	sort!(relDim_df,sort(intCol(relDim_df)))
	relC_arr = unique([c_int,getDescendants(c_int,anyM.sets[:C])...])
	cRes_tup = anyM.cInfo[c_int] |> (x -> (Ts_dis = x.tsDis, R_dis = x.rDis, C = anyM.sets[:C].nodes[c_int].lvl))

	# XXX add demand and size it
	if :out in signVar
		dem_df = matchSetParameter(relDim_df,anyM.parts.bal.par[:dem],anyM.sets,newCol = :value)
		dem_df[!,:value] = dem_df[!,:value] .* getResize(dem_df,anyM.sets[:Ts],anyM.supTs) .* -1
		dem_df[!,:variable] .= "demand"
		filter!(x -> abs(x.value) > minVal, dem_df)
		allData_dic[:out] = vcat(allData_dic[:out],select!(dem_df,Not(:C)))
	end

	# XXX adds all technology related variables
	cBalRes_tup = anyM.cInfo[c_int] |> (x -> (x.tsDis, x.rDis))
	relType_tup = map(x -> x in signVar ? (x == :in ? (:use, :stExtIn) : (:gen,:stExtOut)) : tuple(),(:in,:out)) |> (x -> tuple(vcat(collect.(x)...)...))

	for c in relC_arr
		# gets technologies relevant for respective filterCarrier
		relTech_arr = getRelTech(c,anyM.parts.tech,anyM.sets[:C])

		if isempty(relTech_arr) continue end

		for x in relTech_arr

			# gets resolution and adjusts add_df in case of an agggregated technology
			add_df = select(filter(r -> r.C == c,anyM.parts.tech[x[1]].var[x[2]]),[:Ts_disSup,:Ts_dis,:R_dis,:var])
			tRes_tup = anyM.parts.tech[x[1]].disAgg ? (cRes_tup[1], anyM.parts.tech[x[1]].balLvl.exp[2]) : (cRes_tup[1], cRes_tup[2])
			checkTechReso!(tRes_tup,cBalRes_tup,add_df,anyM.sets)

			# filter values based on filter function and minimum value reported
			add_df = combine(groupby(add_df,[:Ts_disSup,:Ts_dis,:R_dis]), :var => (x -> sum(x)) => :var)
			filter!(filterFunc,add_df)
            if isempty(add_df) continue end
			add_df[!,:value] = value.(add_df[!,:var]) .* (x[2] in (:use,:stExtIn) ? -1.0 : 1.0)
			add_df[!,:variable] .= string(x[2],"; ", x[1])
			filter!(x -> abs(x.value) > minVal, add_df)

			# add to dictionary of dataframe for in or out
			sign_sym = x[2] in (:use,:stExtIn) ? :out : :in
			allData_dic[sign_sym] = vcat(allData_dic[sign_sym] ,select(add_df,Not(:var)))
		end
	end

	# XXX add import and export variables
    if :exc in keys(anyM.parts.exc.var)
		exc_df = filterCarrier(anyM.parts.exc.var[:exc],relC_arr)
		if :out in signVar
			excFrom_df = combine(groupby(filter(filterFunc,rename(copy(exc_df),:R_from => :R_dis)), [:Ts_disSup,:Ts_dis,:R_dis]), :var => (x -> value(sum(x)) * -1) => :value)
			excFrom_df[!,:variable] .= :export
			filter!(x -> abs(x.value) > minVal, excFrom_df)
			if !isempty(excFrom_df)
				allData_dic[:out] = vcat(allData_dic[:out],excFrom_df)
			end
		end

		if :in in signVar
			addLoss_df = rename(getExcLosses(convertExcCol(exc_df),anyM.parts.exc.par,anyM.sets),:R_b => :R_dis)
			excTo_df = combine(x -> (value = value(dot(x.var,(1 .- x.loss))),),groupby(filter(filterFunc,addLoss_df), [:Ts_disSup,:Ts_dis,:R_dis]))
			excTo_df[!,:variable] .= :import
			filter!(x -> abs(x.value) > minVal, excTo_df)
			if !isempty(excTo_df)
				allData_dic[:in] = vcat(allData_dic[:in],excTo_df)
			end
		end
	end

	# XXX add trade
	agg_arr = [:Ts_dis, :R_dis, :C]
	if !isempty(anyM.parts.trd.var)
		for trd in intersect(keys(anyM.parts.trd.var),(:trdBuy,:trdSell))
			trdVar_df = copy(relDim_df)
			trdVar_df[!,:value] = value.(filterCarrier(anyM.parts.trd.var[trd],relC_arr) |> (x -> aggUniVar(x,relDim_df,agg_arr,cRes_tup,anyM.sets))) .* (trd == :trdBuy ? 1.0 : -1.0)
			trdVar_df[!,:variable] .= trd
			filter!(x -> abs(x.value) > minVal, trdVar_df)
			sign_sym = :trdBuy == trd ? :in : :out
			allData_dic[sign_sym] = vcat(allData_dic[sign_sym],select(trdVar_df,Not(:C)))
		end
	end

	# XXX add curtailment
	if :crt in keys(anyM.parts.bal.var)
		crt_df = copy(relDim_df)
		crt_df[!,:value] = value.(filterCarrier(anyM.parts.bal.var[:crt],relC_arr) |> (x -> aggUniVar(x,crt_df,agg_arr, cRes_tup,anyM.sets))) .* -1.0
		crt_df[!,:variable] .= :crt
		filter!(x -> abs(x.value) > minVal, crt_df)
		allData_dic[:out] = vcat(allData_dic[:out],select(crt_df,Not(:C)))
	end

	# XXX add losted load
	if :lss in keys(anyM.parts.bal.var)
		lss_df = copy(relDim_df)
		lss_df[!,:value] = value.(filterCarrier(anyM.parts.bal.var[:lss],relC_arr) |> (x -> aggUniVar(x,lss_df,agg_arr, cRes_tup,anyM.sets)))
		lss_df[!,:variable] .= :lss
		filter!(x -> abs(x.value) > minVal, lss_df)
		allData_dic[:in] = vcat(allData_dic[:in],select(lss_df,Not(:C)))
	end

	# XXX unstack data and write to csv
	if mergeVar
		# merges in and out files and writes to same csv file
		data_df = vcat(values(allData_dic)...)

		if unstck && !isempty(data_df)
			data_df[!,:variable] = CategoricalArray(string.(data_df[!,:variable]))
			data_df = unstack(data_df,:variable,:value)
		end

		if :csv in rtnOpt || :csvDf in rtnOpt
			csvData_df = printObject(data_df,anyM, fileName = string("timeSeries_",car_sym,), rtnDf = rtnOpt)
		end

		if :raw in rtnOpt
			CSV.write("$(anyM.options.outDir)/$(string("timeSeries_",car_sym,))_$(anyM.options.outStamp).csv", data_df)
		end
	else
		# loops over different signs and writes to different csv files
		for signItr in signVar
			data_df = allData_dic[signItr]
			if unstck && !isempty(data_df)
				data_df[!,:variable] = CategoricalArray(string.(data_df[!,:variable]))
				data_df = unstack(data_df,:variable,:value)
			end

			if :csv in rtnOpt || :csvDf in rtnOpt
				csvData_df = printObject(data_df,anyM, fileName = string("timeSeries_",car_sym,"_",signItr), rtnDf = rtnOpt)
			end

			if :raw in rtnOpt
				CSV.write("$(anyM.options.outDir)/$(string("timeSeries_",car_sym,"_",signItr))_$(anyM.options.outStamp).csv", data_df)
			end
		end
	end

	# return dataframes based on specified inputs
	if :rawDf in rtnOpt && :csvDf in rtnOpt
		return data_df, csvData_df
	else
		if :rawDf in rtnOpt return data_df end
		if :csvDf in rtnOpt return csvData_df end
	end
end

# XXX write dual values for constraint dataframe
"""
```julia
printDuals(print_df::DataFrame, model_object::anyModel)
```

Writes duals of a constraint DataFrame to a `.csv` file in readable format (strings instead of ids). See [Individual elements](@ref).
"""
function printDuals(cns_df::DataFrame,anyM::anyModel;filterFunc::Function = x -> true, fileName::String = "", rtnOpt::Tuple{Vararg{Symbol,N} where N} = (:csv,))

    if !(:cns in namesSym(cns_df)) error("No constraint column found!") end
    cns_df = copy(filter(filterFunc,cns_df))
    cns_df[!,:dual] = dual.(cns_df[!,:cns])

	if :csv in rtnOpt || :csvDf in rtnOpt
    	csvData_df = printObject(select(cns_df,Not(:cns)),anyM;fileName = string("dual",fileName != "" ? "_" : "",fileName), rtnDf = rtnOpt)
	end

	if :rawDf in rtnOpt
		CSV.write("$(anyM.options.outDir)/$(string("dual",fileName != "" ? "_" : "",fileName))_$(anyM.options.outStamp).csv", data_df)
	end

	# return dataframes based on specified inputs
	if :rawDf in rtnOpt && :csvDf in rtnOpt
		return select(cns_df,Not(:cns)), csvData_df
	else
		if :rawDf in rtnOpt return data_df end
		if :csvDf in rtnOpt return csvData_df end
	end
end

# </editor-fold>

# <editor-fold desc="plotting tools"

# XXX plots tree graph for input set
"""
```julia
plotTree(tree_sym::Symbol, model_object::anyModel)
```

Plots the hierarchical tree of nodes for the set specified by `tree_sym`. See [Node trees](@ref).

"""
function plotTree(tree_sym::Symbol, anyM::anyModel; plotSize::Tuple{Float64,Float64} = (8.0,4.5), fontSize::Int = 12, useColor::Bool = true, wide::Array{Float64,1} = fill(1.0,30))

    netw = pyimport("networkx")
    plt = pyimport("matplotlib.pyplot")
    PyCall.fixqtpath()

    # <editor-fold desc="initialize variables"
    treeName_dic = Dict(:region => :R,:timestep => :Ts,:carrier => :C,:technology => :Te)

    # convert tree object into a data frame
    tree_obj = anyM.sets[treeName_dic[tree_sym]]
    data_arr = filter(x -> x.idx != 0,collect(values(tree_obj.nodes))) |> (y -> map(x -> getfield.(y,x),(:idx,:val,:lvl,:down,:subIdx)))
    tree_df = DataFrame(idx = data_arr[1], val = data_arr[2], lvl =  data_arr[3], down = data_arr[4], subIdx = data_arr[5], up =map(x -> tree_obj.up[x],data_arr[1]))

    # sets options
    col_dic = Dict(:region => (0.133, 0.545, 0.133),:timestep => (0.251,0.388,0.847),:carrier => (0.584, 0.345, 0.698),:technology => (0.796,0.235,0.2))
    # </editor-fold>

    # <editor-fold desc="computes positon of nodes"
    # adds a new dummy top node
    push!(tree_df,(0,"",0,tree_obj.nodes[0].down ,0,1))
    nodes_int = nrow(tree_df)
    idxPos_dic = Dict(zip(tree_df[:,:idx], 1:(nodes_int)))

    # create vertical position and labels from input tree
    locY_arr = float(tree_df[!,:lvl]) .+ 1.2

    # horizontal position is computed in a two step process
    locX_arr = zeros(Float64, nodes_int)

    # first step, filter all nodes at end of a respective branch and sort them correctly
    lowLvl_df = tree_df[isempty.(tree_df[!,:down]),:]
    lowLvl_df = lowLvl_df[map(y -> findall(x -> x == y, lowLvl_df[:,:idx])[1],sortSiblings(convert(Array{Int64,1},lowLvl_df[:,:idx]),tree_obj)),:]

    # sets distance from next node on the left depending on if they are part of the same subtree
    for (idx2, lowNode) in Iterators.drop(enumerate(eachrow(lowLvl_df)),1)
        if lowNode[:up] == lowLvl_df[idx2-1,:up] distance_fl = wide[lowNode[:lvl]] else distance_fl = 1 end
        locX_arr[idxPos_dic[lowNode[:idx]]] = locX_arr[idxPos_dic[lowLvl_df[idx2-1,:idx]]] + distance_fl
    end

    # second step, remaining horizontal nodes are placed in the middle of their children
    highLvl_df = tree_df[false .== isempty.(tree_df[!,:down]),:]
	highLvl_df = highLvl_df[map(y -> findall(x -> x == y, highLvl_df[:,:idx])[1],sortSiblings(convert(Array{Int64,1},highLvl_df[:,:idx]),tree_obj)),:]

    for highNode in reverse(eachrow(highLvl_df))
        locX_arr[idxPos_dic[highNode[:idx]]] = Statistics.mean(locX_arr[map(x -> idxPos_dic[x],highNode.down)])
    end

    locX_arr[end] = Statistics.mean(locX_arr[map(x -> idxPos_dic[x],tree_df[findall(tree_df[:,:lvl] .== 1),:idx])])
    locY_arr = abs.(locY_arr .- maximum(locY_arr))

    # compute dictionary of final node positions
    pos_dic = Dict(x => (locX_arr[x]/maximum(locX_arr),locY_arr[x]/maximum(locY_arr)) for x in 1:nodes_int)
    posIdx_dic = collect(idxPos_dic) |> (z -> Dict(Pair.(getindex.(z,2),getindex.(z,1))))
    # </editor-fold>

	# <editor-fold desc="determine node colors and labels"
	name_dic = anyM.graInfo.names

	label_dic = Dict(x[1] => x[2] == "" ? "" : name_dic[x[2]] for x in enumerate(tree_df[!,:val]))

	if useColor
		col_arr = [col_dic[tree_sym]]
	else
		col_arr = getNodeColors(collect(1:nodes_int),label_dic,anyM)
	end

	# </editor-fold>

    # <editor-fold desc="draw final tree"
    # draw single nodes
    edges_arr = Array{Tuple{Int,Int},1}()

    for rowTree in eachrow(tree_df)[1:end-1]
      # 0 node in tree_df becomes last node in graph, because there is 0 node within the plots
      if rowTree[:up] == 0 pare_int = nodes_int else pare_int = idxPos_dic[rowTree[:up]] end
      push!(edges_arr, (idxPos_dic[rowTree[:idx]], pare_int))
    end

    # draw graph object
    plt.clf()
    graph_obj = netw.Graph()

    netw.draw_networkx_nodes(graph_obj, pos_dic; nodelist = collect(1:nodes_int), node_color = col_arr)
    netw.draw_networkx_edges(graph_obj, pos_dic; edgelist = edges_arr)
    posLabOff_dic = netw.draw_networkx_labels(graph_obj, pos_dic, font_family = "arial", font_size = fontSize, labels = label_dic)

	figure = plt.gcf()
	figure.set_size_inches(plotSize[1],plotSize[2])

    r = figure.canvas.get_renderer()
    trans = plt.gca().transData.inverted()
    for x in collect(posLabOff_dic)
        down_boo = isempty(tree_obj.nodes[posIdx_dic[x[1]]].down)
        bb = x[2].get_window_extent(renderer=r)
        bbdata = bb.transformed(trans)
		# computes offset of label for leaves and non-leaves by first moving according to size auf letters itself (bbdata) and then by size of the node
		# (node-size in pixel is devided by dpi and plot size to get relative offset)
        offset_arr = [down_boo ? 0.0 : (bbdata.width/2.0 + (150/plotSize[1]/600)), down_boo ? (-bbdata.height/2.0 - 150/plotSize[2]/600) : 0.0]
        x[2].set_position([x[2]."_x" + offset_arr[1],x[2]."_y" + offset_arr[2]])
        x[2].set_clip_on(false)
    end

    # size plot and save
    plt.axis("off")
    plt.savefig("$(anyM.options.outDir)/$(tree_sym)_$(anyM.options.outStamp)", dpi = 600, bbox_inches="tight")
    # </editor-fold>
end

"""
```julia
plotEnergyFlow(plotType::Symbol, model_object::anyModel)
```

Plots the energy flow in a model. Set `plotType` to `:graph` for a qualitative node graph or to `:sankey` for a quantitative Sankey diagram. See [Energy flow](@ref).

"""
plotEnergyFlow(plotType::Symbol,anyM::anyModel; kwargs...) = plotEnergyFlow(Val{plotType}(),anyM::anyModel; kwargs...)

# XXX plot qualitative energy flow graph (applies python modules networkx and matplotlib via PyCall package)
function plotEnergyFlow(objGrp::Val{:graph},anyM::anyModel; plotSize::Tuple{Number,Number} = (16.0,9.0), fontSize::Int = 12, replot::Bool = true, scaDist::Number = 0.5, maxIter::Int = 5000, initTemp::Number = 2.0, useTeColor::Bool = false)

    # XXX import python function
    netw = pyimport("networkx")
    plt = pyimport("matplotlib.pyplot")
    PyCall.fixqtpath()

    # <editor-fold desc="create graph and map edges"

    graph_obj = netw.DiGraph()
    flowGrap_obj = anyM.graInfo.graph

    edges_arr =  vcat(collect.(flowGrap_obj.edgeC),collect.(flowGrap_obj.edgeTe))
    for x in edges_arr
        graph_obj.add_edge(x[1],x[2])
    end

    # </editor-fold>

    # <editor-fold desc="obtain and order graph properties (colors, names, etc.)"

	# get carriers that should be plotted, because they are connected with a technology
	relNodeC1_arr = filter(x -> x[2] in vcat(getindex.(flowGrap_obj.edgeTe,1),getindex.(flowGrap_obj.edgeTe,2)), collect(flowGrap_obj.nodeC))
	# get carriers that shold be plotted, because they are connected with another carrier that should be plotted
	relNodeC2_arr = filter(x -> any(map(y -> x[2] in y && !isempty(intersect(getindex.(relNodeC1_arr,2),y)) , collect.(flowGrap_obj.edgeC))), collect(flowGrap_obj.nodeC))

	# maps node id to node names
    idToC_arr = map(x -> x[2] => anyM.sets[:C].nodes[x[1]].val, filter(y -> y[2] in union(edges_arr...), intersect(flowGrap_obj.nodeC, union(relNodeC1_arr,relNodeC2_arr))))
    idToTe_arr  = map(x -> x[2] => anyM.sets[:Te].nodes[x[1]].val, filter(y -> y[2] in union(edges_arr...), collect(flowGrap_obj.nodeTe)))
    idToName_dic = Dict(vcat(idToC_arr,idToTe_arr))

    # obtain colors of nodes
    ordC_arr = intersect(unique(vcat(edges_arr...)), getindex.(idToC_arr,1))
	ordTe_arr = intersect(unique(vcat(edges_arr...)), getindex.(idToTe_arr,1))
    nodeC_arr = getNodeColors(ordC_arr,idToName_dic,anyM)
	nodeTe_arr = useTeColor ? getNodeColors(ordTe_arr,idToName_dic,anyM) : [(0.85,0.85,0.85)]

	nodesCnt_int = length(idToName_dic)

	# converts edges to sparse matrix for flowLayout function
	id_arr = vcat(getindex.(idToC_arr,1), getindex.(idToTe_arr,1))
	edges_mat = convert(Array{Int64,2},zeros(nodesCnt_int,nodesCnt_int))
	foreach(x -> edges_mat[findall(id_arr .== x[1])[1],findall(id_arr .== x[2])[1]] = 1, filter(x -> x[1] in id_arr && x[2] in id_arr,edges_arr))
	edges_smat = SparseArrays.sparse(edges_mat)

    # compute position of nodes
    if replot || !(isdefined(flowGrap_obj,:nodePos))
        pos_dic = flowLayout(nodesCnt_int,edges_smat; scaDist = scaDist, maxIter = maxIter, initTemp = initTemp)
		flowGrap_obj.nodePos = Dict(id_arr[x] => pos_dic[x] for x in keys(pos_dic))
    end

    # separate into edges between technologies and carriers and between carriers, then get respective colors
    cEdges_arr = filter(x -> x[1] in ordC_arr && x[2] in ordC_arr, collect(graph_obj.edges))
    edgeColC_arr = map(x -> anyM.graInfo.colors[idToName_dic[x[1]]], cEdges_arr)

    teEdges_arr = filter(x -> x[1] in ordTe_arr || x[2] in ordTe_arr, collect(graph_obj.edges))
    edgeColTe_arr = map(x -> x[1] in ordC_arr ? anyM.graInfo.colors[idToName_dic[x[1]]] : anyM.graInfo.colors[idToName_dic[x[2]]], teEdges_arr)

    # </editor-fold>

    # <editor-fold desc="draw and save graph with python"

    # plot final graph object
    plt.clf()

    netw.draw_networkx_nodes(graph_obj, flowGrap_obj.nodePos, nodelist = ordC_arr, node_shape="s", node_size = 300, node_color = nodeC_arr)
    netw.draw_networkx_nodes(graph_obj, flowGrap_obj.nodePos, nodelist = ordTe_arr, node_shape="o", node_size = 185,node_color = nodeTe_arr)

    netw.draw_networkx_edges(graph_obj, flowGrap_obj.nodePos, edgelist = cEdges_arr, edge_color = edgeColC_arr, arrowsize  = 16.2, width = 1.62)
    netw.draw_networkx_edges(graph_obj, flowGrap_obj.nodePos, edgelist = teEdges_arr, edge_color = edgeColTe_arr)

    posLabC_dic = netw.draw_networkx_labels(graph_obj, flowGrap_obj.nodePos, font_size = fontSize, labels = Dict(y[1] => anyM.graInfo.names[y[2]] for y in filter(x -> x[1] in ordC_arr,idToName_dic)), font_weight = "bold", font_family = "arial")
    posLabTe_dic = netw.draw_networkx_labels(graph_obj, flowGrap_obj.nodePos, font_size = fontSize, font_family = "arial", labels = Dict(y[1] => anyM.graInfo.names[y[2]] for y in filter(x -> !(x[1] in ordC_arr),idToName_dic)))

    # adjusts position of carrier labels so that they are right from node, uses code provided by ImportanceOfBeingErnest from here https://stackoverflow.com/questions/43894987/networkx-node-labels-relative-position
	figure = plt.gcf()
	figure.set_size_inches(plotSize[1],plotSize[2])

    r = figure.canvas.get_renderer()
    trans = plt.gca().transData.inverted()
    for x in vcat(collect(posLabC_dic),collect(posLabTe_dic))
		cNode_boo = x[1] in ordC_arr
        bb = x[2].get_window_extent(renderer=r)
        bbdata = bb.transformed(trans)
		# computes offset of label for leaves and non-leaves by first moving according to size auf letters itself (bbdata) and then by size of the node
		# (node-size in pixel is devided by dpi and plot size to get relative offset)
		offset_arr = [cNode_boo ? (bbdata.width/2.0 + (500/plotSize[1]/600)) : 0.0, cNode_boo ? 0.0 : (bbdata.height/2.0 + 200/plotSize[2]/600)]
		x[2].set_position([x[2]."_x" + offset_arr[1],x[2]."_y" + offset_arr[2]])
        x[2].set_clip_on(false)
    end

    plt.axis("off")

    # size plot and save
    plt.savefig("$(anyM.options.outDir)/energyFlowGraph_$(anyM.options.outStamp)", dpi = 600)

    # </editor-fold>
end

# XXX plot quantitative energy flow sankey diagramm (applies python module plotly via PyCall package)
function plotEnergyFlow(objGrp::Val{:sankey},anyM::anyModel; plotSize::Tuple{Number,Number} = (16.0,9.0), minVal::Float64 = 0.1, filterFunc::Function = x -> true, dropDown::Tuple{Vararg{Symbol,N} where N} = (:region,:timestep), rmvNode::Tuple{Vararg{String,N} where N} = tuple(), useTeColor = true)
    plt = pyimport("plotly")
    flowGrap_obj = anyM.graInfo.graph

    # <editor-fold desc="initialize data"

    if !isempty(setdiff(dropDown,[:region,:timestep]))
    error("dropDown only accepts array :region and :timestep as content")
    end

    # get mappings to create buttons of dropdown menue
    drop_dic = Dict(:region => :R_dis, :timestep => :Ts_disSup)
    dropDim_arr = collect(map(x -> drop_dic[x], dropDown))

    # get summarised data and filter dispatch variables
    data_df = reportResults(:summary,anyM,rtnOpt = (:rawDf,))
    filter!(x -> x.variable in (:demand,:gen,:use,:stIn,:stOut,:trdBuy,:trdSell,:demand,:import,:export,:lss,:crt),data_df)

    # filter non relevant entries
    filter!(x -> abs(x.value) > minVal, data_df)
    filter!(filterFunc, data_df)

    # create dictionaries for nodes that are neither technology nor carrier
    othNode_dic = maximum(values(flowGrap_obj.nodeTe)) |> (z -> Dict((x[2].C,x[2].variable) => x[1] + z for x in enumerate(eachrow(unique(filter(x -> x.Te == 0,data_df)[!,[:variable,:C]])))))
    othNodeId_dic = collect(othNode_dic) |> (z -> Dict(Pair.(getindex.(z,2),getindex.(z,1))))

    # </editor-fold>

    # <editor-fold desc="prepare labels and colors"

    # prepare name and color assignment
    names_dic = anyM.graInfo.names
    revNames_dic = collect(names_dic) |> (z -> Dict(Pair.(getindex.(z,2),getindex.(z,1))))
    col_dic = anyM.graInfo.colors

    sortTe_arr = getindex.(sort(collect(flowGrap_obj.nodeTe),by = x -> x[2]),1)
    cColor_dic = Dict(x => anyM.sets[:C].nodes[x].val |> (z -> z in keys(col_dic) ? col_dic[z] : (names_dic[z] in keys(col_dic) ? col_dic[col_dic[z]] : (0.85,0.85,0.85))) for x in sort(collect(keys(flowGrap_obj.nodeC))))

    # create array of node labels
    cLabel_arr = map(x -> names_dic[anyM.sets[:C].nodes[x].val],sort(collect(keys(flowGrap_obj.nodeC))))
    teLabel_arr = map(x -> names_dic[anyM.sets[:Te].nodes[x].val],sortTe_arr)
    othLabel_arr = map(x -> names_dic[String(othNodeId_dic[x][2])],sort(collect(keys(othNodeId_dic))))
    nodeLabel_arr = vcat(cLabel_arr, teLabel_arr, othLabel_arr)
    revNodelLabel_arr = map(x -> revNames_dic[x],nodeLabel_arr)

    # create array of node colors
    cColor_arr = map(x -> anyM.sets[:C].nodes[x].val |> (z -> z in keys(col_dic) ? col_dic[z] : (names_dic[z] in keys(col_dic) ? col_dic[names_dic[z]] : (0.85,0.85,0.85))),sort(collect(keys(flowGrap_obj.nodeC))))
    teColor_arr = map(x -> anyM.sets[:Te].nodes[x].val |> (z -> useTeColor && z in keys(col_dic) ? col_dic[z] : (useTeColor && names_dic[z] in keys(col_dic) ? col_dic[names_dic[z]] : (0.85,0.85,0.85))),sortTe_arr)
    othColor_arr = map(x -> anyM.sets[:C].nodes[othNodeId_dic[x][1]].val |> (z -> z in keys(col_dic) ? col_dic[z] : (names_dic[z] in keys(col_dic) ? col_dic[names_dic[z]] : (0.85,0.85,0.85))),sort(collect(keys(othNodeId_dic))))
    nodeColor_arr = vcat(map(x -> replace.(string.("rgb",string.(map(z -> z .* 255.0,x)))," " => ""),[cColor_arr, teColor_arr, othColor_arr])...)

    dropData_arr = Array{Dict{Symbol,Any},1}()

    # </editor-fold>

    # XXX loop over potential buttons in dropdown menue
    for drop in eachrow(unique(data_df[!,dropDim_arr]))
    # <editor-fold desc="filter data and create flow array"

    dropData_df = copy(data_df)
    if :region in dropDown subR_arr = [drop.R_dis, getDescendants(drop.R_dis,anyM.sets[:R],true)...] end
    for d in dropDown
      filter!(x -> d == :region ? x.R_dis in subR_arr : x.Ts_disSup == drop.Ts_disSup, dropData_df)
    end

    flow_arr = Array{Tuple,1}()

    # write flows reported in data summary
    for x in eachrow(dropData_df)
      a = Array{Any,1}(undef,3)

      # technology related entries
      if x.variable in (:demand,:export,:trdSell,:crt)
        a[1] = flowGrap_obj.nodeC[x.C]
        a[2] = othNode_dic[(x.C,x.variable)]
      elseif x.variable in (:import,:trdBuy,:lss)
        a[1] = othNode_dic[(x.C,x.variable)]
        a[2] = flowGrap_obj.nodeC[x.C]
      elseif x.variable in (:gen,:stOut)

    	  if x.Te in keys(flowGrap_obj.nodeTe) # if technology is not directly part of the graph, use its smallest parent that its
    		  a[1] = flowGrap_obj.nodeTe[x.Te]
    	  else
    		  a[1] = flowGrap_obj.nodeTe[minimum(intersect(keys(flowGrap_obj.nodeTe),getAncestors(x.Te,anyM.sets[:Te],:int)))]
    	  end

    	  a[2] = flowGrap_obj.nodeC[x.C]
      else
        a[1] = flowGrap_obj.nodeC[x.C]

    	if x.Te in keys(flowGrap_obj.nodeTe)
    		a[2] = flowGrap_obj.nodeTe[x.Te]
    	else
    		a[2] = flowGrap_obj.nodeTe[minimum(intersect(keys(flowGrap_obj.nodeTe),getAncestors(x.Te,anyM.sets[:Te],:int)))]
    	end
      end

      a[3] = abs(x.value)

      push!(flow_arr,tuple(a...))
    end

    # create flows connecting different carriers
    idToC_dic = Dict(map(x -> x[2] => x[1], collect(flowGrap_obj.nodeC)))
    for x in filter(x -> anyM.sets[:C].up[x] != 0,intersect(union(getindex.(flow_arr,1),getindex.(flow_arr,2)),values(flowGrap_obj.nodeC)))
      a = Array{Any,1}(undef,3)
      a[1] = flowGrap_obj.nodeC[x]
      a[2] = flowGrap_obj.nodeC[anyM.sets[:C].up[x]]
      a[3] = (getindex.(filter(y -> y[2] == x,flow_arr),3) |> (z -> isempty(z) ? 0.0 : sum(z))) - (getindex.(filter(y -> y[1] == x,flow_arr),3) |> (z -> isempty(z) ? 0.0 : sum(z)))
      push!(flow_arr,tuple(a...))
    end

    # merges flows for different regions that connect the same nodes
    flow_arr = map(unique(map(x -> x[1:2],flow_arr))) do fl
      allFl = filter(y -> y[1:2] == fl[1:2],flow_arr)
      return (allFl[1][1],allFl[1][2],sum(getindex.(allFl,3)))
    end

    # removes nodes accoring function input provided
    for rmv in rmvNode
      # splits remove expression by semicolon and searches for first part
      rmvStr_arr = split(rmv,"; ")
      relNodes_arr = findall(nodeLabel_arr .== rmvStr_arr[1])
      if isempty(relNodes_arr) relNodes_arr = findall(revNodelLabel_arr .== rmvStr_arr[1]) end
      if isempty(relNodes_arr) continue end

      if length(rmvStr_arr) == 2 # if rmv contains two strings separated by a semicolon, the second one should relate to a carrier, carrier is searched for and all related flows are removed
        relC_arr = findall(nodeLabel_arr .== rmvStr_arr[2])
        if isempty(relNodes_arr) relC_arr = findall(revNodelLabel_arr .== rmvStr_arr[2]) end

        if isempty(relC_arr)
            produceMessage(anyM.options,anyM.report, 1," - Remove string contained a carrier not found in graph, check for typos: "*rmv)
            continue
        else
            c_int = relC_arr[1]
        end

        filter!(x -> !((x[1] in relNodes_arr || x[2] in relNodes_arr) && (x[1] == c_int || x[2] == c_int)),flow_arr)
      elseif length(rmvStr_arr) > 2
        error("one remove string contained more then one semicolon, this is not supported")
      else # if rmv only contains one string, only nodes where in- and outgoing flow are equal or only one of both exists
        out_tup = filter(x -> x[1] == relNodes_arr[1],flow_arr)
        in_tup = filter(x -> x[2] == relNodes_arr[1],flow_arr)

        if length(out_tup) == 1 && length(in_tup) == 1 && out_tup[1][3] == in_tup[1][3] # in- and outgoing are the same
          filter!(x -> !(x in (out_tup[1],in_tup[1])),flow_arr)
          push!(flow_arr,(in_tup[1][1],out_tup[1][2],in_tup[1][3]))
        elseif length(out_tup) == 0 # only ingoing flows
          filter!(x -> !(x in in_tup),flow_arr)
        elseif length(in_tup) == 0 # only outgoing flows
          filter!(x -> !(x in out_tup),flow_arr)
        end
      end
    end

    # </editor-fold>

    # <editor-fold desc="create dictionaries for later plotting"

    # collect data for drop in a dictionary

    linkColor_arr = map(x -> collect(x[1] in keys(cColor_dic) ? cColor_dic[x[1]] : cColor_dic[x[2]]) |>
    	(z -> replace(string("rgba",string(tuple([255.0 .*z..., (x[1] in keys(cColor_dic) && x[2] in keys(cColor_dic) ? 0.8 : 0.5)]...)))," " => "")), flow_arr)
    link_dic = Dict(:source => getindex.(flow_arr,1) .- 1, :target => getindex.(flow_arr,2) .- 1, :value => getindex.(flow_arr,3), :color => linkColor_arr)

    fullData_arr = [Dict(:link => link_dic, :node => Dict(:label => nodeLabel_arr, :color => nodeColor_arr))]

    # pushes dictionary to overall array
    label_str = string("<b>",join(map(y -> anyM.sets[Symbol(split(String(y),"_")[1])].nodes[drop[y]].val,dropDim_arr),", "),"</b>")
    push!(dropData_arr,Dict(:args => fullData_arr, :label => label_str, :method => "restyle"))

    # </editor-fold>
    end

    # <editor-fold desc="create various dictionaries to define format and create plot"

    menues_dic =[Dict(:buttons => dropData_arr, :direction => "down", :pad => Dict(:l => 10, :t => 10), :font => Dict(:size => 16, :family => "Arial"), :showactive => true, :x => 0.01, :xanchor => "center", :y => 1.1, :yanchor => "middle")]
    data_dic = Dict(:type => "sankey", :orientation => "h", :valueformat => ".0f", :textfont => Dict(:family => "Arial"), :node => Dict(:pad => 8, :thickness => 36, :line => Dict(:color => "white",:width => 0.01), :hoverinfo => "skip"))
    layout_dic = Dict(:width => 125*plotSize[1], :height => 125*plotSize[2], :updatemenus => menues_dic, :font => Dict(:size => 32, :family => "Arial"))

    fig = Dict(:data => [data_dic], :layout => layout_dic)
    plt.offline.plot(fig, filename="$(anyM.options.outDir)/energyFlowSankey_$(join(string.(dropDown),"_"))_$(anyM.options.outStamp).html")

    # </editor-fold>
end

# XXX define postions of nodes in energy flow graph
# function is mostly taken from [GraphPlot.jl](https://github.com/JuliaGraphs/GraphPlot.jl), who again reference the following source [IainNZ](https://github.com/IainNZ)'s [GraphLayout.jl](https://github.com/IainNZ/GraphLayout.jl)
function flowLayout(nodesCnt_int::Int,edges_smat::SparseMatrixCSC{Int64,Int64}, locsX_arr::Array{Float64,1} = 2*rand(nodesCnt_int).-1.0, locsY_arr::Array{Float64,1} = 2*rand(nodesCnt_int).-1.0; scaDist::Number = 0.5, maxIter::Int=5000, initTemp::Number=2.0)

    # optimal distance bewteen vertices
    k = scaDist * sqrt(4.0 / nodesCnt_int)
    k² = k * k

    # store forces and apply at end of iteration all at once
    force_x = zeros(nodesCnt_int)
    force_y = zeros(nodesCnt_int)

    # iterate maxIter times
    @inbounds for iter = 1:maxIter
        # Calculate forces
        for i = 1:nodesCnt_int
            force_vec_x = 0.0
            force_vec_y = 0.0
            for j = 1:nodesCnt_int
                i == j && continue
                d_x = locsX_arr[j] - locsX_arr[i]
                d_y = locsY_arr[j] - locsY_arr[i]
                dist²  = (d_x * d_x) + (d_y * d_y)
                dist = sqrt(dist²)

                if !( iszero(edges_smat[i,j]) && iszero(edges_smat[j,i]) )
                    # Attractive + repulsive force
                    # F_d = dist² / k - k² / dist # original FR algorithm
                    F_d = dist / k - k² / dist²
                else
                    # Just repulsive
                    # F_d = -k² / dist  # original FR algorithm
                    F_d = -k² / dist²
                end
                force_vec_x += F_d*d_x
                force_vec_y += F_d*d_y
            end
            force_x[i] = force_vec_x
            force_y[i] = force_vec_y
        end
        # Cool down
        temp = initTemp / iter
        # Now apply them, but limit to temperature
        for i = 1:nodesCnt_int
            fx = force_x[i]
            fy = force_y[i]
            force_mag  = sqrt((fx * fx) + (fy * fy))
            scale      = min(force_mag, temp) / force_mag
            locsX_arr[i] += force_x[i] * scale
            locsY_arr[i] += force_y[i] * scale
        end
    end

    # Scale to unit square
    min_x, max_x = minimum(locsX_arr), maximum(locsX_arr)
    min_y, max_y = minimum(locsY_arr), maximum(locsY_arr)
    function scaler(z, a, b)
        2.0*((z - a)/(b - a)) - 1.0
    end
    map!(z -> scaler(z, min_x, max_x), locsX_arr, locsX_arr)
    map!(z -> scaler(z, min_y, max_y), locsY_arr, locsY_arr)

    # converts positions into dictionary
    pos_dic = Dict(z => [locsX_arr[z]*16/9,locsY_arr[z]] for z in 1:nodesCnt_int)

    return pos_dic
end

# XXX returns array of colors for input nodes, which labels can be found in label_dic
function getNodeColors(node_arr::Array{Int,1}, label_dic::Dict{Int64,String},anyM::anyModel)
	revName_dic = collect(anyM.graInfo.names) |> (z -> Dict(Pair.(getindex.(z,2),getindex.(z,1))))
	col_dic = anyM.graInfo.colors

	color_arr = map(node_arr) do x
		str = label_dic[x]
		if str in keys(col_dic) # label is key in color dictionary
			return col_dic[str]
		elseif str in keys(anyM.graInfo.names) && anyM.graInfo.names[str] in keys(col_dic) # internal name is key in dictionary, but label was external
			return col_dic[anyM.graInfo.names[str]]
		elseif str in keys(revName_dic) && revName_dic[str] in keys(col_dic)  # external name is key in dictionary, but label was internal
			return col_dic[revName_dic[str]]
		else # default color
			return (0.85,0.85,0.85)
		end
	end

	return color_arr
end

# XXX move a node after positions were created within energy flow graph
"""
```julia
moveNode!(model_object::anyModel, newPos_arr::Union{Array{Tuple{String,Array{Float64,1}},1},Tuple{String,Array{Float64,1}}})
```

Moves a node within the current layout of the node graph created with `plotEnergyFlow`. See [Energy flow](@ref).

"""
function moveNode!(anyM::anyModel,newPos_arr::Union{Array{Tuple{String,Array{Float64,1}},1},Tuple{String,Array{Float64,1}}})

    flowGrap_obj = anyM.graInfo.graph

    if !isdefined(flowGrap_obj,:nodePos)
        error("Initial positions are not yet defined. Run 'plotEnergyFlow' first.")
    end

    # gets assignment between node ids and names
    edges_arr =  vcat(collect.(flowGrap_obj.edgeC),collect.(flowGrap_obj.edgeTe))

    cToId_arr = map(x -> anyM.sets[:C].nodes[x[1]].val => x[2], filter(y -> y[2] in union(edges_arr...), collect(flowGrap_obj.nodeC)))
    teToId_arr  = map(x -> anyM.sets[:Te].nodes[x[1]].val => x[2], filter(y -> y[2] in union(edges_arr...), collect(flowGrap_obj.nodeTe)))
    nameToId_dic = Dict(vcat(teToId_arr,cToId_arr))

    # if input is just a single tuple and not an array convert to array
    if typeof(newPos_arr) == Tuple{String,Array{Float64,1}}
        newPos_arr = [newPos_arr]
    end

    switchNames_dic = Dict(map(x -> x[2] => x[1],collect(anyM.graInfo.names)))

    # loops overa array of moved notes
    for newPos in newPos_arr
        # get id of node depending on whether it is an orignial name or name just used in plot
        if newPos[1] in keys(nameToId_dic)
            x = nameToId_dic[newPos[1]]
        elseif newPos[1] in values(anyM.graInfo.names)
            x = nameToId_dic[switchNames_dic[newPos[1]]]
        else
            error("Node name not recognized!")
        end
        # actually adjust node position
        flowGrap_obj.nodePos[x] = [flowGrap_obj.nodePos[x][1] + newPos[2][1]*2, flowGrap_obj.nodePos[x][2] + newPos[2][2]*2]
    end
end

# </editor-fold>

# XXX dummy function just do provide a docstring for printIIS (docstring in printIIS wont work, because read-in is conditional)
"""
```julia
printIIS(model_object::anyModel)
```

Uses Gurobi's computeIIS function to determine the constraints of the optimization problem that cause infeasibility.
"""
function printIIS(anyM::anyModel,d::Int)
end

max_stars_repo_name: leonardgoeke/anyMOD

max_issues_repo_name: leonardgoeke/anyMOD

max_forks_repo_name: leonardgoeke/anyMOD


text:
# Supported domains: Real, Boolean, Ordinal, Periodic, Count

# The purpose of domains is to be able to impute over different possible values of `a` regardless of
# the loss that was used in the GLRM. The reason for doing this is to evaluate the performance of GLRMS.
# For instance, let's say we use PCA (QuadLoss losses) to model a binary data frame (not the best idea).
# In order to override the standard imputation with `impute(QuadLoss(), u)`, which assumes imputation over the reals,
# we can use `impute(BoolDomain(), QuadLoss(), u)` and see which of {-1,1} is best. The reason we want to be able to
# do this is to compare a baseline model (e.g. PCA) with a more logical model using heterogenous losses,
# yet still give each model the same amount of information regarding how imputation should be done.

# In order to accomplish this we define a series of domains that tell imputation methods
# what values the data can take. The imputation methods are defined in impute_and_err.jl

# Domains should be assigned to each column of the data and are not part of the low-rank model itself.
# They serve as a way to evaluate the performance of the low-rank model.

export Domain, # the abstract type
	   RealDomain, BoolDomain, OrdinalDomain, PeriodicDomain, CountDomain, CategoricalDomain, # the domains
	   copy

abstract type Domain end

########################################## REALS ##########################################
# Real data can take values from ℜ
struct RealDomain<:Domain
end

########################################## BOOLS ##########################################
# Boolean data should take values from {true, false}
struct BoolDomain<:Domain
end

########################################## ORDINALS ##########################################
# Ordinal data should take integer values ranging from `min` to `max`
struct OrdinalDomain<:Domain
	min::Int
	max::Int
	function OrdinalDomain(min, max)
		if max - min < 2
			@warn("The ordinal variable you've created is degenerate: it has only two levels. Consider using a Boolean variable instead; ordinal loss functions may have unexpected behavior on a degenerate ordinal domain.")
		end
		return new(min, max)
	end
end

########################################## ORDINALS ##########################################
# Categorical data should take integer values ranging from 1 to `max`
struct CategoricalDomain<:Domain
	min::Int
	max::Int
end
CategoricalDomain(m::Int) = CategoricalDomain(1,m)

########################################## PERIODIC ##########################################
# Periodic data can take values from ℜ, but given a period T, we should have error_metric(a,a+T) = 0
struct PeriodicDomain<:Domain
	T::Float64 # the period
end

########################################## COUNTS ##########################################
# Count data can take values over ℕ, which we approximate as {0, 1, 2 ... `max_count`}
struct CountDomain<:Domain
	max_count::Int # the biggest possible count
end

ext: jl


text:
"""
Simulates a SFFM defined by `model` until the `StoppingTime` has occured,
given the `InitialCondition` on (φ(0),X(0),Y(0)).

    SimSFFM(
        model::Model,
        StoppingTime::Function,
        InitCondition::NamedTuple{(:φ, :X, :Y)},
    )

# Arguments
- `model`: A Model object
- `StoppingTime`: A function which takes the value of the process at the current
    time and at the time of the last jump of the phase process, as well as the
    `model` object.
    i.e. `StoppingTime(;model,SFFM,SFFM0)` where `SFFM` and `SFFM0` are tuples
    with keys `(:t::Float64, :φ::Int, :X::Float64, :Y::Float64, :n::Int)` which
    are the value of the SFFM at the current time, and time of the previous jump
    of the phase process, repsectively. The `StoppingTime` must return a
    `NamedTuple{(:Ind, :SFFM)}` type where `:Ind` is a `:Bool` value stating
    whether the stopping time has occured or not and `:SFFM` is a tuple in the
    same form as the input `SFFM` but which contains the value of the SFFM at
    the stopping time.
- `InitCondition`: `NamedTuple` with keys `(:φ, :X, :Y)`, `InitCondition.φ` is a
    vector of length `M` of initial states for the phase, `InitCondition.X` is a
    vector of length `M` of initial states for the X-level, `InitCondition.Y` is
    a vector of length `M` of initial states for the Y-level. `M` is the number
    of simulations to be done.

# Output
- a tuple with keys
    - `t::Array{Float64,1}` a vector of length `M` containing the values of
        `t` at the `StoppingTime`.
    - `φ::Array{Float64,1}` a vector of length `M` containing the values of
        `φ` at the `StoppingTime`.
    - `X::Array{Float64,1}` a vector of length `M` containing the values of
        `X` at the `StoppingTime`.
    - `Y::Array{Float64,1}` a vector of length `M` containing the values of
        `Y` at the `StoppingTime`.
    - `n::Array{Float64,1}` a vector of length `M` containing the number of
        transitions of `φ` at the `StoppingTime`
"""
function DiscretisedFluidQueues.simulate(
    ffq::FluidFluidQueue{<:DiscretisedFluidQueue},
    StoppingTime::Function,
    InitCondition::NamedTuple{(:φ, :X, :Y)},
    rng::Random.AbstractRNG=Random.default_rng(),
)
    model = ffq.dq.model
    # find transition probabilities of phase process
    d = LinearAlgebra.diag(model.T)
    P = (model.T - LinearAlgebra.diagm(0 => d)) ./ -d
    # convenient to sample transitions with
    CumP = cumsum(P, dims = 2)

    # transition probabilities upon hitting a boundary
    # lower
    idx₋ = findall(rates(model).<0.0)
    CumP_lwr = zeros(size(P))
    CumP_lwr[idx₋,:] = cumsum(model.P_lwr, dims = 2)
    # upper
    idx₊ = findall(rates(model).>0.0)
    CumP_upr = zeros(size(P))
    CumP_upr[idx₊,:] = cumsum(model.P_upr, dims = 2)

    # rates for exponential holding times in each phase
    Λ = LinearAlgebra.diag(model.T)

    # containers to store the sims
    M = length(InitCondition.φ)
    tSims = Array{Float64,1}(undef, M)
    φSims = Array{Float64,1}(undef, M)
    XSims = Array{Float64,1}(undef, M)
    YSims = Array{Float64,1}(undef, M)
    nSims = Array{Float64,1}(undef, M)

    # simulate M realisations
    for m = 1:M
        SFM0 = (t=0.0, φ=InitCondition.φ[m], X=InitCondition.X[m], n=0.0)
        Y0 = InitCondition.Y[m]
        # check if initial condition is invalid and return NaNs if it is
        if !(0.0 <= SFM0.X <= model.b) || !(0.0 <= Y0) || !in(SFM0.φ, 1:n_phases(model))
            (tSims[m], φSims[m], XSims[m], nSims[m]) =
                (t = NaN, φ = NaN, X = NaN, n = NaN)
            YSims[m] = NaN
        else
            cell_idx0 = initialise_cell_idx(ffq,SFM0)
            while true
                # generate random exponential holding time via inverse CDF method
                S = log(rand(rng)) / Λ[SFM0.φ]
                X, φ, t = DiscretisedFluidQueues.Update(model, SFM0, S, CumP, CumP_lwr, CumP_upr, rng)
                cell_idx = find_cell_idx(ffq,SFM0.φ,cell_idx0,X) # find the cell in which X resides
                SFM = (t=t, φ=φ, X=X, n=SFM0.n+1)
                Y = UpdateYt(ffq, SFM, SFM0, Y0, cell_idx, cell_idx0)
                τ = StoppingTime(ffq, SFM, SFM0, Y, Y0)
                if τ.Ind
                    (tSims[m], φSims[m], XSims[m], nSims[m]) = τ.SFM
                    YSims[m] = τ.Y
                    break
                end
                cell_idx0 = cell_idx
                SFM0 = SFM
                Y0 = Y
            end
        end
    end
    return (t = tSims, φ = φSims, X = XSims, n = nSims), YSims
end

function initialise_cell_idx(ffq,SFM0)
    cell_idx0 = findfirst(ffq.dq.mesh.nodes .> SFM0.X)
    if (cell_idx0===nothing)
        if (rates(ffq.dq,SFM0.φ)<0.0)
            cell_idx0 = length(ffq.dq.mesh.nodes)-1
        else
            cell_idx0 = length(ffq.dq.mesh.nodes)
        end
    elseif (rates(ffq.dq,SFM0.φ)<=0.0)&&(SFM0.X==0.0)
        cell_idx0 = 0
    else
        cell_idx0 -= 1
    end
    return cell_idx0
end

function find_cell_idx(ffq,φ0,cell_idx0,X)
    # cell_idx = -1
    if (X==0.0)
        cell_idx = 0
    elseif X==ffq.dq.model.b
        cell_idx = length(ffq.dq.mesh.nodes)
    elseif (rates(ffq.dq,φ0)>0.0)
        cell_idx = NaN
        for n in max(1,cell_idx0):length(ffq.dq.mesh.nodes)
            if ffq.dq.mesh.nodes[n]>X
                cell_idx = n-1
                break
            end
        end
    elseif (rates(ffq.dq,φ0)<0.0)
        cell_idx = NaN
        for n in min(length(ffq.dq.mesh.nodes)-1,cell_idx0):-1:1
            if ffq.dq.mesh.nodes[n]<X
                cell_idx = n
                break
            end
        end
    else
        cell_idx = cell_idx0
    end
    return cell_idx
end

"""
Returns ``Y(t+S)`` given ``Y(t)``.

    UpdateYt(
        model::Model,
        SFM0::NamedTuple{(:t, :φ, :X, :n)},
        S::Real,
    )

# Arguments
- `model`: a Model object
- `SFFM0::NamedTuple` containing at least the keys `:X` giving the value of
    ``X(t)`` at the current time, and `:Y` giving the value of ``Y(t)`` at the
    current time, and `:φ` giving the value of `φ(t)`` at the current time.
- `S::Real`: an elapsed amount of time to evaluate ``X`` at, i.e. ``X(t+S)``.
"""
function UpdateYt(
    ffq::FluidFluidQueue,
    SFM::NamedTuple{(:t, :φ, :X, :n)},
    SFM0::NamedTuple{(:t, :φ, :X, :n)},
    Y0::Float64,
    cell_idx::Int, cell_idx0::Int,
)
    # given the position of a SFM, SFM0 at time t, a time step of size s, find the
    # position of Y at time t+s
    speedX = rates(ffq.dq,SFM0.φ)
    if (speedX==0.0)||((speedX<=0.0)&&(SFM0.X==0.0))||((speedX>=0.0)&&(SFM0.X==ffq.dq.model.b))
        # determine if X is at a boundary
        Y = Y0 + (SFM.t-SFM0.t) * getrates(ffq,SFM0.φ,cell_idx0)
    else
        # determine if X has hit a boundary (did not start at the same boundary)
        if (SFM.X==0.0)||(SFM.X==ffq.dq.model.b)
            if speedX>0.0 # must be upper boundary in this case
                cells_traversed_fully = cell_idx0+1:cell_idx-1
                time_spent_in_cells_traversed = [
                    (ffq.dq.mesh.nodes[cell_idx0+1]-SFM0.X)./speedX; # time to leave cell_idx0
                    [Δ(ffq.dq.mesh,k) for k in cells_traversed_fully]./speedX;
                    #(SFM.X-ffq.dq.mesh.nodes[cell_idx])./speedX; # time since last change of cell
                ]
                Y = Y0 + dot(getrates(ffq,SFM0.φ,cell_idx0:cell_idx-1),time_spent_in_cells_traversed)
            elseif speedX<0.0 # must be lower boundary in this case
                cells_traversed_fully = cell_idx+1:cell_idx0-1
                time_spent_in_cells_traversed = -[
                    #(ffq.dq.mesh.nodes[cell_idx+1]-SFM.X)./speedX; # time since last change of cell
                    [Δ(ffq.dq.mesh,k) for k in cells_traversed_fully]./speedX;
                    (SFM0.X-ffq.dq.mesh.nodes[cell_idx0])./speedX; # time to leave first cell
                ]
                Y = Y0 + dot(getrates(ffq,SFM0.φ,cell_idx+1:cell_idx0),time_spent_in_cells_traversed)
            end
        else
            if speedX>0.0
                if cell_idx0==cell_idx
                    time_spent_in_cells_traversed = (SFM.X-SFM0.X)./speedX
                    Y = Y0 + dot(getrates(ffq,SFM0.φ,cell_idx0),time_spent_in_cells_traversed)
                else
                    cells_traversed_fully = cell_idx0+1:cell_idx-1
                    time_spent_in_cells_traversed = [
                        (ffq.dq.mesh.nodes[cell_idx0+1]-SFM0.X)./speedX; # time to leave cell_idx0
                        [Δ(ffq.dq.mesh,k) for k in cells_traversed_fully]./speedX;
                        (SFM.X-ffq.dq.mesh.nodes[cell_idx])./speedX; # time since last change of cell
                    ]
                    Y = Y0 + dot(getrates(ffq,SFM0.φ,cell_idx0:cell_idx),time_spent_in_cells_traversed)
                end
            elseif speedX<0.0
                if cell_idx0==cell_idx
                    time_spent_in_cells_traversed = -(SFM0.X-SFM.X)./speedX
                    Y = Y0 + dot(getrates(ffq,SFM0.φ,cell_idx0),time_spent_in_cells_traversed)
                else
                    cells_traversed_fully = cell_idx+1:cell_idx0-1
                    time_spent_in_cells_traversed = -[
                        (ffq.dq.mesh.nodes[cell_idx+1]-SFM.X)./speedX; # time since last change of cell
                        [Δ(ffq.dq.mesh,k) for k in cells_traversed_fully]./speedX;
                        (SFM0.X-ffq.dq.mesh.nodes[cell_idx0])./speedX; # time to leave first cell
                    ]
                    Y = Y0 + dot(getrates(ffq,SFM0.φ,cell_idx:cell_idx0),time_spent_in_cells_traversed)
                end
            end
        end
    end
    return Y
end

"""
Constructs the `StoppingTime` which is the first exit of the process ``Y(t)``
from the interval ``[u,v]``. ASSUMES ``Y(t)`` is monotonic between jumps.

    first_exit_y( u::Real, v::Real)

# Arguments
- `u`: a lower boundary
- `v`: an upper boundary

# Output
- `FirstExitYFun`: a function with one method
    - `FirstExitYFun(
        model::Model,
        SFFM::NamedTuple{(:t, :φ, :X, :Y, :n)},
        SFFM0::NamedTuple{(:t, :φ, :X, :Y, :n)},
    )`: a stopping time for a SFFM
"""
function first_exit_y( u::Real, v::Real)
    # SFFM Method
    function first_exit_yFun(
        ffq::FluidFluidQueue,
        SFM::NamedTuple{(:t, :φ, :X, :n)},
        SFM0::NamedTuple{(:t, :φ, :X, :n)},
        Y::Float64,
        Y0::Float64,
    )
        Ind = ((Y < u) || (Y > v))
        if Ind
            idx = [Y < u; Y > v]
            boundaryHit = only([u;v][idx])
            cell_idx0 = initialise_cell_idx(ffq,SFM0)
            SFM_path(t) = (t=SFM0.t+t,φ=SFM0.φ,X=DiscretisedFluidQueues.UpdateX(ffq.dq.model,SFM0,t),n=SFM0.n)
            function YFun(t)
                SFMt = SFM_path(t)
                cell_idx = find_cell_idx(ffq,SFM0.φ,cell_idx0,SFMt.X)
                return UpdateYt(ffq, SFMt, SFM0, Y0, cell_idx, cell_idx0) - boundaryHit
            end
            S = SFM.t - SFM0.t
            tstar = fzero(YFun, 1e-64, S)
            X = DiscretisedFluidQueues.UpdateX(ffq.dq.model, SFM0, tstar)
            t = SFM0.t + tstar
            Y = YFun(tstar)+boundaryHit
            SFM = (t, SFM0.φ, X, SFM0.n)
        end
        return (Ind = Ind, SFM = SFM, Y=Y)
    end
    return first_exit_yFun
end

function n_jumps_y(N::Int)
    # SFFM Method
    function n_jumps_yFun(
        ffq::FluidFluidQueue,
        SFM::NamedTuple{(:t, :φ, :X, :n)},
        SFM0::NamedTuple{(:t, :φ, :X, :n)},
        Y::Float64,
        Y0::Float64,
    )
        Ind = SFM.n>=N
        return (Ind = Ind, SFM = SFM, Y=Y)
    end
    return n_jumps_yFun
end

function fixed_time_y(T::Real)
    # SFFM Method
    function fixed_time_yFun(
        ffq::FluidFluidQueue,
        SFM::NamedTuple{(:t, :φ, :X, :n)},
        SFM0::NamedTuple{(:t, :φ, :X, :n)},
        Y::Float64,
        Y0::Float64,
    )
        Ind = SFM.t>=T
        if Ind
            S = T-SFM0.t
            X = DiscretisedFluidQueues.UpdateX(ffq.dq.model,SFM0,S)
            cell_idx = find_cell_idx(ffq,SFM0.φ,cell_idx0,X) # find the cell in which X resides
            SFM = (t=T, φ=SFM0.φ, X=X, n=SFM0.n)
            Y = UpdateYt(ffq, SFM, SFM0, Y0, cell_idx, cell_idx0)
        end
        return (Ind = Ind, SFM = SFM, Y=Y)
    end
    return fixed_time_yFun
end

"""
Finds zero of `f` using the bisection method on the interval `[a,b]` with
error `err`.

    fzero( f::Function, a::Real, b::Real; err::Float64 = 1e-8)

"""
function fzero( f::Function, a::Real, b::Real; err::Float64 = 1e-6)
    # finds zeros of f using the bisection method
    c = a + (b - a) / 2.0
    while a < c < b
        fc = f(c)
        if (abs(fc) < err)&&((b-a)<err)
            break
        end
        p = f(a) * fc
        if p < 0.0
            a, b = a, c
        elseif p > 0.0
            a, b = c, b
        elseif fc == 0.0
            break
        else
            a = a+sqrt(eps())
        end
        c = a + (b - a) / 2.0
    end
    return c
end

function brute_force_zero(f::Function, a::Real, b::Real; err::Float64 = 1e-6)
    x = range(a,b;step=err)
    fx = abs.(f.(x))
    val, ind = findmin(fx)
    (val>err)&&(@warn "no values with |f(x)| sufficiently small, returning the closest: |f(x)| = f("*string(x[ind])*") = "*string(fx[ind]))
    return x[ind]
end


text:
using Documenter
using Literate
using RigidBodyDynamics, RigidBodyDynamics.OdeIntegrators

gendir = joinpath(@__DIR__, "src", "generated")
tutorialpages = String[]

let
    rm(gendir, recursive=true, force=true)
    mkdir(gendir)
    exampledir = joinpath(@__DIR__, "..", "examples")
    excludedirs = String[]
    excludefiles = String[]
    if VERSION < v"1.1.0"
        push!(excludefiles, "Symbolic double pendulum.jl")
    end
    for subdir in readdir(exampledir)
        subdir in excludedirs && continue
        root = joinpath(exampledir, subdir)
        isdir(root) || continue
        for file in readdir(root)
            file in excludefiles && continue
            name, ext = splitext(file)
            lowercase(ext) == ".jl" || continue
            outputdir = joinpath(gendir, subdir)
            cp(root, outputdir)
            preprocess = function (str)
                str = replace(str, "@__DIR__" => "\"$(relpath(root, outputdir))\"")
                str = replace(str, "# PREAMBLE" =>
                """
                # This example is also available as a Jupyter notebook that can be run locally
                # The notebook can be found in the `examples` directory of the package.
                # If the notebooks are missing, you may need to run `using Pkg; Pkg.build()`.
                """)
                return str
            end
            postprocess = function(str)
                str = replace(str, "PKG_SETUP" =>
                """
                ```@setup $name
                import Pkg
                let
                    original_stdout = stdout
                    out_rd, out_wr = redirect_stdout()
                    @async read(out_rd, String)
                    try
                        Pkg.activate("$(relpath(root, outputdir))")
                        Pkg.instantiate()
                    finally
                        redirect_stdout(original_stdout)
                        close(out_wr)
                    end
                end
                ```
                """)
                return str
            end

            absfile = joinpath(root, file)
            tutorialpage = Literate.markdown(absfile, outputdir; preprocess=preprocess, postprocess=postprocess)
            push!(tutorialpages, relpath(tutorialpage, joinpath(@__DIR__, "src")))
        end
    end
end

makedocs(
    modules = [RigidBodyDynamics, RigidBodyDynamics.OdeIntegrators],
    root = @__DIR__,
    checkdocs = :exports,
    sitename ="RigidBodyDynamics.jl",
    authors = "Twan Koolen and contributors.",
    pages = [
        "Home" => "index.md",
        "Tutorials" => tutorialpages,
        "Library" => [
            "Spatial vector algebra" => "spatial.md",
            "Joints" => "joints.md",
            "Rigid bodies" => "rigidbody.md",
            "Mechanism" => "mechanism.md",
            "MechanismState" => "mechanismstate.md",
            "Kinematics/dynamics algorithms" => "algorithms.md",
            "Custom collection types" => "customcollections.md",
            "Cache types" => "caches.md",
            "Simulation" => "simulation.md",
            "URDF parsing and writing" => "urdf.md",
        ],
        "Benchmarks" => "benchmarks.md",
      ],
    format = Documenter.HTML(prettyurls = parse(Bool, get(ENV, "CI", "false")))
)

deploydocs(
    repo = "github.com/JuliaRobotics/RigidBodyDynamics.jl.git"
)


text:
module PowerSimulations

#################################################################################
# Exports

#Core Exports
export PowerSimulationsModel
export PowerResults

#Base Modeling Exports
export CustomModel
export EconomicDispatch
export UnitCommitment

#Device formulation Export

#Network Relevant Exports
export AbstractDCPowerModel
export CopperPlatePowerModel

#Functions
export buildmodel!
export simulatemodel

#################################################################################
# Imports
using JuMP
using TimeSeries
using PowerSystems
import PowerModels
import InfrastructureModels
using GLPK
using MathOptInterface
#using Clp
#using Cbc
using Ipopt
using DataFrames
using LinearAlgebra
using LinearAlgebra.BLAS
using AxisArrays
using Dates

#################################################################################
# Type Alias From other Packages
const PM = PowerModels
const IM = InfrastructureModels
const NetworkModel = PM.AbstractPowerFormulation
const PS = PowerSimulations
const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities
const PTDFArray = AxisArrays.AxisArray{Float64,2,Array{Float64,2},Tuple{AxisArrays.Axis{:branches,Array{String,1}},AxisArrays.Axis{:buses,Array{String,1}}}}

#Type Alias for JuMP containers
const JumpVariable = JuMP.JuMPArray{JuMP.VariableRef,2,Tuple{Array{String,1},UnitRange{Int64}}}
const JumpExpressionMatrix = Matrix{<:JuMP.GenericAffExpr}
const JumpAffineExpressionArray = Array{JuMP.GenericAffExpr{Float64,JuMP.VariableRef},2}
const BalanceNamedTuple = NamedTuple{(:var_active, :var_reactive, :timeseries_active, :timeseries_reactive),Tuple{JumpAffineExpressionArray, UJ, Array{Float64,2}, UF}} where {UJ <: Union{Nothing,JumpAffineExpressionArray}, UF <: Union{Nothing, Array{Float64,2}}}

#Type Alias for Unions
const fix_resource = Union{PowerSystems.RenewableFix, PowerSystems.HydroFix}

#################################################################################
# Includes

#Abstract Models
include("network_models/networks.jl")
include("service_models/services.jl")

#base and core
include("core/abstract_models.jl")
#include("core/dynamic_model.jl")
include("base/instantiate_routines.jl")
include("base/model_constructors.jl")
#include("base/solve_routines.jl")
include("base/simulation_routines.jl")

#utils
include("utils/undef_check.jl")
include("utils/cost_addition.jl")
include("utils/timeseries_injections.jl")
include("utils/device_injections.jl")
include("utils/device_retreval.jl")
include("utils/pm_translator.jl")

#Device Modeling components
include("device_models/common.jl")
include("device_models/renewable_generation.jl")
include("device_models/thermal_generation.jl")
include("device_models/electric_loads.jl")
include("device_models/branches.jl")
#include("device_models/storage.jl")
#include("device_models/hydro_generation.jl")
include("service_models/reserves.jl")

#Network related components
include("network_models/copperplate_balance.jl")
include("network_models/nodal_balance.jl")
include("network_models/powermodels_balance.jl")

#Device constructors
include("component_constructors/thermalgeneration_constructor.jl")
include("component_constructors/branch_constructor.jl")
include("component_constructors/renewablegeneration_constructor.jl")
include("component_constructors/load_constructor.jl")
include("component_constructors/services_constructor.jl")

#Network constructors
include("component_constructors/network_constructor.jl")

#PowerModels
#include("power_models/economic_dispatch.jl")

#Utils

end


text:
export NMU

"""
    NMU(in::Int, out::Int; init=rand)

Neural multiplication unit. Can represent multiplications between inputs.
Weights are clipped to [0,1].

As introduced in in https://openreview.net/pdf?id=H1gNOeHKPS
"""
struct NMU{T<:AbstractMatrix}
    W::T
end

NMU(in::Int, out::Int; init=rand) = NMU(init(out,in))

weights(m::NMU) = min.(max.(m.W, 0), 1)

function (m::NMU)(x::AbstractMatrix)
    W  = weights(m)
    Wr = reshape(W, size(W)..., 1)
    xr = reshape(x, 1, size(x)...)
    z  = Wr .* xr .+ 1 .- W
    dropdims(prod(z, dims=2), dims=2)
end

function (m::NMU)(x::AbstractVector)
    z = m(reshape(x,:,1))
    dropdims(z,dims=2)
end

Flux.@functor NMU

Base.show(io::IO, l::NMU) = print(io,"NMU(in=$(size(l.W,2)), out=$(size(l.W,1)))")


text:
using LogicCircuits
using LinearAlgebra
using Plots
import Pkg
Pkg.activate(".")
using ProbabilisticCircuits

const n_exp = 20        # Number of experiments
const ess = 1.0         # Equivalent sample size
const ess_max = 100.0   # Maximum perturbation for SA
const shuffling = false # Shuffling data

lik_spn = zeros(n_exp)
likl_cspn = zeros(n_exp)
liku_cspn = zeros(n_exp)
lik_psdd = zeros(n_exp)
likl_csdd = zeros(n_exp)
liku_csdd = zeros(n_exp)

# https://github.com/UCLA-StarAI/Density-Estimation-Datasets
for db in ["nltcs"] #accidents"] #"nltcs"] #twenty_dataset_names #

    @show db

    # Training data set
    dataSet=dataset(twenty_datasets(db); do_shuffle=shuffling, batch_size=-1)
    dataTrain = dataSet.train

    # Learning PSDD and CSDD
    psdd = learn_struct_prob_circuit(dataTrain; pseudocount = ess)[1]
    csdd = learn_struct_credal_circuit(WXData(dataTrain), ess)[1]

    # Parsing the number of variables in the db
    matrixValid = convert.(Int8,dataSet.valid.x) # Bool2Int
    n = size(matrixValid)[2] # Number of variables in the db

    # Computing and plotting inferences
    x = 1:n
    p = log_proba(psdd, PlainXData(convert.(Int8, -ones(n,n) + 2*I)))
    l = log_marginal_lower(csdd, PlainXData(convert.(Int8, -ones(n,n) + 2*I)))
    u = log_marginal_upper(csdd, PlainXData(convert.(Int8, -ones(n,n) + 2*I)))
    f = log.(sum(matrixValid,dims=1)./size(matrixValid)[1])[1,:]
    plot(x,exp.(p), seriestype = :scatter, label="PSDD", lw=1)
    plot!(x,exp.(u), seriestype = :scatter, label="CSDD (upper)", lw=1)
    plot!(x,exp.(l), seriestype = :scatter, label="CSDD (lower)", lw=1)
    plot!(x,exp.(f), seriestype = :scatter, label="frequence", lw=1)
    savefig("./examples/marginals-$db.png")

    # Preparing data set
    dataValid = XData(matrixValid)
    matrixTrain_ = convert.(Int8, dataSet.train.x)

    # Loop over n_exp experiments with increasing training set size
    for k = 1:n_exp

        # Fist k slices (out of n_exp) of the training data
        matrixTrain_k = matrixTrain_[1:trunc(Int, (num_examples(dataTrain)/n_exp)*k),:]
        dataTrain_k= PlainXData(convert.(Bool,matrixTrain_k))

        # Learning the circuits (PSDD, SPN and their credal counterparts)
        psdd_k = learn_struct_prob_circuit(dataTrain_k; pseudocount = ess)[1]
        csdd_k =learn_struct_credal_circuit(WXData(dataTrain_k), ess)[1]
        spn_k = learn_probabilistic_circuit(dataTrain_k; pseudocount = ess)
        cspn_k =learn_credal_circuit(WXData(dataTrain_k), ess)

        # Computing db likelihood for each model
        likl_csdd[k] = sum(log_marginal_lower(csdd_k, dataValid))
        liku_csdd[k] = sum(log_marginal_upper(csdd_k, dataValid))
        lik_psdd[k] = sum(log_proba(psdd_k, dataValid))
        likl_cspn[k] = sum(log_marginal_lower(cspn_k, dataValid))
        liku_cspn[k] = sum(log_marginal_upper(cspn_k, dataValid))
        lik_spn[k] = sum(log_proba(spn_k, dataValid))

    end

    # Plotting the likelihoods
    x = 1:n_exp
    plot(x, lik_psdd, xticks=0:n_exp, ribbon =(lik_psdd-likl_csdd,liku_csdd-lik_psdd), title=uppercase(db), xlabel="Training set size",
    label ="PSDD", fillalpha=0.5, ylabel="Log-likelihood", lw=1)
    plot!(x, lik_spn,ribbon =(lik_spn-likl_cspn,liku_cspn-lik_spn),fillalpha=0.5,label ="SPN", lw=1)
    savefig("./examples/loglik-$db.png")

    # Learning credal circuits with increasing ESS
    for k = 1:n_exp
        ess2 = ess_max * k/n_exp
        csdd_k = learn_struct_credal_circuit(WXData(dataTrain), ess2)[1]
        cspn_k =learn_credal_circuit(WXData(dataTrain), ess2)
        likl_csdd[k] = sum(log_marginal_lower(csdd_k, dataValid))
        liku_csdd[k] = sum(log_marginal_upper(csdd_k, dataValid))
        likl_cspn[k] = sum(log_marginal_lower(cspn_k, dataValid))
        liku_cspn[k] = sum(log_marginal_upper(cspn_k, dataValid))
    end

    # Plotting
    plot(x, likl_csdd, xticks=0:n_exp, title=uppercase(db), xlabel="Perturbation level (s)",
    label ="CSDD (lower)", ylabel="Log-likelihood", lw=1)
    plot!(x, liku_csdd, xticks=0:n_exp, title=uppercase(db),label ="CSDD (upper)", lw=1)
    plot!(x, likl_cspn, xticks=0:n_exp, title=uppercase(db),label ="CSPN (lower)", lw=1)
    plot!(x, liku_cspn, xticks=0:n_exp, title=uppercase(db),label ="CSPN (upper)", lw=1)
    savefig("./examples/sensitivity-$db.png")

end

max_stars_repo_name: IDSIA/cjuice

max_issues_repo_name: IDSIA/cjuice

max_forks_repo_name: IDSIA/cjuice


text:
using Test

@testset "conversion" begin
    @test convert(Point, (2, 3)) === Point(2, 3)
    @test convert(Point, (2.0, 3)) === Point(2.0, 3.0)
end

@testset "broadcast" begin
    @testset for T in (Vec, Point)
        x = [T(2, 3), T(7, 3)]

        @test [T(4, 9), T(14, 9)] == x .* T(2, 3)
        @test [T(4, 6), T(9, 6)] == x .+ T(2, 3)
        @test [T(0, 0), T(5, 0)] == x .- T(2, 3)
    end
end


text:
module Kohonen

using Distributions
using Gadfly
using RDatasets

typealias Value Float64
typealias Class Matrix{Value}

################################################################################
##  Datatypes
################################################################################

abstract Metric
abstract Manhattan <: Metric
abstract Chebyshev <: Metric
abstract Linear <: Metric

type Map{Metric}
    nodes :: Matrix{Float64}
    dimension :: Tuple
end

################################################################################
##  Plotting Map
################################################################################

function drawOrthoganal(nodes)
    layers = []

    # vertical lines
    for i in eachindex(nodes[1,:,1])
        line = layer(x = nodes[1,i,:], y = nodes[2,i,:], Geom.path,
            Theme(default_color = colorant"lightgrey", line_width = 0.5px))
        push!(layers, line)
    end

    # horizontal lines
    for i in eachindex(nodes[1,1,:])
        line = layer(x = nodes[1,:,i], y = nodes[2,:,i], Geom.path,
            Theme(default_color = colorant"lightgrey", line_width = 0.5px))
        push!(layers, line)
    end

    return layers
end
function drawDiagonal(nodes)
    layers = []

    # diagnal
    height = size(nodes)[2]
    width = size(nodes)[3]
    for h in 2:(height + width - 2)
        stripe = [ (h - w + 1, w) for w in 1:width ]
        stripe = filter((p) -> p[1] > 0 && p[1] <= height, stripe)
        # @show stripe
        xs = map((p) -> nodes[1,p[1],p[2]], stripe)
        ys = map((p) -> nodes[2,p[1],p[2]], stripe)
        line = layer(x = xs, y = ys, Geom.path,
            Theme(default_color = colorant"lightgrey", line_width = 0.5px))
        push!(layers, line)
    end

    for h in (1 - width):(height - 3)
        stripe = [ (h + w + 1, w) for w in 1:width ]
        stripe = filter((p) -> p[1] > 0 && p[1] <= height, stripe)
        # @show stripe
        xs = map((p) -> nodes[1,p[1],p[2]], stripe)
        ys = map((p) -> nodes[2,p[1],p[2]], stripe)
        line = layer(x = xs, y = ys, Geom.path,
            Theme(default_color = colorant"lightgrey", line_width = 0.5px))
        push!(layers, line)
    end
    return layers
end

function plotMap(m :: Map{Linear}, n = 10000, input = [])
    stack = []
    if length(input) > 0
        push!(stack, layer(x = input[1,:], y = input[2,:], Geom.point,
            Theme(default_color = colorant"coral", default_point_size = 1px)))
    end
    push!(stack, layer(x = m.nodes[1,:], y = m.nodes[2,:], Geom.point))
    push!(stack, layer(x = m.nodes[1,:], y = m.nodes[2,:], Geom.path,
        Theme(default_color = colorant"lightgrey", line_width = 0.5px)))
    return plot(stack...,
        Theme(background_color = colorant"white"),
        Coord.cartesian(fixed = true),
        Scale.x_continuous(minvalue = 0, maxvalue = 1),
        Scale.y_continuous(minvalue = 0, maxvalue = 1),
        Guide.title("n = $n"))
end

function plotMap(m :: Map{Manhattan}, n = 10000, input = [])
    nodes = reshape(m.nodes, (2, m.dimension[1], m.dimension[2]))

    stack = []
    if length(input) > 0
        push!(stack, layer(x = input[1,:], y = input[2,:], Geom.point,
            Theme(default_color = colorant"coral", default_point_size = 1px)))
    end
    push!(stack, layer(x = nodes[1,:,:], y = nodes[2,:,:], Geom.point))
    stack = [stack; drawOrthoganal(nodes)]
    return plot(stack...,
        Theme(background_color = colorant"white"),
        Coord.cartesian(fixed = true),
        Scale.x_continuous(minvalue = 0, maxvalue = 1),
        Scale.y_continuous(minvalue = 0, maxvalue = 1),
        Guide.title("n = $n"))
end

function plotMap(m :: Map{Chebyshev}, n = 10000, input = [])
    nodes = reshape(m.nodes, (2, m.dimension[1], m.dimension[2]))
    stack = []
    if length(input) > 0
        push!(stack, layer(x = input[1,:], y = input[2,:], Geom.point,
            Theme(default_color = colorant"coral", default_point_size = 1px)))
    end
    push!(stack, layer(x = nodes[1,:,:], y = nodes[2,:,:], Geom.point))
    stack = [stack; drawOrthoganal(nodes)]
    stack = [stack; drawDiagonal(nodes)]
    return plot(stack...,
        Theme(background_color = colorant"white"),
        Coord.cartesian(fixed = true),
        Scale.x_continuous(minvalue = 0, maxvalue = 1),
        Scale.y_continuous(minvalue = 0, maxvalue = 1),
        Guide.title("n = $n"))
end

################################################################################
##  Metrics
################################################################################

function metric(m :: Map{Linear}, i, j)
    return j - i
end

function metric(m :: Map{Manhattan}, i, j)
    width = m.dimension[2]
    return abs(div(i,width) - div(j,width)) + abs(i%width - j%width)
end

function metric(m :: Map{Chebyshev}, i, j)
    width = m.dimension[2]
    return max(abs(div(i,width) - div(j,width)), abs(i%width - j%width))
end

################################################################################
##  Best Matching Unit
################################################################################

distanceSq = (w, v) -> ((w - v)' * (w - v))[1]

function findBMU(m, v)
    leastDist = distanceSq(m.nodes[:,1,1], v)
    leastNode = 1
    for i in eachindex(m.nodes[1,:])
        dist = distanceSq(m.nodes[:,i], v)
        if dist < leastDist
            leastDist = dist
            leastNode = i
        end
    end
    return leastNode
end

################################################################################
##  Update the Map
################################################################################

radius(m::Map{Linear}) = div(m.dimension[1], 2)
radius(m::Map{Manhattan}) = div(max(m.dimension[1], m.dimension[2]), 2)
radius(m::Map{Chebyshev}) = div(max(m.dimension[1], m.dimension[2]), 2)

function update(m, vector, time, totalTime, R, G)
    nodes = m.nodes
    bmu = findBMU(m, vector)
    # decreasing influence radius
    influenceRadius = radius(m) * exp(-time/R)

    if influenceRadius >= 1
        for i in eachindex(nodes[1,:])
            dist = metric(m, bmu, i)
            if dist <= influenceRadius
                gain = exp(-time/G) * exp(-dist/G)
                nodes[:,i] = nodes[:,i] + gain * (vector - nodes[:,i])
            end
        end
    else    # adjust BMU only
        # rate = L * exp(-time/totalTime)
        nodes[:,bmu] = nodes[:,bmu] + exp(-time/G) * (vector - nodes[:,bmu])
    end
end

function trainAndPlot(timestamp, input, m, R, G)
    stack = []
    total = div(length(input), 2)
    # at 0
    push!(stack, plotMap(m, 0, input))

    current = 0
    for t in timestamp
        iterations = t - current
        for i in 1:iterations
            update(m, input[:,i + current], i + current, total, R, G)
        end
        push!(stack, plotMap(m, t))
        current = t
    end
    return stack
end

l = Map{Linear}(rand(2, 100), (100,))
lStack = trainAndPlot([25, 100, 500, 1000, 5000], rand(2, 5000), l, 100, 5000)
draw(PNG("linear.png", 1200px, 1800px), vstack(hstack(lStack[1:2]...), hstack(lStack[3:4]...), hstack(lStack[5:6]...)))

m = Map{Manhattan}(rand(2, 10 * 10), (10, 10))
mStack = trainAndPlot([25, 100, 500, 1000, 5000], rand(2, 5000), m, 3000, 1000)
draw(PNG("manhattan.png", 1200px, 1800px), vstack(hstack(mStack[1:2]...), hstack(mStack[3:4]...), hstack(mStack[5:6]...)))

c = Map{Chebyshev}(rand(2, 10 * 10), (10, 10))
cStack = trainAndPlot([25, 100, 500, 1000, 5000], rand(2, 5000), c, 3000, 1000)
draw(PNG("chebyshev.png", 1200px, 1800px), vstack(hstack(cStack[1:2]...), hstack(cStack[3:4]...), hstack(cStack[5:6]...)))

################################################################################
##  Plotting Cities
################################################################################

function ringMetric(tour, i, j)
    tourLength = div(length(tour), 3)
    return min((i - j + tourLength)%tourLength, (j - i + tourLength)%tourLength)
end

function plotTour(cities, tour)
    path = layer(x = tour[1,:], y = tour[2,:], Geom.path,
        Theme(default_color = colorant"lightgrey", line_width = 2px))
    bridge = layer(x = [tour[1,1], tour[1, end]],
        y = [tour[2,1], tour[2, end]], Geom.path,
        Theme(default_color = colorant"lightgrey", line_width = 2px))
    vertices = layer(x = tour[1,:], y = tour[2,:], Geom.point,
            Theme(default_point_size = 3px))
    cities = layer(x = cities[1,:], y = cities[2,:], Geom.point,
            Theme(default_color = colorant"coral", default_point_size = 5px))
    return plot([vertices; cities; path; bridge],
        Theme(background_color = colorant"white"),
        Coord.cartesian(fixed = true),
        Scale.x_continuous(minvalue = 0, maxvalue = 1),
        Scale.y_continuous(minvalue = 0, maxvalue = 1),
        Guide.title("TSP"))
end

function findBMV(tour, city)
    leastDist = distanceSq(tour[1:2,1], city)
    leastNode = 1
    for i in eachindex(tour[1,:])
        dist = distanceSq(tour[1:2,i], city)
        if dist < leastDist
            leastDist = dist
            leastNode = i
        end
    end
    return leastNode
end

function insertVertex(tour, i, newVertex)
    tourLength = div(length(tour), 3)
    firstHalf = tour[:,1:i]
    secondHalf = tour[:,i+1:tourLength]
    return [firstHalf newVertex secondHalf]
end

# 1 iteration
function updateTour(cities, tour, iteration, G, R)
    for c in eachindex(cities[1,:])
        city = cities[:,c]
        bmv = findBMV(tour, city)

        if tour[3,bmv] != iteration     # not visited
            # update the vertex and it's neighbors
            for i in eachindex(tour[1,:])
                dist = ringMetric(tour, bmv, i)
                gain = 0.707106781 * exp(-(dist^2)/G^2)
                tour[1:2,i] = tour[1:2,i] + gain * (city - tour[1:2,i])
            end

            # mark the best matching vertex visited
            tour[3,bmv] = iteration
        else
            # copy the best matching vertex, in case we need to insert it into the
            # tour as a new vertex
            newVertex = [tour[1:2,bmv]; iteration - 1]

            # update the vertex and it's neighbors
            for i in eachindex(tour[1,:])
                dist = ringMetric(tour, bmv, i)
                gain = 0.707106781 * exp(-(dist^2)/G^2)
                tour[1:2,i] = tour[1:2,i] + gain * (city - tour[1:2,i])
            end

            # mark the best matching vertex visited
            tour[3,bmv] = iteration

            tour = insertVertex(tour, bmv, newVertex)
        end
    end

    # delete old nodes
    suvivorIndices = []
    for v in eachindex(tour[1,:])
        if iteration - tour[3,v] < 3
            push!(suvivorIndices, v)
        end
    end
    tour = tour[:,suvivorIndices]
    return tour
end

function tourLength(tour)
    l = sqrt(distanceSq(tour[1:2,end], tour[1:2,1]))
    for v in eachindex(tour[1,1:end-1])
        l = l + sqrt(distanceSq(tour[1:2,v], tour[1:2,v+1]))
    end
    return l
end

function trainAndPlotHistogram(rate, iteration)
    cities = rand(2, 30)
    tour = [0.0, 0.0, 0]
    lengths = []
    # begin with 1 node at (0, 0), visited 0 times
    for i in 1:iteration
        if i % 100 == 0
            @show i
        end
        tour = [0.0, 0.0, 0]
        permCities = cities[:,randcycle(30)]
        G = 1
        R = 1/rate
        for j in 1:(rate + 5)
            tour = updateTour(permCities, tour, j, G, R)
            G = (1 - R) * G
        end
        push!(lengths,tourLength(tour))
    end

    minL = round(minimum(lengths), 3)
    maxL = round(maximum(lengths), 3)
    meanL = round(mean(lengths), 3)
    histo = plot(x=map((x) -> floor(x * 5)/5 , lengths), Geom.histogram,
        Theme(background_color = colorant"white"),
        Guide.title("$iteration tries, α = $(1/rate)\naverage: $meanL, lower bould: $minL, upper bould: $maxL"))
    draw(PNG("tour-$rate.png", 600px, 600px), plotTour(cities, tour))
    draw(PNG("histo-$rate.png", 600px, 300px), histo)
end

trainAndPlotHistogram(1000, 1000)

end

max_stars_repo_name: banacorn/neuron

max_issues_repo_name: banacorn/neuron

max_forks_repo_name: banacorn/neuron


text:
import Pkg

Pkg.activate(temp=true)

Pkg.add(Pkg.PackageSpec(name="Plots", version="1.0.0"))

Pkg.status()


text:
module RulesBIFMMarginalsTest

using Test
using ReactiveMP
import ReactiveMP: @test_marginalrules

@testset "marginalrules:BIFM" begin

    @testset ":in_zprev_znext (m_out::MultivariateNormalDistributionsFamily, m_in::MultivariateNormalDistributionsFamily, m_zprev::ProdFinal{<:MultivariateNormalDistributionsFamily}, m_znext::MultivariateNormalDistributionsFamily)" begin

        meta = BIFMMeta([2.0 0; 0 1], # A
                        [3.0 0; 0 2], # B
                        [4.0 0; 0 3], # C
                        [5.0 0; 0 4], # H
                        [6.0 0; 0 5], # BHBt
                        [1.0, 2],     # ξz
                        [7.0 0; 0 6], # Λz
                        [2.0, 3],     # ξztilde
                        [8.0 0; 0 7], # Λztilde
                        [3.0, 4],     # μu
                        [9.0 0; 0 8]  # Σu
                        )

        @test_marginalrules [ with_float_conversions = true, float32_atol = 1e-2, bigfloat_atol = 1e-8 ] BIFM(:in_zprev_znext) [
            (
                input  = (m_out = MvNormalMeanPrecision([1,2], [2 0; 0 1]),     m_in = MvNormalMeanPrecision([1,2], [1 0; 0 2]), m_zprev = ProdFinal(MvNormalMeanPrecision([1,2], [1 0; 0 2])), m_znext = MvNormalMeanPrecision([1,2], [1 0; 0 3]), meta = meta),
                output = MvNormalWeightedMeanPrecision([28.0, 28.0, 15.0, 13.0], [298.0 0.0 198.0 0.0; 0.0 50.0 0.0 24.0; 198.0 0.0 101.0 0.0; 0.0 24.0 0.0 7.0])
            ),
            (
                input  = (m_out = MvNormalMeanPrecision([3,4], [6 0; 0 1]),     m_in = MvNormalMeanPrecision([8,2], [2 0; 0 2]), m_zprev = ProdFinal(MvNormalMeanPrecision([1,6], [2 0; 0 2])), m_znext = MvNormalMeanPrecision([1,2], [1 0; 0 3]), meta = meta),
                output = MvNormalWeightedMeanPrecision([235.0, 40.0, 144.0, 27.0], [875.0 0.0 582.0 0.0; 0.0 50.0 0.0 24.0; 582.0 0.0 358.0 0.0; 0.0 24.0 0.0 7.0])
            ),
            (
                input  = (m_out = MvNormalMeanPrecision([5,6], [2 0; 0 5]),     m_in = MvNormalMeanPrecision([1,9], [1 0; 0 1]), m_zprev = ProdFinal(MvNormalMeanPrecision([6,2], [1 0; 0 1])), m_znext = MvNormalMeanPrecision([1,2], [1 0; 0 3]), meta = meta),
                output = MvNormalWeightedMeanPrecision([124.0, 201.0, 84.0, 95.0], [298.0 0.0 198.0 0.0; 0.0 193.0 0.0 96.0; 198.0 0.0 101.0 0.0; 0.0 96.0 0.0 42.0])
            ),

        ]

    end
end

end


text:
module testExamples
using LinearAlgebra
using Test

srcDir = dirname(dirname(@__FILE__))*"\\src\\Examples\\"
include(srcDir *"Examples_sep.jl")
include(srcDir *"Examples_ent.jl")
include(srcDir *"Examples.jl")

using .Examples_sep
using .Examples_ent
using .Examples

@testset "Test the separable examples" begin
    ρ_sep  = Examples_sep.get_sep_example()
end

@testset "Test the entangeled examples" begin
    ρ_ent  = Examples_ent.get_ent_example()
end

@testset "Test the  example tools" begin
    ρ = Examples.get_examples()
    dfρ = Examples.get_example_overview()
    println(dfρ)
end

end

max_stars_repo_name: JAndriesJ/sep-rank-bounding

max_issues_repo_name: JAndriesJ/sep-rank-bounding

max_forks_repo_name: JAndriesJ/sep-rank-bounding


text:
# A counter type

type Accumulator{T, V<:Number}
	map::Dict{T,V}
end

## constructors

Accumulator{T,V<:Number}(::Type{T}, ::Type{V}) = Accumulator{T,V}((T=>V)[])
counter(T::Type) = Accumulator(T,Int)

Accumulator{T,V<:Number}(dct::Dict{T,V}) = Accumulator{T,V}(copy(dct))
counter{T}(dct::Dict{T,Int}) = Accumulator{T,Int}(copy(dct))

function counter{T}(seq::AbstractArray{T})
	ct = counter(T)
	for x in seq
		push!(ct, x)
	end
	return ct
end

copy{T,V<:Number}(ct::Accumulator{T,V}) = Accumulator{T,V}(copy(ct.map))

length(a::Accumulator) = length(a.map)

## retrieval

getindex{T,V}(ct::Accumulator{T,V}, x::T) = get(ct.map, x, zero(V))

haskey{T,V}(ct::Accumulator{T,V}, x::T) = haskey(ct.map, x)

keys(ct::Accumulator) = keys(ct.map)

## iteration

start(ct::Accumulator) = start(ct.map)
next(ct::Accumulator, state) = next(ct.map, state)
done(ct::Accumulator, state) = done(ct.map, state)

# manipulation

push!{T,V<:Number}(ct::Accumulator{T,V}, x::T, a::V) = (ct.map[x] = ct[x] + a)
push!{T,V<:Number,V2<:Number}(ct::Accumulator{T,V}, x::T, a::V2) = push!(ct, x, convert(V,a))
push!{T,V<:Number}(ct::Accumulator{T,V}, x::T) = push!(ct, x, one(V))

function push!{T,V<:Number,V2<:Number}(ct::Accumulator{T,V}, r::Accumulator{T,V2})
	for (x::T, v::V2) in r
		push!(ct, x, v)
	end
	ct
end

pop!{T,V<:Number}(ct::Accumulator{T,V}, x::T) = pop!(ct.map, x)

merge{T,V<:Number}(ct1::Accumulator{T,V}, ct2::Accumulator{T,V}) = push!(copy(ct1), ct2)