text:
"""
    FilteringProblem(state_model, obs_model)

Specify a filtering problem for a hidden state model `state_model` and observation model `obs_model`.

Constraints: `state_model` and `obs_model` must be chosen such that
* `state_model` is of type `HiddenStateModel{S1, T1}`,
* `obs_model` is of type `ObservationModel{S1, S2, T2}`,
where `T1<:TimeType`, `T2 <:TimeType`, and `S1`, `S2` are arbitrary.
This ensures that the filtering problem can be evaluated, i.e. the hidden states are of appropriate type to generate observations.
"""
struct FilteringProblem{S1, S2, T1, T2, M1, M2} <: AbstractFilteringProblem{S1, S2, T1, T2}
    state_model::M1
    obs_model::M2
    function FilteringProblem(mod1::HiddenStateModel{S1, T1}, mod2::ObservationModel{S1, S2, T2}) where {S1, S2, T1<:TimeType, T2<:TimeType}
        return new{S1, S2, T1, T2, typeof(mod1), typeof(mod2)}(mod1, mod2)
    end
end

# mandatory methods
state_model(problem::FilteringProblem)      = problem.state_model
obs_model(problem::FilteringProblem)        = problem.obs_model

# pretty printing
function Base.show(io::IO, ::MIME"text/plain", problem::FilteringProblem{S1, S2, T1, T2, M1, M2}) where {S1, S2, T1, T2, M1, M2}
    print(io, T1, " - ", T2," filtering problem
    hidden state model:                     ", problem.state_model,"
    observation model:                      ", problem.obs_model,"
    hidden state type:                      ", S1,"
    observation type:                       ", S2)
end

function Base.show(io::IO, problem::FilteringProblem{S1, S2, T1, T2, M1, M2}) where {S1, S2, T1, T2, M1, M2}
    print(io, T1, " - ", T2," with ", M1, " and ", M2)
end

max_stars_repo_name: simsurace/FeedbackParticleFilters

max_issues_repo_name: simsurace/FeedbackParticleFilters

max_forks_repo_name: simsurace/FeedbackParticleFilters


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Paula B. Staudt
#* $Id$
#*--------------------------------------------------------------------
type thermosyphon
	thermosyphon()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of Components"
			)),
			press_delta(Dict{Symbol,Any}(
				:Brief=>"Pressure Drop in the reboiler",
				:Symbol=>"\\Delta P"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Flow Constant"
			)),
			DanaReal(Dict{Symbol,Any}(
				:Brief=>"Flow Constant",
				:Hidden=>true,
				:Unit=>"mol^3/(kg*m^2)"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Liquid inlet stream",
				:PosX=>0.44,
				:PosY=>1,
				:Symbol=>"_{inL}"
			)),
			streamPH(Dict{Symbol,Any}(
				:Brief=>"Vapour outlet stream",
				:PosX=>0,
				:PosY=>0.09,
				:Symbol=>"_{outV}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Heat supplied",
				:PosX=>1,
				:PosY=>0.77,
				:Symbol=>"Q_{in}",
				:Protected=>true
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Temperature  Indicator of Reboiler",
				:Protected=>true,
				:PosX=>1,
				:PosY=>0.57
			)),
			control_signal(Dict{Symbol,Any}(
				:Brief=>"Pressure Indicator of Reboiler",
				:Protected=>true,
				:PosX=>1,
				:PosY=>0.35
			)),
			[
				:(InletLiquid.F = OutletVapour.F),
				:(InletLiquid.z = OutletVapour.z),
				:(InletLiquid.F*InletLiquid.h + InletQ = OutletVapour.F*OutletVapour.h),
				:(OutletVapour.P = InletLiquid.P - Pdrop),
				:(TI * "K" = OutletVapour.T),
				:(PI * "atm" = OutletVapour.P),
				:(OutletVapour.F^3 = FlowConstant*k*InletQ),
			],
			[
				"Molar Flow Balance","Molar Composition Balance","Energy Balance","Pressure Drop","Temperature indicator","Pressure indicator","Flow through the thermosyphon reboiler",
			],
			[:PP,:NComp,:Pdrop,:FlowConstant,:k,],
			[:InletLiquid,:OutletVapour,:InletQ,:TI,:PI,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Pdrop::press_delta
	FlowConstant::DanaReal
	k::DanaReal
	InletLiquid::stream
	OutletVapour::streamPH
	InletQ::power
	TI::control_signal
	PI::control_signal
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export thermosyphon
function set(in::thermosyphon)
	k = 1*"mol^3/(kg*m^2)"

end
function setEquationFlow(in::thermosyphon)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
end
function atributes(in::thermosyphon,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/Thermosyphon"
	fields[:Brief]="Model of a  Steady State reboiler thermosyphon."
	fields[:Info]="== ASSUMPTIONS ==
* perfect mixing of both phases;
* no thermodynamics equilibrium;

== SET ==
* the pressure drop in the reboiler;
* the FlowConstant that relates the Flow through the reboiler and the heat duty
**	Flow^3 = FlowConstant*InletQ

== SPECIFY ==
* the InletLiquid stream;
* the InletQ (the model requires an energy stream, also you can use a controller for setting the heat duty using the heat_flow model)
OR the outlet temperature (OutletVapour.T);

== OPTIONAL ==
* the reboiler model has two control ports
** TI OutletVapour Temperature Indicator;
** PI OutletVapour Pressure Indicator;
"
	drive!(fields,_)
	return fields
end
thermosyphon(_::Dict{Symbol,Any})=begin
	newModel=thermosyphon()
	newModel.attributes=atributes(newModel,_)
	newModel
end

ext: jl


text:
_parenttfm(t::Transform) = t
_parenttfm(t::BufferedThreadsafe) = _parenttfm(t.buffereds[1])
_parenttfm(t::Buffered) = _parenttfm(t.tfm)

_copyrec(t::Tuple) = Tuple(_copyrec(x) for x in t)
_copyrec(x) = copy(x)

_copyrec!(dstt::Tuple, srct::Tuple) = Tuple(_copyrec!(dst, src) for (dst, src) in zip(dstt, srct))
_copyrec!(dst, src) = copy!(dst, src)


text:
# XXX maps their different levels to carriers
function createCarrierMapping!(setData_dic::Dict,anyM::anyModel)
    # determines number of specified carrier levels
    lvlCar_arr = [Symbol("carrier_",i) for i in 1:anyM.sets[:C].height]

    # tuple of columns with dispatch and expansion resolutions
    resCol_tup =  (:timestep_dispatch, :timestep_expansion, :region_dispatch, :region_expansion)
    resLongShort_tup = Dict(:timestep_dispatch => :lvlTsDis, :timestep_expansion => :lvlTsExp, :region_dispatch => :lvlRDis, :region_expansion => :lvlRExp)

    anyM.cInfo = Dict{Int,NamedTuple{(:tsDis,:tsExp,:rDis,:rExp,:eq),Tuple{Int,Int,Int,Int,Bool}}}()

    # loops over rows in carrier file and writes specific resolutions
    for row in eachrow(setData_dic[:C])
    	# gets tuple of strings for respective carrier
    	car_tup = tuple(map(x -> row[x],lvlCar_arr)...)
    	car_int = lookupTupleTree(car_tup,anyM.sets[:C])[1]
    	# gets resolution value and writes them if they can be parsed to numbers
    	resVal_dic = Dict(resLongShort_tup[x] => row[x] for x in resCol_tup)

		# check, if carrier got an equality constraint or not
		if :carrier_equality in namesSym(row)
			if !(row[:carrier_equality] in ("no","yes"))
				push!(anyM.report,(2,"carrier mapping","","column carrier_equality can only contain keywords 'yes' or 'no'"))
				continue
			else
				eq_boo = row[:carrier_equality] == "yes" ? true : false
			end
		else
			eq_boo = false
		end

    	# check if level values can be converted to integers
    	if any(map(x -> tryparse(Int,x), values(resVal_dic)) .== nothing)
    		push!(anyM.report,(2,"carrier mapping","","no resolutions written for $(createFullString(car_int,anyM.sets[:C])), provide as integer, carrier was skipped"))
    		continue
    	end
    	res_dic = Dict(resLongShort_tup[x] => parse(Int,row[x]) for x in resCol_tup)
    	# writes levels after consistency check
    	if res_dic[:lvlTsDis] < res_dic[:lvlTsExp]
    		push!(anyM.report,(3,"carrier mapping","","temporal resolution of expansion can not be more detailed than for dispatch for $(createFullString(car_int,anyM.sets[:C]))"))
    		continue
    	elseif res_dic[:lvlRDis] > res_dic[:lvlRExp]
    		push!(anyM.report,(3,"carrier mapping","","spatial resolution of expansion must be at least as detailed as dispatch for $(createFullString(car_int,anyM.sets[:C]))"))
    		continue
    	else
    		anyM.cInfo[car_int] = (tsDis = res_dic[:lvlTsDis],tsExp = res_dic[:lvlTsExp],rDis = res_dic[:lvlRDis],rExp = res_dic[:lvlRExp], eq = eq_boo)
    	end
    end

	if any(getindex.(anyM.report,1) .== 3) print(getElapsed(anyM.options.startTime)); errorTest(anyM.report,anyM.options) end

	 #  loops over all carriers and check consistency of resolutions and tries to inherit a resolution where none was defined, cannot be carried if above they have already been errors detected
    for c in filter(x -> x != 0, keys(anyM.sets[:C].nodes))
    	anyM.cInfo = evaluateReso(c,anyM.sets[:C],anyM.cInfo,anyM.report)
    end

    if minimum(map(x -> getfield(x,:tsDis),values(anyM.cInfo))) < maximum(map(x -> getfield(x,:tsExp),values(anyM.cInfo)))
    	push!(anyM.report,(3,"carrier mapping","","one temporal dispatch resoultion is more detailed than one of the temporal expansion resolutions"))
    end

    produceMessage(anyM.options,anyM.report, 3," - Created mapping for carriers")
end

# XXX checks carrier for errors in resolution or derive resolution from lower carriers
function evaluateReso(startIdx_int::Int,car_tree::Tree,cInfo_dic::Dict{Int,NamedTuple{(:tsDis,:tsExp,:rDis,:rExp,:eq),Tuple{Int,Int,Int,Int,Bool}}},report::Array{Tuple,1})
	# extracts all children and all columns related to resolution
	carName_str = createFullString(startIdx_int,car_tree)
	allChildIdx_arr = getDescendants(startIdx_int,car_tree)

	# all entries need to have a resoultion => otherwise evaluateReso on them
	for noResoIdx = setdiff(allChildIdx_arr,collect(keys(cInfo_dic)))
		cInfo_dic = evaluateReso(noResoIdx,car_tree,cInfo_dic,report)
	end

	# tries to inherit resolutions from children, if no data exists yet
	if !haskey(cInfo_dic,startIdx_int)
		if isempty(allChildIdx_arr)
			push!(report,(3,"carrier mapping","","carrier $(carName_str) got no resolution and could not inherit from children either"))
			return cInfo_dic
		else
			newReso_dic = Dict(y => minimum([getfield(cInfo_dic[x],y) for x in allChildIdx_arr]) for y in (:tsDis,:tsExp,:rDis,:rExp))
			cInfo_dic[startIdx_int] = (tsDis = newReso_dic[:tsDis],tsExp = newReso_dic[:tsExp],rDis = newReso_dic[:rDis],rExp = newReso_dic[:rExp], eq = false)
			push!(report,(1,"carrier mapping","","carrier $(carName_str) inherited resolution from children"))
			return cInfo_dic
		end
	# checks if existing resolution is flawed
	else
		for childIdx in allChildIdx_arr
			# check if any children got a smaller resolution value
			if any(map(x -> getfield(cInfo_dic[startIdx_int],x) > getfield(cInfo_dic[childIdx],x),(:tsDis,:tsExp,:rDis,:rExp)))
				push!(report,(3,"carrier mapping","","carrier $(carName_str) got a resolution more detailed than its childrens'"))
			end
		end
	end

	return cInfo_dic
end

# XXX maps information about timesteps used
function createTimestepMapping!(anyM::anyModel)
    # XXX writes the superordinate dispatch level, the timesteps on this level and scaling factor for timesteps depending on the respective superordinate dispatch timestep and the level
    supTsLvl_int = maximum(map(x -> getfield(x,:tsExp),values(anyM.cInfo)))

	if anyM.options.supTsLvl != 0
		if minimum(map(x -> getfield(x,:tsDis),values(anyM.cInfo))) >= anyM.options.supTsLvl
			supTsLvl_int = anyM.options.supTsLvl
			push!(anyM.report,(1,"timestep mapping","","superordinate dispatch level provided via options was used"))
		else
			push!(anyM.report,(2,"timestep mapping","","superordinate dispatch level provided via options could not be used, because it was more detailed than at least one dispatch level provided"))
		end
	end

    supTs_tup = tuple(sort(getfield.(filter(x -> x.lvl == supTsLvl_int,collect(values(anyM.sets[:Ts].nodes))),:idx))...)
    scaSupTs_dic = Dict((x[1],x[2]) => (1/anyM.options.redStep)*8760/length(getDescendants(x[1],anyM.sets[:Ts],false,x[2])) for x in Iterators.product(supTs_tup,filter(x -> x >= supTsLvl_int,1:anyM.sets[:Ts].height)))
    anyM.supTs = (lvl = supTsLvl_int, step = supTs_tup, sca = scaSupTs_dic)

    if length(anyM.supTs.step) > 50
		push!(anyM.report,(2,"timestep mapping","","problem specification resulted in more than 50 superordinate dispatch timesteps, this looks faulty"))
	end

	produceMessage(anyM.options,anyM.report, 3," - Created mapping for time steps")
end

# XXX writes basic information for each technology
function createTechInfo!(tSym::Symbol, setData_dic::Dict,anyM::anyModel)

    part = anyM.parts.tech[tSym]
    t_int = techInt(tSym,anyM.sets[:Te])
    lvlTech_arr = Symbol.(:technology_,1:anyM.sets[:Te].height)

    # tuple of columns with input, output and stored carriers
	typeStr_dic = Dict(:carrier_conversion_in => "conversion input", :carrier_conversion_out => "conversion output", :carrier_stored_in => "storage", :carrier_stored_out => "storage")
    carCol_tup =  (:carrier_conversion_in, :carrier_conversion_out, :carrier_stored_in, :carrier_stored_out)

	# maps carrier strings to their id
	nameC_dic = Dict(collect(values(anyM.sets[:C].nodes)) |> (y -> Pair.(getfield.(y,:val),getfield.(y,:idx))))

    # maps selected strings of tech types to integers
    typeStringInt_dic = Dict("stock" => 0, "mature" => 1,"emerging" => 2)

    # gets datarow for respective technology
	row_df = anyM.sets[:Te].nodes[t_int].val |> (z -> filter(x -> any(map(y -> z == x[y],lvlTech_arr)) ,setData_dic[:Te])[1,:])

    # XXX writes carrier info
    # gets string array of carriers for input, output and stored, looks up respective ids afterwards and writes to mapping file
    carStrArr_dic = Dict(y => y in namesSym(row_df) ? split(replace(row_df[y]," " => ""),";") |> (z -> filter(x -> !isempty(x),z)) : String[] for y in carCol_tup)
	carId_dic = Dict(z => tuple(map(x -> getDicEmpty(nameC_dic,x),carStrArr_dic[z])...) for z in keys(carStrArr_dic))

	for x in filter(x -> Int[] in carId_dic[x], collectKeys(keys(carId_dic)))
		push!(anyM.report,(3,"technology mapping","carrier","$(typeStr_dic[x]) carrier of technology $(string(tSym)) not entered correctly"))
		carId_dic[x] = tuple(filter(y -> y != Int[],collect(carId_dic[x]))...)
	end

	# avoid storage of carriers that are balanced on superordinate dispatch level (e.g. if gas is balanced yearly, there is no need for gas storage)
	for type in (:carrier_stored_out, :carrier_stored_in)
		for c in carId_dic[type]
			if anyM.supTs.lvl == anyM.cInfo[c].tsDis
				carId_dic[type] = tuple(filter(x -> x != c,collect(carId_dic[type]))...)
				push!(anyM.report,(2,"technology mapping","carrier","carrier $(createFullString(c,anyM.sets[:C])) of technology $(string(tSym)) cannot be stored, because carrier is balanced on superordinate dispatch level"))
			end
		end
	end

    # writes all relevant type of dispatch variables and respective carrier
    carGrp_ntup = (use = carId_dic[:carrier_conversion_in], gen = carId_dic[:carrier_conversion_out], stExtIn = carId_dic[:carrier_stored_in], stExtOut = carId_dic[:carrier_stored_out],
                         stIntIn = tuple(intersect(carId_dic[:carrier_conversion_out],carId_dic[:carrier_stored_out])...), stIntOut = tuple(intersect(carId_dic[:carrier_conversion_in],carId_dic[:carrier_stored_in])...))

	if :carrier_stored_active in namesSym(row_df)
		actStStr_arr = split(replace(row_df[:carrier_stored_active]," " => ""),";") |> (z -> filter(x -> !isempty(x),z))
		actSt_tup = tuple(map(x -> getDicEmpty(nameC_dic,x),actStStr_arr)...)
	else
		actSt_tup = tuple()
	end
	part.actSt = actSt_tup

	# report on suspicious looking carrier constellations
    if isempty(union(carId_dic[:carrier_conversion_out],carId_dic[:carrier_stored_out])) push!(anyM.report,(2,"technology mapping","carrier","technology $(string(tSym)) has no output")) end

    if !isempty(setdiff(carId_dic[:carrier_stored_in],union(carGrp_ntup.stIntOut,carGrp_ntup.stExtOut))) && !isempty(carId_dic[:carrier_stored_in])
        push!(anyM.report,(2,"technology mapping","carrier","some carrier of technology $(string(tSym)) can be charged but not discharged"))
    end

    if !isempty(setdiff(carId_dic[:carrier_stored_out],union(carGrp_ntup.stIntIn,carGrp_ntup.stExtIn))) && !isempty(carId_dic[:carrier_stored_out])
        push!(anyM.report,(2,"technology mapping","carrier","some carrier of technology $(string(tSym)) can be discharged but not charged"))
    end

	for c in part.actSt
		if !(c in vcat(map(x -> vcat(getDescendants(x,anyM.sets[:C],true)...,x), union(carGrp_ntup.stExtIn,carGrp_ntup.stExtOut))...))
			push!(anyM.report,(3,"technology mapping","carrier","$(createFullString(c,anyM.sets[:C])) for active storage of technology $(string(tSym)) is not stored or a descendant of a stored carrier"))
		end
	end

    part.carrier = filter(x -> getfield(carGrp_ntup,x) != tuple(),collectKeys(keys(carGrp_ntup))) |> (y -> NamedTuple{Tuple(y)}(map(x -> getfield(carGrp_ntup,x), y)) )

    # detects if any in or out carrier is a parent of another in or out carrier, removes carrier in these cases and reports on it
    for type in (:carrier_conversion_in, :carrier_conversion_out)
        relCar_tup = carId_dic[type]
        inherCar_tup = relCar_tup[findall(map(x -> !(isempty(filter(z -> z != x,intersect(getDescendants(x,anyM.sets[:C],true),relCar_tup)))),relCar_tup))]
        if !isempty(inherCar_tup)
            for inher in inherCar_tup
                push!(anyM.report,(3,"technology mapping","carrier","for technology $(string(tSym)) the $(typeStr_dic[type]) carrier $(createFullString(inher,anyM.sets[:C])) is a parent of another $(typeStr_dic[type]) carrier, this is not supported"))
            end
        end
    end

    # XXX writes technology type
    # finds technology type and tries to convert to an integer
	if :technology_type in namesSym(row_df)
    	type_str = row_df[:technology_type]
	else
		type_str = "mature"
	end
    if !haskey(typeStringInt_dic,String(type_str))
        push!(anyM.report,(3,"technology mapping","type","unknown technology type $type_str used, allowed are: $(join(keys(typeStringInt_dic),", "))"))
        return
    end
    part.type = Symbol(type_str)

    # XXX writes modes of technology
    if :mode in namesSym(row_df) && length(anyM.sets[:M].nodes) > 1
        part.modes = tuple(collect(lookupTupleTree(tuple(string(x),),anyM.sets[:M],1)[1] for x in filter(x -> x != "",split(replace(row_df[:mode]," " => ""),";")))...)
    else
        part.modes = tuple()
    end

    # XXX determines resolution of expansion
    # determines carrier based expansion resolutions
	cEx_boo = true
    if isempty(vcat(collect.(values(carGrp_ntup))...))
		push!(anyM.report,(2,"technology mapping","carrier","for technology $(string(tSym)) no carriers were provided"))
		cEx_boo = false
	end

	tsExp_int = cEx_boo ? maximum(map(y -> getfield(anyM.cInfo[y],:tsExp), vcat(collect.(values(carGrp_ntup))...))) : 0
	rExp_int = cEx_boo ? maximum(map(y -> getfield(anyM.cInfo[y],:rExp), vcat(collect.(values(carGrp_ntup))...))) : 0

	# check if carrier based temporal resolution is overwritten by a technology specifc value
	if cEx_boo && :timestep_expansion in namesSym(row_df)
		tsExpSpc_int = tryparse(Int,row_df[:timestep_expansion])

		if !isnothing(tsExpSpc_int)
			if tsExpSpc_int > anyM.supTs.lvl
				push!(anyM.report,(2,"technology mapping","expansion level","specific temporal expansion level provided for $(string(tSym)) is below superordinate dispatch level and therefore could not be used"))
			else
				push!(anyM.report,(1,"technology mapping","expansion level","specific temporal expansion level provided for $(string(tSym)) was used instead of a carrier based value"))
				tsExp_int = tsExpSpc_int
			end
		end
	end

	# check if carrier based spatial resolution is overwritten by a technology specifc value
	if cEx_boo && :region_expansion in namesSym(row_df)
		rExpSpc_int = tryparse(Int,row_df[:region_expansion])

		if !isnothing(rExpSpc_int)
			if rExpSpc_int < rExp_int
				push!(anyM.report,(2,"technology mapping","expansion level","specific spatial expansion level provided for $(string(tSym)) is less detailed than default value obtained from carriers and therefore could not be used"))
			elseif rExpSpc_int == rExp_int
				push!(anyM.report,(1,"technology mapping","expansion level","specific spatial expansion level provided for $(string(tSym)) is equal to default value obtained from carriers"))
			else
				push!(anyM.report,(1,"technology mapping","expansion level","specific spatial expansion level provided for $(string(tSym)) was used instead of a carrier based value"))
				rExp_int = rExpSpc_int
			end
		end
	end

    expLvl_tup = (tsExp_int,rExp_int)

	# XXX checks if dispatch variables should be disaggregated by expansion regions
	rExpOrg_int = cEx_boo ? maximum(map(y -> getfield(anyM.cInfo[y],:rDis), vcat(collect.(values(carGrp_ntup))...))) : 0

	if :region_disaggregate in namesSym(row_df) && rExp_int > rExpOrg_int # relies on information in explicit column, if disaggregation is possible and column exists
		daggR_str = row_df[:region_disaggregate]
		if daggR_str == "yes"
			disAgg_boo = true
		elseif daggR_str == "no"
			disAgg_boo = false
		else
			push!(anyM.report,(3,"technology mapping","spatial aggregation","unknown keyword $type_str used to control spatial aggregation, please use 'yes' or 'no'"))
			return
		end
	elseif rExp_int > rExpOrg_int # disaggregate by default, if it makes sense
		disAgg_boo = true
	else
		disAgg_boo = false
	end
	part.disAgg = disAgg_boo

    # XXX determines reference resolution for conversion (takes into account "region_disaggregate" by using spatial expansion instead of dispatch level if set to yes)
    if !isempty(union(carGrp_ntup.use,carGrp_ntup.gen))
		refTs_int = minimum([minimum([getproperty(anyM.cInfo[x],:tsDis) for x in getproperty(carGrp_ntup,z)]) for z in intersect(keys(part.carrier),(:gen, :use))])
		refR_int = disAgg_boo ? rExp_int : minimum([minimum([getproperty(anyM.cInfo[x],:rDis) for x in getproperty(carGrp_ntup,z)]) for z in intersect(keys(part.carrier),(:gen, :use))])
        refLvl_tup = (refTs_int, refR_int)
    else
        refLvl_tup = nothing
    end

    part.balLvl = (exp = expLvl_tup, ref = refLvl_tup)

	produceMessage(anyM.options,anyM.report, 3," - Created mapping for technology $(string(tSym))")
end

# XXX maps capacity constraints to technology
function createCapaRestrMap!(tSym::Symbol,anyM::anyModel)

    part = anyM.parts.tech[tSym]

    capaDispRestr_arr = Array{Tuple{String,Array{Int,1},Int,Int},1}()
    # extract tech info
    carGrp_ntup = part.carrier
    balLvl_ntup = part.balLvl
    disAgg_boo  = part.disAgg

    # XXX writes dimension of capacity restrictions for conversion part (even if there are no inputs)
    for side in intersect((:use,:gen),keys(carGrp_ntup))
        # get respective carrier and their reference level
        carDis_tup = map(getfield(carGrp_ntup,side)) do x
                carRow_ntup = anyM.cInfo[x]
                return x, carRow_ntup.tsDis, disAgg_boo ? balLvl_ntup.exp[2] : carRow_ntup.rDis
        end

        carConstr_arr = Tuple{Array{Int,1},Int,Int}[]

        # writes all relevant combinations by going from finest resolution up, separately for temporal and spatial (2 and 3)
        for j = [2,3]
            # sorts descinding by j-th column and ascending by other column
            carDisSort_arr = sort(collect(carDis_tup), by = x -> x[j], rev=true)
            carIt_arr =	map(1:length(carDis_tup)) do x
                if j == 2 (sort([carDisSort_arr[y][1] for y in 1:x]), carDisSort_arr[x][2], minimum([carDisSort_arr[y][3] for y in 1:x]))
                else (sort([carDisSort_arr[y][1] for y in 1:x]), minimum([carDisSort_arr[y][2] for y in 1:x]), carDisSort_arr[x][3]) end
            end
            # filters entries that exceed the reference level or are not below the reference level, if already a constraint on the reference level exists from the previous iteration
			if side == :use && isempty(setdiff((:use,:gen),keys(carGrp_ntup)))
				carIt_arr = carIt_arr[findall(x -> j == 2 ? x[2] > balLvl_ntup.ref[1] : x[3] > balLvl_ntup.ref[2],carIt_arr)]
			else
            	carIt_arr = carIt_arr[findall(x -> j == 2 ? x[2] >= balLvl_ntup.ref[1] : x[3] >= balLvl_ntup.ref[2],carIt_arr)]
			end
            push!(carConstr_arr, carIt_arr...)
        end

        #  identifies and addresses "crossings" (one carrier is temporal more  but spatially less detailed than another one) by converting into new restriction mappings
        carConstrUni_arr = unique(carConstr_arr)
        cross_arr = filter(x -> (any(map(z -> (x[2] > z[2] && x[3] < z[3]) || (x[3] > z[3] && x[2] < z[2]), carConstrUni_arr))),carConstrUni_arr)
        if !isempty(cross_arr)
            newEntry_tup = (union(getindex.(cross_arr,1)...), minimum(getindex.(cross_arr,2)),minimum(getindex.(cross_arr,3)))
            carConstrUni_arr = unique(vcat(newEntry_tup,carConstrUni_arr))
        end

        # filter redundant and "dominated" combinations (less or the same carriers, but not more temporal or spatial detail)
        carConstrUni_arr2 = map(i -> map(x -> x[i],carConstrUni_arr),1:3)
        carConFilt_arr = filter(carConstrUni_arr) do x
                        !(any((BitArray(issubset(x[1],y) for y in carConstrUni_arr2[1]) 		.&
                        (((x[2] .<= carConstrUni_arr2[2]) .& (x[3] .< carConstrUni_arr2[3]))	.|
                        ((x[2] .< carConstrUni_arr2[2]) .& (x[3] .<= carConstrUni_arr2[3])) 	.|
                        ((x[2] .<= carConstrUni_arr2[2]) .& (x[3] .<= carConstrUni_arr2[3]))))	.& BitArray(map(y -> y != x,carConstrUni_arr)))) end

        carConFilt_arr2 = map(i -> map(x -> x[i],carConFilt_arr),1:3)
		typeCapa_sym = side == :use ? "in" : "out"

        # adds necessary capacity restrictions below reference level
        map(x -> push!(capaDispRestr_arr,(typeCapa_sym, carConFilt_arr2[1][x], carConFilt_arr2[2][x], carConFilt_arr2[3][x])),1:length(carConFilt_arr))
    end

    # XXX writes dimension of capacity restrictions for storage
    stInVar_arr, stOutVar_arr = [intersect(x,keys(carGrp_ntup)) for x in ((:stExtIn,:stIntIn),(:stExtOut,:stIntOut))]
    if !isempty(stInVar_arr) || !isempty(stOutVar_arr)
        allCar_arr = unique(vcat(collect.([getproperty(carGrp_ntup,y) for y in union(stInVar_arr,stOutVar_arr)])...))
        for x in allCar_arr
            carRow_ntup = anyM.cInfo[x]
            # storage on carrier level, but at least on reference level, if region is disaggregated balance on expansion (which is at least lower)
            tsLvl_int = balLvl_ntup.ref != nothing ? max(balLvl_ntup.ref[1], carRow_ntup.tsDis) : carRow_ntup.tsDis
            rLvl_int  = disAgg_boo ? balLvl_ntup.exp[2] : carRow_ntup.rDis |> (z -> balLvl_ntup.ref != nothing ? max(balLvl_ntup.ref[2], z) : z)
            if !isempty(stInVar_arr) push!(capaDispRestr_arr,("stIn", [x], tsLvl_int, rLvl_int)) end
            if !isempty(stOutVar_arr) push!(capaDispRestr_arr,("stOut", [x], tsLvl_int, rLvl_int)) end
            push!(capaDispRestr_arr,("stSize", [x], tsLvl_int, rLvl_int))
        end
    end

    part.capaRestr = isempty(capaDispRestr_arr) ? DataFrame() : rename(DataFrame(capaDispRestr_arr), :1 => :cnstrType, :2 => :car, :3 => :lvlTs, :4 => :lvlR)
end

max_stars_repo_name: leonardgoeke/anyMOD

max_issues_repo_name: leonardgoeke/anyMOD

max_forks_repo_name: leonardgoeke/anyMOD


text:
using Base.Test
using DataFrames

import Postgres
P = Postgres
naval = P.Types.naval
TESTDB = "julia_test"

function first_conn(;args...)
    conn = nothing
    try
        conn = P.connect(P.PostgresServer; args...)
    catch err PostgresError
        if (ismatch(r"does not exist$", err.msg))
            run(`createdb julia_test`)
            conn = P.connect(P.PostgresServer; args...)
        else
            throw(err)
        end
    end
    conn
end

function setup_db(curs::P.PostgresCursor)
    queries = [
        """drop table if exists newtable;""",

        """drop type if exists enum_test cascade;
        create type enum_test as enum ('happy', 'sad');""",

        """drop domain if exists domain_test cascade;
        create domain domain_test as int
            check(value >= 0 and value <= 10);""",

        """select setseed(0);
        drop table if exists test;
        create table test as select
             random() as x1
            ,null::domain_test as x2
            ,null::enum_test as f1
            ,null::float8 as y
            from generate_series(1, 1000);
        select setseed(0);
        update test set x2 = random() * 10;
        select setseed(0);
        update test set f1 = (case
            when random() > .7 then 'happy'
            else 'sad'
        end)::enum_test;
        select setseed(0);
        update test set y =
            ((1/4. * x1) + (1/40. * x2) +
                (case
                    when f1='happy' then .7
                    else .3
                end * 1/4)
                + random()*1/4.);

        """
    ]
    for q in queries
        P.query(curs, q)
    end
end

function do_plsql(curs::P.PostgresCursor, cmd::AbstractString)
    query(curs, """
    do
    \$\$begin
    $cmd;
    end\$\$;""")
end
suppress = IOBuffer()

#basic conection
@test_throws P.PostgresError first_conn(db=TESTDB, host="/dev/null/")
conn = first_conn(db=TESTDB, host="localhost")
#try other ways to connect
conn = connect(P.PostgresServer, "postgresql://localhost/julia_test")
conn = connect(P.PostgresServer, "", TESTDB, "localhost", "", "")
version = versioninfo(conn)
@test version[:protocol] == v"3.0.0"
print(suppress, conn)
@test P.status(conn) == :ok
@test isopen(conn)
curs = P.cursor(conn)
print(suppress, curs)
@test P.query(curs, "select 1")[1][1] == 1

setup_db(curs)
#transactions
P.query(curs, "drop table if exists xxx")

P.begin_!(curs)
P.query(curs, "create table xxx (a int); select * from xxx;")
P.rollback!(curs)
@test_throws P.Results.PostgresServerError P.query(curs, "select * from xxx")
P.commit!(curs)

# close connnection so the connection will find the new user defined types.
close(curs)
print(suppress, curs)
close(conn)
@test P.status(conn) == :not_connected
@test_throws P.PostgresError P.query(curs, "select 1")

conn = first_conn(db=TESTDB, host="localhost")
curs = P.cursor(conn)

#round trip types
for t in values(P.Types.base_types)
    print(suppress, t)
    # does not exists in PG
    if t.name == :jlunknown
        continue
    end
    start = repr(P.Types.PostgresValue(naval(t)))
    val = P.query(curs, "select $start")[1][1]
    @test typeof(naval(t)) == typeof(val)
    @test naval(t) == val
end

#extended types
types = [v for v in values(conn.pgtypes)]
enum_test = filter(x->x.name==:enum_test, types)[1]
@test enum_test.enumvals == Set(UTF8String["sad","happy"])
domain_test = filter(x->x.name==:domain_test, types)[1]
print(suppress, domain_test)
print(suppress, enum_test)

#basic query
df = P.query(curs, "select * from test")
@test size(df) == (1000,4)
@test eltype(df[1]) == Float64
@test eltype(df[2]) == Int
@test eltype(df[3]) == UTF8String
df[1, 1] = NA
df[1, 2] = NA
df[1, 3] = NA
P.copyto(curs, df, "test")
P.copyto(curs, df, "newtable", true)

#escaping
hi ="1;select 'powned'"
P.escape_value(conn, "stuff=$hi")

#result interface

res = P.execute(curs, "select 1, null::int, 'HI'::text, 1.2::float8
            from generate_series(1, 100)")
print(suppress, res)
@test !isempty(res)
@test size(res) == (100, 4)
@test size(res[:, 1]) == (100,)
@test length(res[1, :]) == 4
@test length([r for r in res]) == 100
@test length(P.Results.row(res, 1)) == 4
@test length(P.Results.column(res, 1)) == 100
@test_throws BoundsError res[0, :]
@test_throws BoundsError res[:, 0]
@test_throws BoundsError res[0, 0]
@test_throws BoundsError res[:, 5]
@test_throws BoundsError res[101, :]
@test_throws BoundsError res[101, 5]
for i in 1:length(res.types)
    t = res.types[i]
    tt = typeof(naval(t))
    @test isa(res[1,i], Nullable{tt})
end
@test !isnull(res[1, 1])
@test isnull(res[1, 2])

#iteration
curs = P.cursor(conn)
streamed = P.cursor(conn, 10)
function count_iters(curs)
    last = nothing
    for (i, r) in enumerate(curs)
        last = i
    end
    last
end

P.execute(streamed, "select 1 from generate_series(1, 23)")
@test count_iters(streamed) == 4
@test count_iters(streamed) == nothing

P.execute(curs, "select 1 from generate_series(1, 23)")
@test count_iters(curs) == 1
@test count_iters(curs) == nothing

close(curs)
close(streamed)
close(conn)
# this will fail if PG still thinks were connected
#FIXME one connection still open
#run(`dropdb julia_test`)

max_stars_repo_name: NCarson/Postgres

max_issues_repo_name: NCarson/Postgres

max_forks_repo_name: NCarson/Postgres


text:
"""
    ODEProblem(system::GXBeam.System, assembly, tspan; kwargs...)

Construct a `ODEProblem` for the system of nonlinear beams
contained in `assembly` which may be used with the DifferentialEquations package.

Keyword Arguments:
 - `prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}()`:
        A dictionary with keys corresponding to the points at
        which prescribed conditions are applied and elements of type
        [`PrescribedConditions`](@ref) which describe the prescribed conditions
        at those points.  If time varying, this input may be provided as a
        function of time.
 - `distributed_loads = Dict{Int,DistributedLoads{Float64}}()`: A dictionary
        with keys corresponding to the elements to which distributed loads are
        applied and elements of type [`DistributedLoads`](@ref) which describe
        the distributed loads on those elements.  If time varying, this input may
        be provided as a function of time.
 - `point_masses = Dict{Int,Vector{PointMass{Float64}}}()`: A dictionary with keys
        corresponding to the points at which point masses are attached and values
        containing vectors of objects of type [`PointMass`](@ref) which describe
        the point masses attached at those points.  If time varying, this input may
        be provided as a function of time.
 - `gravity`: Gravity vector. If time varying, this input may be provided as a
        function of time.
 - `origin = zeros(3)`: Global frame origin vector. If time varying, this input
        may be provided as a function of time.
 - `linear_velocity = zeros(3)`: Global frame linear velocity vector. If time
        varying, this vector may be provided as a function of time.
 - `angular_velocity = zeros(3)`: Global frame angular velocity vector. If time
        varying, this vector may be provided as a function of time.
 - `linear_acceleration = zeros(3)`: Global frame linear acceleration vector. If time
       varying, this vector may be provided as a function of time.
 - `angular_acceleration = zeros(3)`: Global frame angular acceleration vector. If time
       varying, this vector may be provided as a function of time.
"""
function SciMLBase.ODEProblem(system::System, assembly, tspan;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads = Dict{Int,DistributedLoads{Float64}}(),
    point_masses = Dict{Int,Vector{PointMass{Float64}}}(),
    gravity = (@SVector zeros(3)),
    origin = (@SVector zeros(3)),
    linear_velocity = (@SVector zeros(3)),
    angular_velocity = (@SVector zeros(3)),
    linear_acceleration = (@SVector zeros(3)),
    angular_acceleration = (@SVector zeros(3)),
    )

    # create ODEFunction
    func = SciMLBase.ODEFunction(system, assembly)

    # use initial state from `system`
    u0 = copy(system.x)

    # set parameters
    p = (prescribed_conditions, distributed_loads, point_masses, gravity, origin,
        linear_velocity, angular_velocity, linear_acceleration, angular_acceleration)

    return SciMLBase.ODEProblem{true}(func, u0, tspan, p)
end

"""
    ODEFunction(system::GXBeam.System, assembly)

Construct a `ODEFunction` for the system of nonlinear beams
contained in `assembly` which may be used with the DifferentialEquations package.

The parameters associated with the resulting ODEFunction are defined by the tuple
`(prescribed_conditions, distributed_loads, point_masses, origin, linear_velocity,
angular_velocity, linear_acceleration, angular_acceleration)` where each parameter is defined as follows:
 - `prescribed_conditions`: A dictionary with keys corresponding to the points at
        which prescribed conditions are applied and elements of type
        [`PrescribedConditions`](@ref) which describe the prescribed conditions
        at those points.  If time varying, this input may be provided as a
        function of time.
 - `distributed_loads`: A dictionary with keys corresponding to the elements to
        which distributed loads are applied and elements of type
        [`DistributedLoads`](@ref) which describe the distributed loads at those
        points.  If time varying, this input may be provided as a function of
        time.
 - `point_masses = Dict{Int,Vector{PointMass{Float64}}}()`: A dictionary with keys
        corresponding to the points at which point masses are attached and values
        containing vectors of objects of type [`PointMass`](@ref) which describe
        the point masses attached at those points.  If time varying, this input may
        be provided as a function of time.
 - `gravity`: Gravity vector. If time varying, this input may be provided as a
        function of time.
 - `origin`: Global frame origin vector. If time varying, this input
        may be provided as a function of time.
 - `linear_velocity`: Global frame linear velocity vector. If time
        varying, this vector may be provided as a function of time.
 - `angular_velocity`: Global frame angular velocity vector. If time
        varying, this vector may be provided as a function of time.
 - `linear_acceleration = zeros(3)`: Global frame linear acceleration vector. If time
       varying, this vector may be provided as a function of time.
 - `angular_acceleration = zeros(3)`: Global frame angular acceleration vector. If time
       varying, this vector may be provided as a function of time.
"""
function SciMLBase.ODEFunction(system::System, assembly)

    # check to make sure the system isn't static
    @assert !system.static

    # unpack system pointers
    irow_point = system.irow_point
    irow_elem = system.irow_elem
    irow_elem1 = system.irow_elem1
    irow_elem2 = system.irow_elem2
    icol_point = system.icol_point
    icol_elem = system.icol_elem

    # unpack scaling parameters
    force_scaling = system.force_scaling

    # DAE function
    f = function(resid, u, p, t)

        # get current parameters
        prescribed_conditions = typeof(p[1]) <: AbstractDict ? p[1] : p[1](t)
        distributed_loads = typeof(p[2]) <: AbstractDict ? p[2] : p[2](t)
        point_masses = typeof(p[3]) <: AbstractDict ? p[3] : p[3](t)
        gvec = typeof(p[4]) <: AbstractVector ? SVector{3}(p[4]) : SVector{3}(p[4](t))
        x0 = typeof(p[5]) <: AbstractVector ? SVector{3}(p[5]) : SVector{3}(p[5](t))
        v0 = typeof(p[6]) <: AbstractVector ? SVector{3}(p[6]) : SVector{3}(p[6](t))
        ω0 = typeof(p[7]) <: AbstractVector ? SVector{3}(p[7]) : SVector{3}(p[7](t))
        a0 = typeof(p[8]) <: AbstractVector ? SVector{3}(p[8]) : SVector{3}(p[8](t))
        α0 = typeof(p[9]) <: AbstractVector ? SVector{3}(p[9]) : SVector{3}(p[9](t))

        # calculate residual
        steady_state_system_residual!(resid, u, assembly, prescribed_conditions,
            distributed_loads, point_masses, gvec, force_scaling, irow_point, irow_elem, irow_elem1,
            irow_elem2, icol_point, icol_elem, x0, v0, ω0, a0, α0)

        return resid
    end

    update_mass_matrix! = function(M, u, p, t)

       point_masses = typeof(p[3]) <: AbstractDict ? p[3] : p[3](t)

        # zero out all mass matrix entries
        M .= 0.0

        # calculate mass matrix
        system_mass_matrix!(M, u, assembly, point_masses, force_scaling,
            irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

        return M
    end

    mass_matrix = SciMLBase.DiffEqArrayOperator(copy(system.M), update_func = update_mass_matrix!)

    # jacobian function with respect to states/state rates
    jac = function(J, u, p, t)

        # zero out all jacobian entries
        J .= 0.0

        # get current parameters
        prescribed_conditions = typeof(p[1]) <: AbstractDict ? p[1] : p[1](t)
        distributed_loads = typeof(p[2]) <: AbstractDict ? p[2] : p[2](t)
        point_masses = typeof(p[3]) <: AbstractDict ? p[3] : p[3](t)
        gvec = typeof(p[4]) <: AbstractVector ? SVector{3}(p[4]) : SVector{3}(p[4](t))
        x0 = typeof(p[5]) <: AbstractVector ? SVector{3}(p[5]) : SVector{3}(p[5](t))
        v0 = typeof(p[6]) <: AbstractVector ? SVector{3}(p[6]) : SVector{3}(p[6](t))
        ω0 = typeof(p[7]) <: AbstractVector ? SVector{3}(p[7]) : SVector{3}(p[7](t))
        a0 = typeof(p[8]) <: AbstractVector ? SVector{3}(p[8]) : SVector{3}(p[8](t))
        α0 = typeof(p[9]) <: AbstractVector ? SVector{3}(p[9]) : SVector{3}(p[9](t))

        # calculate jacobian
        steady_state_system_jacobian!(J, u, assembly, prescribed_conditions,
            distributed_loads, point_masses, gvec, force_scaling, irow_point, irow_elem, irow_elem1,
            irow_elem2, icol_point, icol_elem, x0, v0, ω0, a0, α0)

        return J
    end

    # sparsity structure
    sparsity = get_sparsity(system, assembly)

    # jacobian prototype (use dense since sparse isn't working)
    jac_prototype = collect(system.K)

    # TODO: figure out how to use a sparse matrix here.
    # It's failing with a singular exception during the LU factorization.

    return SciMLBase.ODEFunction{true,true}(f; mass_matrix, jac, jac_prototype, sparsity)
end

"""
    DAEProblem(system::GXBeam.System, assembly, tspan; kwargs...)

Construct a `DAEProblem` for the system of nonlinear beams
contained in `assembly` which may be used with the DifferentialEquations package.

A consistent set of initial conditions may be obtained prior to constructing the
DAEProblem using [`initial_condition_analysis!`](@ref) or by constructing a
DAEProblem after a time domain analysis.

Keyword Arguments:
 - `prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}()`:
        A dictionary with keys corresponding to the points at
        which prescribed conditions are applied and elements of type
        [`PrescribedConditions`](@ref) which describe the prescribed conditions
        at those points.  If time varying, this input may be provided as a
        function of time.
 - `distributed_loads = Dict{Int,DistributedLoads{Float64}}()`: A dictionary
        with keys corresponding to the elements to which distributed loads are
        applied and elements of type [`DistributedLoads`](@ref) which describe
        the distributed loads on those elements.  If time varying, this input may
        be provided as a function of time.
 - `point_masses = Dict{Int,Vector{PointMass{Float64}}}()`: A dictionary with keys
        corresponding to the points at which point masses are attached and values
        containing vectors of objects of type [`PointMass`](@ref) which describe
        the point masses attached at those points.  If time varying, this input may
        be provided as a function of time.
 - `gravity = zeros(3)`: Gravity vector. If time varying, this input may be provided as a
        function of time.
 - `origin = zeros(3)`: Global frame origin vector. If time varying, this input
        may be provided as a function of time.
 - `linear_velocity = zeros(3)`: Global frame linear velocity vector. If time
        varying, this vector may be provided as a function of time.
 - `angular_velocity = zeros(3)`: Global frame angular velocity vector. If time
        varying, this vector may be provided as a function of time.
 - `linear_acceleration = zeros(3)`: Global frame linear acceleration vector. If time
       varying, this vector may be provided as a function of time.
 - `angular_acceleration = zeros(3)`: Global frame angular acceleration vector. If time
       varying, this vector may be provided as a function of time.
"""
function SciMLBase.DAEProblem(system::System, assembly, tspan;
    prescribed_conditions = Dict{Int,PrescribedConditions{Float64}}(),
    distributed_loads = Dict{Int,DistributedLoads{Float64}}(),
    point_masses = Dict{Int,Vector{PointMass{Float64}}}(),
    gravity = (@SVector zeros(3)),
    origin = (@SVector zeros(3)),
    linear_velocity = (@SVector zeros(3)),
    angular_velocity = (@SVector zeros(3)),
    linear_acceleration = (@SVector zeros(3)),
    angular_acceleration = (@SVector zeros(3)),
    )

    # create SciMLBase.DAEFunction
    func = SciMLBase.DAEFunction(system, assembly)

    # use initial state from `system`
    u0 = copy(system.x)

    # use initial state rates from `system`
    du0 = zero(u0)
    for (ielem, icol) in enumerate(system.icol_elem)
        du0[icol:icol+2] = system.udot[ielem]
        du0[icol+3:icol+5] = system.θdot[ielem]
        du0[icol+12:icol+14] = system.Vdot[ielem]
        du0[icol+15:icol+17] = system.Ωdot[ielem]
    end

    # set parameters
    p = (prescribed_conditions, distributed_loads, point_masses, gravity, origin,
       linear_velocity, angular_velocity, linear_acceleration, angular_acceleration)

    # get differential variables
    differential_vars = get_differential_vars(system)

    return SciMLBase.DAEProblem{true}(func, du0, u0, tspan, p; differential_vars)
end

"""
    DAEFunction(system::GXBeam.System, assembly)

Construct a `DAEFunction` for the system of nonlinear beams
contained in `assembly` which may be used with the DifferentialEquations package.

The parameters associated with the resulting SciMLBase.DAEFunction are defined by the tuple
`(prescribed_conditions, distributed_loads, point_masses, origin, linear_velocity,
angular_velocity, linear_acceleration, angular_acceleration)`
where each parameter is defined as follows:
 - `prescribed_conditions`: A dictionary with keys corresponding to the points at
        which prescribed conditions are applied and elements of type
        [`PrescribedConditions`](@ref) which describe the prescribed conditions
        at those points.  If time varying, this input may be provided as a
        function of time.
 - `distributed_loads`: A dictionary with keys corresponding to the elements to
        which distributed loads are applied and elements of type [`DistributedLoads`](@ref)
        which describe the distributed loads on those elements.  If time varying, this
        input may be provided as a function of time.
 - `point_masses = Dict{Int,Vector{PointMass{Float64}}}()`: A dictionary with keys
        corresponding to the points at which point masses are attached and values
        containing vectors of objects of type [`PointMass`](@ref) which describe
        the point masses attached at those points.  If time varying, this input may
        be provided as a function of time.
 - `gravity`: Gravity vector. If time varying, this input may be provided as a
        function of time.
 - `origin`: Global frame origin vector. If time varying, this input
        may be provided as a function of time.
 - `linear_velocity`: Global frame linear velocity vector. If time
        varying, this vector may be provided as a function of time.
 - `angular_velocity`: Global frame angular velocity vector. If time
        varying, this vector may be provided as a function of time.
 - `linear_acceleration = zeros(3)`: Global frame linear acceleration vector. If time
       varying, this vector may be provided as a function of time.
 - `angular_acceleration = zeros(3)`: Global frame angular acceleration vector. If time
       varying, this vector may be provided as a function of time.
"""
function SciMLBase.DAEFunction(system::System, assembly)

    # check to make sure the system isn't static
    @assert !system.static

    # unpack system pointers
    irow_point = system.irow_point
    irow_elem = system.irow_elem
    irow_elem1 = system.irow_elem1
    irow_elem2 = system.irow_elem2
    icol_point = system.icol_point
    icol_elem = system.icol_elem

    # unpack scaling parameters
    force_scaling = system.force_scaling

    # DAE function
    f = function(resid, du, u, p, t)

        # get current parameters
        prescribed_conditions = typeof(p[1]) <: AbstractDict ? p[1] : p[1](t)
        distributed_loads = typeof(p[2]) <: AbstractDict ? p[2] : p[2](t)
        point_masses = typeof(p[3]) <: AbstractDict ? p[3] : p[3](t)
        gvec = typeof(p[4]) <: AbstractVector ? SVector{3}(p[4]) : SVector{3}(p[4](t))
        x0 = typeof(p[5]) <: AbstractVector ? SVector{3}(p[5]) : SVector{3}(p[5](t))
        v0 = typeof(p[6]) <: AbstractVector ? SVector{3}(p[6]) : SVector{3}(p[6](t))
        ω0 = typeof(p[7]) <: AbstractVector ? SVector{3}(p[7]) : SVector{3}(p[7](t))
        a0 = typeof(p[8]) <: AbstractVector ? SVector{3}(p[8]) : SVector{3}(p[8](t))
        α0 = typeof(p[9]) <: AbstractVector ? SVector{3}(p[9]) : SVector{3}(p[9](t))

        # calculate residual
        dynamic_system_residual!(resid, du, u, assembly, prescribed_conditions,
            distributed_loads, point_masses, gvec, force_scaling, irow_point, irow_elem,
            irow_elem1, irow_elem2, icol_point, icol_elem, x0, v0, ω0, a0, α0)

        return resid
    end

    # jacobian function with respect to states/state rates
    jac = function(J, du, u, p, gamma, t)

        # zero out all jacobian entries
        J .= 0.0

        # get current parameters
        prescribed_conditions = typeof(p[1]) <: AbstractDict ? p[1] : p[1](t)
        distributed_loads = typeof(p[2]) <: AbstractDict ? p[2] : p[2](t)
        point_masses = typeof(p[3]) <: AbstractDict ? p[3] : p[3](t)
        gvec = typeof(p[4]) <: AbstractVector ? SVector{3}(p[4]) : SVector{3}(p[4](t))
        x0 = typeof(p[5]) <: AbstractVector ? SVector{3}(p[5]) : SVector{3}(p[5](t))
        v0 = typeof(p[6]) <: AbstractVector ? SVector{3}(p[6]) : SVector{3}(p[6](t))
        ω0 = typeof(p[7]) <: AbstractVector ? SVector{3}(p[7]) : SVector{3}(p[7](t))
        a0 = typeof(p[8]) <: AbstractVector ? SVector{3}(p[8]) : SVector{3}(p[8](t))
        α0 = typeof(p[9]) <: AbstractVector ? SVector{3}(p[9]) : SVector{3}(p[9](t))

        # calculate jacobian
        dynamic_system_jacobian!(J, du, u, assembly, prescribed_conditions,
            distributed_loads, point_masses, gvec, force_scaling, irow_point, irow_elem,
            irow_elem1, irow_elem2, icol_point, icol_elem, x0, v0, ω0, a0, α0)

        # add gamma multiplied by the mass matrix
        system_mass_matrix!(J, gamma, u, assembly, point_masses, force_scaling,
            irow_point, irow_elem, irow_elem1, irow_elem2, icol_point, icol_elem)

        return J
    end

    # sparsity structure
    sparsity = get_sparsity(system, assembly)

    # jacobian prototype (use dense since sparse isn't working)
    # jac_prototype = collect(system.K)

    # TODO: figure out how to use a sparse matrix here.
    # It's failing with a singular exception during the LU factorization.
    # Using `jac_prototype` also causes errors

    return SciMLBase.DAEFunction{true,true}(f; jac, sparsity)
end

function get_differential_vars(system::System)
    differential_vars = fill(false, length(system.x))
    for icol in system.icol_elem
        differential_vars[icol:icol+2] .= true # u (for the beam element)
        differential_vars[icol+3:icol+5] .= true # θ (for the beam element)
        differential_vars[icol+12:icol+14] .= true # V (for the beam element)
        differential_vars[icol+15:icol+17] .= true # Ω (for the beam element)
    end
    return differential_vars
end


text:
# # tests for: 1/N \sum_i (<x_i, Q_ix_i> + <x_i,q_i> + η dist(x_i,B_i) + g(\sum_i x_i)

using Test
using LinearAlgebra
using CIAOAlgorithms
using ProximalOperators
using Base.Iterators: take

T = Float64

n = 2;
N = 3;
F = Vector{Sum}(undef, 0)
L = Vector{T}(undef, 0)
η = N * T(10)
box = IndBox(-T(2), T(2))
d = [T[1.0, 2.0], T[-1.0, 3.0], T[0.0, 10.0]]
for i = 1:N
    f2 = SqrDistL2(box, η) # soft box constraint
    Q = diagm(0 => d[i])
    f1 = Quadratic(Q, ones(T, n))
    push!(F, Sum(f1, f2))
    push!(L, opnorm(Q[i]) + η)
end
g = IndBox(-Inf, ones(T, n))   # sum_i x_i \leq 1
x0 = zeros(T, n)

sum_star = [-5.136781609195401, -0.9333333333333327]

maxit = 1000
tol = 1e-4

@testset "Proshi" begin
    # sweeping 1, 2, 3 for randomined, cyclic and shuffled sampling strategies, respectively.

    ## test the solver
    # basic Proshi
    @testset "basic Proshi" for sweeping in collect(1:3)
        solver = CIAOAlgorithms.Proshi{T}(maxit = maxit, sweeping = sweeping)
        x_proshi, it_proshi = solver(x0, F = F, g = g, L = L, N = N)
        @test norm(sum(x_proshi) - sum_star, Inf) < tol
        @test eltype(x_proshi) == Array{T,1}
    end

    # basic Proshi with minibatch
    vec_ref = [(1, 2), (2, 2), (3, 3)] # different samplings and batch sizes
    @testset "Proshi_minibatch" for (sweeping, batch) in vec_ref
        solver = CIAOAlgorithms.Proshi{T}(
            maxit = maxit,
            sweeping = sweeping,
            minibatch = (true, batch),
        )
        x_proshi, it_proshi = solver(x0, F = F, g = g, L = L, N = N)
        @test norm(sum(x_proshi) - sum_star, Inf) < tol
        @test eltype(x_proshi) == Array{T,1}
    end

    # test with user defined stepsizes
    @testset "γ and L as scalars" begin
        @testset "randomized" begin
            γ = N / maximum(L)
            solver = CIAOAlgorithms.Proshi{T}(maxit = maxit, γ = γ)
            x_proshi, it_proshi = solver(x0, F = F, g = g, L = L, N = N)
            @test norm(sum(x_proshi) - sum_star, Inf) < tol
        end
        @testset "cyclic" begin
            solver = CIAOAlgorithms.Proshi{T}(maxit = maxit)
            x_proshi, it_proshi = solver(x0, F = F, g = g, L = maximum(L), N = N)
            @test norm(sum(x_proshi) - sum_star, Inf) < tol
        end
    end

    ## test the iterator
    @testset "the iterator" for sweeping in collect(1:3)
        solver = CIAOAlgorithms.Proshi{T}(sweeping = sweeping)
        iter = CIAOAlgorithms.iterator(solver, x0, F = F, g = g, L = L, N = N)
        @test iter.x0 === x0

        for state in take(iter, 2)
            @test solution(state) === state.s
            @test eltype(solution(state)) == Array{T,1}
        end
    end
end


text:
module TonieBox
    using HTTP
    using JSON3
    using StructTypes

    function __init__()
        authenticate()
    end

    const AUTH_URL = "https://login.tonies.com/auth/realms/tonies/protocol/openid-connect/token"

    const _access_token = Ref{String}()

    access_token() = if isassigned(_access_token)
        _access_token[]
    else
        error("No access token available. Call `authenticate` first.")
    end

    function authenticate(user, password)
        response = HTTP.post(
            AUTH_URL,
            ["Content-Type" => "application/x-www-form-urlencoded"],
             HTTP.URIs.escapeuri(
                Dict(
                    "grant_type" => "password",
                    "client_id" => "my-tonies",
                    "scope" => "openid",
                    "username" => user,
                    "password" => password
                )
            )
        )
        _access_token[] = JSON3.read(response.body)[:access_token]
        @info "Authentication successful"
        return
    end

    function authenticate()
        @info "Authentication needed."
        print("Username: ")
        user = readline()
        print("Password: ")
        password = readline()
        authenticate(user, password)
    end

    function me()
        HTTP.get(
            "https://api.tonie.cloud/v2/me",
            ["Authorization" => "Bearer $(access_token())"],
        ).body |> JSON3.read
    end

    struct Household
        id::String
        name::String
        image::String
        foreignCreativeTonieContent::Bool
        access::String
        canLeave::Bool
        ownerName::String
    end
    StructTypes.StructType(::Type{Household}) = StructTypes.Struct()

    function households()
        JSON3.read(
            HTTP.get(
                "https://api.tonie.cloud/v2/households",
                ["Authorization" => "Bearer $(access_token())"],
            ).body,
            Vector{Household}
        )
    end

    const _current_household = Ref{Household}()
    function current_household()
        if isassigned(_current_household)
            _current_household[]
        else
            hs = households()
            if length(hs) == 1
                current_household!(only(hs))
            else
                error("No current household set. Can't set household automatically because there are $(length(hs)) households.")
            end
        end
    end

    function current_household!(household::Household)
        _current_household[] = household
    end

    struct DeletedChapter
        title::String
        seconds::Float64
    end
    StructTypes.StructType(::Type{DeletedChapter}) = StructTypes.Struct()

    struct TranscodingError
        reason::String
        deletedChapters::Vector{DeletedChapter}
    end
    StructTypes.StructType(::Type{TranscodingError}) = StructTypes.Struct()

    struct Chapter
        id::String
        title::String
        file::String
        seconds::Float64
        transcoding::Bool
    end
    StructTypes.StructType(::Type{Chapter}) = StructTypes.Struct()
    Base.show(io::IO, c::Chapter) = print(io, """Chapter(title: "$(c.title)", seconds: $(c.seconds))""")

    struct CreativeTonie
        id::String
        householdId::String
        name::String
        live::Bool
        private::Bool
        imageUrl::String
        transcodingErrors::Vector{TranscodingError}
        secondsRemaining::Float64
        secondsPresent::Float64
        chaptersRemaining::Int
        chaptersPresent::Int
        transcoding::Bool
        chapters::Vector{Chapter}
    end
    StructTypes.StructType(::Type{CreativeTonie}) = StructTypes.Struct()
    Base.show(io::IO, ct::CreativeTonie) = print(io, """CreativeTonie(name: "$(ct.name)", number of chapters: $(length(ct.chapters)))""")

    function creativetonies(household = current_household())
        JSON3.read(
            HTTP.get(
                "https://api.tonie.cloud/v2/households/$(household.id)/creativetonies",
                ["Authorization" => "Bearer $(access_token())"],
            ).body,
            Vector{CreativeTonie}
        )
    end

    function create_presigned_file_url()
        JSON3.read(
            HTTP.post(
                "https://api.tonie.cloud/v2/file",
                ["Authorization" => "Bearer $(access_token())"],
            ).body
        )
    end

    struct AddChapter
        title::String
        file::String
        origin::String
    end
    StructTypes.StructType(::Type{AddChapter}) = StructTypes.Struct()

    function add_chapter_to_creative_tonie(tonie::CreativeTonie, file, title; origin = "file-julia")
        @assert isfile(file)

        @info "Creating pre-signed file url"
        presigned_file_info = create_presigned_file_url()
        fileId = presigned_file_info.fileId
        request = presigned_file_info.request

        addchapter = AddChapter(
            title,
            fileId,
            origin,
        )

        amazon_url = request.url
        multiparts = collect(pairs(request.fields))

        @info "Uploading file to Amazon S3."
        open(file) do io
            HTTP.post(
                amazon_url,
                [],
                HTTP.Form(
                    vcat(multiparts, "file" => io)
                )
            )
        end

        @info "Adding uploaded file to chapter list. [id=$fileId]"
        HTTP.post(
            "https://api.tonie.cloud/v2/households/$(tonie.householdId)/creativetonies/$(tonie.id)/chapters",
            ["Content-Type" => "application/json", "Authorization" => "Bearer $(access_token())"],
            JSON3.write(addchapter)
        )
        return
    end

    function remove_chapter(creativetonie::CreativeTonie, chapter::Chapter; household = current_household())
        remaining_chapters = filter(!=(chapter), creativetonie.chapters)
        @assert length(remaining_chapters) == length(creativetonie.chapters) - 1
        @info "Deleting chapter '$(chapter.title)' from creative tonie '$(creativetonie.name)'."
        HTTP.patch(
            "https://api.tonie.cloud/v2/households/$(household.id)/creativetonies/$(creativetonie.id)",
            ["Content-Type" => "application/json", "Authorization" => "Bearer $(access_token())"],
            JSON3.write(Dict(
                "chapters" => remaining_chapters
            ))
        )
        @info "Deletion successful."
        return
    end

    function download_mp3(f, url)
        mktempdir() do path
            filetrunk = joinpath(path, "download")
            filename = filetrunk * ".%(ext)s"
            filename_mp3 = filetrunk * ".mp3"
            run(`youtube-dl -o $filename -x --audio-format mp3 $url`)
            @assert isfile(filename_mp3)
            f(filename_mp3)
        end
    end

    function download_mp3_and_add_chapter(creativetonie, url, title)
        download_mp3(url) do filepath
            add_chapter_to_creative_tonie(creativetonie, filepath, title)
        end
    end

    find_chapters(ct::CreativeTonie, r::Regex) = filter(c -> match(r, c.title), ct.chapters)
    find_chapters(ct::CreativeTonie, s::AbstractString) = filter(c -> occursin(s, c.title), ct.chapters)

    function remove_chapters(ct::CreativeTonie, matcher)
        chapters = find_chapters(ct, matcher)
        io = IOBuffer()
        println(io, "This would remove the following chapters:")
        for c in chapters
            println(io, "  - $(c.title)")
        end
        @info String(take!(io))
        println("Type y to remove:")
        s = readline()
        if s == "y"
            for c in chapters
                remove_chapter(ct, c)
            end
        else
            @info "Not removed."
        end
    end
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.concept.type;
#
# import typedb.client.api.TypeDBTransaction;
# import typedb.client.api.concept.thing.Thing;
# import typedb.client.api.concept.type.AttributeType.ValueType;
#
# import javax.annotation.CheckReturnValue;
# import java.util.stream.Stream;
#
# public interface ThingType extends Type {
#
#     @Override
#     @CheckReturnValue
#     default boolean isThingType() {
#         return true;
#     }
#
#     @Override
#     ThingType.Remote asRemote(TypeDBTransaction transaction);
#
#     interface Remote extends Type.Remote, ThingType {
#
#         @Override
#         @CheckReturnValue
#         ThingType getSupertype();
#
#         @Override
#         @CheckReturnValue
#         Stream<? extends ThingType> getSupertypes();
#
#         @Override
#         @CheckReturnValue
#         Stream<? extends ThingType> getSubtypes();
#
#         @CheckReturnValue
#         Stream<? extends Thing> getInstances();
#
#         void setAbstract();
#
#         void unsetAbstract();
#
#         void setPlays(RoleType roleType);
#
#         void setPlays(RoleType roleType, RoleType overriddenType);
#
#         void setOwns(AttributeType attributeType, AttributeType overriddenType, boolean isKey);
#
#         void setOwns(AttributeType attributeType, AttributeType overriddenType);
#
#         void setOwns(AttributeType attributeType, boolean isKey);
#
#         void setOwns(AttributeType attributeType);
#
#         @CheckReturnValue
#         Stream<? extends RoleType> getPlays();
#
#         @CheckReturnValue
#         Stream<? extends AttributeType> getOwns();
#
#         @CheckReturnValue
#         Stream<? extends AttributeType> getOwns(ValueType valueType);
#
#         @CheckReturnValue
#         Stream<? extends AttributeType> getOwns(boolean keysOnly);
#
#         @CheckReturnValue
#         Stream<? extends AttributeType> getOwns(ValueType valueType, boolean keysOnly);
#
#         void unsetPlays(RoleType roleType);
#
#         void unsetOwns(AttributeType attributeType);
#     }
# }


text:
# # Computing Fibonacci Numbers
# The following is an example that everyone likes, computing Fibonacci number recursively.
using NiLang

@i function rfib(out!, n::T) where T
    @routine begin
        n1 ← zero(T)
        n2 ← zero(T)
        n1 += n - 1
        n2 += n - 2
    end
    if (value(n) <= 2, ~)
        out! += 1
    else
        rfib(out!, n1)
        rfib(out!, n2)
    end
    ~@routine
end

# The time complexity of this recursive algorithm is exponential to input `n`. It is also possible to write a reversible linear time with for loops.
# A slightly non-trivial task is computing the first Fibonacci number that greater or equal to a certain number `z`, where a `while` statement is required.

@i function rfibn(n!, z)
    @safe @assert n! == 0
    out ← 0
    rfib(out, n!)
    @from n! == 0 while out < z
        ~rfib(out, n!)
        n! += 1
        rfib(out, n!)
    end
    ~rfib(out, n!)
    out → 0
end

# In this example, the postcondition `n!=0` in the `while` statement is false before entering the loop, and it becomes true in later iterations. In the reverse program, the `while` statement stops at `n==0`.
# If executed correctly, a user will see the following result.

rfib(0, 10)

# compute which the first Fibonacci number greater than 100.

rfibn(0, 100)

# and uncompute

(~rfibn)(rfibn(0, 100)...)

# This example shows how an addition postcondition provided by the user can help to reverse a control flow without caching controls.


text:
import CUDAdrv
import NVTX

include("kernels.jl")

# NOTE: use with `--profile-from-start off`
NVTX.stop()

memory_initial = CUDAdrv.Mem.free()

const TFSTYLE = length(ARGS) >= 1 ? parse(Bool, ARGS[1]) : false
const DIMS = length(ARGS) >= 2 ? parse(Int, ARGS[2]) : 2048
const ITERATIONS = length(ARGS) >= 3 ? parse(Int, ARGS[3]) : 1
const ARITY = length(ARGS) >= 4 ? parse(Int, ARGS[4]) : 0
const UNIFORM = length(ARGS) >= 5 ? parse(Bool, ARGS[5]) : false

ARITY > 0 && @assert !TFSTYLE
UNIFORM   && @assert ARITY==0
const KERNEL!, INPUTS, DERIVS, BUFFERS = if ARITY == 0
    get_hmlstm_kernel(TFSTYLE, true, DIMS, UNIFORM)
else
    get_arity_scaling_kernel(true, DIMS, ARITY)
end

# to mimick TF, we'll copy back the results to the CPU.
# preallocate a page-locked memory buffer to avoid API call overhead.
# TODO: wrap in CUDAdrv
const buffer = Ref{Ptr{Cvoid}}()
const buffer_size = DIMS*DIMS*sizeof(Float32)
CUDAdrv.@apicall(:cuMemAllocHost, (Ptr{Ptr{Cvoid}}, Csize_t), buffer, buffer_size)

function benchmark(kernel!, inputs, derivs, buffers)
    NVTX.@range "kernel" begin
        kernel!(inputs, derivs, buffers)
        for deriv in DERIVS
            # we download results asynchronously to avoid API overhead
            # TODO: `Base.copyto!(; async)` in CUDAdrv
            @assert buffer_size == sizeof(deriv)
            CUDAdrv.Mem.download!(buffer[], deriv.buf, buffer_size; async=true)
        end
        CUDAdrv.synchronize()
    end
end

# warmup
NVTX.mark("Benchmark: $(42)")
benchmark(KERNEL!, INPUTS, DERIVS, BUFFERS)

memory = CUDAdrv.Mem.free()
NVTX.@activate CUDAdrv.@profile begin
    ccall(:jl_dump_compiles, Cvoid, (Ptr{Cvoid},), STDERR.handle)
    for i in 1:ITERATIONS
        GC.gc()
        benchmark(KERNEL!, INPUTS, DERIVS, BUFFERS)
    end
    NVTX.mark("Used Memory: $(memory_initial-memory)")
    ccall(:jl_dump_compiles, Cvoid, (Ptr{Cvoid},), C_NULL)
end

CUDAdrv.@apicall(:cuMemFreeHost, (Ptr{Cvoid},), buffer[])


text:
using Distributions, DiffEqBayesStan

println("Starting the test")
@test DiffEqBayesStan.stan_string(Bernoulli(1)) == "bernoulli(1.0)"
@test DiffEqBayesStan.stan_string(Bernoulli) == "bernoulli"

@test DiffEqBayesStan.stan_string(Binomial(5,0.3)) == "binomial(5, 0.3)"
@test DiffEqBayesStan.stan_string(Binomial) == "binomial"

@test DiffEqBayesStan.stan_string(BetaBinomial(5,1,1)) == "beta_binomial(5, 1.0, 1.0)"
@test DiffEqBayesStan.stan_string(BetaBinomial) == "beta_binomial"

@test DiffEqBayesStan.stan_string(Hypergeometric(5,5,3)) == "hypergeometric(5, 5, 3)"
@test DiffEqBayesStan.stan_string(Hypergeometric) == "hypergeometric"

@test DiffEqBayesStan.stan_string(NegativeBinomial(5,0.3)) == "neg_binomial(5.0, 0.3)"
@test DiffEqBayesStan.stan_string(NegativeBinomial) == "neg_binomial"

@test DiffEqBayesStan.stan_string(Poisson(5)) == "poisson(5.0)"
@test DiffEqBayesStan.stan_string(Poisson) == "poisson"

@test DiffEqBayesStan.stan_string(Normal(0,1)) == "normal(0.0, 1.0)"
@test DiffEqBayesStan.stan_string(Normal) == "normal"

@test DiffEqBayesStan.stan_string(TDist(5)) == "student_t(5.0,0,1)"
@test DiffEqBayesStan.stan_string(TDist) == "student_t"

@test DiffEqBayesStan.stan_string(Cauchy(0,1)) == "cauchy(0.0, 1.0)"
@test DiffEqBayesStan.stan_string(Cauchy) == "cauchy"

@test DiffEqBayesStan.stan_string(Laplace(0,1)) == "double_exponential(0.0, 1.0)"
@test DiffEqBayesStan.stan_string(Laplace) == "double_exponential"

@test DiffEqBayesStan.stan_string(Distributions.Logistic(0,1)) == "logistic(0.0, 1.0)"
@test DiffEqBayesStan.stan_string(Distributions.Logistic) == "logistic"

@test DiffEqBayesStan.stan_string(Gumbel(0,1)) == "gumbel(0.0, 1.0)"
@test DiffEqBayesStan.stan_string(Gumbel) == "gumbel"

@test DiffEqBayesStan.stan_string(LogNormal(0,1)) == "lognormal(0.0, 1.0)"
@test DiffEqBayesStan.stan_string(LogNormal) == "lognormal"

@test DiffEqBayesStan.stan_string(Chisq(5)) == "chi_square(5.0)"
@test DiffEqBayesStan.stan_string(Chisq) == "chi_square"

@test DiffEqBayesStan.stan_string(Exponential(5)) == "exponential(5.0)"
@test DiffEqBayesStan.stan_string(Exponential) == "exponential"

@test DiffEqBayesStan.stan_string(Gamma(2,3)) == "gamma(2.0, 3.0)"
@test DiffEqBayesStan.stan_string(Gamma) == "gamma"

@test DiffEqBayesStan.stan_string(InverseGamma(2,3)) == "inv_gamma(2.0, 3.0)"
@test DiffEqBayesStan.stan_string(InverseGamma) == "inv_gamma"

@test DiffEqBayesStan.stan_string(Weibull(1,1)) == "weibull(1.0, 1.0)"
@test DiffEqBayesStan.stan_string(Weibull) == "weibull"

@test DiffEqBayesStan.stan_string(Frechet(1,1)) =="frechet(1.0, 1.0)"
@test DiffEqBayesStan.stan_string(Frechet) == "frechet"

@test DiffEqBayesStan.stan_string(Rayleigh(5)) == "rayleigh(5.0)"
@test DiffEqBayesStan.stan_string(Rayleigh) == "rayleigh"

@test DiffEqBayesStan.stan_string(Pareto(2,3)) == "pareto(2.0, 3.0)"
@test DiffEqBayesStan.stan_string(Pareto) == "pareto"

@test DiffEqBayesStan.stan_string(GeneralizedPareto(0,1,2)) == "pareto_type_2(0.0, 1.0, 2.0)"
@test DiffEqBayesStan.stan_string(GeneralizedPareto) == "pareto_type_2"

@test DiffEqBayesStan.stan_string(Beta(3,3)) == "beta(3.0, 3.0)"
@test DiffEqBayesStan.stan_string(Beta) == "beta"

@test DiffEqBayesStan.stan_string(Uniform(1,4)) =="uniform(1.0, 4.0)"
@test DiffEqBayesStan.stan_string(Uniform) == "uniform"

@test DiffEqBayesStan.stan_string(VonMises(0,2)) == "von_mises(0.0, 2.0)"
@test DiffEqBayesStan.stan_string(VonMises) == "von_mises"


text:
function Base.:(*)(x::PolyVar{false}, y::PolyVar{false})
    if x === y
        Monomial{false}([x], [2])
    else
        Monomial{false}([x, y], [1, 1])
    end
end

function multiplyvar(v::Vector{PolyVar{false}}, z::Vector{Int}, x::PolyVar{false})
    i = length(v)
    while i > 0 && z[i] == 0
        i -= 1
    end
    if i > 0 && v[i] == x
        multiplyexistingvar(v, x, i)
    else
        #   ---->
        # \  |\  |\
        #  \ | \ | \
        #   \|  \|  \
        # If v[i] > x, we wait either for a rise (v[i] > v[i-1]) or v[i] < x
        # Otherwise, we first wait for a drop and then wait for the same thing
        ndrop = 0
        if i > 0 && v[i] > x
            droplim1 = 0
            droplim2 = 1
        else
            droplim1 = 1
            droplim2 = 2
        end
        i += 1
        while i <= length(v) && v[i] != x
            if i > 1 && v[i] > v[i - 1]
                ndrop += 1
            end
            if ndrop >= droplim2 || (ndrop >= droplim1 && v[i] < x)
                break
            end
            i += 1
        end

        if i <= length(v) && v[i] == x
            multiplyexistingvar(v, x, i)
        else
            insertvar(v, x, i)
        end
    end
end
function multiplyvar(x::PolyVar{false}, v::Vector{PolyVar{false}}, z::Vector{Int})
    i = 1
    while i <= length(v) && z[i] == 0
        i += 1
    end
    if i <= length(v) && v[i] == x
        multiplyexistingvar(v, x, i)
    else
        #   <----
        # \  |\  |\
        #  \ | \ | \
        #   \|  \|  \
        # If z[i] < x, we wait either for a drop (v[i] < v[i+1]) or v[i] > x
        # Otherwise, we first wait for a drop and then wait for the same thing
        ndrop = 0
        if i <= length(v) && v[i] < x
            droplim1 = 0
            droplim2 = 1
        else
            droplim1 = 1
            droplim2 = 2
        end
        i -= 1
        while i > 0 && v[i] != x
            if i < length(v) && v[i] < v[i + 1]
                ndrop += 1
            end
            if ndrop >= droplim2 || (ndrop >= droplim1 && v[i] > x)
                break
            end
            i -= 1
        end
        if i > 0 && v[i] == x
            multiplyexistingvar(v, x, i)
        else
            insertvar(v, x, i+1)
        end
    end
end
function Base.:(*)(x::PolyVar{false}, y::Monomial{false})
    w, updatez = multiplyvar(x, y.vars, y.z)
    Monomial{false}(w, updatez(y.z))
end
function Base.:(*)(y::Monomial{false}, x::PolyVar{false})
    w, updatez = multiplyvar(y.vars, y.z, x)
    Monomial{false}(w, updatez(y.z))
end

function Base.:(*)(x::Monomial{false}, y::Monomial{false})
    i = findlast(z -> z > 0, x.z)
    if i == nothing || i == 0
        return y
    end
    j = findfirst(z -> z > 0, y.z)
    if j == nothing || j == 0
        return x
    end
    if x.vars[i] == y.vars[j]
        w = [x.vars[1:i]; y.vars[j+1:end]]
        z = [x.z[1:i-1]; x.z[i] + y.z[j]; y.z[j+1:end]]
    else
        w = [x.vars[1:i]; y.vars[j:end]]
        z = [x.z[1:i]; y.z[j:end]]
    end
    return Monomial{false}(w, z)
end

function Base.:(*)(y::MonomialVector{false}, x::DMonomialLike{false})
    MonomialVector{false}([yi * x for yi in y])
end
function Base.:(*)(x::DMonomialLike{false}, y::MonomialVector{false})
    # The order may change
    # Example: y * [x^2, y^2] == [y^3, yx^2]
    MonomialVector{false}([x * yi for yi in y])
end


text:
function solve()
    k = parse(Int, readline())
    if k % 2 == 0 || k % 5 == 0
        -1
    else
        ret = 1
        x = 7
        while x % k ≠ 0
            x = x % k * 10 + 7
            ret += 1
        end
        ret
    end
end

println(solve())

max_stars_repo_name: aishikawa/atcoder-julia

max_issues_repo_name: aishikawa/atcoder-julia

max_forks_repo_name: aishikawa/atcoder-julia


text: # Test routines for the various types I've hacked together


text:
module PoseEstimationMathematicalProgram

using JuMP
using Gurobi
using StaticArrays
using Rotations
using GeometryTypes
using FileIO
using CoordinateTransformations
using MultilinearOpt

export PoseEstimationSolveData, constructModel

type PoseEstimationSolveData
    m
    R
    T
    R_abs
    R_diff
    phi
    alpha
    C
    f
    f_outlier
end

function constructModel(scenePoints, points, faces;
    use_relaxed_form = false, rotation_relaxation = :Fixed, other_solver=nothing, closeness=10.0, R0_in=zeros(3, 3), optPhiMax = 0.1, optBigNumber = 10, kwargs...)
    if other_solver == nothing
        m = Model(solver=GurobiSolver(;kwargs...))
    else
        m = Model(solver=other_solver)
    end

    N_scene = size(scenePoints, 2)
    N_model_verts = size(points, 1)
    N_model_faces = size(faces, 1)

    # Set up transform
    @variable(m, -1 <= R[1:3, 1:3] <= 1 )
    @variable(m, -optBigNumber <= T[1:3] <= optBigNumber )

    R_abs = nothing
    R_diff = nothing
    if rotation_relaxation == :Fixed
        @constraint(m, R .== R0_in)
    elseif rotation_relaxation == :Domain
        # Bound R to have row and column L1 norm <= sqrt(3), since
        # row and column L2 norm = 2.
        # (We can't do lower bound L1 norm >= 1 b/c that's not
        # a convex constraint)
        @variable(m, R_abs[1:3, 1:3] >= 0 )
        @constraint(m, R_abs .>= R)
        @constraint(m, R_abs .>= -R)
        for i in 1:3
            @constraint(m, [1 1 1] * R_abs[:, i] .<= sqrt(3))
            @constraint(m, [1 1 1] * R_abs[i, :] .<= sqrt(3))
        end

        # Bound R to be close to R0, in terms of L1 norm
        @variable(m, R_diff[1:3, 1:3] >= 0 )
        @constraint(m, R_diff .>= (R - R0_in))
        @constraint(m, R_diff .>= -(R - R0_in))
        for i in 1:3
            @constraint(m, dot([1; 1; 1], R_diff[:, i]) <= closeness)
            @constraint(m, dot([1; 1; 1], R_diff[i, :]) <= closeness)
        end
    elseif rotation_relaxation == :RtRBilinear
        # Orthogonality
        @constraint(m, transpose(R) * R .== eye(3, 3))
        #@constraint(m, R * transpose(R) .== eye(3, 3))
        # Forcing right-handedness of coordinates via cross products
        # (which should result in det(R) = +1)
        for i in 0:2
          @constraint(m, cross(R[(i+0)%3+1, :], R[(i+1)%3+1, :]) .== R[(i+2)%3+1, :])
          #@constraint(m, cross(R[:, (i+0)%3+1], R[:, (i+1)%3+1]) .== R[:, (i+2)%3+1])
        end
        relaxbilinear!(m, method=:Logarithmic2D, disc_level=4)
    end

    # Terms required for setting up objective:
    # Slack for distance-from-each-scene-point
    @variable(m, phi[1:N_scene] >= 0)

    # Model-scene affine coefficients
    @variable(m, 0 <= C[1:N_scene, 1:N_model_verts] <= 1.0)

    # Face selection
    if (use_relaxed_form)
        @variable(m, 0 <= f[1:N_scene, 1:N_model_faces] <= 1.0)
        @variable(m, 0 <= f_outlier[1:N_scene] <= 1.0)
    else
        @variable(m, 0 <= f[1:N_scene, 1:N_model_faces] <= 1.0, Bin)
        @variable(m, 0 <= f_outlier[1:N_scene] <= 1.0, Bin)
    end
    for i in 1:N_scene
        @constraint(m, sum(f[i, :]) + f_outlier[i] == 1)
    end

    # Coefficients must be affine coefficients if outlier isn't taken
    for i in 1:N_scene
        @constraint(m, sum(C[i, :]) + f_outlier[i] == 1)
    end

    # Build vertex-membership matrix F
    # where F[i, j] == 1 -->
    #         vertex j is a member of face i
    F = Array{Float64}(size(faces, 1), size(points, 1))
    F .= 0.0
    for i = 1:size(faces, 1)
        for j = faces[i]
            F[i, j] = 1
        end
    end

    # Force members of C to be 0 if their respective
    # face is not active
    for i in 1:N_scene, j in 1:N_model_verts
        @constraint(m, C[i, j] <= dot(F[:, j], f[i, :]))
    end

    # Convert points into modelVerts, which is just a
    # 3xN array of the same data. (Allows vectorization
    # in just a bit.)
    modelVerts = Array{Float64}(3, size(points, 1))
    for i=1:size(points, 1)
        for k=1:3
            modelVerts[k, i] = points[i][k]
        end
    end

    # Create L1-norm-helper-members alpha for computing
    # the elementwise absolute value of the correspondence error
    @variable(m, alpha[1:3, 1:N_scene] >= 0)

    for i in 1:N_scene
        # Phi = Sum of Abs of Error X/Y/Z Components
        @constraint(m, phi[i] == dot([1; 1; 1],  alpha[:, i]))

        l1ErrorPos = R * scenePoints[:, i] + T - modelVerts * C[i, :]
        selector = optBigNumber * (f_outlier[i])
        @constraint(m, alpha[:, i] .>= l1ErrorPos - selector)
        @constraint(m, alpha[:, i] .>= -l1ErrorPos - selector)
    end

    @objective(m, Min, sum(phi)/N_scene + sum(f_outlier)*optPhiMax/N_scene);

    return PoseEstimationSolveData(m, R, T, R_abs, R_diff,
        phi, alpha, C, f, f_outlier,)
end

end #module

max_stars_repo_name: gizatt/PoseEstimationSandbox

max_issues_repo_name: gizatt/PoseEstimationSandbox

max_forks_repo_name: gizatt/PoseEstimationSandbox


text:
"""
    $(SIGNATURES)

Worker utility function u(c).
Not a `ModelObject`. Constructed from [`WorkerSet`](@ref).
"""
mutable struct WorkerUtility{F1 <: AbstractFloat}
    util :: UtilityFunctionsLH.UtilityOneArg
    discFactor :: F1
    "Fixed utility from working (per period)"
    utilFixed :: F1
end

"""
    $(SIGNATURES)

Worker object for one school group.
Most calibrated parameters are shared across school groups
(except wage).
Worker does not directly contain calibrated params (no ParamVector).
Ages are model ages.
"""
mutable struct Worker{F1 <: AbstractFloat}
    # Utility from consumption
    util :: WorkerUtility
    "Experience log wage profile. Zero intercept."
    xp :: Vector{F1}
    # Earnings are `wage * h`
    wage :: F1
    "Gross interest rate"
    R :: F1
    "First model age of retirement"
    retireAge :: TimeInt
    retireDuration :: TimeInt
    retireIncome :: F1
    # Bounds for asset and human capital values
    kMin :: F1
    kMax :: F1
    hMin :: F1
    hMax :: F1
end

# Alias for vector with worker object for each school level
# Can now be indexed with integers or with EdLevels (e.g. HSG)
# const WorkerVector  = VectorBySchool{Worker}

# """
# 	$(SIGNATURES)

# Switches that define a `WorkerSet`.
# """
# Base.@kwdef mutable struct WorkerSetSwitches
#     # Does the model have fixed utility from working by schooling?
#     hasUFixed :: Bool = true
#     retireAge :: TimeInt = 65 - 19
#     retireDuration :: TimeInt = 80 - 65
#     nSchool :: Int
# end

# """
# 	$(SIGNATURES)

# ModelObject with calibrated worker parameters.
# """
# Base.@kwdef mutable struct WorkerSet{F1 <: AbstractFloat} <: ModelObject
#     objId :: ObjectId
#     switches :: WorkerSetSwitches
#     "Log Wage by schooling"
#     wages :: IncreasingVector
#     "Experience log wage profile. Zero intercept."
#     xp :: Vector{F1}
# end

"""
	$(SIGNATURES)

Simulated work histories.
"""
mutable struct WorkHistories{F1 <: AbstractFloat}
    iSchoolV :: Vector{UInt8}
    "Age of work start"
    workStartAgeV :: Vector{TimeInt}
    hWorkStartV :: Vector{F1}
    kWorkStartV :: Vector{F1}
    "Earnings at work start"
    earnWorkStartV :: Vector{F1}
    "Lifetime earnings, discounted to work start"
    ltEarnV :: Vector{F1}
    "Marginal utility of wealth as of work start"
    muWealthV :: Vector{F1}
    "Consumption at work start"
    cWorkStartV :: Vector{F1}
end

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


text:
function build_drone!(vis::Visualizer, model::DoubleIntegratorGame; name::Symbol=:Robot, r_col=0.08, r_cost=0.5, α=1.0, r=0.15)
	r = convert(Float32, r)
    p = model.p

	obj_path = joinpath(module_dir(), "resources", "object", "drone.obj")
	orange_mat, blue_mat, black_mat_col = get_material(;α=0.3)
	orange_mat, blue_mat, black_mat_cost = get_material(;α=0.1)
	obj = MeshFileObject(obj_path)
    for i = 1:p
		setobject!(vis[name]["player$i"]["body"], obj)
		setobject!(vis[name]["player$i"]["collision"], GeometryBasics.Sphere(Point3f0(0.0), r_col), black_mat_col)
		setobject!(vis[name]["player$i"]["cost"], GeometryBasics.Sphere(Point3f0(0.0), r_cost), black_mat_cost)
		settransform!(vis[name]["player$i"]["body"], MeshCat.LinearMap(r*RotX(pi/2)))
    end
    return nothing
end

function set_drone!(vis::Visualizer, model::DoubleIntegratorGame, s::AbstractVector; name::Symbol=:Robot)
    p = model.p
    pz = model.pz
	d = Algames.dim(model)
    for i = 1:p
		x = [s[pz[i][1:d]]; zeros(3-d)]
        settransform!(vis[name]["player$i"], MeshCat.Translation(x...))
    end
    return nothing
end


text:
module MockDofsTests

using Test
using Gridap.Fields
using Gridap.Fields: MockField
using Gridap.ReferenceFEs
using Gridap.ReferenceFEs: MockDofBasis

# np = 4
# p = Point(1,2)
# x = fill(p,np)

# v = 3.0
# d = 2
# f = MockField{d}(v)
# fx = evaluate(f,x)

# b = MockDofBasis(x)
# test_dof(b,f,fx)

# l = 10
# ab = fill(b,l)
# af = fill(f,l)
# ax = fill(x,l)

# afx = evaluate(af,ax)
# abf = evaluate(ab,af)
# @test abf == afx

end # module


text:
include("./Tools.jl")
includeall("./Calculators")


text:
"""
"""
module PredictMD_MLJ # begin module PredictMD_MLJ

import ..MLJ
import ..MLJBase
import ..MLJModels

end # end module PredictMD_MLJ


text:
function DataFrame(ta::TimeArray; colnames=Symbol[], timestamp=:Date)
    if length(colnames) != 0
        colnames_str = [string(s) for s in colnames]
        ta = ta[colnames_str...]
    end
    df = DataFrame(hcat(ta.timestamp, ta.values))
    colnames = [Symbol(s) for s in ta.colnames]
    colnames = vcat(timestamp, colnames)
    names!(df, colnames)
    df
end


text:
#==========================================================
This routine initializes an optimization problem to recover
a CNOT gate on a coupled 2-qubit system with 2 energy
levels on each oscillator (with 1 guard state on one and
2 guard states on the other). The drift Hamiltonian in
the rotating frame is
    H_0 = - 0.5*ξ_a(a^†a^†aa)
          - 0.5*ξ_b(b^†b^†bb)
          - ξ_{ab}(a^†ab^†b),
where a,b are the annihilation operators for each qubit.
Here the control Hamiltonian in the rotating frame
includes the usual symmetric and anti-symmetric terms
H_{sym,1} = p_1(t)(a + a^†),  H_{asym,1} = q_1(t)(a - a^†),
H_{sym,2} = p_2(t)(b + b^†),  H_{asym,2} = q_2(t)(b - b^†).
The problem parameters for this example are,
            ω_a    =  2π × 4.10595   Grad/s,
            ξ_a    =  2π × 2(0.1099) Grad/s,
            ω_b    =  2π × 4.81526   Grad/s,
            ξ_b    =  2π × 2(0.1126) Grad/s,
            ξ_{ab} =  2π × 0.1       Grad/s,
We use Bsplines with carrier waves with frequencies
0, ξ_a, 2ξ_a Grad/s for each oscillator.
==========================================================#
using LinearAlgebra
using Ipopt
using Base.Threads
using Random
using DelimitedFiles
using Printf
using FFTW
using Plots
using SparseArrays

Base.show(io::IO, f::Float64) = @printf(io, "%20.13e", f)

using Juqbox

# function my_initial_cond(Ntot, N, Ne, Ng)
#     Ident = Matrix{Float64}(I, Ntot, Ntot)
#     U0 = Ident[1:Ntot, 1:N] # Rectangular subset of identity
#     # adjust the initial guess
#     if Ng[1] + Ng[2] > 0
#         Nt = Ne + Ng
#         # build up a basis for the essential states
#         col = 0
#         m = 0
#         for k2 in 1:Nt[2]
#             for k1 in 1:Nt[1]
#                 m += 1
#                 # is this a guard level?
#                 guard = (k1 > Ne[1]) || (k2 > Ne[2])
#                 if !guard
#                     col += 1
#                     U0[:,col] = Ident[:,m]
#                 end # if ! guard
#             end # for
#         end # for
#     end # if
#     return U0
# end

Ne1 = 2 # essential energy levels per oscillator
Ne2 = 2
Ng1 = 0 # 0 # Osc-1, number of guard states
Ng2 = 0 # 0 # Osc-2, number of guard states

Ne = [Ne1, Ne2]
Ng = [Ng1, Ng2]

N = Ne1*Ne2; # Total number of nonpenalized energy levels
Ntot = (Ne1+Ng1)*(Ne2+Ng2)
Nguard = Ntot - N # total number of guard states

Nt1 = Ne1 + Ng1
Nt2 = Ne2 + Ng2

Tmax = 100.0 # Duration of gate

# frequencies (in GHz, will be multiplied by 2*pi to get angular frequencies in the Hamiltonian matrix)
fa = 4.10595    # official
fb = 4.81526   # official
rot_freq = [fa, fb] # rotational frequencies
x1 = 2* 0.1099  # official
x2 = 2* 0.1126   # official
x12 = 0.1 # Artificially large to allow coupling. Actual value: 1e-6

# construct the lowering and raising matricies: amat, bmat, cmat
# and the system Hamiltonian: H0
#

# Note: The ket psi = ji> = e_j kron e_i.
# We order the elements in the vector psi such that i varies the fastest with i in [1,Nt1] and j in [1,Nt2]
# The matrix amat = I kron a1 acts on alpha in psi = beta kron alpha
# The matrix bmat = a2 kron I acts on beta in psi = beta kron alpha
a1 = Array(Bidiagonal(zeros(Nt1),sqrt.(collect(1:Nt1-1)),:U))
a2 = Array(Bidiagonal(zeros(Nt2),sqrt.(collect(1:Nt2-1)),:U))

I1 = Array{Float64, 2}(I, Nt1, Nt1)
I2 = Array{Float64, 2}(I, Nt2, Nt2)

# create the a, a^\dag, b and b^\dag vectors
amat = kron(I2, a1)
bmat = kron(a2, I1)

adag = Array(transpose(amat))
bdag = Array(transpose(bmat))

# number ops
num1 = Diagonal(collect(0:Nt1-1))
num2 = Diagonal(collect(0:Nt2-1))

# number operators
N1 = Diagonal(kron(I2, num1) )
N2 = Diagonal(kron(num2, I1) )

# System Hamiltonian
H0 = -2*pi*( x1/2*(N1*N1-N1) + x2/2*(N2*N2-N2) + x12*(N1*N2) )

# max coefficients, rotating frame
amax = 0.02 # max amplitude ctrl func for Hamiltonian #1
bmax = 0.025 # max amplitude ctrl func for Hamiltonian #2
maxpar = [amax, bmax]

# estimate max magnitude of eigenvalue
K1 =  H0 +
    (amax.*(amat +  amat') + 1im*amax.*(amat -  amat') +
     bmax.*(bmat + bmat') + 1im*bmax.*(bmat - bmat'))
lamb = eigvals(K1)
maxeig = maximum(abs.(lamb))

# Estimate time step
Pmin = 40 # should be 20 or higher
samplerate1 = maxeig*Pmin/(2*pi)
nsteps = ceil(Int64,Tmax*samplerate1)

println("Number of time steps = ", nsteps)

# package the lowering and raising matrices together into an one-dimensional array of two-dimensional arrays
# Here we choose dense or sparse representation
use_sparse = true
# use_sparse = false

# dense matrices run faster, but take more memory
Hsym_ops=[Array(amat+adag), Array(bmat+bdag)]
Hanti_ops=[Array(amat-adag), Array(bmat-bdag)]
H0 = Array(H0)

Nfreq = 2 # number of carrier frequencies

Ncoupled = length(Hsym_ops)
om = zeros(Ncoupled, Nfreq) # Allocate space for the carrier wave frequencies

@assert(Nfreq==1 || Nfreq==2 || Nfreq==3)
if Nfreq == 2
    om[1:Ncoupled,2] .= -2.0*pi*x12 # coupling freq for both ctrl funcs (re/im)
elseif Nfreq == 3
    om[1,2] = -2.0*pi*x1 # 1st ctrl, re
    om[2,2] = -2.0*pi*x2 # 2nd ctrl, re
    om[1:Ncoupled,3] .= -2.0*pi*x12 # coupling freq for both ctrl funcs (re/im)
end

println("Carrier frequencies 1st ctrl Hamiltonian [GHz]: ", om[1,:]./(2*pi))
println("Carrier frequencies 2nd ctrl Hamiltonian [GHz]: ", om[2,:]./(2*pi))

#U0 = my_initial_cond(Ntot, N, Ne, Ng)
U0 = initial_cond(Ne, Ng)

# SWAP target for the essential levels
gate_swap =  zeros(ComplexF64, N, N)
gate_swap[1,1] = 1.0
gate_swap[2,3] = 1.0
gate_swap[3,2] = 1.0
gate_swap[4,4] = 1.0

utarget = U0 * gate_swap

# rotation matrices
omega1, omega2 = Juqbox.setup_rotmatrices(Ne, Ng, rot_freq)

# Compute Ra*Rb*utarget
rot1 = Diagonal(exp.(im*omega1*Tmax))
rot2 = Diagonal(exp.(im*omega2*Tmax))

vtarget = rot1*rot2*utarget # target in the rotating frame

# assemble problem description for the optimization
params = Juqbox.objparams(Ne, Ng, Tmax, nsteps, Uinit=U0, Utarget=vtarget, Cfreq=om, Rfreq=rot_freq,
                          Hconst=H0, Hsym_ops=Hsym_ops, Hanti_ops=Hanti_ops, use_sparse=use_sparse)

# initial parameter guess
startFromScratch = true # false
startFile = "swap2-pcof-opt.dat"

# dimensions for the parameter vector
D1 = 10 # number of B-spline coeff per oscillator, freq and sin/cos

nCoeff = 2*Ncoupled*Nfreq*D1 # Total number of parameters.

Random.seed!(2456)
if startFromScratch
  pcof0 = amax*0.01 * rand(nCoeff)
  println("*** Starting from pcof with random amplitude ", amax*0.01)
else
  pcof0 = vec(readdlm(startFile))
  println("*** Starting from B-spline coefficients in file: ", startFile)

  nCoeff = length(pcof0)
  D1 = div(nCoeff, 2*Ncoupled*Nfreq)  # number of B-spline coeff per control function
end

samplerate = 32 # for plotting
casename = "cnot2" # for constructing file names

# min and max B-spline coefficient values
useBarrier = true
minCoeff, maxCoeff = Juqbox.assign_thresholds(params,D1,maxpar)
println("Number of min coeff: ", length(minCoeff), "Max Coeff: ", length(maxCoeff))

maxIter = 150 # 0 #250 #50 # optional argument
lbfgsMax = 250 # optional argument

println("*** Settings ***")
# output run information
println("Frequencies: fa = ", fa, " fb = ", fb)
println("Coefficients in the Hamiltonian: x1 = ", x1, " x2 = ", x2, " x12 = ", x12)
println("Essential states in osc = ", [Ne1, Ne2], " Guard states in osc = ", [Ng1, Ng2])
println("Total number of states, Ntot = ", Ntot, " Total number of guard states, Nguard = ", Nguard)
println("Number of B-spline parameters per spline = ", D1, " Total number of parameters = ", nCoeff)
println("Max parameter amplitudes: maxpar = ", maxpar)
println("Tikhonov regularization tik0 (L2) = ", params.tik0)
if use_sparse
    println("Using a sparse representation of the Hamiltonian matrices")
else
    println("Using a dense representation of the Hamiltonian matrices")
end

params.save_pcof_hist = true

# Allocate all working arrays
wa = Juqbox.Working_Arrays(params, nCoeff)
prob = Juqbox.setup_ipopt_problem(params, wa, nCoeff, minCoeff, maxCoeff, maxIter=maxIter, lbfgsMax=lbfgsMax, startFromScratch=startFromScratch)

# uncomment to run the gradient checker for the initial pcof
# addOption( prob, "derivative_test", "first-order"); # for testing the gradient
# addOption(prob, "print_level", 0); # for testing the gradient

println("Initial coefficient vector stored in 'pcof0'")