text:
#-----------------------------------------------------------------------------------------------------------------------------
#
#   Interaction matrix INFORMATION
#   (returns two matrices that indicate the number of periodic and non-periodic connections between respective sites)
#
#   Parameters are
#   - lattice: The complete lattice of which the interaction matrix should be constructed
#   or
#   - unitcell: The unitcell of which the interaction matrix should be constructed
#
#-----------------------------------------------------------------------------------------------------------------------------
function getInteractionMatrixInformation(lattice::Lattice)
    # matrices indicating the number of connections
    con_periodic =  zeros(Int64, size(lattice.positions,1),size(lattice.positions,1))
    con_inside   =  zeros(Int64, size(lattice.positions,1),size(lattice.positions,1))
    # go through all connections of the lattice
    for c in lattice.connections
        # get the indices
        index_from  = Int(c[1])
        index_to    = Int(c[2])
        strength    = c[3]
        wrap        = c[4]
        # check if inside or periodic and add to the counter
        if sum([abs(el) for el in wrap]) == 0
            con_inside[index_from, index_to] = con_inside[index_from, index_to] + 1
        else
            con_periodic[index_from, index_to] = con_periodic[index_from, index_to] + 1
        end
    end
    # return the matrices
    return con_inside, con_periodic
end
function getInteractionMatrixInformation(unitcell::Unitcell)
    # matrices indicating the number of connections
    con_periodic =  zeros(Int64, size(unitcell.basis,1),size(unitcell.basis,1))
    con_inside   =  zeros(Int64, size(unitcell.basis,1),size(unitcell.basis,1))
    # go through all connections of the lattice
    for c in unitcell.connections
        # get the indices
        index_from  = Int(c[1])
        index_to    = Int(c[2])
        strength    = c[3]
        wrap        = c[4]
        # check if inside or periodic and add to the counter
        if sum([abs(el) for el in wrap]) == 0
            con_inside[index_from, index_to] = con_inside[index_from, index_to] + 1
        else
            con_periodic[index_from, index_to] = con_periodic[index_from, index_to] + 1
        end
    end
    # return the matrices
    return con_inside, con_periodic
end

# print the interaction matrix information in compressed LaTeX format
function printInteractionMatrixInformation(connections_inside, connections_periodic)
    # get the dimension of the matrices
    len = size(connections_periodic,1)
    # starting lines
    println("\\begin{equation*}")
    println("   \\mathbf{H} = \\bordermatrix{")
    print("       ")
    for i in 1:len
        print("& $(i) ")
    end
    println("\\cr")
    # iterate over all lines of the matrix
    for j in 1:len
        # start line by printing the current index
        print("       $(j) ")
        # go through all respective numbers
        for i in 1:len
            # check if periodic or internal connections
            if connections_inside[i,j] == 0 && connections_periodic[i,j] == 0
                print("& 0 ")
            elseif connections_inside[i,j] > 0 && connections_periodic[i,j] == 0
                #print("& {\\color{blue} 0},{\\color{red} $(connections_inside[i,j])} ")
                print("& {\\color{blue} $(connections_inside[i,j])} ")
            elseif connections_inside[i,j] == 0 && connections_periodic[i,j] > 0
                #print("& {\\color{blue} $(connections_periodic[i,j])},{\\color{red} 0} ")
                print("& {\\color{red} $(connections_periodic[i,j])} ")
            else
                print("& {\\color{red} $(connections_periodic[i,j])},{\\color{blue} $(connections_inside[i,j])} ")
            end
        end
        # end the line by printing a linebreak
        println("\\cr")
    end
    # end lines
    println("   }")
    println("\\end{equation*}")
end
function printInteractionMatrixInformation(lattice::Lattice)
    println("Interaction matrix for lattice \"$(lattice.filename)\"")
    connections_inside, connections_periodic = getInteractionMatrixInformation(lattice)
    printInteractionMatrixInformation(connections_inside, connections_periodic)
end
function printInteractionMatrixInformation(unitcell::Unitcell)
    println("Interaction matrix for unitcell \"$(unitcell.filename)\"")
    connections_inside, connections_periodic = getInteractionMatrixInformation(unitcell)
    printInteractionMatrixInformation(connections_inside, connections_periodic)
end
export printInteractionMatrixInformation

# function to show the lattice
function showLattice(
		lattice::Lattice;
		conversion = 160,
		site_radius=25,
        site_labels="OFF",
		bond_thickness::Int64=8,
		visualize_periodic=false,
		colorcode_sites = Dict(0 => [255,255,255], 1 => [255,255,255]),
		colorcode_bonds = Dict("0" => [0,0,0], "1.0" => [0,0,0]),
        colorcode_bonds_automation::String = "OFF",
        background_color = [200,200,200],
        forcePyPlot=false
    )
    # if forcing PyPlot, do this and exit
    if forcePyPlot
        showLatticePyPlot(
            lattice,
            conversion = conversion,
            site_radius=site_radius,
            site_labels=site_labels,
            bond_thickness=bond_thickness,
            visualize_periodic=visualize_periodic,
            colorcode_sites = colorcode_sites,
            colorcode_bonds = colorcode_bonds,
            colorcode_bonds_automation=colorcode_bonds_automation,
            background_color = background_color
        )
        return
    end
    # CHECK IF MAYAVI IS PRESENT
    try
        @pyimport mayavi.mlab as mlab
        MAYAVI_AVAILABLE = true
        try
        showLatticeMayavi(
            lattice,
            conversion = conversion,
            site_radius=site_radius,
            site_labels=site_labels,
            bond_thickness=bond_thickness,
            visualize_periodic=visualize_periodic,
            colorcode_sites = colorcode_sites,
            colorcode_bonds = colorcode_bonds,
            colorcode_bonds_automation=colorcode_bonds_automation,
            background_color = background_color
        )
        catch x
        end
        return
    catch error
        if isa(error, PyCall.PyError)
            #println("PyError occured when importing:")
            #println(error)
            MAYAVI_AVAILABLE = false
            #println("Using PyPlot now")
            try
                showLatticePyPlot(
                    lattice,
                    conversion = conversion,
                    site_radius=site_radius,
                    site_labels=site_labels,
                    bond_thickness=bond_thickness,
                    visualize_periodic=visualize_periodic,
                    colorcode_sites = colorcode_sites,
                    colorcode_bonds = colorcode_bonds,
                    colorcode_bonds_automation=colorcode_bonds_automation,
                    background_color = background_color
                )
            catch error2
                println("Error occured when plotting:")
                println(error2)
            end
        else
            println("Error occured when loading:")
            println(error)
        end
    end

end
export showLattice

# show Lattice function to plot in pyplot
function showLatticePyPlot(
		lattice::Lattice;
		conversion = 160,
		site_radius=25,
        site_labels="OFF",
		bond_thickness::Int64=8,
		visualize_periodic=false,
		colorcode_sites = Dict(0 => [255,255,255], 1 => [255,255,255]),
		colorcode_bonds = Dict("0" => [0,0,0], "1.0" => [0,0,0]),
        colorcode_bonds_automation::String = "OFF",
        background_color = [200,200,200]
    )

    # maybe overwrite dictonary
    if colorcode_bonds_automation == "GREY"
        # construct the list of interaction strengths
        cs_list = getConnectionStrengthList(lattice)
        # get the color code list
        cc_list = getGreySequence(size(cs_list, 1))
        # put in a new dictonary
        colorcode_bonds = Dict()
        # insert all pairs
        for i in 1:size(cc_list, 1)
            colorcode_bonds[string(cs_list[i])] = cc_list[i]
        end
    elseif colorcode_bonds_automation == "COLOR"
        # construct the list of interaction strengths
        cs_list = getConnectionStrengthList(lattice)
        # get the color code list
        cc_list = getColorSequence(size(cs_list, 1))
        # put in a new dictonary
        colorcode_bonds = Dict()
        # insert all pairs
        for i in 1:size(cc_list, 1)
            colorcode_bonds[string(cs_list[i])] = cc_list[i]
        end
    end

	# repair color dictonary
	colorcode_bonds["0"] = get(colorcode_bonds, "0", [0,0,0])
	colorcode_sites[0] = get(colorcode_sites, 0, [255,255,255])

    # PLOTTING WITH PYPLOT
    # create a figure
    fig = figure(figsize=(10,10), facecolor=background_color ./ 255.0)
    # add a 3d projection subplot
    ax = fig[:add_subplot](111, projection="3d", axisbg=background_color ./ 255.0)
    # define a function to plot a sphere for every site
    function plotSite(site,radius,color; detail=10)
        x = site[1]
        y = site[2]
        if length(site) == 3
            z = site[3]
        else
            z = 0.0
        end
        # Make data
        n = detail
        u = linspace(0,2*pi,n);
        v = linspace(0,pi,n);
        px = cos(u) * sin(v)';
        py = sin(u) * sin(v)';
        pz = ones(n) * cos(v)';
        px = (px.*radius) .+ x
        py = (py.*radius) .+ y
        pz = (pz.*radius) .+ z
        # Plot the surface
        surf(px,py,pz, rstride=1, cstride=1, linewidth=0, antialiased=true, color=color)
    end
    # define a function to plot a bond as a tube
    function plotBond(from, to, radius, color; detail_ring=8, detail_length=3)
        # check to be 3d
        if length(from) == 2
            from = [from[1], from[2], 0.0]
        end
        if length(to) == 2
            to = [to[1], to[2], 0.0]
        end
        # https://de.mathworks.com/matlabcentral/answers/151235-tube-plot-with-x-y-coordinates-and-radius?requestedDomain=www.mathworks.com
        # plot bond
        vec_u = to .- from;
        t = nullspace(vec_u')';
        vec_v = t[1,:]
        vec_w = t[2,:]
        # normalize
        vec_v = radius.* vec_v ./ norm(vec_v)
        vec_w = radius.* vec_w ./ norm(vec_w)
        #println("v=$(vec_v)    w=$(vec_w)")
        #[S,T] = linspace(0,1,detail)*linspace(0,2*pi,detail)';
        S = linspace(0,1,detail_length)
        T = linspace(0,2*pi,detail_ring)
        px = zeros(detail_length, detail_ring)
        py = zeros(detail_length, detail_ring)
        pz = zeros(detail_length, detail_ring)
        for i in 1:detail_length
        for j in 1:detail_ring
            #P = repmat(from',m*n,1) .+ S.*vec_u .+ cos(T).*vec_v .+ sin(T).*vec_w;
            px[i,j] = from[1] + S[i]*vec_u[1] + cos(T[j])*vec_v[1] + sin(T[j])*vec_w[1];
            py[i,j] = from[2] + S[i]*vec_u[2] + cos(T[j])*vec_v[2] + sin(T[j])*vec_w[2];
            pz[i,j] = from[3] + S[i]*vec_u[3] + cos(T[j])*vec_v[3] + sin(T[j])*vec_w[3];
        end
        end
        surf(px,py,pz, rstride=1, cstride=1, linewidth=0, antialiased=true, color=color)
    end

    # plot all bonds
    for c in lattice.connections
        # skip if wrong directioin
        if c[1] < c[2] || (!visualize_periodic && sum([abs(el) for el in c[4]]) != 0)
            continue
        end
        # get data
        from = lattice.positions[Int(c[1])] .* conversion
        to = lattice.positions[Int(c[2])] .* conversion
        color = get(colorcode_bonds, string(c[3]), colorcode_bonds["0"]) ./ 255.0
        # plot
        plotBond(from, to, bond_thickness, color)
    end
    # plot all sites
    for index in 1:size(lattice.positions,1)
        p = lattice.positions[index]
        if length(p) == 2
            p = [p[1], p[2], 0.0]
        end
        p = p.*conversion
        plotSite(p, site_radius, get(colorcode_sites, lattice.positions_indices[index], colorcode_sites[0])./255.0)
    end
    # plot all site labels
    for index in 1:size(lattice.positions,1)
        p = lattice.positions[index]
        if length(p) == 2
            p = [p[1], p[2], 0.0]
        end
        p = p.*conversion
        # add label
        if site_labels == "LATTICE INDEX"
            ax[:text](p[1]+site_radius, p[2]+site_radius, p[3]+site_radius,  "$(index)", zorder=1, size=site_radius, color="k")
        elseif site_labels == "POSITION INDEX"
            ax[:text](p[1]+site_radius, p[2]+site_radius, p[3]+site_radius,  "$(lattice.positions_indices[index])", zorder=1, size=site_radius, color="k")
        end
    end

    # turn off everthing but the plot
    ax[:set_axis_off]()
    # get the current limits
    x_limits = ax[:get_xlim3d]()
    y_limits = ax[:get_ylim3d]()
    z_limits = ax[:get_zlim3d]()
    # build the current ranges
    x_range = abs(x_limits[2] - x_limits[1])
    x_middle = (x_limits[2] + x_limits[1]) / 2.0
    y_range = abs(y_limits[2] - y_limits[1])
    y_middle = (y_limits[2] + y_limits[1]) / 2.0
    z_range = abs(z_limits[2] - z_limits[1])
    z_middle = (z_limits[2] + z_limits[1]) / 2.0
    # The plot bounding box is a sphere in the sense of the infinity
    # norm, hence I call half the max range the plot radius.
    plot_radius = 0.5*maximum([x_range, y_range, z_range])
    ax[:set_xlim3d](x_middle - plot_radius, x_middle + plot_radius)
    ax[:set_ylim3d](y_middle - plot_radius, y_middle + plot_radius)
    ax[:set_zlim3d](z_middle - plot_radius, z_middle + plot_radius)
    # tighten layout
    tight_layout()
    # show
    show()
end
# show Lattice function to plot in mayavi
function showLatticeMayavi(
		lattice::Lattice;
		conversion = 160,
		site_radius=25,
        site_labels="OFF",
		bond_thickness::Int64=8,
		visualize_periodic=false,
		colorcode_sites = Dict(0 => [255,255,255], 1 => [255,255,255]),
		colorcode_bonds = Dict("0" => [0,0,0], "1.0" => [0,0,0]),
        colorcode_bonds_automation::String = "OFF",
        background_color = [200,200,200]
    )

    # CHECK IF MAYAVI IS PRESENT
    @pyimport mayavi.mlab as mlab

    # maybe overwrite dictonary
    if colorcode_bonds_automation == "GREY"
        # construct the list of interaction strengths
        cs_list = getConnectionStrengthList(lattice)
        # get the color code list
        cc_list = getGreySequence(size(cs_list, 1))
        # put in a new dictonary
        colorcode_bonds = Dict()
        # insert all pairs
        for i in 1:size(cc_list, 1)
            colorcode_bonds[string(cs_list[i])] = cc_list[i]
        end
    elseif colorcode_bonds_automation == "COLOR"
        # construct the list of interaction strengths
        cs_list = getConnectionStrengthList(lattice)
        # get the color code list
        cc_list = getColorSequence(size(cs_list, 1))
        # put in a new dictonary
        colorcode_bonds = Dict()
        # insert all pairs
        for i in 1:size(cc_list, 1)
            colorcode_bonds[string(cs_list[i])] = cc_list[i]
        end
    end

	# repair color dictonary
	colorcode_bonds["0"] = get(colorcode_bonds, "0", [0,0,0])
	colorcode_sites[0] = get(colorcode_sites, 0, [255,255,255])

    # PLOTTING WITH MAYAVI
    # create a figure
    mlab.figure(bgcolor=(background_color[1], background_color[2], background_color[3]))
    # define a function to plot a sphere for every site
    function plotSite(site,radius,color; detail=10)
        x = site[1]
        y = site[2]
        if length(site) == 3
            z = site[3]
        else
            z = 0.0
        end
        mlab.points3d([x],[y],[z],scale_factor=radius, color=(color[1], color[2], color[3]))
    end
    # define a function to plot a bond as a tube
    function plotBond(from, to, radius, color; detail_ring=8, detail_length=3)
        # check to be 3d
        if length(from) == 2
            from = [from[1], from[2], 0.0]
        end
        if length(to) == 2
            to = [to[1], to[2], 0.0]
        end
        mlab.plot3d([from[1], to[1]], [from[2], to[2]], [from[3], to[3]], color=(color[1], color[2], color[3]), tube_radius=radius)
    end

    # plot all bonds
    for c in lattice.connections
        # skip if wrong directioin
        if c[1] < c[2] || (!visualize_periodic && sum([abs(el) for el in c[4]]) != 0)
            continue
        end
        # get data
        from = lattice.positions[Int(c[1])] .* conversion
        to = lattice.positions[Int(c[2])] .* conversion
        color = get(colorcode_bonds, string(c[3]), colorcode_bonds["0"]) ./ 255.0
        # plot
        plotBond(from, to, bond_thickness, color)
    end
    # plot all sites
    for index in 1:size(lattice.positions,1)
        p = lattice.positions[index]
        if length(p) == 2
            p = [p[1], p[2], 0.0]
        end
        p = p.*conversion
        plotSite(p, site_radius, get(colorcode_sites, lattice.positions_indices[index], colorcode_sites[0])./255.0)
    end
    # plot all site labels
    for index in 1:size(lattice.positions,1)
        p = lattice.positions[index]
        if length(p) == 2
            p = [p[1], p[2], 0.0]
        end
        p = p.*conversion
        # add label
        if site_labels == "LATTICE INDEX"
            mlab.text3d(p[1]+site_radius, p[2]+site_radius, p[3]+site_radius, "$(index)", scale=site_radius, color=(0,0,0))
        elseif site_labels == "POSITION INDEX"
            mlab.text3d(p[1]+site_radius, p[2]+site_radius, p[3]+site_radius, "$(lattice.positions_indices[index])", scale=site_radius, color=(0,0,0))
        end
    end

    # show plot
    mlab.show()

    # return the figure
    return fig
end

#-----------------------------------------------------------------------------------------------------------------------------
#
#   METHODS FOR JOB SCHEDULING
#
#-----------------------------------------------------------------------------------------------------------------------------

#-----------------------------------------------------------------------------------------------------------------------------
#
#   FUNCTION TO CREATE A JOB DIRECTORY FOR SCHEDULING
#
#   USAGE: Give a jobdirectory name and a main program to run as well as a default input file and the function creates
#   a job directory containing all input files as well as two scripts to run the program in parallel or single core to
#   calculate all input files.
#   Usage should be in the following steps:
#   1) provide an input file template of the form: <parameter name>\t<default value>
#   2) provide a main program which can be called by a command similar to "PROGRAM_NAME INPUTFILE.txt"
#   3) call this function to create job directory and all input files
#   4) Enter job directory and run one of the two julia execution scripts
#       - Single core execution one after another: "julia script_single_execution.jl"
#       - MPI based multi-core execution: "mpirun -x PATH --host <host1, host2, host3,...> --map-by node -np <number of processes> nice -n19 julia script_MPI_execution.jl"
#   5) collect output files
#
#
#   Parameters (necessary):
#   - directory_path:       The name of the directory to be created
#   - input_file_template:  The input file template from which parameters are extracted
#   - main_program_file:    The file which contains the main program to be run
#
#   Parameters (optional):
#   - parameters:           A dictonary of parameter values to be used
#   - paramters_use_ARGS:   Bool to determine if the ARGS string is searched for parameters
#   - input_file_seperator: The seperator used for input file seperation
#   - main_program_execution_command:   The command to launch the program PROGRAM with a given input file INPUTFILE. Default: julia PROGRAM INPUTFILE
#   - files_to_copy:        A list of file that are copied to the destination folder as well
#
#-----------------------------------------------------------------------------------------------------------------------------
function createJobDirectory(
    directory_path::String,
    input_file_template_name::String,
    main_program_file::String
    ;
    parameters::Dict=Dict(),
    paramters_use_ARGS::Bool=true,
    input_file_seperator::String=":\t",
    main_program_execution_command::String="julia PROGRAM INPUTFILE",
    files_to_copy::Array{String,1}=String[]
        )

    # -----------------
    # START INFORMATION
    # -----------------

    # correct the directory input
    if !endswith(directory_path, "/")
        directory_path = "$(directory_path)/"
    end
    # Start: Print the directory path
    println("Job directory \"$(directory_path)\" will be created")
    # Print the input template filename
    println("Input template is located at \"$(input_file_template_name)\"\n")

    # ---------------------------------------------------
    # OBTAIN THE INPUT FILE TEMPLATE AND EXTRACT KEYWORDS
    # ---------------------------------------------------

    # Open the template file and read all lines
    input_file_template         = open(input_file_template_name)
    input_file_template_lines   = readlines(input_file_template)
    close(input_file_template)

    # Create a new list for all parameters
    parameters_RAW              = Array[]
    parameters_passed_values    = Array[]

    # push in all relevant lines
    for line in input_file_template_lines
        line = strip(line)
        if startswith(line, "#") || line==""
            continue
        end
        push!(parameters_RAW, split(line, "\t"))
    end

    # print all parameters that are found
    println("Found $(size(parameters_RAW, 1)) parameters in template")
    for p in parameters_RAW
        if length(p) == 1
            println("- \"$(p[1])\" without default value")
        elseif length(p) == 2
            println("- \"$(p[1])\" with default value $(p[2])")
        end
    end
    println("")

    # ---------------------------------------------------
    # OBTAIN THE CORRECT VALUES TO BE USED
    # ---------------------------------------------------

    # Iterate over all accessible parameters
    for p in parameters_RAW

        # list of values to use
        values_to_use   = []

        # OPTION 1, check in Command line arguments
        if paramters_use_ARGS
            # look through all command line arguments
            for arg in ARGS
                # check if it fulfils the structure
                if startswith(arg, "$(p[1])=")
                    # get the string representation
                    arg_value = replace(arg, "$(p[1])=", "")
                    # push to the list of values to use
                    push!(values_to_use, arg_value)
                    # break the loop
                    break
                end
            end
        end

        # OPTION 2, check in dictonary
        if length(values_to_use)==0 && p[1] in keys(parameters)
            # get the string representation
            dict_value = parameters[p[1]]
            # push to the list of values to use
            push!(values_to_use, dict_value)
        end

        # OPTION 3, check default values of template
        if length(values_to_use)==0 && length(p) >= 2
            # get the string representation
            default_value = p[2]
            # push to the list of values to use
            push!(values_to_use, default_value)
        end

        # PARSING OF VALUES

        # create a list with parsed entries
        values_to_use_parsed = []

        # iterate over all values to be used
        for v in values_to_use
            # determine which type of value
            if typeof(v) != String && typeof(v) != SubString{String}
                # already parsed, push into array
                push!(values_to_use_parsed, v)
            elseif startswith(v, "[") && endswith(v, "]")
                # list of format [el1, el2, el3, ...]
                v = split(strip(v[2:end-1]), ",")
                for i in 1:length(v)
                    v[i] = strip(v[i])
                    push!(values_to_use_parsed, parse(Float64,v[i]))
                end
            elseif startswith(v, "linspace(") && endswith(v, ")")
                # list via linspace(start, end, steps)
                v = split(v[10:end-1], ",")
                for i in 1:length(v)
                    v[i] = strip(v[i])
                end
                v = linspace(parse(Float64,v[1]),parse(Float64,v[2]),parse(Float64,v[3]))
                for i in 1:length(v)
                    push!(values_to_use_parsed, v[i])
                end
            else
                # it is probably either a number or a string, try parsing to float
                if isnull(tryparse(Float64,v))
                    # not a float, just push in list as it is
                    push!(values_to_use_parsed, v)
                else
                    # push parsed into list
                    push!(values_to_use_parsed, parse(Float64,v))
                end
            end
        end

        # append the list of paramters to the global list
        push!(parameters_passed_values, values_to_use_parsed)

    end

    # print all values that are found
    println("Found values for parameters:")
    for (p_index,p) in enumerate(parameters_RAW)
        println("- \"$(p[1])\", $(length(parameters_passed_values[p_index])) values: $(parameters_passed_values[p_index])")
    end
    println("")

    # -----------
    # CREATE JOBS
    # -----------

    # create a new list of dictonaries which are the jobs
    joblist = []
    # add a dummy job without information
    push!(joblist, Dict())

    # iterate over all parameters
    for p in 1:length(parameters_RAW)

        # find out if single value or multiple values
        if length(parameters_passed_values[p]) == 1
            # only one value, all jobs get this value
            for j in joblist
                j["$(parameters_RAW[p][1])"] = parameters_passed_values[p][1]
            end
        elseif length(parameters_passed_values[p]) >= 1
            # more than one value, every job gets copied for every value
            joblist_backup = joblist
            joblist = []
            for j in joblist_backup
                for i in 1:length(parameters_passed_values[p])
                    j_tmp = copy(j)
                    j_tmp["$(parameters_RAW[p][1])"] = parameters_passed_values[p][i]
                    push!(joblist, j_tmp)
                end
            end
        else
            # PROBLEM: NO VALUE FOR PARAMETER, exlude from job
            println("DANGER: Excluding parameter $(parameters_RAW[p][1]) from jobs because it has no value")
        end

    end

    # print the success
    println("Job building is done, $(length(joblist)) jobs are constructed\n")

    # check if only one job without parameters
    if length(joblist) == 1 && length(keys(joblist[1])) == 0
        println("Some mistake occured, only job has no parameters. Aborting now...\n")
        return
    end

    # --------------------
    # CREATE THE DIRECTORY
    # --------------------

    # check if the directory exists
    if isdir(directory_path)
        # Directory exists
        println("Directory \"$(directory_path)\" already exists")
        print("Delete and reconstruct? [y/n]: ")
        answer = strip(readline()[end-1:end])
        while !in(answer, ["y", "n"])
            print("Delete and reconstruct? [y/n]: ")
            answer = strip(readline()[end-1:end])
        end
        # depending on answer
        if answer == "y"
            # ask for recursive deletion
            print("Delete recursively? [y/n]: ")
            answer = strip(readline()[end-1:end])
            while !in(answer, ["y", "n"])
                print("Delete and reconstruct? [y/n]: ")
                answer = strip(readline()[end-1:end])
            end
            if answer == "y"
                # delete the complete directory
                rm(directory_path, recursive=true)
            else
                # delete the complete directory
                rm(directory_path, recursive=false)
            end
            # print
            println("Deleted already existing directory $(directory_path)")
            # create the directory
            mkdir(directory_path)
            println("Recreated directory $(directory_path)\n")

        else
            # abort the program
            println("Not deleting directory, aborting program now...")
            return
        end
    else
        # Directory does not exist, create it
        mkdir(directory_path)
        println("Created directory $(directory_path)\n")
    end

    # ----------------------
    # CREATE THE INPUT FILES
    # ----------------------

    # iterate over all jobs
    for i in 1:length(joblist)
        # dump job i to file input_i.txt
        jobfilename = "input_$(i).txt"
        f = open("$(directory_path)$(jobfilename)", "w")
        # write entries of the dictonary into the file
        for k in keys(joblist[i])
            write(f, "$(k)$(input_file_seperator)$(joblist[i][k])\n")
        end
        # close the file
        close(f)
        # set the correct input filename in the job
        joblist[i]["input_filename"] = jobfilename
    end
    # print the success
    println("Created all $(length(joblist)) jobfiles\n")

    # ------------------------------
    # COPY / CREATE THE SCRIPT FILES
    # ------------------------------

    # Copy main program file
    cp(main_program_file, "$(directory_path)$(main_program_file)")
    println("Copied main program file \"$(main_program_file)\"")

    # Create single execution julia script
    script_single_exe = open("$(directory_path)script_single_execution.jl", "w")
    # write all lines
    for i in 1:length(joblist)
        # compose a line
        line = "run(`$(main_program_execution_command)`)\n"
        # change the program part
        line = replace(line, "PROGRAM", main_program_file)
        # change the input part
        line = replace(line, "INPUTFILE", joblist[i]["input_filename"])
        # set the correct input filename in the job
        write(script_single_exe, line)
    end
    # close the script
    close(script_single_exe)
    println("Created single execution script \"script_single_execution.jl\"")

    # Create single execution julia script
    script_MPI_exe = open("$(directory_path)script_MPI_execution.jl", "w")

    # write all lines, starting with imports
    write(script_MPI_exe, "# use with:\n")
    write(script_MPI_exe, "# mpirun -x PATH --host <host1,host2,...> --map-by node -np <number of processes> nice -n19 julia script_MPI_execution.jl\n")
    write(script_MPI_exe, "using MPI\n")
    write(script_MPI_exe, "MPI.Init()\n")
    # get the number of the MPI process
    write(script_MPI_exe, "commworld = MPI.COMM_WORLD\n")
    write(script_MPI_exe, "commrank  = MPI.Comm_rank(commworld)\n")
    write(script_MPI_exe, "commsize  = MPI.Comm_size(commworld)\n")
    write(script_MPI_exe, "max_i     = $(length(joblist))\n")
    write(script_MPI_exe, "init_i    = commrank+1\n")
    write(script_MPI_exe, "for i in init_i:commsize:max_i\n")
    # compose main line
    line = "\trun(pipeline(`$(main_program_execution_command)`, stdout=\"bash_stdout_INPUTFILE\", stderr=\"bash_stderr_INPUTFILE\"))\n\tprintln(\"Done input file \\\"INPUTFILE\\\"\")\n"
    # change the program part
    line = replace(line, "PROGRAM", main_program_file)
    # change the input part
    line = replace(line, "INPUTFILE", "input_\$(i).txt")
    # set the correct input filename in the job
    write(script_MPI_exe, line)
    write(script_MPI_exe, "end\n")
    # Finalize MPI script
    write(script_MPI_exe, "MPI.Finalize()\n")
    # close the script
    close(script_MPI_exe)
    println("Created MPI execution script \"script_MPI_execution.jl\"")

    # Copy additional files
    for fn in files_to_copy
        cp(fn, "$(directory_path)$(fn)")
        println("Copied additional file \"$(fn)\"")
    end

end
export createJobDirectory


text:
using Distributed
using ThreadSafeDicts

using Test

Threads.nthreads() == 1 && println("Hmm, you only have one thread allocated for Julia.")

function testThreadSafeDicts()
    dict = ThreadSafeDict()
    dict = ThreadSafeDict{String,Int64}()
    @test isempty(dict)
    dict["ten"] = 10
    dict["twenty"] = 20
    delete!(dict, "ten")
    @test pop!(dict) == ("twenty" => 20)
    @test isempty(dict)
    dict = ThreadSafeDict{String,Int64}(["a" => 0, "b" => 1, "c" => 1, "d" => 2,
        "e" => 3, "f" => 5, "g" => 8, "h" => 13, "i" => 21, "j" => 34, "extra" => -1])

    @test dict["b"] == 1 && dict["e"] == 3
    @Threads.threads for k in collect(keys(dict))
        dict[k] *= 100
    end

    x = get!(dict, "another", 77)
    @test x == 77
    @test dict["e"] == 300
    @test get(dict, "e", nothing) == 300
    dict["extra"] = 0
    y = dict["extra"]
    @test y == 0
    @test haskey(dict, "a") == true

    empty!(dict)

    @sync begin
        @Threads.threads for i in 1:1000
            dict[string(i)] = i
        end
    end

    @test ((x, y) = iterate(dict)) != nothing
    @test iterate(dict, y) != nothing

    @test length(dict.d) == 1000
    empty!(dict)
    @sync begin
        @Threads.threads for i in 1:1000
            sleep(rand() / 100)
            dict["number"] = i
        end
    end
    @test dict["number"] >= 1000 / Threads.nthreads()
    println("testing printing: ", dict)
end

testThreadSafeDicts()

ext: jl


text:
# ### axis flip

#md # ![FILE_NAME.png](images/FILE_NAME.png)

plot(x->x, 0:0.01:2π,
    proj=:polar,
    xflip=true,
    yflip=true,
    ## flip=true # x,y-flip
)


text:
using Test, JSONTables, Tables, JSON3

@testset "JSONTables" begin

cjson = replace(replace("""{
"a": [1,2,3],
"b": [4.1, null, 6.3],
"c": ["7", "8", null]
}""", " "=>""), "\n"=>"")

rjson = replace(replace("""[
{"a": 1, "b": 4.1, "c": "7"},
{"a": 2, "b": null, "c": "8"},
{"a": 3, "b": 6.3, "c": null}
]""", " "=>""), "\n"=>"")

ctable = (a=[1,2,3], b=[4.1, missing, 6.3], c=["7", "8", missing])
rtable = Tables.rowtable(ctable)

cjtable = JSONTables.jsontable(cjson)
rjtable = JSONTables.jsontable(rjson)

@test Tables.istable(typeof(cjtable))
@test Tables.istable(typeof(rjtable))

@test Tables.columnaccess(typeof(cjtable))
@test !Tables.rowaccess(typeof(cjtable))
@test Tables.columns(cjtable) === cjtable
@test propertynames(cjtable) == (:a, :b, :c)
@test getproperty(cjtable, :a) == [1, 2, 3]

@test !Tables.columnaccess(typeof(rjtable))
@test Tables.rowaccess(typeof(rjtable))
@test Tables.rows(rjtable) === rjtable
@test Base.IteratorSize(typeof(rjtable)) == Base.HasLength()
@test length(rjtable) == 3
@test Base.IteratorEltype(typeof(rjtable)) == Base.HasEltype()
@test eltype(rjtable) == Any

@test isequal(Tables.columntable(cjtable), ctable)
@test isequal(Tables.columntable(rjtable), ctable)
@test isequal(Tables.rowtable(cjtable), rtable)
@test isequal(Tables.rowtable(rjtable), rtable)

@test JSONTables.objecttable(ctable) == cjson
@test JSONTables.objecttable(rtable) == cjson
@test JSONTables.arraytable(ctable) == rjson
@test JSONTables.arraytable(rtable) == rjson

io = IOBuffer()
JSONTables.objecttable(io, ctable)
@test String(take!(io)) == cjson

JSONTables.arraytable(io, rtable)
@test String(take!(io)) == rjson

# #7
text = """{
        "color_scheme": "Packages/Color Scheme - Default/Mariana.sublime-color-scheme",
        "dictionary": "Packages/Language - English/en_US.dic",
        "draw_white_space": "all",
        "font_face": "monospace 821",
        "font_size": "10",
        "theme": "Adaptive.sublime-theme"
}"""
@test_throws ArgumentError JSONTables.jsontable(text)

nonhomogenous = """
[
    {"a": 1, "b": 2, "c": 3},
    {"b": 4, "c": null, "d": 5},
    {"a": 6, "d": 7},
    {"a": 8, "b": 9, "c": 10, "d": null},
    {"d": 11, "c": 10, "b": 9, "a": 8}
]
"""

jt = JSONTables.jsontable(nonhomogenous)
@test Tables.schema(jt).names == (:a, :b, :c, :d)
ct = Tables.columntable(jt)
@test isequal(ct.a, [1, missing, 6, 8, 8])
@test isequal(ct.d, [missing, 5, 7, missing, 11])

end


text:
import Base: rand,
             ∈,
             isempty

export Line,
       an_element

"""
    Line{N<:Real, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}

Type that represents a line in 2D of the form ``a⋅x = b`` (i.e., a special case
of a `Hyperplane`).

### Fields

- `a` -- normal direction (non-zero)
- `b` -- constraint

### Examples

The line ``y = -x + 1``:

```jldoctest
julia> Line([1., 1.], 1.)
Line{Float64,Array{Float64,1}}([1.0, 1.0], 1.0)
```
"""
struct Line{N<:Real, VN<:AbstractVector{N}} <: AbstractPolyhedron{N}
    a::VN
    b::N

    # default constructor with length constraint
    function Line(a::VN, b::N) where {N<:Real, VN<:AbstractVector{N}}
        @assert length(a) == 2 "lines must be two-dimensional"
        @assert !iszero(a) "a line needs a non-zero normal vector"
        return new{N, VN}(a, b)
    end
end

isoperationtype(::Type{<:Line}) = false
isconvextype(::Type{<:Line}) = true

# constructor from a LinearConstraint
Line(c::LinearConstraint{N}) where {N<:Real} = Line(c.a, c.b)

"""
    Line(p::AbstractVector{N}, q::AbstractVector{N}) where {N<:Real}

Constructor a line give two points.

### Input

- `p` -- point in 2D
- `q` -- another point in 2D

### Output

The line which passes through `p` and `q`.

### Algorithm

Given two points ``p = (x₁, y₁)`` and ``q = (x₂, y₂)`` the line that passes through these
two points is

```math
ℓ:~~y - y₁ = \\dfrac{(y₂ - y₁)}{(x₂ - x₁)} ⋅ (x-x₁).
```
The particular case ``x₂ = x₁`` defines a line parallel to the ``y``-axis (vertical line).
"""
function Line(p::AbstractVector{N}, q::AbstractVector{N}) where {N<:Real}
    x₁, y₁ = p[1], p[2]
    x₂, y₂ = q[1], q[2]

    if x₁ == x₂  # line is vertical
        a = [one(N), zero(N)]
        b = x₁
        return LazySets.Line(a, b)
    end

    k = (y₁ - y₂)/(x₂ - x₁)
    a = [k, one(N)]
    b = y₁ + k * x₁
    return Line(a, b)
end

# --- polyhedron interface functions ---

"""
    constraints_list(L::Line{N}) where {N<:Real}

Return the list of constraints of a line.

### Input

- `L` -- line

### Output

A list containing two half-spaces.
"""
function constraints_list(L::Line{N}) where {N<:Real}
    return _constraints_list_hyperplane(L.a, L.b)
end

# --- LazySet interface functions ---

"""
    dim(L::Line)

Return the ambient dimension of a line.

### Input

- `L` -- line

### Output

The ambient dimension of the line, which is 2.
"""
function dim(L::Line)
    return 2
end

"""
    σ(d::AbstractVector{N}, L::Line{N}) where {N<:Real}

Return the support vector of a line in a given direction.

### Input

- `d` -- direction
- `L` -- line

### Output

The support vector in the given direction, which is defined the same way as for
the more general `Hyperplane`.
"""
function σ(d::AbstractVector{N}, L::Line{N}) where {N<:Real}
    return σ(d, Hyperplane(L.a, L.b))
end

"""
    isbounded(L::Line)

Determine whether a line is bounded.

### Input

- `L` -- line

### Output

`false`.
"""
function isbounded(::Line)
    return false
end

"""
    isuniversal(L::Line{N}, [witness]::Bool=false) where {N<:Real}

Check whether a line is universal.

### Input

- `P`       -- line
- `witness` -- (optional, default: `false`) compute a witness if activated

### Output

* If `witness` option is deactivated: `false`
* If `witness` option is activated: `(false, v)` where ``v ∉ P``

### Algorithm

Witness production falls back to `isuniversal(::Hyperplane)`.
"""
function isuniversal(L::Line{N}, witness::Bool=false) where {N<:Real}
    if witness
        return isuniversal(Hyperplane(L.a, L.b), true)
    else
        return false
    end
end

"""
    an_element(L::Line{N}) where {N<:Real}

Return some element of a line.

### Input

- `L` -- line

### Output

An element on the line.

### Algorithm

If the ``b`` value of the line is zero, the result is the origin.
Otherwise the result is some ``x = [x1, x2]`` such that ``a·[x1, x2] = b``.
We first find out in which dimension ``a`` is nonzero, say, dimension 1, and
then choose ``x1 = 1`` and accordingly ``x2 = \\frac{b - a1}{a2}``.
"""
function an_element(L::Line{N}) where {N<:Real}
    if L.b == zero(N)
        return zeros(N, 2)
    end
    i = L.a[1] == zero(N) ? 2 : 1
    x = Vector{N}(undef, 2)
    x[3-i] = one(N)
    x[i] = (L.b - L.a[3-i]) / L.a[i]
    return x
end

"""
    ∈(x::AbstractVector{N}, L::Line{N}) where {N<:Real}

Check whether a given point is contained in a line.

### Input

- `x` -- point/vector
- `L` -- line

### Output

`true` iff `x ∈ L`.

### Algorithm

The point ``x`` belongs to the line if and only if ``a⋅x = b`` holds.
"""
function ∈(x::AbstractVector{N}, L::Line{N}) where {N<:Real}
    @assert length(x) == dim(L)
    return dot(L.a, x) == L.b
end

"""
    rand(::Type{Line}; [N]::Type{<:Real}=Float64, [dim]::Int=2,
         [rng]::AbstractRNG=GLOBAL_RNG, [seed]::Union{Int, Nothing}=nothing)

Create a random line.

### Input

- `Line` -- type for dispatch
- `N`    -- (optional, default: `Float64`) numeric type
- `dim`  -- (optional, default: 2) dimension
- `rng`  -- (optional, default: `GLOBAL_RNG`) random number generator
- `seed` -- (optional, default: `nothing`) seed for reseeding

### Output

A random line.

### Algorithm

All numbers are normally distributed with mean 0 and standard deviation 1.
Additionally, the constraint `a` is nonzero.
"""
function rand(::Type{Line};
              N::Type{<:Real}=Float64,
              dim::Int=2,
              rng::AbstractRNG=GLOBAL_RNG,
              seed::Union{Int, Nothing}=nothing)
    @assert dim == 2 "cannot create a random Line of dimension $dim"
    rng = reseed(rng, seed)
    a = randn(rng, N, dim)
    while iszero(a)
        a = randn(rng, N, dim)
    end
    b = randn(rng, N)
    return Line(a, b)
end

"""
    isempty(L::Line)

Return if a line is empty or not.

### Input

- `L` -- line

### Output

`false`.
"""
function isempty(L::Line)
    return false
end

"""
    constrained_dimensions(L::Line{N}) where {N<:Real}

Return the indices in which a line is constrained.

### Input

- `L` -- line

### Output

A vector of ascending indices `i` such that the line is constrained in dimension
`i`.

### Examples

A line with constraint ``x1 = 0`` is constrained in dimension 1 only.
"""
function constrained_dimensions(L::Line{N}) where {N<:Real}
    return nonzero_indices(L.a)
end

function _linear_map_hrep_helper(M::AbstractMatrix{N}, P::Line{N},
                                 algo::AbstractLinearMapAlgorithm) where {N<:Real}
    constraints = _linear_map_hrep(M, P, algo)
    if length(constraints) == 2
        # assuming these constraints define a line
        c = first(constraints)
        return Line(c.a, c.b)
    elseif isempty(constraints)
        return Universe{N}(size(M, 1))
    else
        error("unexpected number of $(length(constraints)) constraints")
    end
end

"""
    translate(L::Line{N}, v::AbstractVector{N}; share::Bool=false
             ) where {N<:Real}

Translate (i.e., shift) a line by a given vector.

### Input

- `L`     -- line
- `v`     -- translation vector
- `share` -- (optional, default: `false`) flag for sharing unmodified parts of
             the original set representation

### Output

A translated line.

### Notes

The normal vector of the line (vector `a` in `a⋅x = b`) is shared with the
original line if `share == true`.

### Algorithm

A line ``a⋅x = b`` is transformed to the line ``a⋅x = b + a⋅v``.
In other words, we add the dot product ``a⋅v`` to ``b``.
"""
function translate(L::Line{N}, v::AbstractVector{N}; share::Bool=false
                  ) where {N<:Real}
    @assert length(v) == dim(L) "cannot translate a $(dim(L))-dimensional " *
                                "set by a $(length(v))-dimensional vector"
    a = share ? L.a : copy(L.a)
    b = L.b + dot(L.a, v)
    return Line(a, b)
end


text:
using StrTables, LaTeX_Entities

using Test

const def = LaTeX_Entities.default

# Test the functions lookupname, matches, longestmatches, completions
# Check that characters from all 3 tables (BMP, non-BMP, 2 character) are tested

@testset "LaTeX_Entities" begin

@testset "lookupname" begin
    @test lookupname(def, SubString("My name is Spock", 12)) == ""
    @test lookupname(def, "foobar")   == ""
    @test lookupname(def, "dagger")   == "†" # \u2020
    #@test lookupname(def, "mscrl")    == "𝓁" # \U1f4c1
    @test lookupname(def, "c_l")       == "𝓁" # \U1f4c1
    @test lookupname(def, "nleqslant") == "⩽̸" # \u2a7d\u338
end

@testset "matches" begin
    @test isempty(matches(def, ""))
    @test isempty(matches(def, "\U1f596"))
    @test isempty(matches(def, SubString("My name is \U1f596", 12)))
    for (chrs, exp) in (("√", ["sqrt", "surd"]),
                        #("𝓁", ["mscrl"]),
                        ("𝓁", ["c_l"]),
                        ("⩽̸", ["nleqslant"]))
        res = matches(def, chrs)
        @test length(res) >= length(exp)
        @test intersect(res, exp) == exp
    end
end

@testset "longestmatches" begin
    @test isempty(longestmatches(def, "\U1f596 abcd"))
    @test isempty(longestmatches(def, SubString("My name is \U1f596", 12)))
    for (chrs, exp) in (("√abcd", ["sqrt", "surd"]),
                        #("𝓁abcd", ["mscrl"]),
                        ("𝓁abcd", ["c_l"]),
                        ("⩽̸abcd", ["nleqslant"]))
        res = longestmatches(def, chrs)
        @test length(res) >= length(exp)
        @test intersect(res, exp) == exp
    end
end

@testset "completions" begin
    @test isempty(completions(def, "ScottPaulJones"))
    @test isempty(completions(def, SubString("My name is Scott", 12)))
    for (chrs, exp) in (("A", ["AA", "AE", "Alpha"]),
                        #("mtt", ["mtta", "mttthree", "mttzero"]),
                        ("varp", ["varperspcorrespond", "varphi", "varpi"]),
                        ("nleq", ["nleq", "nleqslant"]))
        res = completions(def, chrs)
        @test length(res) >= length(exp)
        @test intersect(res, exp) == exp
    end
end
end


text:
# https://juliagraphics.github.io/Gtk.jl/latest/manual/listtreeview/#List-Store-1

using Gtk

# List Store

# Lets start with a very simple example: A table with three columns
# representing the name, the age and the gender of a person. Each column must
# have a specific type. Here, we chose to represent the gender using a boolean
# value where true represents female and false represents male. We thus
# initialize the list store using

ls = GtkListStore(String, Int, Bool)

# Now we will the store with data

push!(ls,("Peter",20,false))
push!(ls,("Paul",30,false))
push!(ls,("Mary",25,true))

# If we want so insert the data at a specific position we can use the insert
# function

insert!(ls, 2, ("Susanne", 35, true))

# You can use ls like a matrix like container. Calling length and size will
# give you

@show length(ls)
@show size(ls)

# Specific element can be be accessed using

@show ls[1,1]
@show ls[1,1]

# List View

# Now we actually want to display our data. To this end we create a tree view
# object

tv = GtkTreeView(GtkTreeModel(ls))

# Then we need specific renderers for each of the columns. Usually you will
# only need a text renderer, but in our example we want to display the boolean
# value using a checkbox.

rTxt = GtkCellRendererText()
rTog = GtkCellRendererToggle()

# Finally we create for each column a TreeViewColumn object

c1 = GtkTreeViewColumn("Name", rTxt, Dict([("text",0)]))
c2 = GtkTreeViewColumn("Age", rTxt, Dict([("text",1)]))
c3 = GtkTreeViewColumn("Female", rTog, Dict([("active",2)]))

# We need to push these column description objects to the tree view

push!(tv, c1, c2, c3)

# Then we can display the tree view widget in a window

win = GtkWindow(tv, "List View")
showall(win)

# If you prefer that the columns are resizable by the user call

for c in [c1, c2, c3]
    GAccessor.resizable(c, true)
end

# Sorting

# We next want to make the tree view sortable

for (i,c) in enumerate([c1,c2,c3])
  GAccessor.sort_column_id(c,i-1)
end

# If you now click on one of the column headers, the data will be
# sorted with respect to the selected column. You can even make
# the columns reorderable

for (i,c) in enumerate([c1, c2, c3])
    #GAccessor.reorderable(c, i)
    # gtk_tree_view_column_set_reorderable
    GAccessor.reorderable(c, true)
end

# Selection
# Usually the interesting bit of a list will be the entry being selected. This
# is done using an additional GtkTreeSelection object that can be retrieved by

selection = GAccessor.selection(tv)

# One either have single selection or multiple selections. We toggle this by
# calling

#selection = GAccessor.mode(selection,Gtk.GConstants.GtkSelectionMode.MULTIPLE) # not supported yet
# selected_rows(selection) crashes
selection = GAccessor.mode(selection,Gtk.GConstants.GtkSelectionMode.SINGLE)

# We will stick with single selection for now and want to know the index of the selected item

@show ls[selected(selection),1]

#"Pete"
# Since it can happen that no item has been selected at all, it is a good idea
# to put this into an if statement

if hasselection(selection)
  # do something with selected(selection)
end
# Sometimes you want to invoke an action of an item is selected. This can be done by

signal_connect(selection, "changed") do widget
  if hasselection(selection)
    currentIt = selected(selection)

    # now you can to something with the selected item
    println("Name: ", ls[currentIt,1], " Age: ", ls[currentIt,1])
  end
end
# Another useful signal is "row-activated" that will be triggered by a double click of the user.

# Note: getting multiple selections still not implemented


text:
# Tue 23 Oct 2012 07:10:59 PM EDT
#
# NeHe Tut 2 - Draw a triangle and a square
#
# Q - quit

# load necessary GLUT/OpenGL routines

global OpenGLver="1.0"
using OpenGL
using GLUT

# intialize variables

global window

width  = 640
height = 480

# function to init OpenGL context

function initGL(w::Integer,h::Integer)
    glViewport(0,0,w,h)
    glClearColor(0.0, 0.0, 0.0, 0.0)
    glClearDepth(1.0)
    glDepthFunc(GL_LESS)
    glEnable(GL_DEPTH_TEST)
    glShadeModel(GL_SMOOTH)

    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    gluPerspective(45.0,w/h,0.1,100.0)

    glMatrixMode(GL_MODELVIEW)
end

# prepare Julia equivalents of C callbacks that are typically used in GLUT code

function ReSizeGLScene(w::Int32,h::Int32)
    if h == 0
        h = 1
    end

    glViewport(0,0,w,h)

    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    gluPerspective(45.0,w/h,0.1,100.0)

    glMatrixMode(GL_MODELVIEW)

    return nothing
end

_ReSizeGLScene = cfunction(ReSizeGLScene, Void, (Int32, Int32))

function DrawGLScene()
    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glLoadIdentity()

    glColor(1.0,1.0,1.0)

    glTranslate(-1.5,0.0,-6.0)

    glBegin(GL_POLYGON)
      glVertex(0.0,1.0,0.0)
      glVertex(1.0,-1.0,0.0)
      glVertex(-1.0,-1.0,0.0)
    glEnd()

    glTranslate(3.0,0.0,0.0)

    glBegin(GL_QUADS)
        glVertex(-1.0,1.0,0.0)
        glVertex(1.0,1.0,0.0)
        glVertex(1.0,-1.0,0.0)
        glVertex(-1.0,-1.0,0.0)
    glEnd()

    glutSwapBuffers()

    return nothing
end

_DrawGLScene = cfunction(DrawGLScene, Void, ())

function keyPressed(the_key::Char,x::Int32,y::Int32)
    if the_key == int('q')
        glutDestroyWindow(window)
    end

    return nothing # keyPressed returns "void" in C. this is a workaround for Julia's "automatically return the value of the last expression in a function" behavior.
end

_keyPressed = cfunction(keyPressed, Void, (Char, Int32, Int32))

# run GLUT routines

glutInit()
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH)
glutInitWindowSize(width, height)
glutInitWindowPosition(0, 0)

window = glutCreateWindow("NeHe Tut 2")

glutDisplayFunc(_DrawGLScene)
glutFullScreen()

glutIdleFunc(_DrawGLScene)
glutReshapeFunc(_ReSizeGLScene)
glutKeyboardFunc(_keyPressed)

initGL(width, height)

glutMainLoop()


text:
"""
    rdm!(rbm, data)

Trains the RBM on data using contrastive divergence with randomly initialized chains.
See http://arxiv.org/abs/2105.13889.
"""
function rdm!(rbm::RBM, data::AbstractArray;
    batchsize = 1,
    epochs = 1,
    optim = ADAM(),
    wts = nothing,
    steps::Int = 1,
)
    @assert size(data) == (size(visible(rbm))..., size(data)[end])
    @assert isnothing(wts) || _nobs(data) == _nobs(wts)

    stats = suffstats(rbm, data; wts)

    for epoch in 1:epochs
        batches = minibatches(data, wts; batchsize = batchsize)
        for (vd, wd) in batches
            # fantasy particles, initialized randomly
            vm = sample_from_inputs(visible(rbm), Falses(size(visible(rbm)))..., batchsize)
            vm = sample_v_from_v(rbm, vm; steps = steps)
            ∂ = ∂contrastive_divergence(rbm, vd, vm; wd = wd, wm = wd, stats)
            update!(rbm, update!(∂, rbm, optim))
        end
    end
    return rbm
end


text:
using ValueOrientedRiskManagementInsurance

using Distributions
using DataFrames

include("S2NonLife_Input.jl")

#################################################################
println("Start S2NonLife.jl")

lobs = Array(NLLob, 0)

for lob in select_lob_names
  push!(lobs, NLLob(df_lobs[df_lobs[:name] .== lob, :],
                    convert(Array, β_raw[lob]),
                    rf,
                    corr_prem_res, df_lobs_prem_risk))
end

## premium reserve risk #########################################
scr_prem_res = 3 * premrestotalsd(lobs, corr_lob)[1]
scr_lapse = 0.0

## catastrophe risk #############################################
## proportional reinsurance seems to be ignored for cat scr

## natural catastrophe risk -------------------------------------
scr_cat_nat = 0.0   ## no natural catastrophe risk in portfolio

## man-made catastrophe risk ------------------------------------
## only fire risk and liability risk in "man-made" portfolio

## fire
scr_cat_fire = sum(cat_fire_is) ## no risk management mechanisms

## liability
prem_liab_vec = Array(Float64, nrow(df_cat_liability))
scr_cat_liab_vec = Array(Float64, nrow(df_cat_liability))
cat_liab_grp = Array(Int, 0)
for lob in lobs
  if lob.name == :liability
    prem_liab_vec = lob.prem_gross_cy * cat_liability_mix
    scr_cat_liab_vec =
      convert(Array,
              prem_liab_vec .* df_cat_liability[:,:f_liab])
    for i = 1:length(prem_liab_vec)
      if prem_liab_vec[i] > 0
        push!(cat_liab_grp, i)
      end
    end
  end
end
scr_cat_liab =
  sqrt(scr_cat_liab_vec' * corr_liab * scr_cat_liab_vec)[1,1]

## total man-made catastrophe risk ------------------------------
scr_cat_man_made_vec = [scr_cat_fire, scr_cat_liab]
scr_cat_man_made =
  sqrt(scr_cat_man_made_vec ⋅ scr_cat_man_made_vec)

## other catastrophe risk ---------------------------------------
scr_cat_other = 0.0

## catastrophe risk for non-prop reinsurance --------------------
scr_cat_nprop = 0.0  ## XYZ does not act as reinsurer

## total catastrophe risk ---------------------------------------
scr_cat = sqrt((scr_cat_nat + scr_cat_nprop)^2 +
                 scr_cat_man_made^2 +
                 scr_cat_other^2)

## Total SCR ####################################################
scr_nl_vec = [scr_prem_res, scr_lapse, scr_cat]
scr_nl = sqrt(scr_nl_vec ⋅ (corr_scr * scr_nl_vec))

## main results #################################################
println("scr_prem_risk        :  $(round(scr_prem_res, 2))")
println("scr_lapse            :  $(round(scr_lapse, 2))")
println("scr_cat_man_made_fire:  $(round(scr_cat_fire, 2))")
println("scr_cat_man_made_liab:  $(round(scr_cat_liab, 2))")
println("scr_cat_man_made     :  $(round(scr_cat_man_made, 2))")
println("scr_cat_other        :  $(round(scr_cat_other, 2))")
println("scr_cat_nprob        :  $(round(scr_cat_nprop, 2))")
println("scr_cat              :  $(round(scr_cat, 2))")
println("scr_total            :  $(round(scr_nl, 2))")

println("End S2NonLife.jl")


text:
using Pkg; Pkg.instantiate(); Pkg.activate(".")
using CSV
using DataFrames
A = Array{Float64}(undef, 12, 5)
open("data/day03.txt") do f

    # line_number
  line = 1

  # read till end of file
  while ! eof(f)

     # read a new / next line for every iteration
     s = readline(f)

     for i in 1:5
        println("$s")
        A[line,i] = parse(Int, SubString(s, i, i))
     end

     line += 1
     println("$line . $s")
  end

end

show(A)

gamma_one = 0
gamma_zero = 0

for i in size(A)[2]
    if A[i,1] > 0
        gamma_one += 1
    else
        gamma_zero += 1
    end
end

gamma_one
gamma_zero


text:
module KROME

using Reexport: @reexport
using CBinding: 𝐣𝐥

include("libkrome.jl")

"""
    examples_dir()

Return the directory where the examples provided with KROME.jl are located. If KROME is
installed as a regular package (with `]add KROME`), these files are read-only and should *not* be
modified. To find out which files are available, use, e.g., `readdir`:

# Examples
```julia
julia> readdir(examples_dir())
2-element Array{String,1}:
 "av-slab-benchmark"
 "test_hello"
```
"""
examples_dir() = joinpath(pathof(KROME) |> dirname |> dirname, "examples")

end


text:
using Oceananigans.Operators: Δy_qᶠᶜᶜ, Δx_qᶜᶠᶜ, Δx_qᶠᶜᶜ, ℑxyᶠᶠᵃ, ℑxzᶠᵃᶠ, ℑyzᵃᶠᶠ

# Interface for "conditional fluxes" (see ImmersedBoundaries module)
@inline _viscous_flux_ux(args...) = viscous_flux_ux(args...)
@inline _viscous_flux_uy(args...) = viscous_flux_uy(args...)
@inline _viscous_flux_uz(args...) = viscous_flux_uz(args...)
@inline _viscous_flux_vx(args...) = viscous_flux_vx(args...)
@inline _viscous_flux_vy(args...) = viscous_flux_vy(args...)
@inline _viscous_flux_vz(args...) = viscous_flux_vz(args...)
@inline _viscous_flux_wx(args...) = viscous_flux_wx(args...)
@inline _viscous_flux_wy(args...) = viscous_flux_wy(args...)
@inline _viscous_flux_wz(args...) = viscous_flux_wz(args...)

@inline _diffusive_flux_x(args...) = diffusive_flux_x(args...)
@inline _diffusive_flux_y(args...) = diffusive_flux_y(args...)
@inline _diffusive_flux_z(args...) = diffusive_flux_z(args...)

#####
##### Viscous flux divergences
#####

@inline function ∂ⱼ_τ₁ⱼ(i, j, k, grid, closure::AbstractTurbulenceClosure, args...)
    disc = time_discretization(closure)
    return 1/Vᶠᶜᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Ax_qᶜᶜᶜ, _viscous_flux_ux, disc, closure, args...) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_qᶠᶠᶜ, _viscous_flux_uy, disc, closure, args...) +
                                    δzᵃᵃᶜ(i, j, k, grid, Az_qᶠᶜᶠ, _viscous_flux_uz, disc, closure, args...))
end

@inline function ∂ⱼ_τ₂ⱼ(i, j, k, grid, closure::AbstractTurbulenceClosure, args...)
    disc = time_discretization(closure)
    return 1/Vᶜᶠᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶠᶜ, _viscous_flux_vx, disc, closure, args...) +
                                    δyᵃᶠᵃ(i, j, k, grid, Ay_qᶜᶜᶜ, _viscous_flux_vy, disc, closure, args...) +
                                    δzᵃᵃᶜ(i, j, k, grid, Az_qᶜᶠᶠ, _viscous_flux_vz, disc, closure, args...))
end

@inline function ∂ⱼ_τ₃ⱼ(i, j, k, grid, closure::AbstractTurbulenceClosure, args...)
    disc = time_discretization(closure)
    return 1/Vᶜᶜᶠ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶜᶠ, _viscous_flux_wx, disc, closure, args...) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_qᶜᶠᶠ, _viscous_flux_wy, disc, closure, args...) +
                                    δzᵃᵃᶠ(i, j, k, grid, Az_qᶜᶜᶜ, _viscous_flux_wz, disc, closure, args...))
end

@inline function ∇_dot_qᶜ(i, j, k, grid, closure::AbstractTurbulenceClosure, tracer_index, args...)
    disc = time_discretization(closure)
    return 1/Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶜᶜ, _diffusive_flux_x, disc, closure, tracer_index, args...) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_qᶜᶠᶜ, _diffusive_flux_y, disc, closure, tracer_index, args...) +
                                    δzᵃᵃᶜ(i, j, k, grid, Az_qᶜᶜᶠ, _diffusive_flux_z, disc, closure, tracer_index, args...))
end

#####
##### Products of viscosity and stress, divergence, vorticity
#####

@inline ν_σᶜᶜᶜ(i, j, k, grid, closure, K, clock, σᶜᶜᶜ, args...) = νᶜᶜᶜ(i, j, k, grid, closure, K, clock) * σᶜᶜᶜ(i, j, k, grid, args...)
@inline ν_σᶠᶠᶜ(i, j, k, grid, closure, K, clock, σᶠᶠᶜ, args...) = νᶠᶠᶜ(i, j, k, grid, closure, K, clock) * σᶠᶠᶜ(i, j, k, grid, args...)
@inline ν_σᶠᶜᶠ(i, j, k, grid, closure, K, clock, σᶠᶜᶠ, args...) = νᶠᶜᶠ(i, j, k, grid, closure, K, clock) * σᶠᶜᶠ(i, j, k, grid, args...)
@inline ν_σᶜᶠᶠ(i, j, k, grid, closure, K, clock, σᶜᶠᶠ, args...) = νᶜᶠᶠ(i, j, k, grid, closure, K, clock) * σᶜᶠᶠ(i, j, k, grid, args...)

@inline ν_δᶜᶜᶜ(i, j, k, grid, closure, K, clock, u, v) = νᶜᶜᶜ(i, j, k, grid, closure, K, clock) * div_xyᶜᶜᶜ(i, j, k, grid, u, v)
@inline ν_ζᶠᶠᶜ(i, j, k, grid, closure, K, clock, u, v) = νᶠᶠᶜ(i, j, k, grid, closure, K, clock) * ζ₃ᶠᶠᶜ(i, j, k, grid, u, v)

@inline κ_σᶠᶜᶜ(i, j, k, grid, closure, K, id, clock, σᶠᶜᶜ, args...) = κᶠᶜᶜ(i, j, k, grid, closure, K, id, clock) * σᶠᶜᶜ(i, j, k, grid, args...)
@inline κ_σᶜᶠᶜ(i, j, k, grid, closure, K, id, clock, σᶜᶠᶜ, args...) = κᶜᶠᶜ(i, j, k, grid, closure, K, id, clock) * σᶜᶠᶜ(i, j, k, grid, args...)
@inline κ_σᶜᶜᶠ(i, j, k, grid, closure, K, id, clock, σᶜᶜᶠ, args...) = κᶜᶜᶠ(i, j, k, grid, closure, K, id, clock) * σᶜᶜᶠ(i, j, k, grid, args...)

#####
##### Viscosity "extractors"
#####

# Number
@inline νᶜᶜᶜ(i, j, k, grid, clock, loc, ν::Number) = ν
@inline νᶠᶜᶠ(i, j, k, grid, clock, loc, ν::Number) = ν
@inline νᶜᶠᶠ(i, j, k, grid, clock, loc, ν::Number) = ν
@inline νᶠᶠᶜ(i, j, k, grid, clock, loc, ν::Number) = ν

@inline κᶠᶜᶜ(i, j, k, grid, clock, loc, κ::Number) = κ
@inline κᶜᶠᶜ(i, j, k, grid, clock, loc, κ::Number) = κ
@inline κᶜᶜᶠ(i, j, k, grid, clock, loc, κ::Number) = κ

# Array / Field at `Center, Center, Center`
const Lᶜᶜᶜ = Tuple{Center, Center, Center}
@inline νᶜᶜᶜ(i, j, k, grid, clock, ::Lᶜᶜᶜ, ν::AbstractArray) = @inbounds ν[i, j, k]
@inline νᶠᶜᶠ(i, j, k, grid, clock, ::Lᶜᶜᶜ, ν::AbstractArray) = ℑxzᶠᵃᶠ(i, j, k, grid, ν)
@inline νᶜᶠᶠ(i, j, k, grid, clock, ::Lᶜᶜᶜ, ν::AbstractArray) = ℑyzᵃᶠᶠ(i, j, k, grid, ν)
@inline νᶠᶠᶜ(i, j, k, grid, clock, ::Lᶜᶜᶜ, ν::AbstractArray) = ℑxyᶠᶠᵃ(i, j, k, grid, ν)

@inline κᶠᶜᶜ(i, j, k, grid, clock, ::Lᶜᶜᶜ, κ::AbstractArray) = ℑxᶠᵃᵃ(i, j, k, grid, κ)
@inline κᶜᶠᶜ(i, j, k, grid, clock, ::Lᶜᶜᶜ, κ::AbstractArray) = ℑyᵃᶠᵃ(i, j, k, grid, κ)
@inline κᶜᶜᶠ(i, j, k, grid, clock, ::Lᶜᶜᶜ, κ::AbstractArray) = ℑzᵃᵃᶠ(i, j, k, grid, κ)

# Array / Field at `Center, Center, Face`
const Lᶜᶜᶠ = Tuple{Center, Center, Face}
@inline νᶜᶜᶜ(i, j, k, grid, clock, ::Lᶜᶜᶠ, ν::AbstractArray) = ℑzᵃᵃᶜ(i, j, k, grid, ν)
@inline νᶠᶜᶠ(i, j, k, grid, clock, ::Lᶜᶜᶠ, ν::AbstractArray) = ℑxᶠᵃᵃ(i, j, k, grid, ν)
@inline νᶜᶠᶠ(i, j, k, grid, clock, ::Lᶜᶜᶠ, ν::AbstractArray) = ℑyᵃᶠᵃ(i, j, k, grid, ν)
@inline νᶠᶠᶜ(i, j, k, grid, clock, ::Lᶜᶜᶠ, ν::AbstractArray) = ℑxyzᶠᶠᶜ(i, j, k, grid, ν)

@inline κᶠᶜᶜ(i, j, k, grid, clock, ::Lᶜᶜᶠ, κ::AbstractArray) = ℑxzᶠᵃᶠ(i, j, k, grid, κ)
@inline κᶜᶠᶜ(i, j, k, grid, clock, ::Lᶜᶜᶠ, κ::AbstractArray) = ℑyzᵃᶠᶠ(i, j, k, grid, κ)
@inline κᶜᶜᶠ(i, j, k, grid, clock, ::Lᶜᶜᶠ, κ::AbstractArray) = @inbounds κ[i, j, k]

# Function

const c = Center()
const f = Face()

@inline νᶜᶜᶜ(i, j, k, grid, clock, loc, ν::F) where F<:Function = ν(node(c, c, c, i, j, k, grid)..., clock.time)
@inline νᶠᶜᶠ(i, j, k, grid, clock, loc, ν::F) where F<:Function = ν(node(f, c, f, i, j, k, grid)..., clock.time)
@inline νᶜᶠᶠ(i, j, k, grid, clock, loc, ν::F) where F<:Function = ν(node(c, f, f, i, j, k, grid)..., clock.time)
@inline νᶠᶠᶜ(i, j, k, grid, clock, loc, ν::F) where F<:Function = ν(node(f, f, c, i, j, k, grid)..., clock.time)

@inline κᶠᶜᶜ(i, j, k, grid, clock, loc, κ::F) where F<:Function = κ(node(f, c, c, i, j, k, grid)..., clock.time)
@inline κᶜᶠᶜ(i, j, k, grid, clock, loc, κ::F) where F<:Function = κ(node(c, f, c, i, j, k, grid)..., clock.time)
@inline κᶜᶜᶠ(i, j, k, grid, clock, loc, κ::F) where F<:Function = κ(node(c, c, f, i, j, k, grid)..., clock.time)

# "DiscreteDiffusionFunction"
@inline νᶜᶜᶜ(i, j, k, grid, clock, loc, ν::DiscreteDiffusionFunction) = ν.func(i, j, k, grid, c, c, c)
@inline νᶠᶜᶠ(i, j, k, grid, clock, loc, ν::DiscreteDiffusionFunction) = ν.func(i, j, k, grid, f, c, f)
@inline νᶜᶠᶠ(i, j, k, grid, clock, loc, ν::DiscreteDiffusionFunction) = ν.func(i, j, k, grid, c, f, f)
@inline νᶠᶠᶜ(i, j, k, grid, clock, loc, ν::DiscreteDiffusionFunction) = ν.func(i, j, k, grid, f, f, c)

@inline κᶠᶜᶜ(i, j, k, grid, clock, loc, κ::DiscreteDiffusionFunction) = κ.func(i, j, k, grid, f, c, c)
@inline κᶜᶠᶜ(i, j, k, grid, clock, loc, κ::DiscreteDiffusionFunction) = κ.func(i, j, k, grid, c, f, c)
@inline κᶜᶜᶠ(i, j, k, grid, clock, loc, κ::DiscreteDiffusionFunction) = κ.func(i, j, k, grid, c, c, f)


text:
using FermiCG
using Printf
using Test
using LinearAlgebra
using Profile
using HDF5
using Random
using PyCall
using Arpack
using StaticArrays

#@testset "tpsci" begin
    atoms = []

    r = 1
    a = 1
    push!(atoms,Atom(1,"H", [0, 0*a, 0*r]))
    push!(atoms,Atom(2,"H", [0, 0*a, 1*r]))
    push!(atoms,Atom(3,"H", [0, 1*a, 2*r]))
    push!(atoms,Atom(4,"H", [0, 1*a, 3*r]))
    push!(atoms,Atom(5,"H", [0, 2*a, 4*r]))
    push!(atoms,Atom(6,"H", [0, 2*a, 5*r]))
    push!(atoms,Atom(7,"H", [0, 3*a, 6*r]))
    push!(atoms,Atom(8,"H", [0, 3*a, 7*r]))
    push!(atoms,Atom(9,"H", [0, 4*a, 8*r]))
    push!(atoms,Atom(10,"H",[0, 4*a, 9*r]))
    push!(atoms,Atom(11,"H",[0, 5*a, 10*r]))
    push!(atoms,Atom(12,"H",[0, 5*a, 11*r]))

    clusters    = [(1:2),(3:4),(5:6),(7:8),(9:12)]
    init_fspace = [(1,1),(1,1),(1,1),(1,1),(2,2)]
    clusters    = [(1:4),(5:8),(9:12)]
    init_fspace = [(2,2),(2,2),(2,2)]
    clusters    = [(1:4),(5:8),(9:10),(11:12)]
    init_fspace = [(2,2),(2,2),(1,1),(1,1)]
    na = 6
    nb = 6

    basis = "sto-3g"
    mol     = Molecule(0,1,atoms,basis)

    nroots = 4

    # get integrals
    mf = FermiCG.pyscf_do_scf(mol)
    nbas = size(mf.mo_coeff)[1]
    ints = FermiCG.pyscf_build_ints(mol,mf.mo_coeff, zeros(nbas,nbas));

    # localize orbitals
    C = mf.mo_coeff
    Cl = FermiCG.localize(mf.mo_coeff,"lowdin",mf)
    FermiCG.pyscf_write_molden(mol,Cl,filename="lowdin.molden")
    S = FermiCG.get_ovlp(mf)
    U =  C' * S * Cl
    println(" Rotate Integrals")
    flush(stdout)
    ints = FermiCG.orbital_rotation(ints,U)
    println(" done.")
    flush(stdout)

    #
    # define clusters
    clusters = [Cluster(i,collect(clusters[i])) for i = 1:length(clusters)]
    display(clusters)

    rdm1 = zeros(size(ints.h1))

    e_cmf, U, Da, Db  = FermiCG.cmf_oo(ints, clusters, init_fspace, rdm1,
                                       max_iter_oo=40, verbose=0, gconv=1e-6, method="bfgs")
    FermiCG.pyscf_write_molden(mol,Cl*U,filename="cmf.molden")
    ints = FermiCG.orbital_rotation(ints,U)

    e_ref = e_cmf - ints.h0

    max_roots = 100
    # build Hamiltonian, cluster_basis and cluster ops
    #display(Da)
    #cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=2, max_roots=max_roots)
    cluster_bases = FermiCG.compute_cluster_eigenbasis(ints, clusters, verbose=0, max_roots=max_roots,
                                                       init_fspace=init_fspace, rdm1a=Da, rdm1b=Db)
    clustered_ham = FermiCG.extract_ClusteredTerms(ints, clusters)
    cluster_ops = FermiCG.compute_cluster_ops(cluster_bases, ints);

    p_spaces = Vector{FermiCG.ClusterSubspace}()
    q_spaces = Vector{FermiCG.ClusterSubspace}()

    # define p spaces
    for ci in clusters
        tss = FermiCG.ClusterSubspace(ci)
        tss[init_fspace[ci.idx]] = 1:1
        push!(p_spaces, tss)
    end

    # define q spaces
    for tssp in p_spaces
        tss = FermiCG.get_ortho_compliment(tssp, cluster_bases[tssp.cluster.idx])
        push!(q_spaces, tss)
    end

    println(" ================= Cluster P Spaces ===================")
    display.(p_spaces)
    println(" ================= Cluster Q Spaces ===================")
    display.(q_spaces)

    FermiCG.add_cmf_operators!(cluster_ops, cluster_bases, ints, Da, Db, verbose=0);

    ci_vector = FermiCG.ClusteredState(clusters, R=nroots)

    ref_fock = FermiCG.FockConfig(init_fspace)
    FermiCG.add_fockconfig!(ci_vector, ref_fock)

    #1 e hops
    if false
        focks1e = []
        for ci in clusters
            for cj in clusters
                ci != cj || continue
                tmp = [ref_fock.config...]
                tmp[ci.idx] = (tmp[ci.idx][1] + 1, tmp[ci.idx][2])
                tmp[cj.idx] = (tmp[cj.idx][1] - 1, tmp[cj.idx][2])
                push!(focks1e, FermiCG.FockConfig(tmp))
                tmp = [ref_fock.config...]
                tmp[ci.idx] = (tmp[ci.idx][1], tmp[ci.idx][2] + 1)
                tmp[cj.idx] = (tmp[cj.idx][1], tmp[cj.idx][2] - 1)
                push!(focks1e, FermiCG.FockConfig(tmp))
            end
        end

        for fock in focks1e
            FermiCG.add_fockconfig!(ci_vector, fock)
        end
    end
    # 1excitons
    if true
        focks1e = []

        for ci in clusters
            config = [1 for i in 1:length(clusters)]
            for i in 1:size(cluster_bases[ci.idx][ref_fock[ci.idx]],2)
            #for i in 1:2

                config[ci.idx] = i
                ci_vector[ref_fock][FermiCG.ClusterConfig(config)] = MVector{nroots}(zeros(nroots))
            end
        end
    end

    e0, e2, v0, v1 = FermiCG.tpsci_ci(ci_vector, cluster_ops, clustered_ham,
                                        thresh_cipsi=1e-3, thresh_foi=1e-6, thresh_asci=1e-2);

    ref = [-18.32983268907018,
           -18.05273493473686,
           -18.02737336416847,
           -17.995438987323247]

    @test isapprox(abs.(ref), abs.(e0), atol=1e-6)

#end

max_stars_repo_name: nmayhall-vt/FermiCG

max_issues_repo_name: nmayhall-vt/FermiCG

max_forks_repo_name: nmayhall-vt/FermiCG


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ d7cb8fd0-8ff7-11eb-3acf-c960087ff109
begin
	using PlutoUI
	using NativeSVG
end

# ╔═╡ 0bb21f1a-8ff8-11eb-3385-718fb337b749
md"""# Multiple Dispatch

In Julia you have the ability to write code that can operate on different types. This is called “generic programming.”
In this chapter I will discuss the use of type declarations in Julia, and I will introduce methods that offer ways to implement different behavior for a function depending on the types of their arguments. This is called “multiple dispatch.”
"""

# ╔═╡ 1f325d60-8ff8-11eb-1ebe-1f0a41b7149a
md"""## Type Declarations

The `::` operator attaches *type annotations* to expressions and variables, indicating what types they should have:

```julia
julia> (1 + 2) :: Float64
ERROR: TypeError: in typeassert, expected Float64, got Int64
julia> (1 + 2) :: Int64
3
```

This helps to confirm that your program works the way you expect.

The `::` operator can also be appended to the lefthand side of an assignment, or included as part of a declaration:

```julia
julia> function returnfloat()
	       x::Float64 = 100
	       x
       end
returnfloat (generic function with 1 method)
julia> x = returnfloat()
100.0
julia> typeof(x)
Float64
```

The variable `x` is always of type `Float64` and the value is converted to a floating point if needed.

A type annotation can also be attached to the header of a function definition:

```julia
function sinc(x)::Float64
	if x === 0
		return 1
	end
	sin(x)/(x)
end
```

The return value of `sinc` is always converted to type `Float64`.
The default behavior in Julia when types are omitted is to allow values to be of any
type (`Any`).
"""

# ╔═╡ 82ea73c4-8ff8-11eb-2e25-bde15dc079ad
md"""# Methods

In “Time”, we defined a struct named `MyTime` and you wrote a function named `printtime`:

```julia
using Printf

struct MyTime
	hour :: Int64
	minute :: Int64
	second :: Int64
end

function printtime(time)
	@printf("%02d:%02d:%02d", time.hour, time.minute, time.second)
end
```

As you can see, type declarations can (and should, for performance reasons) be added to the fields in a struct definition.

To call this function, we have to pass a `MyTime` object as an argument:

```julia
julia> start = MyTime(9, 45, 0)
MyTime(9, 45, 0)
julia> printtime(start)
09:45:00
```

To add a *method* to the function `printtime` that only accepts a `MyTime` object as an argument, all we have to do is append `::` followed by `MyTime` to the argument time in the function definition:

```julia
function printtime(time::MyTime)
	@printf("%02d:%02d:%02d", time.hour, time.minute, time.second)
end
```

A method is a function definition with a specific signature: `printtime` has one argument of type `MyTime`.

Calling the function `printtime` with a `MyTime` object yields the same result as before:

```julia
julia> printtime(start)
09:45:00
```

We can now redefine the first method without the :: type annotation, allowing an argument of any type:

```julia
function printtime(time)
	println("I don't know how to print the argument time.")
end
```

If you call the function `printtime` with an object that isn’t a `MyTime` object, you now get:

```julia
julia> printtime(150)
I don't know how to print the argument time.
```
"""

# ╔═╡ f34a98f4-8ff8-11eb-3cd3-f170eae0c30d
md"""#### Exercise 17-1

Rewrite `timetoint` and `inttotime` (from “Prototyping Versus Planning”) to specify their arguments.
"""

# ╔═╡ 0aa3ccbe-8ff9-11eb-3ac1-e343e9bb813b
md"""## Additional Examples

Here’s a version of `increment` (from “Modifiers”) rewritten to specify its arguments:

```julia
function increment(time::MyTime, seconds::Int64)
	seconds += timetoint(time)
	inttotime(seconds)
end
```

Note that now it is a pure function, not a modifier. Here’s how you would invoke `increment`:

```julia
julia> start = MyTime(9, 45, 0)
MyTime(9, 45, 0)
julia> increment(start, 1337)
MyTime(10, 7, 17)
```

If you put the arguments in the wrong order, you get an error:

```julia
julia> increment(1337, start)
ERROR: MethodError: no method matching increment(::Int64, ::MyTime)
```

The signature of the method is `increment(time::MyTime, seconds::Int64)`, not `increment(seconds::Int64, time::MyTime)`.

Rewriting `isafter` to act only on `MyTime` objects is as easy:

```julia
function isafter(t1::MyTime, t2::MyTime)
	(t1.hour, t1.minute, t1.second) > (t2.hour, t2.minute, t2.second)
end
```

By the way, optional arguments are implemented as syntax for multiple method definitions. For example, this definition:

```julia
function f(a=1, b=2)
	a + 2b
end
```

translates to the following three methods:

```julia
f(a,b) = a + 2b
f(a) = f(a, 2)
f() = f(1, 2)
```

These expressions are valid Julia method definitions. This is shorthand notation for defining functions/methods.
"""

# ╔═╡ a2e88438-8ff9-11eb-02b8-3f084a3f441d
md"""## Constructors

A constructor is a special function that is called to create an object. The *default constructor* methods of `MyTime`, which take all fields as parameters, have the following signatures:

```julia
MyTime(hour, minute, second)
MyTime(hour::Int64, minute::Int64, second::Int64)
```

We can also add our own *outer constructor* methods:

```julia
function MyTime(time::MyTime)
	MyTime(time.hour, time.minute, time.second)
end
```

This method is called a *copy constructor* because the new `MyTim` object is a copy of its argument.

To enforce invariants, we need *inner constructor* methods:

```julia
struct MyTime
	hour :: Int64
	minute :: Int64
	second :: Int64
	function MyTime(hour::Int64=0, minute::Int64=0, second::Int64=0)
		@assert(0 ≤ minute < 60, "Minute is not between 0 and 60.")
		@assert(0 ≤ second < 60, "Second is not between 0 and 60.")
		new(hour, minute, second)
	end
end
```

The struct `MyTime` now has four inner constructor methods:

```julia
MyTime()
MyTime(hour::Int64)
MyTime(hour::Int64, minute::Int64)
MyTime(hour::Int64, minute::Int64, second::Int64)
```

An inner constructor method is always defined inside the block of a type declaration, and it has access to a special function called new that creates objects of the newly declared type.

!!! danger
    The default constructor is not available if any inner constructor is defined. You have to write explicitly all the inner constructors you need.

A second method without arguments of the local function new exists:

```julia
struct MyTime
	hour :: Int
	minute :: Int
	second :: Int
	function MyTime(hour::Int64=0, minute::Int64=0, second::Int64=0)
		@assert(0 ≤ minute < 60, "Minute is between 0 and 60.")
		@assert(0 ≤ second < 60, "Second is between 0 and 60.")
		time = new()
		time.hour = hour
		time.minute = minute
		time.second = second
		time
	end
end
```

This allows us to construct recursive data structures—i.e., structs where one of the fields is the struct itself. In this case the struct has to be mutable because its fields are modified after instantiation.
"""

# ╔═╡ 286fca94-8ffa-11eb-1b9b-4d36482bd021
md"""## `show`

`show` is a special function that returns a string representation of an object. For exam‐ ple, here is a `show` method for `MyTime` objects:

```julia
using Printf

function Base.show(io::IO, time::MyTime)
	@printf(io, "%02d:%02d:%02d", time.hour, time.minute, time.second)
end
```

The prefix `Base` is needed because we want to add a new method to the `Base.show` function.

When you print an object, Julia invokes the `show` function:

```julia
julia> time = MyTime(9, 45)
09:45:00
```

When I write a new composite type, I almost always start by writing an outer constructor, which makes it easier to instantiate objects, and a show method, which is useful for debugging.
"""

# ╔═╡ 6021d5b0-8ffa-11eb-053a-85008d0a7f1a
md"""#### Exercise 17-2

Write an outer constructor method for the `Point` class that takes `x` and `y` as optional parameters and assigns them to the corresponding fields.
"""

# ╔═╡ 757c2abc-8ffa-11eb-1b8c-edb1f8938d82
md"""## Operator Overloading

By defining operator methods, you can specify the behavior of operators on programmer-defined types. For example, if you define a method named `+` with two `MyTime` arguments, you can use the `+` operator on `MyTime` objects.

Here is what the definition might look like:

```julia
import Base.+

function +(t1::MyTime, t2::MyTime)
	seconds = timetoint(t1) + timetoint(t2)
	inttotime(seconds)
end
```

The import statement adds the `+` operator to the local scope so that methods can be added.

And here is how you could use it:

```julia
julia> start = MyTime(9, 45)
09:45:00
julia> duration = MyTime(1, 35, 0)
01:35:00
julia> start + duration
11:20:00
```

When you apply the `+` operator to `MyTime` objects, Julia invokes the newly added method. When the REPL shows the result, Julia invokes show. So, there is a lot happening behind the scenes!

Adding to the behavior of an operator so that it works with programmer-defined types is called *operator overloading*.
"""

# ╔═╡ d92ea206-8ffa-11eb-1b9a-5ba6198ccb5f
md"""## Multiple Dispatch

In the previous section we added two `MyTime` objects, but you also might want to add an integer to a `MyTime` object:

```julia
function +(time::MyTime, seconds::Int64)
	increment(time, seconds)
end
```

Here is an example that uses the `+` operator with a `MyTime` object and an integer:

```julia
julia> start = MyTime(9, 45)
09:45:00
julia> start + 1337
10:07:17
```

Addition is a commutative operator, so we have to add another method:

```julia
function +(seconds::Int64, time::MyTime)
	time + seconds
end
```

And we get the same result:

```julia
julia> 1337 + start
10:07:17
```

The *dispatch* mechanism determines which method to execute when a function is called. Julia allows the dispatch process to choose which of a function’s methods to call based on the number of arguments given, and on the types of all of the function’s arguments. Using all of a function’s arguments to choose which method should be invoked is known as *multiple dispatch*.
"""

# ╔═╡ 20fdca9e-8ffb-11eb-3066-59f9a3689bf8
md"""#### Exercise 17-3

Write `+` methods for Point objects:

* If both operands are `Point` objects, the method should return a new `Point` object whose `x` coordinate is the sum of the `x` coordinates of the operands, and likewise for the `y` coordinates.
* If the first or the second operand is a tuple, the method should add the first ele‐ ment of the tuple to the `x` coordinate and the second element to the `y` coordinate, and return a new `Point` object with the result.
"""

# ╔═╡ 5c319622-8ffb-11eb-1f0f-3d222d4b7408
md"""## Generic Programming

Multiple dispatch is useful when it is necessary, but (fortunately) it is not always necessary. Often you can avoid it by writing functions that work correctly for arguments with different types. This is known as *generic programming*.

Many of the functions we wrote for strings also work for other sequence types. For example, in “Dictionaries as Collections of Counters” we used `histogram` to count the number of times each letter appears in a word:

```julia
function histogram(s)
	d = Dict()
	for c in s
		if c ∉ keys(d)
			d[c] = 1
		else
			d[c] += 1
		end
	end
	d
end
```

This function also works for lists, tuples, and even dictionaries, as long as the elements of `s` are hashable so they can be used as keys in `d`:

```julia
julia> t = ("spam", "egg", "spam", "spam", "bacon", "spam")
("spam", "egg", "spam", "spam", "bacon", "spam")
julia> histogram(t)
Dict{Any,Any} with 3 entries:
  "bacon" => 1
  "spam"  => 4
  "egg"   => 1
```

Functions that work with several types are called *polymorphic*. Polymorphism can facilitate code reuse.

For example, the built-in function sum, which adds the elements of a sequence, works as long as the elements of the sequence support addition.

Since a `+` method is provided for `MyTime` objects, they work with `sum`:

```julia
julia> t1 = MyTime(1, 7, 2)
01:07:02
julia> t2 = MyTime(1, 5, 8)
01:05:08
julia> t3 = MyTime(1, 5, 0)
01:05:00
julia> sum((t1, t2, t3))
03:17:10
```

In general, if all of the operations inside a function work with a given type, the function works with that type.

The best kind of polymorphism is the unintentional kind, where you discover that a function you already wrote can be applied to a type you never planned for.
"""

# ╔═╡ d1255e32-8ffb-11eb-336f-17da44827693
md"""## Interface and Implementation

One of the goals of multiple dispatch is to make software more maintainable, which means that you can keep the program working when other parts of the system change, and modify the program to meet new requirements.

A design principle that helps achieve that goal is to keep interfaces separate from implementations. This means that the methods having an argument annotated with a type should not depend on how the fields of that type are represented.

For example, in this chapter we developed a struct that represents a time of day. Methods having an argument annotated with this type include `timetoint`, `isafter`, and `+`.

We could implement those methods in several ways. The details of the implementation depend on how we represent `MyTime`. In this chapter, the fields of a `MyTime` object are `hour`, `minute`, and `second`.

As an alternative, we could replace these fields with a single integer representing the number of seconds since midnight. This implementation would make some functions, like `isafter`, easier to write, but other functions harder.

After you deploy a new type, you might discover a better implementation. If other parts of the program are using your type, it might be time-consuming and errorprone to change the interface.

But if you designed the interface carefully, you can change the implementation without changing the interface, which means that other parts of the program don’t have to change.
"""

# ╔═╡ a563a33e-8ffc-11eb-1bba-23967190b2a7
md"""## Debugging

Calling a function with the correct arguments can be difficult when more than one method for the function is specified. To help with this, Julia allows us to introspect the signatures of the methods of a function.

To know what methods are available for a given function, you can use the function `methods`:

```julia
julia> methods(printtime)
# 2 methods for generic function "printtime":
[1] printtime(time::MyTime) in Main at REPL[3]:2
[2] printtime(time) in Main at REPL[4]:2
```

In this example, the function `printtime` has two methods: one with a `MyTime` argument and one with an `Any` argument.
"""

# ╔═╡ d187fbf4-8ffc-11eb-2b76-75fa897f44e7
md"""## Glossary

*type annotation*:
The operator `::` followed by a type, indicating that an expression or a variable is of that type.

*method*:
A definition of a possible behavior for a function.

*signature*:
The number and type of the arguments of a method, allowing the dispatch to select the most specific method of a function during the function call.

*constructor*:
A special function called to create an object.

*default constructor*:
An inner constructor that is available when no programmer-defined inner constructors are provided.

*outer constructor*:
A constructor defined outside the type definition to define convenience methods for creating an object.

*copy constructor*:
An outer constructor method of a type with as its only argument an object of the type. It creates a new object that is a copy of the argument.

*inner constructor*:
A constructor defined inside the type definition to enforce invariants or to construct self-referential objects.

*operator overloading*:
Adding to the behavior of an operator like + so it works with a programmer-defined type.

*dispatch*:
The choice of which method to execute when a function is executed.

*multiple dispatch*:
Dispatch based on all of a function’s arguments.

*generic programming*:
Writing code that can work with more than one type.

*polymorphic function*:
A function whose argument(s) can be of several types.
"""

# ╔═╡ 2250f414-8ffd-11eb-3cd1-db7908b3544f
md"""## Exercises

#### Exercise 17-4
Change the fields of `MyTime` to be a single integer representing seconds since midnight. Then modify the methods defined in this chapter to work with the new implementation.
"""

# ╔═╡ 5ce1e624-8ffd-11eb-0ca8-63ad122f224d
md"""#### Exercise 17-5

Write a definition for a type named `Kangaroo` with a field named `pouchcontents` of type `Array` and the following methods:

* A constructor that initializes `pouchcontents` to an empty array
* A method named `putinpouch` that takes a `Kangaroo` object and an object of any
type and adds it to `pouchcontents`
* A `show` method that returns a string representation of the `Kangaroo` object and
the contents of the pouch

Test your code by creating two `Kangaroo` objects, assigning them to variables named `kanga` and `roo`, and then adding `roo` to the contents of `kanga`’s pouch.
"""

# ╔═╡ Cell order:
# ╟─d7cb8fd0-8ff7-11eb-3acf-c960087ff109
# ╟─0bb21f1a-8ff8-11eb-3385-718fb337b749
# ╟─1f325d60-8ff8-11eb-1ebe-1f0a41b7149a
# ╟─82ea73c4-8ff8-11eb-2e25-bde15dc079ad
# ╟─f34a98f4-8ff8-11eb-3cd3-f170eae0c30d
# ╟─0aa3ccbe-8ff9-11eb-3ac1-e343e9bb813b
# ╟─a2e88438-8ff9-11eb-02b8-3f084a3f441d
# ╟─286fca94-8ffa-11eb-1b9b-4d36482bd021
# ╟─6021d5b0-8ffa-11eb-053a-85008d0a7f1a
# ╟─757c2abc-8ffa-11eb-1b8c-edb1f8938d82
# ╟─d92ea206-8ffa-11eb-1b9a-5ba6198ccb5f
# ╟─20fdca9e-8ffb-11eb-3066-59f9a3689bf8
# ╟─5c319622-8ffb-11eb-1f0f-3d222d4b7408
# ╟─d1255e32-8ffb-11eb-336f-17da44827693
# ╟─a563a33e-8ffc-11eb-1bba-23967190b2a7
# ╟─d187fbf4-8ffc-11eb-2b76-75fa897f44e7
# ╟─2250f414-8ffd-11eb-3cd1-db7908b3544f
# ╟─5ce1e624-8ffd-11eb-0ca8-63ad122f224d


text:
using QuantumSpatialSearch
using BenchmarkTools

N = 10
psi0 = 1/sqrt(N) * ones(ComplexF64, N)
maxdt = 0.01
println("Precompiling")
spatialsearch(psi0, star_graph_Ad;
                    maxdt=maxdt, # maximum allowed time for dt
                    gamma=1. / N,
                    posW=1,
                    time=2*sqrt(N),
                    mu=0.01,
                    coupling=1.,
                    noiseStrength=1.,
                    noiseRealizations=2,
                    dysonOrder=4);

println("Running benchmark")
res = @benchmark spatialsearch(psi0, star_graph_Ad;
                    maxdt=maxdt, # maximum allowed time for dt
                    gamma=1. / N,
                    posW=1,
                    time=2*sqrt(N),
                    mu=0.01,
                    coupling=1.,
                    noiseStrength=1.,
                    noiseRealizations=100,
                    dysonOrder=4)

max_stars_repo_name: matteoacrossi/QuantumSpatialSearch

max_issues_repo_name: matteoacrossi/QuantumSpatialSearch

max_forks_repo_name: matteoacrossi/QuantumSpatialSearch


text:
module TestConstraintSemiToBinary

using Test

using MathOptInterface
const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
    return
end

include("../utilities.jl")

function test_SemiToBinary()
    mock = MOI.Utilities.MockOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
    )
    config = MOI.Test.Config()
    bridged_mock = MOI.Bridges.Constraint.SemiToBinary{Float64}(mock)
    bridge_type = MOI.Bridges.Constraint.SemiToBinaryBridge{
        Float64,
        MOI.Semiinteger{Float64},
    }
    @test MOI.supports_constraint(
        bridge_type,
        MOI.VariableIndex,
        MOI.Semiinteger{Float64},
    )
    @test MOI.Bridges.Constraint.concrete_bridge_type(
        bridge_type,
        MOI.VariableIndex,
        MOI.Semiinteger{Float64},
    ) == bridge_type
    @test MOI.supports(bridged_mock, MOI.ConstraintPrimalStart(), bridge_type)
    MOI.Test.runtests(
        bridged_mock,
        config,
        include = [
            "test_basic_VariableIndex_Semiinteger",
            "test_basic_VariableIndex_Semicontinuous",
        ],
    )
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 0.0)
            MOI.Utilities.mock_optimize!(mock, [0.0, 0.0, 0.0])
        end,
        (mock::MOI.Utilities.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 2.0)
            MOI.Utilities.mock_optimize!(mock, [2.0, 1.0, 1.0])
        end,
        (mock::MOI.Utilities.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 2.0)
            MOI.Utilities.mock_optimize!(mock, [2.0, 2.0, 1.0])
        end,
        (mock::MOI.Utilities.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 2.5)
            MOI.Utilities.mock_optimize!(mock, [2.5, 2.5, 1.0])
        end,
        (mock::MOI.Utilities.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 3.0)
            MOI.Utilities.mock_optimize!(mock, [3.0, 3.0, 1.0])
        end,
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.set(mock, MOI.TerminationStatus(), MOI.INFEASIBLE),
    )
    MOI.Test.test_linear_Semicontinuous_integration(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VariableIndex,
                MOI.Semicontinuous{Float64},
            }(),
        ),
    )
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        (
            (MOI.VariableIndex, MOI.EqualTo{Float64}, 1),
            (MOI.VariableIndex, MOI.ZeroOne, 0),
            (MOI.VariableIndex, MOI.Integer, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 1),
        ),
    )
    MOI.empty!(bridged_mock)
    MOI.Utilities.set_mock_optimize!(
        mock,
        (mock::MOI.Utilities.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 0.0)
            MOI.Utilities.mock_optimize!(mock, [0.0, 0.0, 0.0])
        end,
        (mock::MOI.Utilities.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 2.0)
            MOI.Utilities.mock_optimize!(mock, [2.0, 1.0, 1.0])
        end,
        (mock::MOI.Utilities.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 2.0)
            MOI.Utilities.mock_optimize!(mock, [2.0, 2.0, 1.0])
        end,
        (mock::MOI.Utilities.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 3.0)
            MOI.Utilities.mock_optimize!(mock, [3.0, 2.5, 1.0])
        end,
        (mock::MOI.Utilities.MockOptimizer) -> begin
            MOI.set(mock, MOI.ObjectiveBound(), 3.0)
            MOI.Utilities.mock_optimize!(mock, [3.0, 3.0, 1.0])
        end,
        (mock::MOI.Utilities.MockOptimizer) ->
            MOI.set(mock, MOI.TerminationStatus(), MOI.INFEASIBLE),
    )
    MOI.Test.test_linear_Semiinteger_integration(bridged_mock, config)
    ci = first(
        MOI.get(
            bridged_mock,
            MOI.ListOfConstraintIndices{
                MOI.VariableIndex,
                MOI.Semiinteger{Float64},
            }(),
        ),
    )

    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), ci) == 3
    new_set = MOI.Semiinteger{Float64}(19.0, 20.0)
    MOI.set(bridged_mock, MOI.ConstraintSet(), ci, new_set)
    @test MOI.get(bridged_mock, MOI.ConstraintSet(), ci) == new_set
    for attr in [MOI.ConstraintPrimalStart()]
        @test MOI.supports(bridged_mock, attr, typeof(ci))
        value = 2.0
        MOI.set(bridged_mock, attr, ci, value)
        @test MOI.get(bridged_mock, attr, ci) ≈ value
    end
    _test_delete_bridge(
        bridged_mock,
        ci,
        2,
        (
            (MOI.VariableIndex, MOI.EqualTo{Float64}, 1),
            (MOI.VariableIndex, MOI.ZeroOne, 0),
            (MOI.VariableIndex, MOI.Integer, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.LessThan{Float64}, 0),
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}, 1),
        ),
    )
    s = """
    variables: x, y
    y == 4.0
    x in Semiinteger(2.0, 3.0)
    minobjective: x
    """
    model = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(model, s)
    sb = """
    variables: x, y, z
    y == 4.0
    z in ZeroOne()
    x in Integer()
    clo: x + -2.0z >= 0.0
    cup: x + -3.0z <= 0.0
    minobjective: x
    """
    modelb = MOI.Utilities.Model{Float64}()
    MOI.Utilities.loadfromstring!(modelb, sb)
    MOI.empty!(bridged_mock)
    @test MOI.is_empty(bridged_mock)
    MOI.Utilities.loadfromstring!(bridged_mock, s)
    MOI.Test.util_test_models_equal(
        bridged_mock,
        model,
        ["x", "y"],
        String[],
        [
            ("x", MOI.Semiinteger{Float64}(2.0, 3.0)),
            ("y", MOI.EqualTo{Float64}(4.0)),
        ],
    )
    # setting names on mock
    ci = first(
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.ZeroOne}(),
        ),
    )
    z = MOI.VariableIndex(ci.value)
    MOI.set(mock, MOI.VariableName(), z, "z")
    ci = first(
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{MOI.VariableIndex,MOI.Integer}(),
        ),
    )
    ci = first(
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Float64},
                MOI.GreaterThan{Float64},
            }(),
        ),
    )
    MOI.set(mock, MOI.ConstraintName(), ci, "clo")
    ci = first(
        MOI.get(
            mock,
            MOI.ListOfConstraintIndices{
                MOI.ScalarAffineFunction{Float64},
                MOI.LessThan{Float64},
            }(),
        ),
    )
    MOI.set(mock, MOI.ConstraintName(), ci, "cup")
    MOI.Test.util_test_models_equal(
        mock,
        modelb,
        ["x", "y", "z"],
        ["clo", "cup"],
        [
            ("x", MOI.Integer()),
            ("y", MOI.EqualTo{Float64}(4.0)),
            ("z", MOI.ZeroOne()),
        ],
    )
    return
end

"""
    test_lower_bound_already_set()

The second call to `add_constraint` is broken because it should throw:
```julia
MOI.LowerBoundAlreadySet{
    MOI.Semicontinuous{Float64},
    MOI.GreaterThan{Float64},
}
```
See MathOptInterface issue #1431.
"""
function test_lower_bound_already_set()
    model = MOI.Utilities.Model{Float64}()
    bridged = MOI.Bridges.Constraint.SemiToBinary{Float64}(model)
    x = MOI.add_variable(bridged)
    MOI.add_constraint(bridged, x, MOI.Semicontinuous(1.0, 2.0))
    @test_broken(
        MOI.add_constraint(bridged, x, MOI.GreaterThan(0.0)) === nothing,
    )
    return
end

end  # module

TestConstraintSemiToBinary.runtests()


text:
using Pkg
Pkg.activate(".")
println("Loading SpineBasedRecordLinkage.jl")
using SpineBasedRecordLinkage
if length(ARGS) != 2
    error("The number of arguments is $(length(ARGS)), should be 2.")
end
directory1 = ARGS[1]
outfile    = ARGS[2]
summarise_linkage_run(directory1, outfile)


text:
# Usage:
#   julia run_one_benchmark.jl [<benchmark> <#runs> <#threads>]

using Statistics
using Serialization
using Printf

gctime(stat) = stat.total_time

const BENCH = isempty(ARGS) ? "benches/append/append.jl" : ARGS[1]
const RUNS = isempty(ARGS) ? 10 : parse(Int, ARGS[2])
const THREADS = isempty(ARGS) ? 1 : parse(Int, ARGS[3])
const JULIAVER = Base.julia_cmd()[1]

path = split(BENCH, "/")
file = pop!(path)
path = join(path, "/")

dir = joinpath(@__DIR__, path)
cd(dir)
@printf("Running Benchmark path = %s file = %s times = %d threads = %d\n", path, file, RUNS, THREADS)
@show BENCH

value = []
times = []
stats = []
for _ in 1:RUNS
   r = open(deserialize, `$JULIAVER --project=. --threads=$THREADS $file SERIALIZE`)
   push!(value, r.value)
   push!(times, r.times)
   push!(stats, r.stats)
end
@printf("run time: %0.0fms min, %0.0fms max %0.0fms median\n",
   minimum(times) / 1_000_000,
   maximum(times) / 1_000_000,
   median(times)  / 1_000_000)

time = map(gctime, stats)
@printf("gc time: %0.0fms min, %0.0fms max, %0.0fms median\n",
   minimum(time)/ 1_000_000,
   maximum(time)/ 1_000_000,
   median(time) / 1_000_000)

max_stars_repo_name: oscardssmith/GCBenchmarks

max_issues_repo_name: oscardssmith/GCBenchmarks

max_forks_repo_name: oscardssmith/GCBenchmarks