text:
immutable SpriteElem{N,T}
    offset::Vec{2,T}
    uv_offset_width::Vec{4,T}
    position::Point{N, T}
    color::RGBA{Float32}
    scale::Vec{2,T}
end
type Text
    offsets
    uv_offset_width
    positions
    colors
    scales
    rotations
    atlas
end
Base.length(rt::Text) = length(rt.positions)
function Base.sizehint!(t::Text, size::Int)
    for elem in (:positions, :colors, :scales, :uv_offset_width, :offsets, :rotations)
        field = getfield(t, elem)
        if isa(field, AbstractArray) # all fields can be scalars
            sizehint!(field, size)
        end
    end
    nothing
end
function Base.resize!(t::Text, size::Int)
    for elem in (:positions, :colors, :scales, :uv_offset_width, :offsets, :rotations)
        field = getfield(t, elem)
        if isa(field, AbstractArray)
            resize!(field, size)
        end
    end
    nothing
end
function Base.setindex!(t::Text, value::NTuple{6}, i::Integer)
    char, offset, position, color, scale, rotation = value
    gs = GLVisualize.glyph_scale!(t.atlas, char, GLVisualize.DEFAULT_FONT_FACE)
    guvw = GLVisualize.glyph_uv_width!(t.atlas, char, GLVisualize.DEFAULT_FONT_FACE)
    t.offsets[i] = offset
    t.uv_offset_width[i] = guvw
    t.positions[i] = position
    t.colors[i] = color
    t.scales[i] = gs .* Vec2f0(scale)
    t.rotations[i] = rotation
    nothing
end

function Base.setindex!(t::Text, value::Char, i::Integer)
    gs = GLVisualize.glyph_scale!(rt.text.atlas, c, GLVisualize.DEFAULT_FONT_FACE)
    guvw = GLVisualize.glyph_uv_width!(rt.text.atlas, c, GLVisualize.DEFAULT_FONT_FACE)
    t.uv_offset_width[i] = guvw
    t.scales[i] = gs
end

function Base.setindex!(t::Text, value::Associative, i::Integer)
    for (k,v) in value
        getfield(t, k)[i] = v
    end
end
function Base.setindex!(t::Text, value::SpriteElem, i::Integer)
    for f=1:5
        getfield(t, f)[i] = getfield(value, f)
    end
end
function Text(robj::RenderObject)
    @materialize atlas, offset, uv_offset_width, position, color, scale, rotation = robj.uniforms
    text = Text(
        offset,
        uv_offset_width,
        position,
        color,
        scale,
        rotation,
        atlas
    )
end
type RichText
    text::Text
    inlined_objects::Vector
    links::Vector
    cursor
    defaults::Dict
end
function RichText(text::Text)
    style_default = Dict(
        :color => RGBA{Float32}(0,0,0,1),
        :scale => Vec2f0(1),
        :font => GLVisualize.DEFAULT_FONT_FACE
    )
    RichText(text, [], [], 0:1, style_default)
end

function Base.sizehint!(t::RichText, size::Int)
    sizehint!(t.text, size)
    nothing
end
function Base.resize!(t::RichText, size::Int)
    resize!(t.text, size)
    nothing
end
Base.length(rt::RichText) = length(rt.text)
make_iter(x) = repeated(x)
make_iter(x::AbstractArray) = x

function get_iter(defaultfunc, dictlike, key)
    make_iter(get(defaultfunc, dictlike, key))
end

function getposition(text, text2, fonts, scales, start_pos)
    GLVisualize.calc_position(text2, start_pos, scales, fonts, text.text.atlas)
end
function getoffsets(text, text2, fonts, scales)
    GLVisualize.calc_offset(text2, scales, fonts, text.text.atlas)
end
font_color(text::RichText) = text.defaults[:color]
font(text::RichText) = text.defaults[:font]
font_size(text::RichText) = text.defaults[:scale]
"""
Inserts text at [`position`] with `style`
"""
function Base.insert!(text::RichText, text2::AbstractString, position=text.cursor, style=text.defaults)
    fonts     = get_iter(()->font(text), style, :font)
    scales    = get_iter(()->font_size(text), style, :scale)
    start_pos = get(style, :start_position, text.text.positions[max(last(text.cursor), 1)])
    positions = get_iter(()->getposition(text, text2, fonts, scales, start_pos), style, :position)
    offsets   = get_iter(()->getoffsets(text, text2, fonts, scales), style, :offset)
    colors    = get_iter(()->font_color(text), style, :color)
    rotation  = get_iter(()->Vec3f0(0,0,1), style, :rotation)

    pos = max(last(text.cursor), 1)
    if pos + length(text2) > length(text)
        resize!(text, pos + length(text2))
    end
    for elem in zip(text2, offsets, positions, colors, scales, rotation)
        text.text[pos] = elem
        pos += 1
    end
    nothing
end

iter_or_array(x) = repeated(x)
iter_or_array(x::Base.Repeated) = x
iter_or_array(x::Array) = x
iter_or_array(x::Vector{Ptr{FreeType.FT_FaceRec}}) = repeated(x)

function calc_position(last_pos, start_pos, atlas, glyph, font, scale)
    advance_x, advance_y = glyph_advance!(atlas, glyph, font, scale)
    if isnewline(glyph)
        return Point2f0(start_pos[1], last_pos[2]-advance_y) #reset to startx
    else
        return last_pos + Point2f0(advance_x, 0)
    end
end
function calc_position(glyphs, start_pos, scales, fonts, atlas)
    positions = zeros(Point2f0, length(glyphs))
    last_pos  = Point2f0(start_pos)
    s, f = iter_or_array(scales), iter_or_array(fonts)
    c1 = first(glyphs)
    for (i, (c2, scale, font)) in enumerate(zip(glyphs, s, f))
        c2 == '\r' && continue # stupid windows!
        positions[i] = last_pos
        last_pos = calc_position(last_pos, start_pos, atlas, c2, font, scale)
        #k = FreeTypeAbstraction.kerning(c1, c2, font, scale)
        #last_pos += k # add kerning
    end
    positions
end
function calc_offset(glyphs, scales, fonts, atlas)
    offsets = fill(Point2f0(0.0), length(glyphs))
    s, f = iter_or_array(scales), iter_or_array(fonts)
    c1 = first(glyphs)
    for (i, (c2, scale, font)) in enumerate(zip(glyphs, s, f))
        c2 == '\r' && continue # stupid windows!
        offsets[i] = Point2f0(glyph_bearing!(atlas, c2, font, scale))
        c1 = c2
    end
    offsets # bearing is the glyph offset
end

isnewline(x) = x == '\n'
# i must be a valid character index
function next_newline(text, i::Integer)
	res = findnext(isnewline, text, i)
	res == 0 ? length(text) : res #jump to end of text if no newline found
end
previous_newline(text, i::Integer) = max(1, findprev(isnewline, text, i)) #jump to start

export previous_newline, next_newline

#=
textextetext\n
texttext<current pos>texttext\n
texttext<finds this pos>text\n
=#
function down_after_newline(text, current_position)
	i 	= current_position
	pnl = previous_newline(text, i)
	nnl = next_newline(text, i)
	nl_distance = i-pnl # distance from previous newline
	min(length(text), nnl+nl_distance)
end
#=
textexte<finds this pos>text\n
texttext<current pos>texttext\n
texttexttext\n
=#
function up_before_newline(text, current_position)
	i 	 = current_position
	pnl  = previous_newline(text, i)
	ppnl = previous_newline(text, max(1, pnl-1))
	nl_distance = i-pnl # distance from previous newline
	min(length(text), ppnl+nl_distance)
end

function next_arrow_selection(dir, text, first_index, last_index)
	(dir == :up)    && return up_before_newline(text, first_index)
	(dir == :down)  && return down_after_newline(text, last_index)
	(dir == :left)  && return max(1, first_index-1)
	(dir == :right) && return min(length(text)+1, last_index+1)
    -1
end
function move_cursor(t0, dir, mouseselection, text_selection, is_movingselection)
	text, selection = text_selection.text, text_selection.selection
	mouseselection0, selection0 = t0
	selection0 = selection
	mouseselection0 != mouseselection && return (mouseselection, mouseselection) # if mouse selection has changed, return the new position
	if selection0 != 0:0
		first_i   = first(selection0) # first is always valid, if its not zero
		# last is not valid, if selection is in between characters
		last_i    = isempty(selection0) ? first(selection0) : last(selection0) # if only single char selected use first, otherwise last position of selection
		nas       = next_arrow_selection(dir, text, first_i, last_i) # -1 for no new arrow selection
        if nas != -1
            if is_movingselection
                return (mouseselection, nas:last_i)
            else
                return (mouseselection, nas:0)
            end
        end
	end
	(mouseselection0, selection0)
end
#
#
# immutable TextCursor
#     selections::Vector{UnitRange}
# end
# type Text{T <: Integer}
#     text::Vector{T}
# end
# function Base.setindex!(text, index::TextCursor, newtext)
#     for selection in index.selections
#
#     end
# end
deletetext(_, text_selection) = deletetext(text_selection)
function deletetext(text_selection)
	selection, text = text_selection.selection, text_selection.text
	offset = 0
	if first(selection) > 0 && last(selection) > 0
		if single_selection(selection)
			splice!(text, last(selection))
			offset = -1
		else
			splice!(text, selection)
		end
		text_selection.selection = max(1, first(selection)+offset):0 # when text gets removed, selection will turn into single selection
	end
	nothing
end
export inserttext

function inserttext(text_selection, to_insert)
	selection, text = text_selection.selection, text_selection.text
	if first(selection) > 0 && first(selection) <= length(text)
		splice!(text, selection, to_insert)
		chars_added = length(to_insert)
		text_selection.selection = (first(selection)+chars_added):0 # when text gets removed, selection will turn into single selection
	end
	nothing
end

function insert_text(text, cursor, new_text)
    text[cursor] = new_text
    push!(text, text)
    push!(cursor, text)
end

function selection_click(index)

    push!(cursor, )
end
function selection_drag(index)
    push!(cursor.selection, index)
    push!(cursor, )
end
function selection_add(cursor, selection)
    push!(cursor.selection, selection)
end
function selection_reset(cursor)
    push!(cursor)
end

function textedit_signals(inputs, background, text)
    @materialize unicodeinput, selection, buttonspressed, arrow_navigation, mousedragdiff_objectid = inputs
    # create object which can globally hold the text and selection
    text_raw    = TextWithSelection(text[:glyphs], 0:0)
    text_edit   = Signal(text_raw)
    shift       = const_lift(in, GLFW.KEY_LEFT_SHIFT, buttonspressed)
    selection   = preserve(const_lift(
        last,
        foldp(
            move_cursor,
            (selection.value, selection.value),
            arrow_navigation, selection,
            text_edit,
            shift
        )
    ))

    is_text(x) = x[2][1] == background.id || x[2][1] == text.id
    selection  = filterwhen(
        const_lift(is_text, mousedragdiff_objectid),
        0:0, selection
    )
    preserve(const_lift(s->(text_edit.value.selection=s), selection)) # is there really no other way?!

    strg_v          = const_lift(==, buttonspressed, Config.shortcuts["paste"])
    strg_c          = const_lift(==, buttonspressed, Config.shortcuts["copy"])
    strg_x          = const_lift(==, buttonspressed, Config.shortcuts["cut"])
    enter_key       = const_lift(==, buttonspressed, [GLFW.KEY_ENTER])
    del             = const_lift(==, buttonspressed, [GLFW.KEY_BACKSPACE])

    enter_insert    = const_lift(insert_enter,   filterwhen(enter_key, true, enter_key))
    clipboard_copy  = const_lift(copyclipboard,  filterwhen(strg_c,    true, strg_v), text_edit)

    delete_text     = const_lift(deletetext,     filterwhen(del,       true, del),    text_edit)
    cut_text        = const_lift(cutclipboard,   filterwhen(strg_x,    true, strg_x), text_edit)

    clipboard_paste = const_lift(clipboardpaste, filterwhen(strg_v,    true, strg_v))

    text_gate       = const_lift(isnotempty, unicodeinput)
    unicode_input   = filterwhen(text_gate, Char['0'], unicodeinput)
    text_to_insert  = merge(clipboard_paste, unicode_input, enter_insert)
    text_to_insert  = const_lift(process_for_gl, text_to_insert)

    text_inserted   = const_lift(inserttext, text_edit, text_to_insert)

    text_updates    = merge(
        const_lift(return_nothing, text_inserted),
        const_lift(return_nothing, clipboard_copy),
        const_lift(return_nothing, delete_text),
        const_lift(return_nothing, cut_text),
        const_lift(return_nothing, selection)
    )
    text_selection_signal = sampleon(text_updates, text_edit)

    selection   = const_lift(x->x.selection,  text_selection_signal)
    text_sig    = const_lift(x->x.text,       text_selection_signal)

    preserve(const_lift(update_positions, text_sig, Signal(text), Signal(background[:style_index])))
    preserve(foldp(visualize_selection, 0:0, selection,    Signal(background[:style_index])))
    const_lift(Compat.String, text_sig), selection
end
# # i must be a valid character index
# function next_newline(text, i::Integer)
# 	res = findnext(isnewline, text, i)
# 	res == 0 ? length(text) : res
# end
# previous_newline(text, i::Integer) = max(1, findprev(isnewline, text, i))
#
# export previous_newline
# export next_newline
#
#
#
# #=
# textextetext\n
# texttext<current pos>texttext\n
# texttext<finds this pos>text\n
# =#
# function down_after_newline(text, current_position)
# 	i 	= current_position
# 	pnl = previous_newline(text, i)
# 	nnl = next_newline(text, i)
# 	nl_distance = i-pnl # distance from previous newline
# 	min(length(text), nnl+nl_distance)
# end
# #=
# textexte<finds this pos>text\n
# texttext<current pos>texttext\n
# texttexttext\n
# =#
# function up_before_newline(text, current_position)
# 	i 	 = current_position
# 	pnl  = previous_newline(text, i)
# 	ppnl = previous_newline(text, max(1, pnl-1))
# 	nl_distance = i-pnl # distance from previous newline
# 	min(length(text), ppnl+nl_distance)
# end
#
# function next_arrow_selection(dir, text, first_index, last_index)
# 	(dir == :up)    && return up_before_newline(text, first_index)
# 	(dir == :down)  && return down_after_newline(text, last_index)
# 	(dir == :left)  && return max(1, first_index-1)
# 	(dir == :right) && return min(length(text)+1, last_index+1)
#     -1
# end
#
# function move_cursor(t0, dir, mouseselection, text_selection, is_movingselection)
# 	text, selection = text_selection.text, text_selection.selection
# 	mouseselection0, selection0 = t0
# 	selection0 = selection
# 	mouseselection0 != mouseselection && return (mouseselection, mouseselection) # if mouse selection has changed, return the new position
# 	if selection0 != 0:0
# 		first_i   = first(selection0) # first is always valid, if its not zero
# 		# last is not valid, if selection is in between characters
# 		last_i    = isempty(selection0) ? first(selection0) : last(selection0) # if only single char selected use first, otherwise last position of selection
# 		nas       = next_arrow_selection(dir, text, first_i, last_i) # -1 for no new arrow selection
#         if nas != -1
#             if is_movingselection
#                 return (mouseselection, nas:last_i)
#             else
#                 return (mouseselection, nas:0)
#             end
#         end
# 	end
# 	(mouseselection0, selection0)
# end
# export move_cursor
#
#
# function visualize_selection(
# 		last_selection::UnitRange{Int},
# 		selection 	  ::UnitRange{Int},
# 		style 		  ::GPUVector{GLSpriteStyle}
# 	)
# 	fl, ll =  first(last_selection), last(last_selection)
# 	if !isempty(last_selection) && fl > 0 && ll > 0 && (fl <= length(style)) && (ll <= length(style))
# 		style[last_selection] = fill(GLSpriteStyle(0,0), length(last_selection))
# 	end
# 	fs, ls =  first(selection), last(selection)
# 	if !isempty(selection) && fs > 0 && ls > 0 && (fs <= length(style)) && (ls <= length(style))
# 		style[selection] = fill(GLSpriteStyle(1,0), length(selection))
# 	end
# 	selection
# end
#
#
# AND(a,b) 	  = a&&b
# isnotempty(x) = !isempty(x)
# return_nothing(x...) = nothing
# export AND
# export isnotempty
# export return_nothing
#
# single_selection(selection::UnitRange) 	= isempty(selection) && first(selection)!=0
# is_textinput_modifiers(buttons::Vector{Int}) = isempty(buttons) || buttons == [GLFW.KEY_LEFT_SHIFT]
#
# function clipboardpaste(_)
# 	clipboard_data = ""
# 	try
# 		clipboard_data = clipboard()
# 	catch e # clipboard throws error when there is no data (WTF)
# 	end
# 	return Compat.String(clipboard_data)
# end
#
# export clipboardpaste
# export copyclipboard
#
# function back2julia(x::GLSprite)
# 	isnewline(x[1]) && return '\n'
# 	ID_TO_CHAR[x[1]]
# end
# function Base.Compat.String(v::GPUVector{GLSprite})
# 	data = gpu_data(v)
# 	Compat.String(join(map(back2julia, data)))
# end
# # const_lift will have a boolean value at first argument position
# copyclipboard(_, text_selection) = copyclipboard(text_selection)
# function copyclipboard(text_selection)
# 	selection, text = text_selection.selection, text_selection.text
# 	if first(selection) > 0
# 		if single_selection(selection) # for single selection we do a sublime style line copy
# 			#i 	= chr2ind(text, first(selection))
# 			i 	= min(length(text), first(selection)) # can be on position behind last character
# 			pnl = previous_newline(text, i)
# 			nnl = next_newline(text, i)
# 			tocopy = text[pnl:nnl]
# 		else # must be range selection
# 			tocopy = text[selection]
# 		end
# 		clipboard(join(map(x->ID_TO_CHAR[x[1]], tocopy)))
# 	end
# 	nothing
# end
# export cutclipboard
# cutclipboard(_, text_selection) = cutclipboard(text_selection)
# function cutclipboard(text_selection)
# 	copyclipboard(text_selection)
# 	deletetext(text_selection)
# 	nothing
# end
# export deletetext
# deletetext(_, text_selection) = deletetext(text_selection)
# function deletetext(text_selection)
# 	selection, text = text_selection.selection, text_selection.text
# 	offset = 0
# 	if first(selection) > 0 && last(selection) > 0
# 		if single_selection(selection)
# 			splice!(text, last(selection))
# 			offset = -1
# 		else
# 			splice!(text, selection)
# 		end
# 		text_selection.selection = max(1, first(selection)+offset):0 # when text gets removed, selection will turn into single selection
# 	end
# 	nothing
# end
# export inserttext
# inserttext(_, text_selection) = inserttext(text_selection)
# function inserttext(text_selection, to_insert)
# 	selection, text = text_selection.selection, text_selection.text
# 	if first(selection) > 0 && first(selection) <= length(text)
# 		splice!(text, selection, to_insert)
# 		chars_added = length(to_insert)
# 		text_selection.selection = (first(selection)+chars_added):0 # when text gets removed, selection will turn into single selection
# 	end
# 	nothing
# end
#
# type TextWithSelection{S } #<: AbstractString}
# 	text::S
# 	selection::UnitRange{Int}
# end
# export TextWithSelection
#
#
#
#
# export visualize_selection


text:
export AtmosFilterPerturbations
export AtmosSpecificFilterPerturbations

struct AtmosFilterPerturbations{M} <: AbstractFilterTarget
    atmos::M
end

vars_state_filtered(target::AtmosFilterPerturbations, FT) =
    vars_state(target.atmos, Prognostic(), FT)

function compute_filter_argument!(
    target::AtmosFilterPerturbations,
    filter_state::Vars,
    state::Vars,
    aux::Vars,
)
    # copy the whole state
    parent(filter_state) .= parent(state)
    # remove reference state
    filter_state.ρ -= aux.ref_state.ρ
    filter_state.energy.ρe -= aux.ref_state.ρe
    if !(moisture_model(target.atmos) isa DryModel)
        filter_state.moisture.ρq_tot -= aux.ref_state.ρq_tot
    end
    if (moisture_model(target.atmos) isa NonEquilMoist)
        filter_state.moisture.ρq_liq -= aux.ref_state.ρq_liq
        filter_state.moisture.ρq_ice -= aux.ref_state.ρq_ice
    end
end
function compute_filter_result!(
    target::AtmosFilterPerturbations,
    state::Vars,
    filter_state::Vars,
    aux::Vars,
)
    # copy the whole filter state
    parent(state) .= parent(filter_state)
    # add reference state
    state.ρ += aux.ref_state.ρ
    state.energy.ρe += aux.ref_state.ρe
    if !(moisture_model(target.atmos) isa DryModel)
        state.moisture.ρq_tot += aux.ref_state.ρq_tot
    end
    if (moisture_model(target.atmos) isa NonEquilMoist)
        filter_state.moisture.ρq_liq += aux.ref_state.ρq_liq
        filter_state.moisture.ρq_ice += aux.ref_state.ρq_ice
    end
end

struct AtmosSpecificFilterPerturbations{M} <: AbstractFilterTarget
    atmos::M
end

vars_state_filtered(target::AtmosSpecificFilterPerturbations, FT) =
    vars_state_filtered(target.atmos, FT)

function compute_filter_argument!(
    target::AtmosSpecificFilterPerturbations,
    filter_state::Vars,
    state::Vars,
    aux::Vars,
)

    ρ_inv = 1 / state.ρ
    ρ_ref_inv = 1 / aux.ref_state.ρ

    # copy the whole state
    parent(filter_state) .= parent(state) * ρ_inv

    # remove reference state
    filter_state.energy.e -= aux.ref_state.ρe * ρ_ref_inv
    if !(moisture_model(target.atmos) isa DryModel)
        filter_state.moisture.q_tot -= aux.ref_state.ρq_tot * ρ_ref_inv
    end
    if (moisture_model(target.atmos) isa NonEquilMoist)
        filter_state.moisture.q_liq -= aux.ref_state.ρq_liq * ρ_ref_inv
        filter_state.moisture.q_ice -= aux.ref_state.ρq_ice * ρ_ref_inv
    end

end

function compute_filter_result!(
    target::AtmosSpecificFilterPerturbations,
    state::Vars,
    filter_state::Vars,
    aux::Vars,
)

    ρ = state.ρ
    ρ_ρ_ref_ratio = ρ / aux.ref_state.ρ

    # copy the whole filter state
    parent(state) .= parent(filter_state) * ρ

    # add reference state
    state.energy.ρe += aux.ref_state.ρe * ρ_ρ_ref_ratio
    if !(moisture_model(target.atmos) isa DryModel)
        state.moisture.ρq_tot += aux.ref_state.ρq_tot * ρ_ρ_ref_ratio
    end
    if (moisture_model(target.atmos) isa NonEquilMoist)
        state.moisture.ρq_liq += aux.ref_state.ρq_liq * ρ_ρ_ref_ratio
        state.moisture.ρq_ice += aux.ref_state.ρq_ice * ρ_ρ_ref_ratio
    end
end

ext: jl


text:
    # Functions and methods for cascade computation

    """
    `Cascade.computeDecayYieldOutcome(outcome::DecayYield.Outcome, linesR::Array{PhotoEmission.Line,1},
                                      linesA::Array{AutoIonization.Line,1}, settings::DecayYield.Settings)`
        ... to compute the flourescence and Auger yields for a single decay yield outcome as specified by the corresponding
            level; an outcome::DecayYield.Outcome is returned in which all physical parameters are now specified for the given
            decay-yield.
    """
    function computeDecayYieldOutcome(outcome::DecayYield.Outcome, linesR::Array{PhotoEmission.Line,1},
                                      linesA::Array{AutoIonization.Line,1}, settings::DecayYield.Settings)
        # Identify the level key of the given level also in the lists of radiative and Auger lines
        level = outcome.level;    levelKey = LevelKey( LevelSymmetry(level.J, level.parity), level.index, level.energy, 0.)
        similarKey = LevelKey();  rateR = 0.;    rateA = 0.;   NoPhotonLines = 0;   NoAugerLines = 0
        for  line in linesR
            compareKey = LevelKey( LevelSymmetry(line.initialLevel.J, line.initialLevel.parity), line.initialLevel.index, line.initialLevel.energy, 0.)
            if   Basics.isSimilar(levelKey, compareKey, 1.0e-3)    println("** compareKey = $compareKey");   similarKey = deepcopy(compareKey)    end
        end
        if   similarKey == LevelKey()    error("No similar level found !")   end

        for  line in linesR
            if  similarKey == LevelKey( LevelSymmetry(line.initialLevel.J, line.initialLevel.parity), line.initialLevel.index, line.initialLevel.energy, 0.)
                rateR = rateR + line.photonRate;   NoPhotonLines = NoPhotonLines + 1
            end
        end
        for  line in linesA
            if  similarKey == LevelKey( LevelSymmetry(line.initialLevel.J, line.initialLevel.parity), line.initialLevel.index, line.initialLevel.energy, 0.)
                rateA = rateA + line.totalRate;   NoAugerLines = NoAugerLines + 1
            end
        end

        omegaR     = rateR / (rateR + rateA);   omegaA = rateA / (rateR + rateA)
        newOutcome = DecayYield.Outcome(level, NoPhotonLines, NoAugerLines, rateR, rateA, omegaR, omegaA)
        return( newOutcome )
    end

    """
    `Cascade.computeTotalAugerRate(level::Cascade.Level)`
        ... computes the total Auger rate of level as given by its daugther levels; a rate::Float64 is returned.
    """
    function computeTotalAugerRate(level::Cascade.Level)
        rate = 0.
        for daugther in level.daugthers
            if  daugther.process != Basics.Auger();      continue   end
            aLine   = daugther.lineSet.linesA[daugther.index]
            rate    = rate + aLine.totalRate
        end
        return( rate )
    end

    """
    `Cascade.computeTotalPhotonRate(level::Cascade.Level)`
        ... computes the total photon (radiative) rate of level as given by its daugther levels; a rate::EmProperty is returned.
    """
    function computeTotalPhotonRate(level::Cascade.Level)
        rate = Basics.EmProperty(0.)
        for daugther in level.daugthers
            if  daugther.process != Basics.Radiative();      continue   end
            rLine   = daugther.lineSet.linesR[daugther.index]
            rate    = rate + rLine.photonRate
        end

        return( rate )
    end

    """
    `Cascade.displayBlocks(stream::IO, blockList::Array{Cascade.Block,1}; sa::String="")`
        ... group & display the blocks of the cascade with same No. of electrons; this blocks are displayed with the
            minimum and maximum energy of each multiplet. The optional sa::String can be used to display some details
            about the given blocks. nothing is returned.
    """
    function displayBlocks(stream::IO, blockList::Array{Cascade.Block,1}; sa::String="")
        #
        nx = 150
        println(stream, "\n* Configuration 'blocks' (multiplets): " * sa * "\n")
        println(stream, "  ", TableStrings.hLine(nx))
        println(stream, "      No.   Configurations                                                                       " *
                        "      No. CSF  ",
                        "      Range of total energies " * TableStrings.inUnits("energy") )
        println(stream, "  ", TableStrings.hLine(nx))
        i = 0
        for  block  in  blockList
            i = i + 1;
            sa = "   " * TableStrings.flushright( 6, string(i); na=2)
            sb = " ";         for conf  in blockList[i].confs   sb = sb * string(conf) * ", "    end
            en = Float64[];   for level in  block.multiplet.levels    push!(en, level.energy)    end
            minEn = minimum(en);   minEn = Defaults.convertUnits("energy: from atomic", minEn)
            maxEn = maximum(en);   maxEn = Defaults.convertUnits("energy: from atomic", maxEn)
            sa = sa * TableStrings.flushleft(87, sb[1:end-2]; na=2)
            sb = "          " * string( length(block.multiplet.levels) )
            sa = sa * sb[end-9:end] * "        "
            sa = sa * TableStrings.flushleft(30, string( round(minEn)) * " ... " * string( round(maxEn)); na=2)
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(nx))

        return( nothing )
    end

    """
    `Cascade.displayLevels(stream::IO, multiplets::Array{Multiplet,1}; sa::String="")`
        ... display on stream the initial configurations as well as the calculated levels for all initial multiplets.
    """
    function displayLevels(stream::IO, multiplets::Array{Multiplet,1}; sa::String="")
        nx = 44
        println(stream, " ")
        println(stream, "* Configurations and levels for all given " * sa * "multiplets of the cascade, relative to the lowest:")
        for  multiplet  in multiplets
            println(stream, "  ")
            confList = Basics.extractNonrelativisticConfigurations(multiplet.levels[1].basis)
            for  conf in confList
                println(stream, "  $conf")
            end
            println(stream, "  ", TableStrings.hLine(nx))
            println(stream, "    Level  J Parity          Energy " * TableStrings.inUnits("energy") )
            println(stream, "  ", TableStrings.hLine(nx))
            for  i = 1:length(multiplet.levels)
                lev = multiplet.levels[i]
                en  = lev.energy - multiplet.levels[1].energy;    en_requested = Defaults.convertUnits("energy: from atomic", en)
                sc  = "   "  * TableStrings.level(i) * "     " * string(LevelSymmetry(lev.J, lev.parity)) * "     "
                @printf(stream, "%s %.15e %s", sc, en_requested, "\n")
            end
            println(stream, "  ", TableStrings.hLine(nx))
        end
        return( nothing )
    end

    """
    `Cascade.displaySteps(stream::IO, steps::Array{Cascade.Step,1}; sa::String="")`
        ... displays all predefined steps in a neat table and supports to delete individual steps from the list.
            sa::String can be used to display details about the given steps
    """
    function displaySteps(stream::IO, steps::Array{Cascade.Step,1}; sa::String="")
       nx = 170
        println(stream, " ")
        println(stream, "* Steps that are defined for the current " * sa * "cascade due to the given approach:")
        println(stream, " ")
        println(stream, "  ", TableStrings.hLine(nx))
        sa = "  "
        sa = sa * TableStrings.center( 9, "Step-No"; na=2)
        sa = sa * TableStrings.flushleft(11, "Process"; na=1)
        sa = sa * TableStrings.flushleft(55, "Initial:  No CSF, configuration(s)"; na=4)
        sa = sa * TableStrings.flushleft(55, "Final:  No CSF, configuration(s)"; na=4)
        sa = sa * TableStrings.flushleft(40, "Energies from ... to in " * TableStrings.inUnits("energy"); na=4)
        println(stream, sa)
        println(stream, "  ", TableStrings.hLine(nx))
        #
        for  i = 1:length(steps)
            sa = " " * TableStrings.flushright( 7, string(i); na=5)
            sa = sa  * TableStrings.flushleft( 11, string(steps[i].process); na=1)
            sb = "";   for conf in steps[i].initialConfigs   sb = sb * string(conf) * ", "    end
            sa = sa  * TableStrings.flushright( 5, string( length(steps[i].initialMultiplet.levels[1].basis.csfs) )*", "; na=0)
            sa = sa  * TableStrings.flushleft( 50, sb[1:end-2]; na=4)
            sb = "";   for conf in steps[i].finalConfigs     sb = sb * string(conf) * ", "    end
            sa = sa  * TableStrings.flushright( 5, string( length(steps[i].finalMultiplet.levels[1].basis.csfs) )*", "; na=0)
            sa = sa  * TableStrings.flushleft( 50, sb[1:end-2]; na=4)
            minEn = 1000.;   maxEn = -1000.;
            for  p = 1:length(steps[i].initialMultiplet.levels),  q = 1:length(steps[i].finalMultiplet.levels)
                minEn = min(minEn, steps[i].initialMultiplet.levels[p].energy - steps[i].finalMultiplet.levels[q].energy)
                maxEn = max(maxEn, steps[i].initialMultiplet.levels[p].energy - steps[i].finalMultiplet.levels[q].energy)
            end
            minEn = Defaults.convertUnits("energy: from atomic", minEn);   maxEn = Defaults.convertUnits("energy: from atomic", maxEn)
            sa = sa * string( round(minEn)) * " ... " * string( round(maxEn))
            println(stream, sa)
        end
        println(stream, "  ", TableStrings.hLine(nx))
    end

    """
    `Cascade.generateBlocks(comp::Cascade.Computation, confs::Array{Configuration,1}, initalOrbitals::Dict{Subshell, Orbital};
                            sa::String="", printout::Bool=true)`
        ... generate all block::Cascade.Block's, that need to be computed for this cascade, and compute also the corresponding multiplets.
            The different cascade approches enables one to realized follow different strategies how these block are selected and computed.
            A blockList::Array{Cascade.Block,1} is returned.
    """
    function generateBlocks(comp::Cascade.Computation, confs::Array{Configuration,1}, initalOrbitals::Dict{Subshell, Orbital};
                            sa::String="", printout::Bool=true)
        blockList = Cascade.Block[]
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        if    comp.approach == AverageSCA()
            if  printout
            println("\n* Generate blocks " * sa)
            println("\n  In the cascade approach $(comp.approach), the following assumptions/simplifications are made: ")
            println("    + orbitals from the initial multiplet are applied throughout; ")
            println("    + all blocks (multiplets) are generated from single-CSF levels and without any configuration mixing even in the SC; ")
            println("    + only E1 dipole transitions are applied in all radiative decay stets; ")
            println("    + for each decay step, a (single) set of continuum orbitals with an configuration averaged energy is applied " *
                           "for all transitions of the same step. \n")
            if  printSummary
            println(iostream, "\n* Generate blocks " * sa)
            println(iostream, "\n* In the cascade approach $(comp.approach), the following assumptions/simplifications are made: ")
            println(iostream, "    + orbitals from the initial multiplet are applied throughout; ")
            println(iostream, "    + all blocks (multiplets) are generated from single-CSF levels and without any configuration mixing even in the SC; ")
            println(iostream, "    + only E1 dipole transitions are applied in all radiative decay stets; ")
            println(iostream, "    + for each decay step, a (single) set of continuum orbitals with an configuration averaged energy is applied " *
                              "for all transitions of the same step. \n")
            end
            end     # printout
            #
            for  confa  in confs
                print("  Multiplet computations for $(string(confa)[1:end]) with $(confa.NoElectrons) electrons ... ")
                if  printSummary   println(iostream, "\n*  Multiplet computations for $(string(confa)[1:end]) with $(confa.NoElectrons) electrons ... ")   end
                multiplet = Basics.perform("computation: mutiplet from orbitals, no CI, CSF diagonal", [confa],  initalOrbitals,
                                           comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                # Shift the total energies of all levels if requested for the StepwiseDecayScheme
                if  typeof(comp.scheme) == StepwiseDecayScheme   &&   haskey(comp.scheme.chargeStateShifts, confa.NoElectrons)
                    energyShift = comp.scheme.chargeStateShifts[confa.NoElectrons]
                    multiplet   = Basics.shiftTotalEnergies(multiplet, energyShift)
                    print("shift all levels by $energyShift [a.u.] ... ")
                end
                push!( blockList, Cascade.Block(confa.NoElectrons, [confa], true, multiplet) )
                println("and $(length(multiplet.levels[1].basis.csfs)) CSF done. ")
            end
        elseif    comp.approach == SCA()
            if printout
            println("\n* Generate blocks " * sa)
            println("\n  In the cascade approach $(comp.approach), the following assumptions/simplifications are made: ")
            println("    + orbitals are generated independently for each multiplet (block); ")
            println("    + configuration interaction is included for each block; ")
            println("    + only E1 dipole transitions are applied in all radiative decay stets; \n")
            if  printSummary
            println(iostream, "\n* Generate blocks " * sa)
            println(iostream, "\n* In the cascade approach $(comp.approach), the following assumptions/simplifications are made: ")
            println(iostream, "    + orbitals are generated independently for each multiplet (block); ")
            println(iostream, "    + configuration interaction is included for each block; ")
            println(iostream, "    + only E1 dipole transitions are applied in all radiative decay stets; \n")
            end
            end     # printout
            #
            i = 0
            for  confa  in confs
                ## i = i + 1;    if   i in [1,2, 4,5,6,7,8,9,10,11,12,13,14]  ||  i > 15   println("  Block $i omitted.");    continue    end
                ## i = i + 1;    if   i < 11  ||  i > 11   println("  Block $i omitted.");    continue    end
                print("  Multiplet computations for $(string(confa)[1:end]) with $(confa.NoElectrons) electrons ... ")
                if  printSummary   print(iostream, "* Multiplet computations for $(string(confa)[1:end]) with $(confa.NoElectrons) electrons ... ")   end
                basis     = Basics.performSCF([confa], comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                multiplet = Basics.performCI(basis, comp.nuclearModel, comp.grid, comp.asfSettings; printout=false)
                push!( blockList, Cascade.Block(confa.NoElectrons, [confa], true, multiplet) )
                println("and $(length(multiplet.levels[1].basis.csfs)) CSF done. ")
            end
        else  error("Unsupported cascade approach.")
        end

        return( blockList )
    end

    """
    `Cascade.generateConfigurationList(multiplets::Array{Multiplet,1}, further::Int64, NoShake::Int64)`
        ... generates all possible (decay) configurations with up to further holes and with NoShake displacements with regard
            to the given multiplets. First, all configuratons are generated for which the hole is either moved 'outwards' or
            is moved and a second 'outer' hole is created; this step is repated further + 2 times to make sure that all relevant
            configurations are met. From the generated list, however, only those configurations are kept eventually with
            up to further holes, when compared to the configurations of the given multiplets. A confList::Array{Configuration,1}
            is returned.
    """
    function generateConfigurationList(multiplets::Array{Multiplet,1}, further::Int64, NoShake::Int64)
        # Determine all (different) configurations from multiplets
        confList = Configuration[]
        for mp  in  multiplets
            cfList = Basics.extractNonrelativisticConfigurations(mp.levels[1].basis)
            for  cf in cfList   if  cf in confList   nothing   else   push!(confList, cf)      end      end
        end
        cList = copy(confList);   initialNoElectrons = multiplets[1].levels[1].basis.NoElectrons
        # First, move and generate new 'outer' hole without displacements
        for  fur = 1:further+1
            newConfList = Configuration[]
            for conf  in cList
                holeList = Basics.determineHoleShells(conf)
                for  holeShell in holeList
                    wa = generateConfigurationsWith1OuterHole(conf,  holeShell);   append!(newConfList, wa)
                    wa = generateConfigurationsWith2OuterHoles(conf, holeShell);   append!(newConfList, wa)
                end
            end
            newConfList = unique(newConfList)
            ##x if  length(newConfList) > 0    newConfList = Basics.excludeDoubles(newConfList)    end
            cList = newConfList
            append!(confList, newConfList)
        end
        # Make sure that only configurations with up to further holes are returned
        newConfList = Configuration[]
        for   conf in confList
            if  conf.NoElectrons + further >= initialNoElectrons   push!(newConfList, conf)    end
        end
        # Add further shake-displacements if appropriate
        ##x newConfList = Basics.excludeDoubles(newConfList)
        newConfList = unique(newConfList)
        return( newConfList )
    end

    """
    `Cascade.generateConfigurationsWith1OuterHole(conf,  holeShell)`
        ... generates all possible (decay) configurations where the hole in holeShell is moved 'outwards'.
            A confList::Array{Configuration,1} is returned.
    """
    function generateConfigurationsWith1OuterHole(conf::Configuration,  holeShell::Shell)
         shList = Basics.generate("shells: ordered list for NR configurations", [conf]);   i0 = 0
         for  i = 1:length(shList)
             if   holeShell == shList[i]    i0 = i;    break    end
         end
         if  i0 == 0   error("stop a")   end
         #
         # Now move the hole 'outwards'
         confList = Configuration[]
         for  i = i0+1:length(shList)
             if  haskey(conf.shells, shList[i])  &&  conf.shells[ shList[i] ] >= 1
                 newshells = copy( conf.shells )
                 newshells[ shList[i] ] = newshells[ shList[i] ] - 1
                 newshells[ holeShell ] = newshells[ holeShell ] + 1
                 push!(confList, Configuration( newshells, conf.NoElectrons ) )
             end
         end
         return( confList )
    end

    """
    `Cascade.generateConfigurationsWith2OuterHoles(conf,  holeShell)`
        ... generates all possible (decay) configurations where the hole in holeShell is moved 'outwards'.
            A confList::Array{Configuration,1} is returned.
    """
    function generateConfigurationsWith2OuterHoles(conf::Configuration,  holeShell::Shell)
         shList = Basics.generate("shells: ordered list for NR configurations", [conf]);   i0 = 0
         for  i = 1:length(shList)
             if   holeShell == shList[i]    i0 = i;    break    end
         end
         if  i0 == 0   error("stop a")   end
         #
         # Now move the hole 'outwards'
         confList = Configuration[]
         for  i = i0+1:length(shList)
             if  haskey(conf.shells, shList[i])  &&  conf.shells[ shList[i] ] >= 2
                 newshells = copy( conf.shells )
                 newshells[ shList[i] ] = newshells[ shList[i] ] - 2
                 newshells[ holeShell ] = newshells[ holeShell ] + 1
                 push!(confList, Configuration( newshells, conf.NoElectrons - 1 ) )
             end
             #
             for  j = i0+1:length(shList)
                 if  i != j   &&   haskey(conf.shells, shList[i])  &&  conf.shells[ shList[i] ] >= 1   &&
                                   haskey(conf.shells, shList[j])  &&  conf.shells[ shList[j] ] >= 1
                     newshells = copy( conf.shells )
                     newshells[ shList[i] ] = newshells[ shList[i] ] - 1
                     newshells[ shList[j] ] = newshells[ shList[j] ] - 1
                     newshells[ holeShell ] = newshells[ holeShell ] + 1
                     push!(confList, Configuration( newshells, conf.NoElectrons - 1 ) )
                 end
             end
         end
         return( confList )
    end

    """
    `Cascade.groupDisplayConfigurationList(Z::Float64, confs::Array{Configuration,1}; sa::String="")`
        ... group & display the configuration list into sublists with the same No. of electrons; this lists are displayed together
            with an estimated total energy. An ordered confList::Array{Configuration,1} is returned with configurations of decreasing
            number of electrons.
    """
    function groupDisplayConfigurationList(Z::Float64, confs::Array{Configuration,1}; sa::String="")
        minNoElectrons = 1000;   maxNoElectrons = 0
        for  conf in confs
            minNoElectrons = min(minNoElectrons, conf.NoElectrons)
            maxNoElectrons = max(maxNoElectrons, conf.NoElectrons)
        end
        #
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        #
        println("\n* Electron configuration(s) used:    " * sa)
        ## @warn "*** Limit to just 4 configurations for each No. of electrons. ***"                       ## delete nxx
        if  printSummary   println(iostream, "\n* Electron configuration(s) used:    " * sa)    end
        confList = Configuration[];   nc = 0
        for  n = maxNoElectrons:-1:minNoElectrons
            nxx = 0                                                                                ## delete nxx
            println("\n  Configuration(s) with $n electrons:")
            if  printSummary   println(iostream, "\n    Configuration(s) with $n electrons:")      end
            nd = 0
            for  conf in confs  nd = max(nd, length("      " * string(conf)))   end
            for  conf in confs
                if n == conf.NoElectrons
                    ## nxx = nxx + 1;    if nxx > 4   break    end                                         ## delete nxx
                    nc = nc + 1
                    push!(confList, conf )
                    if  Z > 36.0    wa = 0.
                    else            wa = Semiempirical.estimate("binding energy", round(Int64, Z), conf);
                                    wa = Defaults.convertUnits("energy: from atomic", wa)
                    end
                    sb = "   av. BE = "  * string( round(-wa) ) * "  " * TableStrings.inUnits("energy")
                    sd = "      " * string(conf) * "                                "
                    println(sd[1:nd+3] * sb * "      ($nc)" )
                    if  printSummary   println(iostream, sd[1:nd+3] * sb * "      ($nc)")      end
                end
            end
        end

        println("\n  A total of $nc configuration have been defined for this " * sa * "cascade, and selected configurations could be " *
                "removed here:  [currently not supported]")
        if  printSummary   println(iostream, "\n* A total of $nc configuration have been defined for this cascade, and selected " *
                                             "configurations could be removed here:  [currently not supported]")      end
        return( confList )
    end

    """
    `Cascade.modifySteps(stepList::Array{Cascade.Step,1})`
        ... allows the user to modify the steps, for instance, by deleting selected steps of the cascade or by modifying the settings of
            one or several steps. A newStepList::Array{Cascade.Step,1} for which the transition data are eventually computed.
    """
    function modifySteps(stepList::Array{Cascade.Step,1})
        #
        newStepList = Cascade.Step[]
        #
        println("\n* Here, modify the individual steps explicitly in the code, if needed, ...... and just do it !!")
        #
        #  Delete individual steps from stepList
        #  if  i in [1,2,5, ...] modify the particular settings, etc.
        for  i = 1:length(stepList)
            step = stepList[i]
            #
            if  i in []
                println("  Modify step $i :")
                newStep = Cascade.Step(step.process, step.settings, step.initialConfs, step.finalConfs, step.initialMultiplet, step.initialMultiplet)
                push!(newStepList, newStep)
            else
                push!(newStepList, step)
            end
        end
        #
        # wa = [1,2,3]
        # delete from list
        #
        println("\n  A total of $(length(newStepList)) steps are still defined in the cascade.")
        printSummary, iostream = Defaults.getDefaults("summary flag/stream")
        if  printSummary   println(iostream, "\n* A total of $(length(newStepList)) steps are still defined in the cascade.")    end

        return( newStepList )
    end


text:
import MechanicalSketch: @import_expand, empty_figure, WI, HE, EM, O, finish
import MechanicalSketch: readsvg, place_image, settext, Quantity, ustrip, unit
import MechanicalSketch: Point, julialogo, PALETTE, color_with_lumin, @svg
import MechanicalSketch: background, color_with_alpha, scale, @layer, setopacity
import MechanicalSketch: rotate, translate
using MechGluecode
import Plots.plot
# Plot dependencies in 'test_functions_37.jl'

if !@isdefined N
    @import_expand(~m, s, °, N)
end
include("test_functions_37.jl")
empty_figure(filename = joinpath(@__DIR__, "test_37.png");
      backgroundcolor = color_with_lumin(PALETTE[5], 90),
      hue = "black");
@layer begin
    scale(4)
    rotate(30°)
    translate(-WI /10, -HE / 8)
    background(color_with_alpha(PALETTE[5], 0))
    julialogo(;bodycolor = color_with_alpha(PALETTE[1], 0.1) )
end;
pt = O + (-WI / 2 + EM, -HE / 2 + 2EM)
settext("<b>Plot recipe with svg style modifications</b>", pt, markup = true)
pt += (0, EM)

x = 0.0s:1.0s:5.0s
y = x .* N ./ s

settext("MechGlueCode has quantity plot recipes:", pt)
pl1 = plot(x, y);
ptul, ptlr = place_image_37(pt, pl1, height = 0.25HE)
pt += (0, 2EM + ptlr[2] - ptul[2])

settext("When guides are specified, postfix units:", pt)
pl2 = plot(x, y, xguide = "Time", yguide = "Force");
ptul, ptlr = place_image_37(pt, pl2, height = 0.25HE)
pt += (0, 2EM + ptlr[2] - ptul[2])

settext("This recipe treats axes independently:", pt)
pl2 = plot(x, y, xguide = "Time");
ptul, ptlr = place_image_37(pt, pl2, height = 0.25HE)
pt += (0, 2EM + ptlr[2] - ptul[2])

pt = O + (0, -HE / 4 + EM)
settext("No recipe change required for 3d:", pt)
n = 100
ts = range(0, stop = 8π, length = n)
x = ts * N .* map(cos, ts)
y = (0.1ts) * N .* map(sin, ts) .|> N
z = (1:n)s
pl3 = plot(x, y, z, zcolor = reverse(z./s),
      m = (10, 0.8, :blues, Plots.stroke(0)),
      leg = false, cbar = true, w = 5,
      zguide = "Time");
plot!(pl3, zeros(n), zeros(n), 1:n, w = 10);
ptul, ptlr = place_image_37(pt, pl3; width = 0.45WI)

finish()


text:
using JudyDicts
using Base.Test

function test_juHS()
    #println("testing JUHS...")
    local ja = JudyDict{Array{Uint8}, Int}()

    local indices = [(b"0001", 11), (b"00021", 121), (b"0005", 15)]

    for i in indices
        ret_ptr = ju_set(ja, i[1], i[2])
        @test C_NULL != ret_ptr
        @test i[2] == unsafe_load(ret_ptr)
    end

    for i in indices
        ret_tuple = ju_get(ja, i[1])
        @test i[2] == ret_tuple[1]
        @test C_NULL != ret_tuple[2]
    end

    @test (indices[1])[2] == ja[(indices[1])[1]]
    ja[(indices[1])[1]] = 2 * (indices[1])[2]
    @test (2 * ((indices[1])[2])) == ja[(indices[1])[1]]

    @test 1 == ju_unset(ja, (indices[1])[1])
    ja[(indices[1])[1]] = 2 * (indices[1])[2]
    @test (2 * ((indices[1])[2])) == ja[(indices[1])[1]]
    ja[(indices[1])[1]] = nothing
    @test nothing == ja[(indices[1])[1]]

    ja[(indices[1])[1]] = (indices[1])[2]
    @test (indices[1])[2] == ja[(indices[1])[1]]
    delete!(ja, (indices[1])[1])
    @test nothing == ja[(indices[1])[1]]

    @test C_NULL == (ju_get(ja, (indices[1])[1]))[2]
end

function test_juSO()
    #println("testing JUSL...")
    local ja = JudyDict{String, ASCIIString}()

    local indices = [("0001", "11"), ("00021", "121"), ("0005", "15")]

    for i in indices
        ja[i[1]] = i[2]
    end

    for i in indices
        @test i[2] == ja[i[1]]
    end

    ret_tuple = (nothing, C_NULL, 0)
    for i in 1:length(indices)
        ret_tuple = (i == 1) ? ju_first(ja) : ju_next(ja)
        @test ((indices[i])[2] == ret_tuple[1]) && (C_NULL != ret_tuple[2]) && ((indices[i])[1] == ret_tuple[3])
    end
    @test C_NULL == (ju_next(ja))[2]
    for i in 1:length(indices)
        ret_tuple = (i == 1) ? ju_last(ja) : ju_prev(ja)
        local rev_idx = length(indices)-i+1
        @test ((indices[rev_idx])[2] == ret_tuple[1]) && (C_NULL != ret_tuple[2]) && ((indices[rev_idx])[1] == ret_tuple[3])
    end
    @test C_NULL == (ju_prev(ja))[2]

    i = 1
    for x in ja
        @test ((indices[i])[2] == x[1]) && ((indices[i])[1] == x[2])
        i+=1
    end

    @test (indices[1])[2] == ja[(indices[1])[1]]
    ja[(indices[1])[1]] = "2" * (indices[1])[2]
    @test ("2" * ((indices[1])[2])) == ja[(indices[1])[1]]

    @test 1 == ju_unset(ja, (indices[1])[1])
    @test C_NULL == (ju_get(ja, (indices[1])[1]))[2]

    ja[(indices[1])[1]] = (indices[1])[2]
    @test (indices[1])[2] == ja[(indices[1])[1]]
    delete!(ja, (indices[1])[1])
    @test nothing == ja[(indices[1])[1]]
end

function test_juSL()
    #println("testing JUSL...")
    local ja = JudyDict{String, Int}()

    local indices = [("0001", 11), ("00021", 121), ("0005", 15)]

    for i in indices
        ret_ptr = ju_set(ja, i[1], i[2])
        @test C_NULL != ret_ptr
        @test i[2] == unsafe_load(ret_ptr)
    end

    for i in indices
        ret_tuple = ju_get(ja, i[1])
        @test i[2] == ret_tuple[1]
        @test C_NULL != ret_tuple[2]
    end

    ret_tuple = (0, C_NULL, 0)
    for i in 1:length(indices)
        ret_tuple = (i == 1) ? ju_first(ja) : ju_next(ja)
        @test ((indices[i])[2] == ret_tuple[1]) && (C_NULL != ret_tuple[2]) && ((indices[i])[1] == ret_tuple[3])
    end
    @test C_NULL == (ju_next(ja))[2]
    for i in 1:length(indices)
        ret_tuple = (i == 1) ? ju_last(ja) : ju_prev(ja)
        local rev_idx = length(indices)-i+1
        @test ((indices[rev_idx])[2] == ret_tuple[1]) && (C_NULL != ret_tuple[2]) && ((indices[rev_idx])[1] == ret_tuple[3])
    end
    @test C_NULL == (ju_prev(ja))[2]

    i = 1
    for x in ja
        @test ((indices[i])[2] == x[1]) && ((indices[i])[1] == x[2])
        i+=1
    end

    @test (indices[1])[2] == ja[(indices[1])[1]]
    ja[(indices[1])[1]] = 2 * (indices[1])[2]
    @test (2 * ((indices[1])[2])) == ja[(indices[1])[1]]

    @test 1 == ju_unset(ja, (indices[1])[1])
    @test C_NULL == (ju_get(ja, (indices[1])[1]))[2]

    ja[(indices[1])[1]] = (indices[1])[2]
    @test (indices[1])[2] == ja[(indices[1])[1]]
    delete!(ja, (indices[1])[1])
    @test nothing == ja[(indices[1])[1]]
end

function test_juLO()
    #println("testing JUL...")
    local ja = JudyDict{Int, ASCIIString}()

    @test 0 == ju_mem_used(ja)

    local indices = [(1, "11"), (5, "15"), (21, "121"), (41, "141"), (51, "151"), (61, "161")]

    for i in indices
        ja[i[1]] = i[2]
    end

    for i in indices
        @test i[2] == ja[i[1]]
    end

    @test length(indices) == ju_count(ja)

    for i in 1:length(indices)
        ret_tuple = ju_by_count(ja, i)
        @test ((indices[i])[2] == ret_tuple[1]) && (C_NULL != ret_tuple[2]) && ((indices[i])[1] == ret_tuple[3])
    end

    @test ju_mem_used(ja) > 0

    ret_tuple = (0, C_NULL, 0)
    for i in 1:length(indices)
        ret_tuple = (i == 1) ? ju_first(ja) : ju_next(ja)
        @test ((indices[i])[2] == ret_tuple[1]) && (C_NULL != ret_tuple[2]) && ((indices[i])[1] == ret_tuple[3])
    end
    @test C_NULL == (ju_next(ja))[2]
    for i in 1:length(indices)
        ret_tuple = (i == 1) ? ju_last(ja) : ju_prev(ja)
        local rev_idx = length(indices)-i+1
        @test ((indices[rev_idx])[2] == ret_tuple[1]) && (C_NULL != ret_tuple[2]) && ((indices[rev_idx])[1] == ret_tuple[3])
    end
    @test C_NULL == (ju_prev(ja))[2]

    i = 1
    for x in ja
        @test ((indices[i])[2] == x[1]) && ((indices[i])[1] == x[2])
        i+=1
    end

    for i in 1:20
        ret_tuple = (i == 1) ? ju_first_empty(ja) : ju_next_empty(ja)
        @test 1 == ret_tuple[1]
    end
    for i in 1:20
        ret = (i == 1) ? ju_last_empty(ja) : ju_prev_empty(ja)
        @test 1 == ret_tuple[1]
    end

    @test (indices[1])[2] == ja[(indices[1])[1]]
    ja[(indices[1])[1]] = "2" * (indices[1])[2]
    @test "2" * (indices[1])[2] == ja[(indices[1])[1]]

    @test 1 == ju_unset(ja, (indices[1])[1])
    @test C_NULL == (ju_get(ja, (indices[1])[1]))[2]

    ja[(indices[1])[1]] = (indices[1])[2]
    @test (indices[1])[2] == ja[(indices[1])[1]]
    delete!(ja, (indices[1])[1])
    @test nothing == ja[(indices[1])[1]]
end

function test_juL()
    #println("testing JUL...")
    local ja = JudyDict{Int, Int}()

    @test 0 == ju_mem_used(ja)

    local indices = [(1, 11), (5, 15), (21, 121), (41, 141), (51, 151), (61, 161)]

    for i in indices
        ret_ptr = ju_set(ja, i[1], i[2])
        @test C_NULL != ret_ptr
        @test i[2] == unsafe_load(ret_ptr)
    end

    for i in indices
        ret_tuple = ju_get(ja, i[1])
        @test i[2] == ret_tuple[1]
        @test C_NULL != ret_tuple[2]
    end

    @test length(indices) == ju_count(ja)

    for i in 1:length(indices)
        ret_tuple = ju_by_count(ja, i)
        @test ((indices[i])[2] == ret_tuple[1]) && (C_NULL != ret_tuple[2]) && ((indices[i])[1] == ret_tuple[3])
    end

    @test ju_mem_used(ja) > 0

    ret_tuple = (0, C_NULL, 0)
    for i in 1:length(indices)
        ret_tuple = (i == 1) ? ju_first(ja) : ju_next(ja)
        @test ((indices[i])[2] == ret_tuple[1]) && (C_NULL != ret_tuple[2]) && ((indices[i])[1] == ret_tuple[3])
    end
    @test C_NULL == (ju_next(ja))[2]
    for i in 1:length(indices)
        ret_tuple = (i == 1) ? ju_last(ja) : ju_prev(ja)
        local rev_idx = length(indices)-i+1
        @test ((indices[rev_idx])[2] == ret_tuple[1]) && (C_NULL != ret_tuple[2]) && ((indices[rev_idx])[1] == ret_tuple[3])
    end
    @test C_NULL == (ju_prev(ja))[2]

    i = 1
    for x in ja
        @test ((indices[i])[2] == x[1]) && ((indices[i])[1] == x[2])
        i+=1
    end

    for i in 1:20
        ret_tuple = (i == 1) ? ju_first_empty(ja) : ju_next_empty(ja)
        @test 1 == ret_tuple[1]
    end
    for i in 1:20
        ret = (i == 1) ? ju_last_empty(ja) : ju_prev_empty(ja)
        @test 1 == ret_tuple[1]
    end

    @test (indices[1])[2] == ja[(indices[1])[1]]
    ja[(indices[1])[1]] = 2 * (indices[1])[2]
    @test (2 * ((indices[1])[2])) == ja[(indices[1])[1]]

    @test 1 == ju_unset(ja, (indices[1])[1])
    @test C_NULL == (ju_get(ja, (indices[1])[1]))[2]

    ja[(indices[1])[1]] = (indices[1])[2]
    @test (indices[1])[2] == ja[(indices[1])[1]]
    delete!(ja, (indices[1])[1])
    @test nothing == ja[(indices[1])[1]]
end

function test_ju1()
    #println("testing JU1...")
    local ja = JudyDict{Int, Bool}()

    @test 0 == ju_mem_used(ja)

    local indices = [1, 5]
    for i in indices
        @test 1 == ju_set(ja, i)
    end
    for i in indices
        @test 1 == ju_get(ja, i)
    end
    @test length(indices) == ju_count(ja)

    for i in 1:length(indices)
        ret_tuple = ju_by_count(ja, i)
        @test (1 == ret_tuple[1]) && (indices[i] == ret_tuple[2])
    end

    @test ju_mem_used(ja) > 0

    local ret = (0, 0)
    for i in 1:length(indices)
        ret = (i == 1) ? ju_first(ja) : ju_next(ja)
        @test (1 == ret[1]) && (indices[i] == ret[2])
    end
    @test 0 == (ju_next(ja))[1]
    for i in 1:length(indices)
        ret = (i == 1) ? ju_last(ja) : ju_prev(ja)
        @test (1 == ret[1]) && (indices[length(indices)-i+1] == ret[2])
    end
    @test 0 == (ju_prev(ja))[1]

    i = 1
    for x in ja
        @test (indices[i] == x[2])
        i+=1
    end

    for i in 1:20
        ret = (i == 1) ? ju_first_empty(ja) : ju_next_empty(ja)
        @test 1 == ret[1]
    end
    for i in 1:20
        ret = (i == 1) ? ju_last_empty(ja) : ju_prev_empty(ja)
        @test 1 == ret[1]
    end

    @test true == ja[indices[1]]
    ja[indices[1]] = false
    @test false == ja[indices[1]]
    ja[indices[1]] = true

    @test 1 == ju_unset(ja, indices[1])
    @test 1 != ju_get(ja, indices[1])

    ja[indices[1]] = true
    @test true == ja[indices[1]]
    delete!(ja, indices[1])
    @test false == ja[indices[1]]
end

test_juHS()
test_juSL()
test_juL()
test_ju1()
test_juLO()
test_juSO()


text:
#_______________________________________________________________________________
# Writting results

if SaveInputs
    CPworkbookResults = CPworkbook
else
    CPworkbookResults = Workbook()
end

    # requirements fulfillement
RequirSheet = getSheet( CPworkbookResults, "Requirements" )
if RequirSheet.ptr === Ptr{Void}( 0 )
    RequirSheet = createSheet( CPworkbookResults, "Requirements" )
end

rowIdx = 0
Row = createRow(RequirSheet, rowIdx)
rowIdx += 1

cellIdx = 1
for i in 1:length(MPObjectives)
    obj = ""
    for j in 1:length(MPObjectives[i].Objectives)
        obj = obj * MPObjectives[i].Objectives[j]
    end
    Cell=createCell(Row, cellIdx); setCellValue(Cell, obj)
    cellIdx += 1
end
for i in 1:NBYears
    Row = createRow(RequirSheet, rowIdx)
    rowIdx += 1
    cellIdx = 0
    Cell=createCell(Row, cellIdx); setCellValue(Cell, i)
    cellIdx += 1
    for j in 1:length(MPObjectives)
        Cell=createCell(Row, cellIdx);
        if IntegerSolution
            setCellValue(Cell, round(Int, RequirementFulfillment[i, j]))
        else
            setCellValue(Cell, RequirementFulfillment[i, j])
        end
        cellIdx += 1
    end
end

    # Initial manpower distribution
InitDistSheet = getSheet( CPworkbookResults, "Initial distribution" )
if InitDistSheet.ptr === Ptr{Void}( 0 )
    InitDistSheet = createSheet( CPworkbookResults, "Initial distribution" )
end
maxCols = 0
rowIdx = 0
Row = createRow(InitDistSheet, rowIdx)
rowIdx += 1

cellIdx = 2
#for i in 1:length(GuyCareerPaths)
#    Cell=createCell(Row, cellIdx); setCellValue(Cell, i)
#    cellIdx += 1
#end

index = 1
for i in 1:length(InitManpower)
    Row = createRow(InitDistSheet, rowIdx)
    rowIdx += 1
    cellIdx = 0
    Cell=createCell(Row, cellIdx); setCellValue(Cell, InitManpower[i].Academiclvl)
    cellIdx += 1
    Cell=createCell(Row, cellIdx); setCellValue(Cell, InitManpower[i].Seniority)
    cellIdx += 1
    for j in 1:length(GuyCareerPaths)
        if j in InitMPPartCP[i]
            if sol.sol[index] != 0
                Cell=createCell(Row, cellIdx); setCellValue(Cell, j); cellIdx += 1
                Cell=createCell(Row, cellIdx); setCellValue(Cell, sol.sol[index]); cellIdx += 1
            end
            index += 1
        end
#        Cell=createCell(Row, cellIdx);

#        if j in InitMPPartCP[i]
#            setCellValue(Cell, sol.sol[index])
#            index += 1
#        else
#            setCellValue(Cell, 0)
    end
    maxCols = max(maxCols, cellIdx)
end
Row = createRow(InitDistSheet, 0)
setCellValue(createCell(Row, 0), "SubPop")
setCellValue(createCell(Row, 1), "Seniority")
index = 2
while index < maxCols
    setCellValue(createCell(Row, index), "CareerPath"); index += 1;
    setCellValue(createCell(Row, index), "Nb assigned"); index += 1;
end

    # Annual Recruitment
RecruitmentSheet = getSheet( CPworkbookResults, "Perforemed Recruitment" )
if RecruitmentSheet.ptr === Ptr{Void}( 0 )
    RecruitmentSheet = createSheet( CPworkbookResults, "Perforemed Recruitment" )
end
maxCols = 0
rowIdx = 0
rowIdx += 1

for i in 1:NBYears
    sum = 0
    Row = createRow(RecruitmentSheet, rowIdx)
    rowIdx += 1
    cellIdx = 0
    Cell=createCell(Row, cellIdx); setCellValue(Cell, i)
    cellIdx += 2
    for j in 1:length(GuyCareerPaths)
        if sol.sol[InitMPDivisionNb + (i - 1)*length(GuyCareerPaths) + j] != 0
            sum += sol.sol[InitMPDivisionNb + (i - 1)*length(GuyCareerPaths) + j]
            Cell=createCell(Row, cellIdx); setCellValue(Cell, j)
            cellIdx += 1
            Cell=createCell(Row, cellIdx); setCellValue(Cell, sol.sol[InitMPDivisionNb + (i - 1)*length(GuyCareerPaths) + j])
            cellIdx += 1
        end
    end
    Cell=createCell(Row, 1); setCellValue(Cell, sum)
    maxCols = max(maxCols, cellIdx - 1)
end

Row = createRow(RecruitmentSheet, 0)
setCellValue(createCell(Row, 0), "Year")
setCellValue(createCell(Row, 1), "Total recruitment")
index = 2
while index < maxCols
    setCellValue(createCell(Row, index), "CareerPath"); index += 1;
    setCellValue(createCell(Row, index), "Nb recruits"); index += 1;
end

    # Deviation Variables

DeviationVariables = getSheet( CPworkbookResults, "DeviationVariables" )
if DeviationVariables.ptr === Ptr{Void}( 0 )
    DeviationVariables = createSheet( CPworkbookResults, "DeviationVariables" )
end

rowIdx = 0
Row = createRow(DeviationVariables, rowIdx)
rowIdx += 1

cellIdx = 0
Cell=createCell(Row, cellIdx); setCellValue(Cell, rowIdx)
cellIdx += 1
for i in 1:(length(GuyCareerPaths) * NBYears)

    Cell=createCell(Row, cellIdx); setCellValue(Cell, sol.sol[InitMPDivisionNb + i])
    cellIdx += 1
    if i%length(GuyCareerPaths) == 0
        Row = createRow(DeviationVariables, rowIdx)
        rowIdx += 1

        cellIdx = 0
        Cell=createCell(Row, cellIdx); setCellValue(Cell, rowIdx)
        cellIdx += 1
    end
end
        #Positive deviation
rowIdx = 0
Row = createRow(DeviationVariables, rowIdx)
rowIdx += 1

cellIdx = 0
Cell=createCell(Row, cellIdx); setCellValue(Cell, rowIdx)
cellIdx += 1
for i in 1:(length(MPObjectives) * NBYears)

    Cell=createCell(Row, cellIdx); setCellValue(Cell, sol.sol[InitMPDivisionNb + AnnualRecDivNb + i])
    cellIdx += 1
    if i%length(MPObjectives) == 0
        Row = createRow(DeviationVariables, rowIdx)
        rowIdx += 1

        cellIdx = 0
        Cell=createCell(Row, cellIdx); setCellValue(Cell, rowIdx)
        cellIdx += 1
    end
end
        # Negative deviation
rowIdx += 2
Row = createRow(DeviationVariables, rowIdx)
rowIdx += 1

cellIdx = 0
Cell=createCell(Row, cellIdx); setCellValue(Cell, rowIdx)
cellIdx += 1
for i in 1:(length(MPObjectives) * NBYears)

    Cell=createCell(Row, cellIdx); setCellValue(Cell, sol.sol[InitMPDivisionNb + AnnualRecDivNb + (length(MPObjectives) * NBYears) + i])
    cellIdx += 1
    if i%length(MPObjectives) == 0
        Row = createRow(DeviationVariables, rowIdx)
        rowIdx += 1

        cellIdx = 0
        Cell=createCell(Row, cellIdx); setCellValue(Cell, rowIdx)
        cellIdx += 1
    end
end

try rm(joinpath( dirname( Base.source_path() ), "..", "..", "Results_$SimulationName.xlsx" )) end
write(joinpath( dirname( Base.source_path() ), "..", "..", "Results_$SimulationName.xlsx" ), CPworkbookResults)
#XlsxFile, CPworkbookResults)#

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
module TextModels

    using WordTokenizers
    using Languages

    using DataDeps
    using Pkg.Artifacts

    using Flux, Tracker
    using Flux: identity, onehot, onecold, @treelike, onehotbatch

    using TextAnalysis
    import TextAnalysis: fit!, predict

    export CRF, viterbi_decode, crf_loss

    export NERTagger, PoSTagger, SentimentAnalyzer

    include("averagePerceptronTagger.jl")
    include("sentiment.jl")

    # CRF
    include("CRF/crf.jl")
    include("CRF/predict.jl")
    include("CRF/crf_utils.jl")
    include("CRF/loss.jl")

    # NER and POS
    include("sequence/ner_datadeps.jl")
    include("sequence/ner.jl")
    include("sequence/pos_datadeps.jl")
    include("sequence/pos.jl")
    include("sequence/sequence_models.jl")

    # ULMFiT
    module ULMFiT
        using ..TextAnalysis
        using DataDeps
        using Flux
        using Tracker
        using BSON
        include("ULMFiT/utils.jl")
        include("ULMFiT/datadeps.jl")
        include("ULMFiT/data_loaders.jl")
        include("ULMFiT/custom_layers.jl")
        include("ULMFiT/pretrain_lm.jl")
        include("ULMFiT/fine_tune_lm.jl")
        include("ULMFiT/train_text_classifier.jl")
    end
    export ULMFiT

    function __init__()
        pos_tagger_datadep_register()
        ner_datadep_register()
        pos_datadep_register()
        ULMFiT.ulmfit_datadep_register()

        global sentiment_model = artifact"sentiment_model"
    end
end


text:
import OmegaCore
"""

```
using OmegaCore
using SoftPredicates
using ReplicaExchange
using InferenceBase

# Model
x = 1 ~ StdNormal{Float64}()
y = 2 ~ StdNormal{Float64}()
z = pw(==ₛ, x, y)
joint_post = @joint(x,y) |ᶜ z

# Target
loglikelihood = logerr ∘ softconstraints(joint_post)
temp_range = [0.001, 0.1, 1, 10, 100]
logposterior(ω) = logenergyexo(ω) + loglikelihood(ω)

struct TemperedLogPosterior{T}
  α::T
end
(at::TemperedLogPosterior)(ω) = withkernel(kseα(at.α)) do
  logposterior(ω)
end

logenergys = TemperedLogPosterior.(temp_range)
randsample(joint_post, 10; alg = Replica, logenergy = logposterior)

sim_chain_keep_n(rng, logenergy, n; state_init) =
  omegarandsample(rng, logenergy, n, MH; state_init = init_state)
```
"""
function OmegaCore.omegarandsample(rng,
                                   logenergy,
                                   n,
                                   ::ReplicaAlg;
                                   num_swaps = 10,
                                   state_init = auto_init(rng, x),
                                   sim_chain_keep_n,
                                   propose_and_logratio = defpropose_and_logratio(x),
                                   logenergys,
                                   kwargs...)
  samples_per_swap = div(n, num_swaps) # zt: might not fit right
  states_init = # do-something

  ωsamples = re(rng,
                logenergys,
                samples_per_swap,
                num_swaps,
                states_init,
                n,
                state_init,
                sim_chain_keep_n,
                kwargs...)
end

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
############################################################################################
### Types
############################################################################################
@doc raw"""
Subtype of `SemImply` that implements the RAM notation with symbolic precomputation.

# Constructor

    RAMSymbolic(;specification,
        vech = false,
        gradient = true,
        hessian = false,
        approximate_hessian = false,
        meanstructure = false,
        kwargs...)

# Arguments
- `specification`: either a `RAMMatrices` or `ParameterTable` object
- `meanstructure::Bool`: does the model have a meanstructure?
- `gradient::Bool`: is gradient-based optimization used
- `hessian::Bool`: is hessian-based optimization used
- `approximate_hessian::Bool`: for hessian based optimization: should the hessian be approximated
- `vech::Bool`: should the half-vectorization of Σ be computed (instead of the full matrix)
    (automatically set to true if any of the loss functions is SemWLS)

# Extended help

## Implementation
Subtype of `SemImply`.

## Interfaces
- `identifier(::RAMSymbolic) `-> Dict containing the parameter labels and their position
- `n_par(::RAMSymbolic)` -> Number of parameters

- `Σ(::RAMSymbolic)` -> model implied covariance matrix
- `μ(::RAMSymbolic)` -> model implied mean vector

Jacobians (only available in gradient! calls)
- `∇Σ(::RAMSymbolic)` -> ``∂vec(Σ)/∂θᵀ``
- `∇μ(::RAMSymbolic)` -> ``∂μ/∂θᵀ``

- `∇Σ_function(::RAMSymbolic)` -> function to overwrite `∇Σ` in place,
    i.e. `∇Σ_function(∇Σ, θ)`. Normally, you do not want to use this but simply
    query `∇Σ(::RAMSymbolic)`.

Hessians
The computation of hessians is more involved, and uses the "chain rule for
hessian matrices".
Therefore, we desribe it at length in the mathematical appendix of the online documentation,
and the relevant interfaces are omitted here.

Additional interfaces
- `has_meanstructure(::RAMSymbolic)` -> `Val{Bool}` does the model have a meanstructure?

## RAM notation
The model implied covariance matrix is computed as
```math
    \Sigma = F(I-A)^{-1}S(I-A)^{-T}F^T
```
and for models with a meanstructure, the model implied means are computed as
```math
    \mu = F(I-A)^{-1}M
```
"""
struct RAMSymbolic{F1, F2, F3, A1, A2, A3, S1, S2, S3, V, V2, F4, A4, F5, A5, D1, B} <: SemImplySymbolic
    Σ_function::F1
    ∇Σ_function::F2
    ∇²Σ_function::F3
    Σ::A1
    ∇Σ::A2
    ∇²Σ::A3
    Σ_symbolic::S1
    ∇Σ_symbolic::S2
    ∇²Σ_symbolic::S3
    n_par::V
    ram_matrices::V2
    μ_function::F4
    μ::A4
    ∇μ_function::F5
    ∇μ::A5
    identifier::D1
    has_meanstructure::B
end

############################################################################################
### Constructors
############################################################################################

function RAMSymbolic(;
        specification,
        loss_types = nothing,
        vech = false,
        gradient = true,
        hessian = false,
        meanstructure = false,
        approximate_hessian = false,
        kwargs...)

    ram_matrices = RAMMatrices(specification)
    identifier = StructuralEquationModels.identifier(ram_matrices)

    n_par = length(ram_matrices.parameters)
    n_var, n_nod = ram_matrices.size_F

    par = (Symbolics.@variables θ[1:n_par])[1]

    A = zeros(Num, n_nod, n_nod)
    S = zeros(Num, n_nod, n_nod)
    !isnothing(ram_matrices.M_ind) ? M = zeros(Num, n_nod) : M = nothing
    F = zeros(ram_matrices.size_F); F[CartesianIndex.(1:n_var, ram_matrices.F_ind)] .= 1.0

    set_RAMConstants!(A, S, M, ram_matrices.constants)
    fill_A_S_M(A, S, M, ram_matrices.A_ind, ram_matrices.S_ind, ram_matrices.M_ind, par)

    A, S, F = sparse(A), sparse(S), sparse(F)

    if !isnothing(loss_types)
        any(loss_types .<: SemWLS) ? vech = true : nothing
    end

    # Σ
    Σ_symbolic = get_Σ_symbolic_RAM(S, A, F; vech = vech)
    #print(Symbolics.build_function(Σ_symbolic)[2])
    Σ_function = Symbolics.build_function(Σ_symbolic, par, expression=Val{false})[2]
    Σ = zeros(size(Σ_symbolic))
    precompile(Σ_function, (typeof(Σ), Vector{Float64}))

    # ∇Σ
    if gradient
        ∇Σ_symbolic = Symbolics.sparsejacobian(vec(Σ_symbolic), [par...])
        ∇Σ_function = Symbolics.build_function(∇Σ_symbolic, par, expression=Val{false})[2]
        constr = findnz(∇Σ_symbolic)
        ∇Σ = sparse(constr[1], constr[2], fill(1.0, nnz(∇Σ_symbolic)), size(∇Σ_symbolic)...)
        precompile(∇Σ_function, (typeof(∇Σ), Vector{Float64}))
    else
        ∇Σ_symbolic = nothing
        ∇Σ_function = nothing
        ∇Σ = nothing
    end

    if hessian & !approximate_hessian
        n_sig = length(Σ_symbolic)
        n_par = size(par, 1)
        ∇²Σ_symbolic_vec = [Symbolics.sparsehessian(σᵢ, [par...]) for σᵢ in vec(Σ_symbolic)]

        @variables J[1:n_sig]
        ∇²Σ_symbolic = zeros(Num, n_par, n_par)
        for i in 1:n_sig
            ∇²Σ_symbolic += J[i]*∇²Σ_symbolic_vec[i]
        end

        ∇²Σ_function = Symbolics.build_function(∇²Σ_symbolic, J, par, expression=Val{false})[2]
        ∇²Σ = zeros(n_par, n_par)
    else
        ∇²Σ_symbolic = nothing
        ∇²Σ_function = nothing
        ∇²Σ = nothing
    end

    # μ
    if meanstructure
        has_meanstructure = Val(true)
        μ_symbolic = get_μ_symbolic_RAM(M, A, F)
        μ_function = Symbolics.build_function(μ_symbolic, par, expression=Val{false})[2]
        μ = zeros(size(μ_symbolic))
        if gradient
            ∇μ_symbolic = Symbolics.jacobian(μ_symbolic, [par...])
            ∇μ_function = Symbolics.build_function(∇μ_symbolic, par, expression=Val{false})[2]
            ∇μ = zeros(size(F, 1), size(par, 1))
        else
            ∇μ_function = nothing
            ∇μ = nothing
        end
    else
        has_meanstructure = Val(false)
        μ_function = nothing
        μ = nothing
        ∇μ_function = nothing
        ∇μ = nothing
    end

    return RAMSymbolic(
        Σ_function,
        ∇Σ_function,
        ∇²Σ_function,
        Σ,
        ∇Σ,
        ∇²Σ,
        Σ_symbolic,
        ∇Σ_symbolic,
        ∇²Σ_symbolic,
        n_par,
        ram_matrices,
        μ_function,
        μ,
        ∇μ_function,
        ∇μ,
        identifier,
        has_meanstructure
    )
end

############################################################################################
### objective, gradient, hessian
############################################################################################

# dispatch on meanstructure
objective!(imply::RAMSymbolic, par, model) =
    objective!(imply, par, model, imply.has_meanstructure)
gradient!(imply::RAMSymbolic, par, model) =
    gradient!(imply, par, model, imply.has_meanstructure)

# objective
function objective!(imply::RAMSymbolic, par, model, has_meanstructure::Val{T}) where T
    imply.Σ_function(imply.Σ, par)
    T && imply.μ_function(imply.μ, par)
end

# gradient
function gradient!(imply::RAMSymbolic, par, model, has_meanstructure::Val{T}) where T
    objective!(imply, par, model, imply.has_meanstructure)
    imply.∇Σ_function(imply.∇Σ, par)
    T && imply.∇μ_function(imply.∇μ, par)
end

# other methods
hessian!(imply::RAMSymbolic, par, model) = gradient!(imply, par, model)
objective_gradient!(imply::RAMSymbolic, par, model) = gradient!(imply, par, model)
objective_hessian!(imply::RAMSymbolic, par, model) = gradient!(imply, par, model)
gradient_hessian!(imply::RAMSymbolic, par, model) = gradient!(imply, par, model)
objective_gradient_hessian!(imply::RAMSymbolic, par, model) = gradient!(imply, par, model)

############################################################################################
### Recommended methods
############################################################################################

identifier(imply::RAMSymbolic) = imply.identifier
n_par(imply::RAMSymbolic) = imply.n_par

function update_observed(imply::RAMSymbolic, observed::SemObserved; kwargs...)
    if n_man(observed) == size(imply.Σ)
        return imply
    else
        return RAMSymbolic(;observed = observed, kwargs...)
    end
end

############################################################################################
### additional methods
############################################################################################

Σ(imply::RAMSymbolic) = imply.Σ
∇Σ(imply::RAMSymbolic) = imply.∇Σ
∇²Σ(imply::RAMSymbolic) = imply.∇²Σ

μ(imply::RAMSymbolic) = imply.μ
∇μ(imply::RAMSymbolic) = imply.∇μ

Σ_function(imply::RAMSymbolic) = imply.Σ_function
∇Σ_function(imply::RAMSymbolic) = imply.∇Σ_function
∇²Σ_function(imply::RAMSymbolic) = imply.∇²Σ_function

has_meanstructure(imply::RAMSymbolic) = imply.has_meanstructure

ram_matrices(imply::RAMSymbolic) = imply.ram_matrices

############################################################################################
### additional functions
############################################################################################

function get_Σ_symbolic_RAM(S, A, F; vech = false)
    invia = neumann_series(A)
    Σ_symbolic = F*invia*S*permutedims(invia)*permutedims(F)
    Σ_symbolic = Array(Σ_symbolic)
    # Σ_symbolic = Symbolics.simplify.(Σ_symbolic)
    Threads.@threads for i in eachindex(Σ_symbolic)
        Σ_symbolic[i] = Symbolics.simplify(Σ_symbolic[i])
    end
    if vech Σ_symbolic = Σ_symbolic[tril(trues(size(F, 1), size(F, 1)))] end
    return Σ_symbolic
end

function get_μ_symbolic_RAM(M, A, F)
    invia = neumann_series(A)
    μ_symbolic = F*invia*M
    μ_symbolic = Array(μ_symbolic)
    Threads.@threads for i in eachindex(μ_symbolic)
        μ_symbolic[i] = Symbolics.simplify(μ_symbolic[i])
    end
    return μ_symbolic
end


text:
"""
Package: ForecastData

art_dir(name)

Retrieve the directory where artifact `name` name is located

"""
function art_path(name)

    art_toml = joinpath(splitdir(@__DIR__)[1], "Artifacts.toml")
    art_hash = artifact_hash(name, art_toml)

    if isnothing(art_hash)
        @info "Not available dataset named: `$(name)`"
        return nothing
    end

    ensure_artifact_installed(name,art_toml)

    # art_file = artifact_path(art_hash) * "/" * name * ".tar.gz"
    # tmp_path = mktempdir("/tmp"; prefix="fdata_", cleanup=true)
    # Tar.extract(IOBuffer(inflate_gzip(art_file)), tmp_path)

    artifact_path(art_hash)

end


text:
const Float = Float64

SVectorF{N} = SVector{N,Float}
SVectorU{N} = SVector{N,UInt}

cross(a, b) = [
    a[2] * b[3] - a[3] * b[2],
    a[3] * b[1] - a[1] * b[3],
    a[1] * b[2] - a[2] * b[1],
]

max_stars_repo_name: lars-frogner/PIC

max_issues_repo_name: lars-frogner/PIC

max_forks_repo_name: lars-frogner/PIC


text:
function intergrid_1d_interp(N, m)

    ic = Int(ceil(m/2.0))
    Minv = stencil1d(m)
    sten = transpose(Minv) * Float64[(0.5)^(j-1)/factorial(j-1) for j=1:m]

    Ih = zeros(Float64, (N,N))
    Is = I + Ih
    for i=1:m
        k = i - ic
        vec = sten[i] * ones(N - abs(k))
        Ih += diagm(k => vec)
    end

    Iout = zeros(Float64, (2*N,N))
    Iout[1:2:end,:] = Is
    Iout[2:2:end,:] = Ih
    return Iout
end

function intergrid_transforms(Nx0, Ny0, mx, my, P0, P1)

    Ix = intergrid_1d_interp(Nx0, mx)
    Iy = intergrid_1d_interp(Ny0, my)
    Ic = kron(sparse(Iy), sparse(Ix))
    Rc = transpose(Ic) / 4.0

    Restr = P0 * Rc * transpose(P1)
    Interp = P1 * Ic * transpose(P0)
    return Restr, Interp
end

function edge_trim_box(walls::Vector{Wall}, deltas::Vector{Float64}, Nx::Int64, Ny::Int64, bottom_left::Vector{Float64}, top_right::Vector{Float64})

    deltaX = top_right[1] - bottom_left[1]
    deltaY = top_right[2] - bottom_left[2]
    dx = deltaX / Nx
    dy = deltaY / Ny

    points = zeros(Float64, (2, Nx*Ny))
    proj_rows = Int32[i for i=1:Nx*Ny]
    proj_cols = zeros(Int32, Nx*Ny)
    proj_vals = ones(Int32, Nx*Ny)
    nanmask = fill(NaN, (Nx, Ny))
    Nk = 0

    for j=1:Ny
        for i=1:Nx

            x = bottom_left[1] + (i-1)*dx
            y = bottom_left[2] + (j-1)*dx
            if inside(x, y, deltas, walls)
                Nk += 1
                points[:,Nk] = [x, y]
                proj_cols[Nk] = i + Nx*(j-1)
                nanmask[i,j] = 1.0
            end
        end
    end

    P = sparse(proj_rows[1:Nk], proj_cols[1:Nk], proj_vals[1:Nk], Nk, Nx*Ny)
    return points[:,1:Nk], P, nanmask, dx, dy
end

function edge_trim_double(points0::Matrix{Float64}, walls::Vector{Wall}, deltas::Vector{Float64}, dx0::Float64, dy0::Float64, Nx0::Int64, Ny0::Int64, bottom_left::Vector{Float64})

    dx = dx0 / 2.0
    dy = dy0 / 2.0
    Nk0 = size(points0)[2]

    points = zeros(Float64, (2, 4*Nk0))
    proj_rows = Int32[i for i=1:4*Nk0]
    proj_cols = zeros(Int32, 4*Nk0)
    proj_vals = ones(Int32, 4*Nk0)
    nanmask = fill(NaN, (2*Nx0, 2*Ny0))
    Nk = 0

    for l=1:Nk0
        for j=1:2
            for i=1:2
                x = points0[1,l] + (i-1)*dx
                y = points0[2,l] + (j-1)*dy
                deltaX = x - bottom_left[1]
                deltaY = y - bottom_left[2]
                ic = Int(round(deltaX/dx)) + 1
                jc = Int(round(deltaY/dy)) + 1

                if inside(x, y, deltas, walls)
                    Nk += 1
                    points[:,Nk] = [x, y]
                    proj_cols[Nk] = ic + 2*Nx0*(jc-1)
                    nanmask[ic, jc] = 1.0
                end

            end
        end
    end

    P = sparse(proj_rows[1:Nk], proj_cols[1:Nk], proj_vals[1:Nk], Nk, 4*Nx0*Ny0)
    return points[:,1:Nk], P, nanmask, dx, dy
end

#
#
#
#
#
#
# function edge_trim_double(Nx0::Int64, Ny0::Int64, dx0::Float64, dy0::Float64, points0::Matrix{Float64}, proj_cols0::Vector{Int32}, deltas::Vector{Float64}, walls::Vector{Wall})
#
#     dx = dx0 / 2.0
#     dy = dy0 / 2.0
#     Nk0 = size(points0)[2]
#
#     points = zeros(Float64, (2, 4*Nk0))
#     proj_rows = Int32[i for i=1:4*Nk0]
#     proj_cols = zeros(Int32, 4*Nk0)
#     proj_vals = ones(Int32, 4*Nk0)
#     nanmask = fill(NaN, (2*Nx0, 2*Ny0))
#     Nk = 0
#
#     for l=1:Nk0
#         for j=1:2
#             for i=1:2
#                 x = points0[1,l] + (i-1)*dx
#                 y = points0[2,l] + (j-1)*dy
#
#                 i0 = rem(proj_cols0[l], Nx0)
#                 j0 = div(proj_cols0[l], Ny0)
#                 i2 = (i0-1)*2 + i
#                 j2 = (j0-1)*2 + j
#
#                 if inside(x, y, deltas, walls)
#                     Nk += 1
#                     points[:,Nk] = [x, y]
#                     proj_cols[Nk] = i2 + 2*Nx0*(j2-1)
#                     nanmask[i2, j2] = 1.0
#                 end
#
#             end
#         end
#     end
#
#     P = sparse(proj_rows[1:Nk], proj_cols[1:Nk], proj_vals[1:Nk], Nk, 4*Nx0*Ny0)
#     return points[:,1:Nk], P, nanmask, dx, dy
# end

struct Grid

    points::Matrix{Float64}
    origin::Vector{Float64}
    projection::SparseMatrixCSC
    inverse_projection::SparseMatrixCSC
    restrictions::Vector{SparseMatrixCSC}
    interpolations::Vector{SparseMatrixCSC}
    nanmask::Matrix{Float64}
    mx::Int64
    my::Int64
    Nx::Int64
    Ny::Int64
    Nk::Int64
    Mxinv::Matrix{Float64}
    Myinv::Matrix{Float64}
    MxyinvT::Matrix{Float64}
    dx::Float64
    dy::Float64

end

function Grid(walls::Vector{Wall}, deltas::Vector{Float64}, corners::Matrix{Float64}, N_levels::Int64, N0::Vector{Int}, m::Vector{Int})

    Nx, Ny = N0
    mx, my = m

    # trim the bounding box, then double the resolutions and compute Restr and Interp
    points, P0, nanmask, dx, dy = edge_trim_box(walls, deltas, Nx, Ny, corners[:,1], corners[:,2])
    points, P1, nanmask, dx, dy = edge_trim_double(points, walls, deltas/2, dx, dy, Nx, Ny, corners[:,1])
    Restr, Interp = intergrid_transforms(Nx, Ny, mx+1, my+1, P0, P1)
    Nx = Nx * 2
    Ny = Ny * 2

    restrictions = SparseMatrixCSC[Restr]
    interpolations = SparseMatrixCSC[Interp]

    for l=3:N_levels

        P0 = P1
        points, P1, nanmask, dx, dy = edge_trim_double(points, walls, deltas/l, dx, dy, Nx, Ny, corners[:,1])
        Restr, Interp = intergrid_transforms(Nx, Ny, mx+1, my+1, P0, P1)
        prepend!(restrictions, [Restr])
        append!(interpolations, [Interp])
        Nx = Nx * 2
        Ny = Ny * 2
    end

    origin = corners[:,1]
    projection = P1
    inverse_projection = transpose(P1)
    Nk = size(points)[2]
    Mxinv = stencil1d(mx)
    Myinv = stencil1d(my)
    MxyinvT = stencil2d(mx, my)

    return Grid(points, origin, projection, inverse_projection, restrictions, interpolations,
                nanmask, mx, my, Nx, Ny, Nk, Mxinv, Myinv, MxyinvT, dx, dy)
end

function Grid(walls::Vector{Wall}, deltas::Vector{Float64}, corners::Matrix{Float64}, N_levels::Int64, N0::Int64, m::Int64)

    return Grid(walls, deltas, corners, N_levels, [N0, N0], [m, m])
end

function vec_to_mesh(vec, grid::Grid)

    vals = grid.inverse_projection * vec

    Nx = grid.Nx
    Ny = grid.Ny
    out = zeros(Float64, (Nx, Ny))

    for j=1:Ny
        out[:,j] = vals[1+(j-1)*Nx:j*Nx]
    end

    return out .* grid.nanmask
end

function f_to_grid(f::Function, grid::Grid)

    Nk = grid.Nk
    vec = zeros(Float64, Nk)

    for k=1:Nk
        x = grid.points[1,k]
        y = grid.points[2,k]
        vec[k] = f(x, y)
    end

    return vec
end

# overloading Base.map() to operate on the points in the grid, using pmap instead of usual map
function grid_map(f::Function, outsize::Tuple{Int64, Int64}, grd::Grid)

    n, m = outsize
    output = SharedArray{Float64}((n*grd.Nk,m))

    # @distributed unfortunately only works if nprocs > 1
    if nprocs() > 1
        @distributed for i=1:grd.Nk
            k = 1 + n*(i-1)
            output[k:k+n-1,:] = f(grd.points[:,i]..., grd)
        end
    else
        for i=1:grd.Nk
            k = 1 + n*(i-1)
            output[k:k+n-1,:] = f(grd.points[:,i]..., grd)
        end
    end

    return output[:,:]
end


text:
# This file is a part of DensityInterface.jl, licensed under the MIT License (MIT).

"""
    abstract type DensityKind end

    DensityKind(object)

Subtypes of `DensityKind` indicate if an `object` *is* a density or if it *has*
a density, in the sense of the `DensityInterface` API, or if is *not*
associated with a density (not compatible with `DensityInterface`).

`DensityKind(object)` returns either `IsDensity()`, `HasDensity()` or
`NoDensity()`.

In addition to the subtypes [`IsDensity`](@ref), [`HasDensity`](@ref) or
[`NoDensity`](@ref), a union `IsOrHasDensity = Union{IsDensity, HasDensity}`
is defined for convenience.

`DensityKind(object) isa IsOrHasDensity` implies that `object` is either a
density itself or can be said to have an associated density. It also implies
that the value of that density at given points can be calculated via
[`logdensityof`](@ref) and [`densityof`](@ref).

`DensityKind(object)` defaults to `NoDensity()` (object is not and does not
have a density). For a type that *is* (directly represents) a density, like a
probability density, define

```julia
@inline DensityKind(::MyDensityType) = IsDensity()
```

For a type that *has* (is associated with) a density in some way, like
a probability distribution has a probability density, define

```julia
@inline DensityKind(::MyDensityType) = HasDensity()
```
"""
abstract type DensityKind end
export DensityKind

@inline DensityKind(object) = NoDensity()

"""
    struct IsDensity <: DensityKind end

As a return value of [`DensityKind(object)`](@ref), indicates that
`object` *is* (represents) a density, like a probability density
object.

See [`DensityKind`](@ref) for details.
"""
struct IsDensity <: DensityKind end
export IsDensity

"""
    struct HasDensity <: DensityKind end

As a return value of [`DensityKind(object)`](@ref), indicates that
`object` *has* a density, like a probability distribution has
a probability density.

See [`DensityKind`](@ref) for details.
"""
struct HasDensity <: DensityKind end
export HasDensity

"""
    struct NoDensity <: DensityKind end

As a return value of [`DensityKind(object)`](@ref), indicates that
`object` *is not* and *does not have* a density, as understood by
`DensityInterface`.

See [`DensityKind`](@ref) for details.
"""
struct NoDensity <: DensityKind end
export NoDensity

"""
    IsOrHasDensity = Union{IsDensity, HasDensity}

As a return value of [`DensityKind(object)`](@ref), indicates that `object`
either *is* or *has* a density, as understood by `DensityInterface`.

See [`DensityKind`](@ref) for details.
"""
const IsOrHasDensity = Union{IsDensity, HasDensity}
export IsOrHasDensity

function _check_is_or_has_density(object)
    DensityKind(object) isa IsOrHasDensity || throw(ArgumentError("Object of type $(typeof(object)) neither is nor has a density"))
end

"""
    logdensityof(object, x)::Real

Compute the logarithmic value of the density `object` (resp. its associated density)
at a given point `x`.

```jldoctest a
julia> DensityKind(object)
IsDensity()

julia> logy = logdensityof(object, x); logy isa Real
true
```

See also [`DensityKind`](@ref) and [`densityof`](@ref).
"""
function logdensityof end
export logdensityof

"""
    logdensityof(object)

Return a function that computes the logarithmic value of the density `object`
(resp. its associated density) at a given point.

```jldoctest a
julia> log_f = logdensityof(object); log_f isa Function
true

julia> log_f(x) == logdensityof(object, x)
true
```

`logdensityof(object)` defaults to `Base.Fix1(logdensityof, object)`, but may be
specialized. If so, [`logfuncdensity`](@ref) will typically have to be
specialized for the return type of `logdensityof` as well.

[`logfuncdensity`](@ref) is the inverse of `logdensityof`, so
`logfuncdensity(log_f)` must be equivalent to `object`.
"""
function logdensityof(object)
    _check_is_or_has_density(object)
    Base.Fix1(logdensityof, object)
end

"""
    densityof(object, x)::Real

Compute the value of the density `object` (resp. its associated density)
at a given point `x`.

```jldoctest a
julia> DensityKind(object)
IsDensity()

julia> densityof(object, x) == exp(logdensityof(object, x))
true
```

`densityof(object, x)` defaults to `exp(logdensityof(object, x))`, but
may be specialized.

See also [`DensityKind`](@ref) and [`densityof`](@ref).
"""
densityof(object, x) = exp(logdensityof(object, x))
export densityof

"""
    densityof(object)

Return a function that computes the value of the density `object`
(resp. its associated density) at a given point.

```jldoctest a
julia> f = densityof(object);

julia> f(x) == densityof(object, x)
true
```

`densityof(object)` defaults to `Base.Fix1(densityof, object)`, but may be specialized.
"""
function densityof(object)
    _check_is_or_has_density(object)
    Base.Fix1(densityof, object)
end

"""
    logfuncdensity(log_f)

Return a `DensityInterface`-compatible density that is defined by a given
log-density function `log_f`:

```jldoctest
julia> object = logfuncdensity(log_f);

julia> DensityKind(object)
IsDensity()

julia> logdensityof(object, x) == log_f(x)
true
```

`logfuncdensity(log_f)` returns an instance of [`DensityInterface.LogFuncDensity`](@ref)
by default, but may be specialized to return something else depending on the
type of `log_f`). If so, [`logdensityof`](@ref) will typically have to be
specialized for the return type of `logfuncdensity` as well.

`logfuncdensity` is the inverse of `logdensityof`, so the following must
hold true:

* `d = logfuncdensity(logdensityof(object))` is equivalent to `object` in
  respect to `logdensityof` and `densityof`. However, `d` may not be equal to
  `object`, especially if `DensityKind(object) == HasDensity()`: `logfuncdensity` always
  creates something that *is* density, never something that just *has*
  a density in some way (like a distribution or a measure in general).
* `logdensityof(logfuncdensity(log_f))` is equivalent (typically equal or even
  identical to) to `log_f`.

See also [`DensityKind`](@ref).
"""
function logfuncdensity end
export logfuncdensity

@inline logfuncdensity(log_f) = LogFuncDensity(log_f)

# For functions stemming from objects that *have* a density, create a new density:
@inline _logfuncdensity_impl(::HasDensity, log_f::Base.Fix1{typeof(logdensityof)}) = LogFuncDensity(log_f)
# For functions stemming from objects that *are* a density, recover original object:
@inline _logfuncdensity_impl(::IsDensity, log_f::Base.Fix1{typeof(logdensityof)}) = log_f.x

@inline logfuncdensity(log_f::Base.Fix1{typeof(logdensityof)}) = _logfuncdensity_impl(DensityKind(log_f.x), log_f)

InverseFunctions.inverse(::typeof(logfuncdensity)) = logdensityof
InverseFunctions.inverse(::typeof(logdensityof)) = logfuncdensity

"""
    struct DensityInterface.LogFuncDensity{F}

Wraps a log-density function `log_f` to make it compatible with `DensityInterface`
interface. Typically, `LogFuncDensity(log_f)` should not be called
directly, [`logfuncdensity`](@ref) should be used instead.
"""
struct LogFuncDensity{F}
    _log_f::F
end
LogFuncDensity

@inline DensityKind(::LogFuncDensity) = IsDensity()

@inline logdensityof(object::LogFuncDensity, x) = object._log_f(x)
@inline logdensityof(object::LogFuncDensity) = object._log_f

@inline densityof(object::LogFuncDensity, x) = exp(object._log_f(x))
@inline densityof(object::LogFuncDensity) = exp ∘ object._log_f

function Base.show(io::IO, object::LogFuncDensity)
    print(io, nameof(typeof(object)), "(")
    show(io, object._log_f)
    print(io, ")")
end

"""
    funcdensity(f)

Return a `DensityInterface`-compatible density that is defined by a given
non-log density function `f`:

```jldoctest
julia> object = funcdensity(f);

julia> DensityKind(object)
IsDensity()

julia> densityof(object, x) == f(x)
true
```

`funcdensity(f)` returns an instance of [`DensityInterface.FuncDensity`](@ref)
by default, but may be specialized to return something else depending on the
type of `f`). If so, [`densityof`](@ref) will typically have to be
specialized for the return type of `funcdensity` as well.

`funcdensity` is the inverse of `densityof`, so the following must
hold true:

* `d = funcdensity(densityof(object))` is equivalent to `object` in
  respect to `logdensityof` and `densityof`. However, `d` may not be equal to
  `object`, especially if `DensityKind(object) == HasDensity()`: `funcdensity` always
  creates something that *is* density, never something that just *has*
  a density in some way (like a distribution or a measure in general).
* `densityof(funcdensity(f))` is equivalent (typically equal or even
  identical to) to `f`.

See also [`DensityKind`](@ref).
"""
function funcdensity end
export funcdensity

@inline funcdensity(f) = FuncDensity(f)

# For functions stemming from objects that *have* a density, create a new density:
@inline _funcdensity_impl(::HasDensity, f::Base.Fix1{typeof(densityof)}) = FuncDensity(f)
# For functions stemming from objects that *are* a density, recover original object:
@inline _funcdensity_impl(::IsDensity, f::Base.Fix1{typeof(densityof)}) = f.x

@inline funcdensity(f::Base.Fix1{typeof(densityof)}) = _funcdensity_impl(DensityKind(f.x), f)

InverseFunctions.inverse(::typeof(funcdensity)) = densityof
InverseFunctions.inverse(::typeof(densityof)) = funcdensity

"""
    struct DensityInterface.FuncDensity{F}

Wraps a non-log density function `f` to make it compatible with
`DensityInterface` interface. Typically, `FuncDensity(f)` should not be
called directly, [`funcdensity`](@ref) should be used instead.
"""
struct FuncDensity{F}
    _f::F
end
FuncDensity

@inline DensityKind(::FuncDensity) = IsDensity()

@inline logdensityof(object::FuncDensity, x) = log(object._f(x))
@inline logdensityof(object::FuncDensity) = log ∘ object._f

@inline densityof(object::FuncDensity, x) = object._f(x)
@inline densityof(object::FuncDensity) = object._f

function Base.show(io::IO, object::FuncDensity)
    print(io, nameof(typeof(object)), "(")
    show(io, object._f)
    print(io, ")")
end


text:
function read_single(filename::String)
    ret = Int64[]
    open(filename) do f
        for line in eachline(f)
            push!(ret, parse(Int64, line))
        end
    end
    return ret
end

function read_node_map(filename::String)
    nm = Dict{Int64,Int64}()
    open(filename) do f
        for (i, line) in enumerate(eachline(f))
            if i == 1; continue; end
            new, orig = [parse(Int64, v) for v in split(line)]
            nm[orig] = new
        end
    end
    return nm
end

function read_labels(filename::String, nm, n::Int64, incl_gender::Bool)
    labels = zeros(Int64, n)
    label_map = Dict{AbstractString,Int64}()

    open(filename) do f
        for line in eachline(f)
            data = split(line)
            orig_id = parse(Int64, data[1])
            # Some nodes do not participate in any edge
            if !haskey(nm, orig_id); continue; end
            node = nm[orig_id]

            # Label is either class or class and gender (note: some genders listed as "UNKNOWN")
            label = data[2]
            if incl_gender
                label = string(label, " ", data[3])
            end
            if !haskey(label_map, label)
                label_map[label] = length(label_map) + 1
            end
            labels[node] = label_map[label]
        end
    end

    return labels, label_map
end

function main()
    for dataset in ["high-school", "primary-school"]
        for incl_gender_in_label in [true, false]
            # Read data
            nverts = read_single("raw-$dataset/contact-$dataset-nverts.txt")
            simplices = read_single("raw-$dataset/contact-$dataset-simplices.txt")
            nm = read_node_map("raw-$dataset/contact-$dataset-nodemap.txt")
            metadata =
                if     dataset == "high-school";    "metadata_2013"
                elseif dataset == "primary-school"; "metadata_primaryschool"
                else                                ""
                end
            labels, label_map = read_labels("raw-$dataset/$metadata.txt", nm, maximum(simplices),
                                            incl_gender_in_label)

            # Get all unique hyperedges
            hedges = Set{Set{Int64}}()
            let curr_ind = 0
	        for nvert in nverts
    	            simplex = simplices[(curr_ind + 1):(curr_ind + nvert)]
                    push!(hedges, Set{Int64}(simplex))
	            curr_ind += nvert
                end
            end

            # setup directory
            dir = string("contact-", dataset, "-classes")
            if incl_gender_in_label
                dir = string(dir, "-gender")
            end
            if !isdir(dir)
                mkdir(dir)
            end

            # hyperedge data
            open("$dir/hyperedges-$dir.txt", "w") do f
                for hedge in hedges
                    write(f, join(sort(collect(hedge)), ','))
                    write(f, "\n")
                end
            end

            # names of the integer labels
            open("$dir/label-names-$dir.txt", "w") do f
                output = sort([(v, k) for (k, v) in label_map])
                for (v, k) in output
                    write(f, "$k\n")
                end
            end

            # labels for each node
            open("$dir/node-labels-$dir.txt", "w") do f
                for l in labels
                    write(f, "$l\n")
                end
            end
        end
    end
end

max_stars_repo_name: PhilChodrow/HypergraphModularity

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: PhilChodrow/HypergraphModularity