text:
struct RenderSettings

    # horizontal lines. if character, repeat.
    toprule::String
    midrule::String
    bottomrule::String

    headerrule::Function  # Function that takes the headerCellStartEnd array and returns a
                            # sting that describes the text line below the header titles
                            # if it's one single character, it's put into a table format and repeated (e.g. "-" for ascii)

    headercolsep::String # separator between columns for header
    colsep::String  # separator between columns
    linestart::String # start of each line string
    linebreak::String   # link break string

    label_fe_yes::String    # what the FE block prints if the FE is present. override with __LABEL_FE_YES__ in 'label' argument
    label_fe_no::String    # what the FE block prints if the FE is not present. override with __LABEL_FE_NO__ in 'label' argument

    label_statistic_n::String # label for number of observations. override with __LABEL_STATISTIC_N__ in 'label' argument
    label_statistic_r2::String # label for R^2. override with __LABEL_STATISTIC_R2__ in 'label' argument
    label_statistic_adjr2::String # label for adjusted R^2. override with __LABEL_STATISTIC_adjr2__
    label_statistic_r2_within::String # label for within-R^2. override with __LABEL_STATISTIC_R2_WITHIN__
    label_statistic_f::String # label for F-Stat. override with __LABEL_STATISTIC_F__ in 'label' argument
    label_statistic_p::String # label for F-test p value. override with __LABEL_STATISTIC_P__
    label_statistic_f_kp::String # label for first-stage F statistic. override with __LABEL_STATISTIC_F_KP__
    label_statistic_p_kp::String # label for first-stage F-stat p value. override with __LABEL_STATISTIC_P_KP__
    label_statistic_dof::String # label for degrees of freedom. override with __LABEL_STATISTIC_DOF__

    label_estimator::String # label for the Estimator block. Override with __LABEL_ESTIMATOR__
    label_estimator_ols::String # label for the Estimator block. Override with __LABEL_ESTIMATOR_OLS__
    label_estimator_iv::String # label for the Estimator block. Override with __LABEL_ESTIMATOR_IV__
    label_estimator_nl::String # label for the Estimator block. Override with __LABEL_ESTIMATOR_NL__

    outfile::String    # file to print output into.
                       # if empty, print to STDOUT.

    encapsulateRegressand::Function     # function that takes a string and
                                        # min and max column index and returns
                                        # a formatted string

    header::Function                    # function that return the header string
    footer::Function                    # function that returns the footer string
                                        # both should take the number of results and align string as arguments

end


text:
#!/usr/bin/env julia
# -*- coding: UTF-8 -*-
# __julia-version__ = 1.7.2
# __author__        = "Martin Scheidt, Max Kannenberg"
# __copyright__     = "2022"
# __license__       = "ISC"

"""
    Settings(file)

Settings is a datastruture for calculation context.
The function Settings() will create a set of settings for the train run calculation.
`file` is optinal may be used to load settings in the YAML format.

# Example
```jldoctest
julia> my_settings = Settings() # will generate default settings
Settings(mass_point, :distance, 20, 3, running_time, :dataframe)
```
"""
function Settings(file="DEFAULT")

    ## default values
    massModel    = :mass_point
    stepVariable = :distance
    stepSize     = 20
    approxLevel  = 3
    outputDetail = :running_time
    outputFormat = :dataframe

    ## load from file
    if file != "DEFAULT"

        ## JSON schema for YAML-file validation
        schema = Schema("""{
            "properties": {
                "massModel": {
                    "description": "type of train model",
                    "type": "string",
                    "enum": [ "mass_point", "homogeneous_strip" ]
                },
                "stepVariable": {
                    "description": "variable of the linear multistep method",
                    "type": "string",
                    "enum": [ "distance", "time", "velocity" ]
                },
                "stepSize": {
                    "description": "step size acording to stepVariable",
                    "type": "number",
                    "exclusiveMinimum": 0
                },
                "outputDetail": {
                    "description": "Selecting the detail of the result",
                    "type": "string",
                    "enum": [ "running_time", "points_of_interest", "driving_course", "everything" ]
                },
                "outputFormat": {
                    "description": "Output format",
                    "type": "string",
                    "enum": [ "dataframe", "dict" ]
                }
            }
        }""")

        settings = YAML.load(open(file))["settings"]

        ## validate the loaded file
        try
            validate(schema, settings)
        catch err
            println("Could not load settings file '$file'.\n Format is not recognized - using default as fall back.")
            settings = Dict()
        end

        ## set the variables in "settings"
        haskey(settings, "massModel")    ? massModel    = Symbol(settings["massModel"])    : nothing
        haskey(settings, "stepVariable") ? stepVariable = Symbol(settings["stepVariable"]) : nothing
        haskey(settings, "stepSize")     ? stepSize     =        settings["stepSize"]      : nothing
        haskey(settings, "approxLevel")  ? approxLevel  =        settings["approxLevel"]   : nothing
        haskey(settings, "outputDetail") ? outputDetail = Symbol(settings["outputDetail"]) : nothing
        haskey(settings, "outputFormat") ? outputFormat = Symbol(settings["outputFormat"]) : nothing
    end

    Settings(massModel, stepVariable, stepSize, approxLevel, outputDetail, outputFormat)

end #function Settings() # outer constructor

"""
    Path(file, type = :YAML)

Path is a datastruture for calculation context.
The function Path() will create a running path for the train.
Supported formats are: railtoolkit/schema (2022.05)

# Example
```jldoctest
julia> my_path = Path("file.yaml") # will generate a path from a YAML file.
Path(variables)
```
"""
function Path(file, type = :YAML)

    ## default values
    name      = ""
    id        = ""
    uuid      = UUIDs.uuid4()
    poi       = []
    sections  = []

    ## process flags
    POI_PRESENT = false

    ## load from file
    if type == :YAML

        data = YAML.load(open(file))
        if data["schema"] != "https://railtoolkit.org/schema/running-path.json"
            error("Could not load path file '$file'.\n
                    YAML format is not recognized.
                    Currently supported: railtoolkit/schema/running-path (2022.05)")
        end
        if data["schema_version"] != "2022.05"
            error("Could not load path file '$file'.\n
                    YAML format is not recognized.
                    Currently supported: railtoolkit/schema/running-path (2022.05)")
        end

        ## JSON schema for YAML-file validation
        railtoolkit_schema = Schema("""{
            "required": [ "schema", "schema_version", "paths" ],
            "properties": {
                "schema": {
                "description": "Identifier of the schema",
                "enum": [ "https://railtoolkit.org/schema/running-path.json" ]
                },
                "schema_version": {
                "description": "Version of the schema",
                "type": "string",
                "pattern": "[2-9][0-9][0-9][0-9].[0-1][0-9]"
                },
                "paths": {
                "type": "array",
                "minItems": 1,
                "items": {
                    "required": [ "name", "id", "characteristic_sections" ],
                    "type": "object",
                    "properties": {
                    "characteristic_sections": {
                        "description": "",
                        "type": "array",
                        "minItems": 2,
                        "uniqueItems": true,
                        "items": {
                        "type": "array",
                        "minItems": 3,
                        "maxItems": 3,
                        "description": "",
                        "prefixItems": [
                            {
                            "description": "milage in meter",
                            "type": "number"
                            },
                            {
                            "description": "speed in kilometers per hour",
                            "type": "number",
                            "exclusiveMinimum": 0
                            },
                            {
                            "description": "resistance in permil",
                            "type": "number"
                            }
                        ]
                        }
                    },
                    "id": {
                        "description": "Identifier of the path",
                        "type": "string"
                    },
                    "name": {
                        "description": "Name of the path",
                        "type": "string"
                    },
                    "points_of_interest": {
                        "description": "",
                        "type": "array",
                        "uniqueItems": true,
                        "items": {
                        "type": "array",
                        "minItems": 3,
                        "maxItems": 3,
                        "description": "",
                        "prefixItems": [
                            { "type": "number" },
                            { "type": "string" },
                            { "enum": [ "front", "rear" ] }
                        ]
                        }
                    },
                    "UUID": {
                        "description": "The unique identifier for a path",
                        "type": "string",
                        "format": "uuid"
                    }
                    }
                }
                }
            }
        }""")

        paths = data["paths"]
        try
            validate(railtoolkit_schema, paths)
        catch err
            error("Could not load path file '$file'.\n
                    YAML format is not recognized.
                    Currently supported: railtoolkit/schema/running-path (2022.05)")
        end
        if length(paths) > 1
            println("WARNING: the loaded file contains more than one path. Using only the first!")
        end
        path = paths[1]

        ## set the variables in "path"
        # required
        name    = path["name"]
        id      = path["id"]
        tmp_sec = path["characteristic_sections"]
        # optional
        haskey(path, "UUID")               ? uuid        = parse(UUID, path["UUID"] ) : nothing
        haskey(path, "points_of_interest") ? POI_PRESENT = true                       : nothing
        haskey(path, "points_of_interest") ? tmp_points  = path["points_of_interest"] : nothing

    else
        error("Unknown file type '$type'")
    end #if type

    ## process characteristic sections
    sort!(tmp_sec, by = x -> x[1])
    for row in 2:length(tmp_sec)
        s_start = tmp_sec[row-1][1]     # first point of the section (in m)
        s_end   = tmp_sec[row][1]       # first point of the next section (in m)
        v_limit = tmp_sec[row-1][2]/3.6 # paths speed limt (in m/s)
        f_Rp    = tmp_sec[row-1][3]     # specific path resistance of the section (in ‰)

        section = Dict(:s_start => s_start,
                        :s_end   => s_end,
                        :v_limit => v_limit,
                        :f_Rp    => f_Rp)
        push!(sections, section)
    end #for row
    # s_start in first entry defines the path's beginning
    # s_end in last entry defines the path's ending

    ## process points of interest
    if POI_PRESENT
        sort!(tmp_points, by = x -> x[1])
        for elem in tmp_points
            station = elem[1]     # first point of the section (in m)
            label   = elem[2] # paths speed limt (in m/s)
            measure = elem[3]     # specific path resistance of the section (in ‰)

            point = Dict(:station => station,
                            :label   => label,
                            :measure => measure)
            push!(poi, point)
        end #for elem
    end #if !isempty(points)

    Path(name, id, uuid, poi, sections)

end #function Path() # outer constructor

"""
    Train(file, type = :YAML)

Train is a datastruture for calculation context.
The function Train() will create a train to use in calculations.
Supported formats for the YAML files are: railtoolkit/schema (2022.05)

# Example
```jldoctest
julia> my_train = Train("file.yaml") # will generate a train from a YAML file.
Train(variables)
```
"""
function Train(file, type = :YAML)

    ## default values
    name          = ""            #
    id            = ""            #
    uuid          = UUIDs.uuid4() #
    length        = 0             # in meter
    m_train_full  = 0             # in kilogram
    m_train_empty = 0             # in kilogram
    m_loco        = 0             # in kilogram
    m_td          = 0             # in kilogram
    m_tc          = 0             # in kilogram
    m_car_full    = 0             # in kilogram
    m_car_empty   = 0             # in kilogram
    ξ_train       = 1.08          # rotation mass factor, source: "Fahrdynamik des Schienenverkehrs" by Wende, 2003, p. 13 for "Zug, überschlägliche Berechnung"
    ξ_loco        = 1.09          # rotation mass factor
    ξ_cars        = 1.06          # rotation mass factor
    transportType = :freight      # "freight" or "passenger" for resistance calculation
    v_limit       = 140           # in m/s (default 504 km/h)
    a_braking     = 0             # in m/s^2, todo: implement as function
    f_Rtd0        = 0             # coefficient for basic resistance due to the traction units driving axles (in ‰)
    f_Rtc0        = 0             # coefficient for basic resistance due to the traction units carring axles (in ‰)
    F_Rt2         = 3000          # coefficient for air resistance of the traction units (in N)
    f_Rw0         = 0             # coefficient for the consists basic resistance (in ‰)
    f_Rw1         = 0             # coefficient for the consists resistance to rolling (in ‰)
    f_Rw2         = 0             # coefficient fo the consistsr air resistance (in ‰)
    F_v_pairs     = []            # [v in m/s, F_T in N]

    ## load from file
    if type == :YAML

        data = YAML.load(open(file))
        if data["schema"] != "https://railtoolkit.org/schema/rolling-stock.json"
            error("Could not load path file '$file'.\n
                    YAML format is not recognized.
                    Currently supported: railtoolkit/schema/rolling-stock (2022.05)")
        end
        if data["schema_version"] != "2022.05"
            error("Could not load path file '$file'.\n
                    YAML format is not recognized.
                    Currently supported: railtoolkit/schema/rolling-stock (2022.05)")
        end

        ## JSON schema for YAML-file validation
        railtoolkit_schema = Schema("""{
            "required": [ "schema", "schema_version" ],
            "anyOf": [
            {"required": [ "trains" ] },
            {"required": [ "vehicles" ] }
            ],
            "properties": {
            "schema": {
                "description": "Identifier of the schema",
                "enum": [ "https://railtoolkit.org/schema/rolling-stock.json" ]
            },
            "schema_version": {
                "description": "Version of the schema",
                "type": "string",
                "pattern": "[2-9][0-9][0-9][0-9].[0-1][0-9]"
            },
            "trains": {
                "type": "array",
                "minItems": 1,
                "items": {
                "required": [ "name", "id", "formation" ],
                "type": "object",
                "properties": {
                    "id": {
                    "description": "Identifier of the train",
                    "type": "string"
                    },
                    "name": {
                    "description": "Name of the train",
                    "type": "string"
                    },
                    "UUID": {
                    "description": "The unique identifier for a train",
                    "type": "string",
                    "format": "uuid"
                    },
                    "formation": {
                    "description": "Collection of vehicles that form the train",
                    "type": "array",
                    "minItems": 1,
                    "uniqueItems": false,
                    "items": {
                        "type": "string"
                    }
                    }
                }
                }
            },
            "vehicles": {
                "type": "array",
                "minItems": 1,
                "items": {
                "required": [ "name", "id", "vehicle_type", "length", "mass" ],
                "type": "object",
                "properties": {
                    "air_resistance": {
                    "description": "coefficient for air resistance in permil",
                    "type": "number",
                    "exclusiveMinimum": 0
                    },
                    "base_resistance": {
                    "description": "coefficient for basic resistance in permil",
                    "type": "number",
                    "exclusiveMinimum": 0
                    },
                    "id": {
                    "description": "Identifier of the vehicle",
                    "type": "string"
                    },
                    "length": {
                    "description": "The length of the vehicle in meter",
                    "type": "number",
                    "exclusiveMinimum": 0
                    },
                    "load_limit": {
                    "description": "The maximum permitted load of the vehicle in metric ton",
                    "type": "number",
                    "exclusiveMinimum": 0
                    },
                    "mass_traction": {
                    "description": "The mass on the powered axles of the vehicle in metric ton",
                    "type": "number",
                    "exclusiveMinimum": 0
                    },
                    "mass": {
                    "description": "The empty mass of the vehicle in metric ton",
                    "type": "number",
                    "exclusiveMinimum": 0
                    },
                    "name": {
                    "description": "Name of the vehicle",
                    "type": "string"
                    },
                    "picture": {
                    "description": "A URI with a picture for humans",
                    "type": "string",
                    "format": "uri"
                    },
                    "power_type": {
                    "description": "Type of propulsion",
                    "enum": [ "diesel", "electric", "steam" ]
                    },
                    "rolling_resistance": {
                    "description": "coefficient for resistance of rolling axles in permil",
                    "type": "number",
                    "exclusiveMinimum": 0
                    },
                    "rotation_mass": {
                    "description": "Factor for rotating mass; >= 1",
                    "type": "number",
                    "minimum": 1
                    },
                    "speed_limit": {
                    "description": "Maximum permitted speed in kilometers per hour",
                    "type": "number",
                    "exclusiveMinimum": 0
                    },
                    "tractive_effort": {
                    "description": "Tractive effort as pairs of speed in kilometers per hour and tractive force in newton",
                    "type": "array",
                    "minItems": 3,
                    "uniqueItems": true,
                    "items": {
                        "type": "array",
                        "minItems": 2,
                        "maxItems": 2,
                        "uniqueItems": true,
                        "items": {
                        "type": "number",
                        "minimum": 0
                        }
                    }
                    },
                    "UUID": {
                    "description": "The unique identifier for a vehicle",
                    "type": "string",
                    "format": "uuid"
                    },
                    "vehicle_type": {
                    "description": "Type of vehicle",
                    "enum": [ "traction unit", "freight", "passenger", "multiple unit" ]
                    }
                }
                }
            }
            }
        }""")

        try
            validate(railtoolkit_schema, data)
        catch err
            error("Could not load path file '$file'.\n
                    YAML format is not recognized.
                    Currently supported: railtoolkit/schema/rolling-stock (2022.05)")
        end

    else
        error("Unknown file type '$type'")
    end #if type

    trains = data["trains"]
    Base.length(trains) > 1 ? println("WARNING: the loaded file contains more than one train. Using only the first!") : nothing
    Base.length(trains) == 0 ? error("No train present in file '$file'") : nothing
    train = trains[1]
    used_vehicles = unique(train["formation"])

    included_vehicles = []
    for vehicle in data["vehicles"]
        push!(included_vehicles,vehicle["id"])
    end

    ## test if all vehicles of the formation are avilable
    for vehicle in used_vehicles
        vehicle ∉ included_vehicles ? error("'$vehicle' is not present in '$file'") : nothing
    end

    ## gather the count of vehicles and usage in the formation
    vehicles = NamedTuple[]
    for vehicle in data["vehicles"]
        if vehicle["id"] in used_vehicles
            n = count(==(vehicle["id"]),train["formation"])
            type = vehicle["vehicle_type"]
            type == "traction unit" || type == "multiple unit" ? propulsion    = true       : propulsion = false
            type == "passenger"     || type == "multiple unit" ? transportType = :passenger : nothing
            push!(vehicles, (data=vehicle, n=n, propulsion=propulsion) )
        end
    end

    ## set the variables in "train"
    name = train["name"]
    id   = train["id"]
    haskey(train, "UUID") ? uuid = parse(UUID, train["UUID"] ) : nothing
    transportType == :freight ? a_braking = -0.225 : a_braking = -0.375  # set a default a_braking value depending on the train type

    ## set the variables for all vehicles
    for vehicle in vehicles
        length           += vehicle.data["length"] * vehicle.n
        m_train_full     += vehicle.data["mass"]   * vehicle.n * 1000 # in kg
        m_train_empty    += vehicle.data["mass"]   * vehicle.n * 1000 # in kg
        haskey(vehicle.data, "load_limit")    ?
            m_train_full += vehicle.data["load_limit"] * vehicle.n * 1000 :  # in kg
            nothing
        haskey(vehicle.data, "speed_limit")   ?
            v_limit > vehicle.data["speed_limit"]/3.6 ? v_limit = vehicle.data["speed_limit"]/3.6 : nothing :
            nothing
    end

    ## divide vehicles in propulsion and non-propulsion
    loco = []
    for i in 1:Base.length(vehicles)
        if vehicles[i].propulsion
            push!(loco, vehicles[i])
            deleteat!(vehicles, i)
        end
    end
    Base.length(loco)  > 1 ? println("WARNING: the loaded file contains more than one traction unit or multiple unit. Using only the first!") : nothing
    loco[1].n     > 1 ? println("WARNING: the loaded file contains more than one traction unit or multiple unit. Using only one!") : nothing
    Base.length(loco) == 0 ? error("No traction unit or multiple unit present in file '$file'") : nothing
    loco = loco[1].data
    cars = vehicles

    ## set the variables for locos
    m_loco= loco["mass"] * 1000
    haskey(loco, "a_braking")          ? a_braking = loco["a_braking"]                : nothing
    haskey(loco, "base_resistance")    ? f_Rtd0 = loco["base_resistance"]             : nothing
    haskey(loco, "rolling_resistance") ? f_Rtc0 = loco["rolling_resistance"]          : nothing
    haskey(loco, "air_resistance")     ? F_Rt2  = loco["air_resistance"] * g * m_loco : nothing
    haskey(loco, "mass_traction")      ? m_td   = loco["mass_traction"] * 1000        : m_td = m_t
    haskey(loco, "rotation_mass")      ? ξ_loco = loco["rotation_mass"]               : nothing
    m_tc = m_loco- m_td
    haskey(loco, "tractive_effort")    ? F_v_pairs = loco["tractive_effort"] : F_v_pairs = [ [0.0, m_td * g * μ],[v_limit*3.6, m_td * g * μ] ]
    F_v_pairs = reduce(hcat,F_v_pairs)'       # convert to matrix
    F_v_pairs[:,1] ./= 3.6                    # convert km/h to m/s
    F_v_pairs = tuple.(eachcol(F_v_pairs)...) # convert each row to tuples

    ## set the variables for cars
    if !isempty(cars)
        resis_base = []
        resis_roll = []
        resis_air  = []
        rotMassFac = []
        for car in cars
            haskey(car.data, "base_resistance")    ?
                append!(resis_base,repeat([car.data["base_resistance"]],car.n))    :
                append!(resis_base,repeat([f_Rw0],car.n))
            haskey(car.data, "rolling_resistance") ?
                append!(resis_roll,repeat([car.data["rolling_resistance"]],car.n)) :
                append!(resis_roll,repeat([f_Rw1],car.n))
            haskey(car.data, "air_resistance")     ?
                append!(resis_air,repeat([car.data["air_resistance"]],car.n))      :
                append!(resis_air, repeat([f_Rw2],car.n))
            haskey(car.data, "rotation_mass") ?
                append!(rotMassFac,repeat([(car.data["rotation_mass"],car.data["mass"])],car.n)) :
                append!(rotMassFac,repeat([(ξ_cars ,car.data["mass"])],car.n))
            m_car_empty += car.data["mass"] * car.n * 1000 # in kg
            m_car_full  += car.data["mass"] * car.n * 1000 # in kg
            haskey(car.data, "load_limit")    ?
                m_car_full += car.data["load_limit"] * car.n * 1000 :  # in kg
                nothing
        end
        f_Rw0   = Statistics.mean(resis_base)
        f_Rw1   = Statistics.mean(resis_roll)
        f_Rw2   = Statistics.mean(resis_air)
        carRotMass = 0
        for elem in rotMassFac
            carRotMass += elem[1]*elem[2] * 1000 # in kg
        end
        ξ_cars  = carRotMass/m_car_empty
        ξ_train = (ξ_loco * m_loco+ carRotMass)/m_train_empty
    else
        ξ_cars  = 0
        ξ_train = ξ_loco
    end

    Train(
        name, id, uuid, length,
        m_train_full, m_td, m_tc, m_car_full,
        ξ_train, ξ_loco, ξ_cars,
        transportType, v_limit,
        a_braking,
        f_Rtd0, f_Rtc0, F_Rt2, f_Rw0, f_Rw1, f_Rw2,
        F_v_pairs
    )

end #function Train() # outer constructor

## create a moving section containing characteristic sections
function createMovingSection(path::Path, v_trainLimit::Real, s_trainLength::Real)
    # this function creates and returns a moving section dependent on the paths attributes

    s_entry = path.sections[1][:s_start]          # first position (in m)
    s_exit = path.sections[end][:s_end]           # last position (in m)
    pathLength = s_exit - s_entry                   # total length (in m)

    CSs=Vector{Dict}()
    s_csStart=s_entry
    csId=1
    for row in 2:length(path.sections)
        previousSection = path.sections[row-1]
        currentSection = path.sections[row]
        speedLimitIsDifferent = min(previousSection[:v_limit], v_trainLimit) != min(currentSection[:v_limit], v_trainLimit)
        pathResistanceIsDifferent = previousSection[:f_Rp] != currentSection[:f_Rp]
        if speedLimitIsDifferent || pathResistanceIsDifferent
        # 03/09 old: if min(previousSection[:v_limit], v_trainLimit) != min(currentSection[:v_limit], v_trainLimit) || previousSection[:f_Rp] != currentSection[:f_Rp]
            push!(CSs, createCharacteristicSection(csId, s_csStart, previousSection, min(previousSection[:v_limit], v_trainLimit), s_trainLength, path))
            s_csStart = currentSection[:s_start]
            csId = csId+1
        end #if
    end #for
    push!(CSs, createCharacteristicSection(csId, s_csStart, path.sections[end], min(path.sections[end][:v_limit], v_trainLimit), s_trainLength, path))

    movingSection= Dict(:id => 1,                       # identifier    # if there is more than one moving section in a later version of this tool the id should not be constant anymore
                        :length => pathLength,          # total length (in m)
                        :s_entry => s_entry,            # first position (in m)
                        :s_exit => s_exit,              # last position (in m)
                        :t => 0.0,                      # total running time (in s)
                        :E => 0.0,                      # total energy consumption (in Ws)
                        :characteristicSections => CSs) # list of containing characteristic sections

    return movingSection
end #function createMovingSection

## create a characteristic section for a path section. A characteristic section is a part of the moving section. It contains behavior sections.
function createCharacteristicSection(id::Integer, s_entry::Real, section::Dict, v_limit::Real, s_trainLength::Real, path::Path)
    # Create and return a characteristic section dependent on the paths attributes
    characteristicSection= Dict(:id => id,                            # identifier
                                :s_entry => s_entry,                    # first position (in m)
                                :s_exit => section[:s_end],             # last position  (in m)
                                :length => section[:s_end] -s_entry,    # total length  (in m)
                                :r_path => section[:f_Rp],              # path resistance (in ‰)
                                :behaviorSections => Dict(),            # list of containing behavior sections
                                :t => 0.0,                              # total running time (in s)
                                :E => 0.0,                              # total energy consumption (in Ws)
                                :v_limit => v_limit,                    # speed limit (in m/s)
                                # initializing :v_entry, :v_peak and :v_exit with :v_limit
                                :v_peak => v_limit,                     # maximum reachable speed (in m/s)
                                :v_entry => v_limit,                    # maximum entry speed (in m/s)
                                :v_exit => v_limit)                     # maximum exit speed (in m/s)

    # list of positions of every point of interest (POI) in this charateristic section for which data points should be calculated
    s_exit = characteristicSection[:s_exit]

    ##TODO: use a tuple with naming
    pointsOfInterest = Tuple[]
    # pointsOfInterest = Real[]
    if !isempty(path.poi)
        for POI in path.poi
            s_poi = POI[:station]
            if POI[:measure] == "rear"
                s_poi -= s_trainLength
            end
            if s_entry < s_poi && s_poi < s_exit
                push!(pointsOfInterest, (s_poi, POI[:label]) )
                # push!(pointsOfInterest, s_poi )
            end
        end
    end
    push!(pointsOfInterest, (s_exit,""))     # s_exit has to be the last POI so that there will always be a POI to campare the current position with
    # push!(pointsOfInterest, s_exit)     # s_exit has to be the last POI so that there will always be a POI to campare the current position with

    merge!(characteristicSection, Dict(:pointsOfInterest => pointsOfInterest))

    return characteristicSection
end #function createCharacteristicSection

"""
a DataPoint is the smallest element of the driving course. One step of the step approach is between two data points
"""
function createDataPoint()
    dataPoint = Dict(
        :i => 0,            # identifier and counter variable of the driving course
        :behavior => "",    # type of behavior section the data point is part of - see createBehaviorSection()
                            # a data point which is the last point of one behavior section and the first point of the next behavior section will be attached to the latter
        :s => 0.0,          # position (in m)
        :Δs => 0.0,         # step size (in m)
        :t => 0.0,          # point in time (in s)
        :Δt => 0.0,         # step size (in s)
        :v => 0.0,          # velocity (in m/s)
        :Δv => 0.0,         # step size (in m/s)
        :a => 0.0,          # acceleration (in m/s^2)
        :W => 0.0,          # mechanical work (in Ws)
        :ΔW => 0.0,         # mechanical work in this step (in Ws)
        :E => 0.0,          # energy consumption (in Ws)
        :ΔE => 0.0,         # energy consumption in this step (in Ws)
        :F_T => 0.0,        # tractive effort (in N)
        :F_R => 0.0,        # resisting force (in N)
        :R_path => 0.0,     # path resistance (in N)
        :R_train => 0.0,    # train resistance (in N)
        :R_traction => 0.0, # traction unit resistance (in N)
        :R_wagons => 0.0,   # set of wagons resistance (in N)
        :label => ""        # a label for importend points
    )
    return dataPoint
end #function createDataPoint

"""
BehaviorSection() TODO!
"""
function createBehaviorSection(type::String, s_entry::Real, v_entry::Real, startingPoint::Integer)
    BS= Dict(
        :type => type,                 # type of behavior section: "breakFree", "clearing", "accelerating", "cruising", "downhillBraking", "diminishing", "coasting", "braking" or "standstill"
        :length => 0.0,                # total length  (in m)
        :s_entry => s_entry,           # first position (in m)
        :s_exit => 0.0,                # last position  (in m)
        :t => 0.0,                     # total running time (in s)
        :E => 0.0,                     # total energy consumption (in Ws)
        :v_entry => v_entry,           # entry speed (in m/s)
        :v_exit => 0.0,                # exit speed (in m/s)
        :dataPoints => [startingPoint] # list of identifiers of the containing data points starting with the initial point
    )
    return BS
end #function createBehaviorSection

ext: jl


text:
function gaussian(windowLength::Integer, sigma::Float64 = 1.0)
  samplingStep = 6/(windowLength-1)
  return [exp(-((x/sigma)^2)) for x in -3:samplingStep:3]
end


text:
module TestMissing

using Test
using RLEVectors

standard_rle = RLEVector(["a", "d", "c", "b"], [2, 4, 6, 8])
missing_ok_rle = RLEVector(Vector{Union{Missing,String}}(["a", "d", "c", "b"]), [2, 4, 6, 8])
with_missing_rle = RLEVector(Vector{Union{Missing,String}}(["a", missing, "c", "b"]), [2, 4, 6, 8])

@testset "Staring with simple vector" begin
    x = standard_rle
    @test disallowmissing(x) == standard_rle
    @test allowmissing(x) == missing_ok_rle
end

@testset "Staring with missings possible but not present" begin
    x = missing_ok_rle
    @test allowmissing(x) == x
    @test disallowmissing(x) == standard_rle
end

@testset "Staring with missing possible" begin
    x = with_missing_rle
    @test_throws MethodError disallowmissing(x)
    x = with_missing_rle
    @test typeof(allowmissing(x)) == typeof(x)
end

end # module


text:
"""
```julia
rate(q::Integer, M::Integer, n::Integer) -> Real
```

Calculate the rate of a code.  That is, how efficient the code is.

Parameters:
  - `q::Integer`: the number of symbols in the code.
  - `M::Integer`: the size/number of elements in the code.
  - `n::Integer`: The word length.

Returns:
  - `Real`: Rate of the code.

---

### Examples

```julia
julia> rate(3, 5, 4) # the rate of the code which has 3 symbols, 5 words in the code, and word length of 4 (e.g., Σ = {A, B, C}, C = {ABBA,CABA,BBBB,CAAB,ACBB})
0.3662433801794817
```
"""
rate(q::T, M::T, n::T) where {T <: Integer} = log(q, M) / n

__spheres(q::T, n::T, r::T) where {T <: Integer} = sum(Integer[((big(q) - 1)^i) * binomial(big(n), big(i)) for i in 0:r])
__sphere_bound(round_func::Function, q::T, n::T, d::T) where {T <: Integer} = round_func((big(q)^n) / __spheres(q, n, d))

"""
```julia
sphere_covering_bound(q::Integer, n::Integer, d::Integer) -> Integer
```

Computes the sphere covering bound of a ``[n, d]_q``-code.

Parameters:
  - `q::Integer`: the number of symbols in the code.
  - `n::Integer`: the word length.
  - `d::Integer`: the distance of the code.

Returns:
  - `Integer`: the sphere covering bound.

---

### Examples

```julia
julia> sphere_covering_bound(5,7,3)
215
```
"""
sphere_covering_bound(q::T, n::T, d::T) where {T <: Integer} = __sphere_bound(ceil, q, n, d - 1)

"""
```julia
sphere_packing_bound(q::Integer, n::Integer, d::Integer) -> Integer
sphere_packing_bound(q::Integer, n::Integer, d::Integer, ::Rounding) -> Real
```

Computes the sphere packing bound of a ``[n, d]_q``-code.  The sphere packing bound is also known as the hamming bound.  You can use `hamming_bound` to compute the same thing.

Parameters:
  - `q::Integer`: the number of symbols in the code.
  - `n::Integer`: the word length.
  - `d::Integer`: the distance of the code.
  - `::Rounding`: use the argument `no_round` in this position to preserve the rounding of the code — which usually by default rounds down.

Returns:
  - `Integer`: the sphere packing bound.

---

### Examples

```julia
julia> sphere_packing_bound(5,7,3)
2693
```
"""
sphere_packing_bound(q::T, n::T, d::T) where T <: Integer =
	__sphere_bound(a -> floor(T, a), q, n, floor(T, (d - 1) / 2))
sphere_packing_bound(q::T, n::T, d::T, ::Rounding) where T <: Integer =
	__sphere_bound(identity, q, n, floor(T, (d - 1) / 2))
hamming_bound(q::T, n::T, d::T) where T <: Integer =
	sphere_packing_bound(q, n, d)
hamming_bound(q::T, n::T, d::T, ::Rounding) where T <: Integer =
	sphere_packing_bound(q, n, d, no_round)

"""
```julia
sphere_packing_bound(q::Integer, n::Integer, d::Integer) -> Real
```

Computes the Singleton bound of a ``[n, d]_q``-code.

Parameters:
  - `q::Integer`: the number of symbols in the code.
  - `n::Integer`: the word length.
  - `d::Integer`: the distance of the code.

Returns:
  - `Real`: the Singleton bound.  Can round down, as it is an equivalent to the Hamming bound in that it is an upper bound.
"""
# promote()
# _T = typeof(T)

singleton_bound(q::T, n::T, d::T) where T <: Integer =
	floor(T, float(big(q))^(big(n) - big(d) + 1))
singleton_bound(q::T, n::T, d::T, ::Rounding) where T <: Integer =
	float(big(q))^(big(n) - big(d) + 1)

gilbert_varshamov_bound(q::T, n::T, d::T) where T <: Integer =
	__sphere_bound(a -> floor(T, a), q, n, d - 1)
gilbert_varshamov_bound(q::T, n::T, d::T, ::Rounding) where T <: Integer =
	__sphere_bound(identity, q, n, d - 1)

function __plotkin_bound_core(round_func::Function, q::T, n::T, d::T) where T <: Integer
	if ! isequal(q, 2)
		throw(error("The Plotkin bound only works for the binary code."))
	end

	if iseven(d) && 2d > n
		return round_func((d) / (2d + 1 - n))
	elseif isodd(d) && 2d + 1 > n
		return round_func((d + 1) / (2d + 1 - n))
	elseif iseven(d)
		return T(4d)
		# return A_2(2d, d) ≤ 4d
	elseif isodd(d)
		return T(4d + 4)
		# return A_2(2d + 1, d) ≤ 4d + 4
	end
end

plotkin_bound(q::T, n::T, d::T) where T <: Integer =
	__plotkin_bound_core(a -> floor(T, a), q, n, d)
plotkin_bound(q::T, n::T, d::T, ::Rounding) where T <: Integer =
	__plotkin_bound_core(identity, q, n, d, no_round)

elias_bassalygo_bound(q::T, n::T, d::T) where T <: Integer =
elias_bassalygo_bound(q::T, n::T, d::T, ::Rounding) where T <: Integer =

function __johnson_bound_core(round_func::Function, q::T, n::T, d::T) where T <: Integer
	if isinteger((d - 1) / 2) # is odd
		t = T((d - 1) / 2)
		__sphere_bound(round_func, q, n, t) # if d = 2t + 1
	elseif isinteger(d / 2)
		t = T(d / 2)
		__sphere_bound(round_func, q, n, t)
	end
end

@doc raw"""
```julia
construct_ham_matrix(r::Int, q::Int) -> Matrix
```

Construct a Hamming parity-check matrix.

Parameters:
  - `r::Int`: number of rows of a parity check matrix.
  - `q:::Int`: The size of the alphabet of the code.

Returns:
  - `Matrix`: The Hamming matrix, denoted as ``\text{Ham}(r, q)``

---

### Examples

```julia
julia> construct_ham_matrix(3,2)
3×7 Array{Int64,2}:
 0  0  0  1  1  1  1
 0  1  1  0  0  1  1
 1  0  1  0  1  0  1
```
"""
function construct_ham_matrix(r::Int, q::Int)
    ncols = Int(floor((q^r - 1) / (q - 1)))
    M = Matrix{Int}(undef, r, ncols)

    for i in 1:ncols
        M[:, i] = reverse(digits(parse(Int, string(i, base = q)), pad = r), dims = 1)
    end

    return M
end

"""
```julia
isperfect(n::Int, k::Int, d::Int, q::Int) -> Bool
```

Checks if a code is perfect.  That is, checks if the number of words in the code is exactly the "Hamming bound", or the "Sphere Packing Bound".

Parameters:
  - `q:::Int`: The size of the alphabet of the code.
  - `n::Int`: The length of the words in the code (block length).
  - `d::Int`: The distance of the code.
  - `k::Int`: The dimension of the code.

Returns:
  - `Bool`: true or false

---

### Examples

```julia
julia> isperfect(11, 6, 5, 3)
true
```
"""
function isperfect(n::T, k::T, d::T, q::T) where T <: Int
	isprimepower(q) || throw(error("Cannot check if the code is perfect with q not a prime power."))
    M = q^k

    isequal(sphere_packing_bound(q, n, d), M) && return true
	return false
end

"""
```julia
ishammingbound(r::Int, q::Int) -> Bool
```

Checks if the code is a perfect code that is of the form of a generalised Hamming code.

Parameters:
  - `r::Int`: number of rows of a parity check matrix.
  - `q::Int`: The size of the alphabet of the code.

Returns:
  - `Bool`: true or false
"""
function ishammingperfect(r::Int, q::Int)
    n = 2^r - 1
    k = n - r
    M = q^k
    d = size(construct_ham_matrix(r, q), 1) # the number of rows of the hamming matrix (which is, by design, linearly independent)
    d = r
    # r is dim of dueal code; dim of code itself is block length minus r
    # println(n)
    # println((q^r - 1) / (q - 1))

    isequal(n, (q^r - 1) / (q - 1)) && \
		isequal(d, 3) && \
		isequal(M, q^(((q^r - 1) / (q - 1)) - r)) && \
        return true
    return false
end

"""
```julia
ishammingperfect(n::Int, k::Int, d::Int, q::Int) -> Bool
ishammingperfect(q::Int, n::Int, d::Int) -> Bool
```

Checks if the code is a perfect code that is of the form of a generalised Hamming code.

Parameters:
  - `q:::Int`: The size of the alphabet of the code.
  - `n::Int`: The length of the words in the code (block length).
  - `d::Int`: The distance of the code.
  - `k::Int`: The dimension of the code.

Returns:
  - `Bool`: true or false

---

### Examples

```julia
julia> isgolayperfect(11, 6, 5, 3) # this is one of golay's perfect codes
true
```
"""
function ishammingperfect(n::T, k::T, d::T, q::T) where T <: Int
    isprimepower(q) || return false

    M = q^k
    r = log(ℯ, ((n * log(ℯ, 1)) / (log(ℯ, 2))) + 1) / log(ℯ, 2)

    if isequal(n, (q^(r - 1)) / (q - 1)) && isequal(d, 3) && isequal(M, q^(((q^r - 1) / (q - 1)) - r))
        return true
    end

    return false
end
function ishammingperfect(q::Int, n::Int, d::Int)
	isprimepower(q) || return false # we are working in finite fields, so q must be a prime power
	d ≠ 3 && return false

	r = 1
	while ((q^r - 1) / (q - 1)) < n
		r = r + 1
	end

	return ifelse(isequal(((q^r - 1) / (q - 1)), n), true, false)
end

"""
```julia
isgolayperfect(n::Int, k::Int, d::Int, q::Int) -> Bool
```

Golay found two perfect codes.  `isgolayperfect` checks if a code of block length n, distance d, alphabet size q, and dimension k, is a perfect code as described by Golay.

Parameters:
  - `n::Int`: The block length of words in the code (e.g., word "abc" has block length 3).
  - `k::Int`: The dimension of the code.
  - `d::Int`: The distance of the code (i.e., the minimum distance between codewords in the code).
  - `q::Int`: An Int that is a prime power.  The modulus of the finite field.

Returns:
  - `Bool`: true or false.

---

### Examples

```julia
julia> isgolayperfect(11, 6, 5, 3) # this is one of golay's perfect codes
true
```
"""
function isgolayperfect(n::T, k::T, d::T, q::T) where T <: Int
	isprimepower(q) ||  false # we are working in finite fields, so q must be a prime power
    M = q^k
    (isequal(q, 2) && isequal(n, 23) && isequal(d, 7) && isequal(M, 2^12)) && return true
    (isequal(q, 3) && isequal(n, 11) && isequal(d, 5) && isequal(M, 3^6)) && return true
    return false
end


text:
using POMDPs, StatsBase, POMDPToolbox, SARSOP, QMDP, AutomotiveDrivingModels, JLD, Parameters, Distributions

include("../../admenv/crosswalk/occluded_crosswalk_env.jl")
include("../../admenv/crosswalk/pomdp_types.jl")
include("../../admenv/crosswalk/spaces.jl")
include("../../admenv/crosswalk/transition.jl")
include("../../admenv/crosswalk/observation.jl")
include("../../admenv/crosswalk/belief.jl")
include("../../admenv/crosswalk/adm_helpers.jl")
include("../../admenv/crosswalk/decomposition.jl")

include("config.jl")
include("ego_control.jl")
include("sensor.jl")
include("policy.jl")
include("state_estimation.jl")
include("constant_pedestrian.jl")
include("simulation.jl")
include("baseline_policy.jl")

config = EvalConfig(time_out = 200,
                    n_ped=2,
                    n_episodes = 3)

params = EnvParams(ped_rate = 0.9)
env = CrosswalkEnv(params)

pomdp = OCPOMDP()
pomdp.collision_cost = float(ARGS[1])

solver = SARSOPSolver(randomization=true, precision = 0.001, timeout = 3600.)

policy = solve(solver, pomdp,
               SARSOP.create_policy(solver, pomdp, "crosswalk$(abs(pomdp.collision_cost)).policy"),
               "model$(round(abs(pomdp.collision_cost))).pomdpx" )

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
export MSet

mutable struct MSet{T} <: AbstractSet{T}
  dict::Dict{T, Int}
  MSet{T}() where {T} = new(Dict{T,Int}())
  MSet{T}(itr) where {T} = union!(new(Dict{T,Int}()), itr)
end

MSet() = MSet{Any}()
MSet(itr) = MSet{eltype(itr)}(itr)

Base.eltype(::Type{MSet{T}}) where {T} = T
Base.similar(s::MSet{T}) where {T} = MSet{T}()
Base.similar(s::MSet, T::Type) = MSet{T}()

#TODO: compact printing, remove trailing , ... the works...
function Base.show(io::IO, s::MSet)
    print(io,"MSet")
    if isempty(s)
        print(io,"{",eltype(s),"}()")
        return
    end
    print(io,"(")
    for (k,v) = s.dict
      if v > 1
        print(io, "$k : $v, ") # ugly
      else
        print(io, "$k, ") # ugly
      end
    end
    print(io,")")
end

Base.isempty(s::MSet) = isempty(s.dict)
Base.length(s::MSet)  = BigInt(sum(values(s.dict)))
Base.IteratorSize(s::MSet) = Base.HasLength()
Base.IteratorEltype(s::MSet) = Base.HasEltype()
Base.eltype(s::MSet{T}) where {T} = T
Base.in(x, s::MSet) = haskey(s.dict, x)
function Base.push!(s::MSet, x, mult::Int = 1)
  add_to_key!(s.dict, x, mult)
  #if haskey(s.dict, x)
  #  s.dict[x] += mult
  #else
  #  s.dict[x] = mult
  #end
end

function Base.pop!(s::MSet, x)
  s.dict[x] -= 1
  if s.dict[x] == 0
    delete!(s.dict, x)
  end
  return x
end
Base.pop!(s::MSet, x, deflt) = x in s ? pop!(s, x) : deflt
Base.pop!(s::MSet) = (val = iterate(s.dict)[1][1]; pop!(s, val))
Base.delete!(s::MSet, x) = (delete!(s.dict, x); s)

Base.copy(s::MSet) = union!(similar(s), s)

function Base.iterate(s::MSet)
  I = iterate(s.dict)
  I === nothing && return I
  return I[1][1], (I[1], I[2], 1)
end

function Base.iterate(s::MSet, state)
  if state[3] < state[1][2]
    return state[1][1], (state[1], state[2], state[3]+1)
  else
    I = iterate(s.dict, state[2])
    I === nothing && return I
    val, st = I
    return (val[1], (val, st, 1))
  end
end

Base.union(s::MSet) = copy(s)
function Base.union(s::MSet, sets::Set...)
    u = MSet{join_eltype(s, sets...)}()
    union!(u,s)
    for t in sets
        union!(u,t)
    end
    return u
end
Base.union!(s::MSet, xs) = (for x=xs; push!(s,x); end; s)
Base.union!(s::MSet, xs::AbstractArray) = (for x=xs; push!(s,x); end; s)

function multiplicities(s::MSet)
  return values(s.dict)
end

function Base.unique(s::MSet)
  return collect(keys(s.dict))
end

############################################
# subsets iterator
############################################

struct MSubSetItr{T}
  b::Array{T, 1}
  m::Array{Int, 1}
  length::Int
end

function subsets(s::MSet{T}) where T
  # subsets are represented by integers in a weird base
  # the distinct elements are b0...bn with mult mi
  # subset (bi, ni) -> sum ni gi where gi = prod (mj+1)
  b = collect(unique(s))
  m = [s.dict[x] for x = b]
  #= not needed for the iterator
  g = [1]
  for i=2:length(b)
    push!(g, g[end]*(m[i]+1))
  end
  =#
  return MSubSetItr{T}(b, m, length(m) == 0 ? 1 : prod(x+1 for x=m))
end

function int_to_elt(M::MSubSetItr{T}, i::Int) where T
  s = MSet{T}()
  for j=1:length(M.b)
    k = i % (M.m[j]+1)
    for l=1:k
      push!(s, M.b[j])
    end
    i = div(i-k, M.m[j]+1)
  end
  return s
end

function Base.iterate(M::MSubSetItr)
  return int_to_elt(M, 0), 0
end

function Base.iterate(M::MSubSetItr, st::Int)
  st += 1
  st >= M.length && return nothing
  return int_to_elt(M, st), st
end

function Base.length(M::MSubSetItr)
  return M.length
end

Base.IteratorSize(M::MSubSetItr) = Base.HasLength()
Base.IteratorEltype(M::MSubSetItr) = Base.HasEltype()
Base.eltype(M::MSubSetItr{T}) where {T} = MSet{T}

function Base.show(io::IO, M::MSubSetItr)
  println(io, "subset iterator of length $(M.length) for $(M.b) with multiplicities $(M.m)")
end

#... to be completed from base/Set.jl ...

#subsets for Set
struct SubSetItr{T}
  b::Array{T, 1}
  length::Int
end

function subsets(s::Set{T}) where T
  # subsets are represented by integers in base 2
  b = collect(unique(s))
  return SubSetItr{T}(b, 2^length(b))
end

function int_to_elt(M::SubSetItr{T}, i::Int) where T
  s = Set{T}()
  j = 1
  while i > 0
    if isodd(i)
      push!(s, M.b[j])
    end
    j += 1
    i = div(i, 2)
  end
  return s
end

function Base.iterate(M::SubSetItr)
  return int_to_elt(M, 0), 0
end

function Base.iterate(M::SubSetItr, st::Int)
  st += 1
  st >= M.length && return nothing
  return int_to_elt(M, st), st
end

function Base.length(M::SubSetItr)
  return M.length
end

Base.IteratorSize(M::SubSetItr) = Base.HasLength()
Base.IteratorEltype(M::SubSetItr) = Base.HasEltype()
Base.eltype(M::SubSetItr{T}) where {T} = Set{T}

function Base.show(io::IO, M::SubSetItr)
  println(io, "subset iterator of length $(M.length) for $(M.b)")
end

#only subsets of a given size

struct SubSetSizeItr{T}
  b::Array{T, 1}
  k::Int #subsets of size k only
  B::Array{Array{Int, 1}, 1}
  length::Int
end

function subsets(s::Set{T}, k::Int) where T
  # subsets are represented by integers in the Combinatorial_number_system
  # https://en.wikipedia.org/wiki/Combinatorial_number_system
  # this iterator could indexed, the other one below not
  # maybe use "The coolest way to generate combinations" instead
  b = collect(unique(s))
  m = Int(binom(length(b), k))
  C = Array{Array{Int, 1}, 1}()
  while k >= 1
    B = Int[]
    i = k-1
    while true
      c = Int(binom(i, k))
      if c < m && length(B) < length(b)
        push!(B, c)
        i += 1
      else
        break
      end
    end
    push!(C, B)
    k -=1
  end

  return SubSetSizeItr{T}(b, length(C), C, m)
end

function int_to_elt(M::SubSetSizeItr{T}, i::Int) where T
  s = Set{T}()
  j = 1
  while j <= length(M.B)
    z = findlast(x -> x <= i, M.B[j])
    push!(s, M.b[z + M.k - j])
    i -= M.B[j][z]
    j += 1
  end
  while length(s) < M.k
    push!(s, M.b[length(s)])
  end

  return s
end

function Base.iterate(M::SubSetSizeItr)
  return int_to_elt(M, 0), 0
end

function Base.iterate(M::SubSetSizeItr, st::Int)
  st += 1
  st >= M.length && return nothing
  return int_to_elt(M, st), st
end

function Base.length(M::SubSetSizeItr)
  return M.length
end

Base.IteratorSize(M::SubSetSizeItr) = Base.HasLength()
Base.IteratorEltype(M::SubSetSizeItr) = Base.HasEltype()
Base.eltype(M::SubSetSizeItr{T}) where {T} = Set{T}

function Base.show(io::IO, M::SubSetSizeItr)
  println(io, "subset iterator of length $(M.length) for $(M.b) and subsets of size $(M.k)")
end


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

function Base.read(io::IO, ::Type{SMPSModel})
    raw = read(io, RawSMPS)
    return SMPSModel(raw)
end
"""
    read(io::IO,
         ::Type{SMPSSampler})

Return an `StochasticModel` from the model definition read from `io` in SMPS format.
"""
function Base.read(io::IO, ::Type{StochasticModel})
    smps::SMPSModel{2} = read(io, SMPSModel)
    return stochastic_model(smps)
end
"""
    read(io::IO,
         ::Type{SMPSSampler})

Return an `SMPSSampler` capable of sampling `SMPSScenario` using the model definition read from `io` in SMPS format.
"""
function Base.read(io::IO, ::Type{SMPSSampler})
    smps::SMPSModel{2} = read(io, SMPSModel)
    return SMPSSampler(smps.raw.sto, smps.stages[2])
end
"""
    read(io::IO,
         ::Type{StochasticProgram};
         num_scenarios::Union{Nothing, Integer} = nothing,
         instantiation::StochasticInstantiation = UnspecifiedInstantiation(),
         optimizer = nothing;
         defer::Bool = false,
         kw...)

Instantiate a two-stage stochastic program using the model definition read from `io` in SMPS format, of size `num_scenarios`. If `num_scenarios = nothing`, instantiate using the full support. Optionally, supply an `optimizer`. If no explicit `instantiation` is provided, the structure is induced by the optimizer. The structure is `Deterministic` by default.
"""
function Base.read(io::IO,
                   ::Type{StochasticProgram};
                   num_scenarios::Union{Nothing, Integer} = nothing,
                   instantiation::StochasticInstantiation = StochasticPrograms.UnspecifiedInstantiation(),
                   optimizer = nothing,
                   defer::Bool = false,
                   direct_model::Bool = false,
                   kw...)
    smps::SMPSModel{2} = read(io, SMPSModel)
    sm = stochastic_model(smps)
    sampler = SMPSSampler(smps.raw.sto, smps.stages[2])
    if num_scenarios != nothing
        return instantiate(sm, sampler, num_scenarios; instantiation = instantiation, optimizer = optimizer, defer = defer, direct_model = direct_model, kw...)
    else
        return instantiate(sm, full_support(sampler); instantiation = instantiation, optimizer = optimizer, defer = defer, direct_model = direct_model, kw...)
    end
end
Base.read(filename::AbstractString, ::Type{StochasticProgram}; kw...) = open(io -> read(io, StochasticProgram; kw...), filename)


text:
using jInv.Mesh;
using MAT;
using jInv
using jInv.InverseSolve
using ShapeReconstructionPaLS
using ParamLevelSet;
using ShapeReconstructionPaLS.Utils;
using ShapeReconstructionPaLS.PointCloud;
using DelimitedFiles
using LinearAlgebra
using Statistics
using Distributed
using SparseArrays
using Random
ENV["MPLBACKEND"] = "Qt4Agg"

function ind2subv(shape, indices)
    CI = CartesianIndices(shape)
    return getindex.(Ref(CI), indices)
end

plotting = false;

if plotting
	using PyPlot;
	using jInvVis;
	close("all");
end

invertDip = true;
invertVis = false;

locateRBFwithGrads = true;

################################################
### Reading the model and pad ##################
################################################
#model = "skull_160"; ## skull
 model = "hand_160";
file = matopen(string(pwd(),"/models/",model,".mat"));

A = read(file, "B");
close(file);
A = convert(Array{Float32,3},A[1:4:end,1:4:end,1:4:end]);
n_inner = collect(size(A));
pad = ceil.(Int64,0.1*n_inner[1]);
n = n_inner .+ 2*pad;

n_tup = tuple(n...)
volume = prod(n);
n1 = n_tup[1];
n2 = n_tup[2];
n3 = n_tup[3];

m = zeros(Float64,n_tup);
mask = zeros(Int8,n_tup);
m[pad+1:end-pad,pad+1:end-pad,pad+1:end-pad] = A; A = 0;
mask[pad+1:end-pad,pad+1:end-pad,pad+1:end-pad] .= 1;
if plotting
	plotModel(m);
	savefig("TrueModel.png");
end
origm = copy(m);
Mesh = getRegularMesh([0.0 5.0 0.0 5.0 0.0 5.0],n);

noiseAnglesDeg = 0.0;
#original was mean(2*Mesh.h);
noiseTrans = mean(0.0*Mesh.h);
noiseAngles = deg2rad(noiseAnglesDeg);
#########################################################################################################
### Dipping Data ########################################################################################
#########################################################################################################
samplingBinning = 1;

Mesh = getRegularMesh([0.0 5.0 0.0 5.0 0.0 5.0],n);
pad = div(pad,samplingBinning);

method = RBF10BasedSimple; methodName="RBF10BasedSimple";

println("Inv mesh:")
println(Mesh.n)

n_tup = tuple(Mesh.n...);
nWorkers = nworkers();
println("nworkers=",nWorkers)
misfun = PCFun; ## least squares

	dipDataFilename = string("pointCloudData",model,"_dataRes",n,1);
	b = [0.0 0 0] #; 0.0 0 0 ]
	theta_phi_dip = [0.0 0.0] #; 0.0 0.0];
	noiseTrans = 0.0*prod(Mesh.h);
	noiseAngles = deg2rad(0.0);
	npc = 1;
	prepareSyntheticPointCloudData(m,Mesh,npc,theta_phi_dip,b,noiseTrans,noiseAngles,dipDataFilename);

	#########################################################################################################
	### Read the dipping data ###############################################################################
	#########################################################################################################
	n,domain,Data,P,Normals,Margin,b,theta_phi_dip,npc = readPointCloudDataFile(dipDataFilename);
	println("initial obtained translation:",b);
        P = Array{Array{Int64,1}}(P);
	writedlm(string("simPCpart2",".txt"),convert(Array{Int64},P[1]));
#       writedlm(string("simPCpart1",".txt"),convert(Array{Int64},P[2]));
#	quit()
	### Set up the dip inversion
	println("type of p:",typeof(P));
	pForDip = getPointCloudParam(Mesh,P,Normals,Margin,theta_phi_dip,b,samplingBinning,nWorkers,method);
	ndips = npc;
	if(locateRBFwithGrads)
		pForDip_MATFree = getPointCloudParam(Mesh,P,Normals,Margin,theta_phi_dip,b,samplingBinning,nWorkers,MATFree);
	end

	### Create Dip param if we use the GRAD of MATFREE to locate new RBFs
	dobsDirect = divideDataToWorkersData(nWorkers,Data);
	println("size of dobs",size(dobsDirect))
	#nnz = filter(x -> x != 0,dobsDirect[1]);
	#println("length of dobsdirect: ",length(dobsDirect[1]),"_",length(nnz));
	Wd_Dip = Array{Array{Float32}}(undef,nWorkers);
	for k=1:nWorkers
		Wd_Dip[k] = ones(Float32,size(dobsDirect[k]));
	end
	println("Wd dip:",size(Wd_Dip))

n_Moves_all = npc;
println("npc = ",npc);

isSimple = 0;
isRBF10 = 0;

###############################################################################################
############## USE RBF dictionary #############################################################
###############################################################################################

	nRBF = 20;

	isRBF10 = (method == RBF10BasedSimple1 || method==RBF10BasedSimple2 || method==RBF10Based);
	isSimple = (method==RBFBasedSimple1 || method==RBF10BasedSimple1 || method==RBFBasedSimple2 || method==RBF10BasedSimple2);

	if isRBF10
		numParamOfRBF = 10;
	else
		numParamOfRBF = 5;
	end
	isSimple = true;
	n_m_simple = numParamOfRBF*nRBF;
	nAll = isSimple ? numParamOfRBF*nRBF : numParamOfRBF*nRBF + n_Moves_all*5;

	modfun = identityMod;

	if isSimple
		modfunForPlotting = (m)->ParamLevelSetModelFunc(Mesh,m;computeJacobian = 0,
			sigma = getDefaultHeavySide() ,bf = 1, numParamOfRBF = numParamOfRBF);
	else
		modfunForPlotting = (m)->ParamLevelSetModelFunc(Mesh,m[1:(numParamOfRBF*nRBF)];computeJacobian = 0,
			sigma = getDefaultHeavySide() ,bf = 1, numParamOfRBF = numParamOfRBF);
	end

	Iact = 1.0;
	IactPlot = 1.0;
	sback = zero(Float32);
	a = -100000000.0;
	b = 10000000.0;
	boundsHigh = zeros(Float32,nAll) .+ b;
	boundsLow = zeros(Float32,nAll) .+ a;
	# m0 - initial guess. mref - reference for regularization.
	mref = zeros(Float32,nAll);
	m0   = zeros(Float64,n_m_simple);
	mCenters = zeros(Float32,n_tup);
	mref[1:numParamOfRBF:n_m_simple] .= ParamLevelSet.centerHeavySide*2.0*rand(nRBF);
	if isRBF10
		# here we initialize B as an identity matrix
		mref[2:10:n_m_simple] .= 0.5;
		mref[5:10:n_m_simple] .= 0.5;
		mref[7:10:n_m_simple] .= 0.5;
		boundsLow[2:10:n_m_simple] .= 0.05;
		boundsLow[5:10:n_m_simple] .= 0.05;
		boundsLow[7:10:n_m_simple] .= 0.05;
	else
		mref[2:5:n_m_simple] .= 1.0;
		boundsLow[2:5:n_m_simple] .= 0.05;
	end
	mref[(numParamOfRBF-2):numParamOfRBF:n_m_simple] .= (Mesh.domain[1] + Mesh.domain[2])/2;
	mref[(numParamOfRBF-1):numParamOfRBF:n_m_simple] .= (Mesh.domain[3] + Mesh.domain[4])/2;
	mref[numParamOfRBF:numParamOfRBF:n_m_simple] .= (Mesh.domain[5] + Mesh.domain[6])/2;

	m0[:] = mref[1:n_m_simple];
	if !isSimple
		m0 = wrapRBFparamAndRotationsTranslations(m0,theta_phi_dip,zeros(size(theta_phi_dip,1),3));
	end

	mref[:] = m0;
	m0[(numParamOfRBF-2):numParamOfRBF:n_m_simple] .+= 0.2*(Mesh.domain[2] - Mesh.domain[1])*randn(nRBF);
	m0[(numParamOfRBF-1):numParamOfRBF:n_m_simple] .+= 0.2*(Mesh.domain[4] - Mesh.domain[3])*randn(nRBF);
	m0[numParamOfRBF:numParamOfRBF:n_m_simple]     .+= 0.2*(Mesh.domain[6] - Mesh.domain[5])*randn(nRBF);

	alpha = 5e-0;

	II = (sparse(1.0I, nAll,nAll));
	regfun = (m, mref, M)->TikhonovReg(m,mref,M,II);
	HesPrec = getEmptyRegularizationPreconditioner();

	if isRBF10
		spd_reg = (m,mref,M) -> RBF_SPD_regularization(m,mref,nRBF);
		regfun 	 		= [regfun;spd_reg];
		alpha 		 	= [alpha;1e-1*alpha]; #was 1e-4 for cube\prism !!!!
		mref 			= [mref zeros(size(mref))];
	end

pMisRFs = getMisfitParam(pForDip, Wd_Dip, dobsDirect, misfun, Iact,sback);

locateRBFwithGrads = true;
if(locateRBFwithGrads)
		mask=0;
		Iact_free = 1.0;
		#### USE A BOUND MODEL or bounds
		a_free = 0.0;
		b_free= 1.0;

		pMis_Free = getMisfitParam(pForDip_MATFree, Wd_Dip, dobsDirect, misfun, Iact_free,sback);
end

cgit  = 50; # not used.
maxit = 10;
pcgTol = 1e-2; #not used

@everywhere function updateAngles(pFor::Array{RemoteChannel},theta_phi)
@sync begin
	@async begin
		for k=1:length(pFor)
			pFor[k] = remotecall_fetch(updateAngles,pFor[k].where,pFor[k],theta_phi);
		end
	end
end
return pFor;
end

@everywhere function updateAngles(pFor::RemoteChannel,theta_phi)
	pForMATFREE  = take!(pFor)
	worker_theta_phi = theta_phi[pForMATFREE.workerSubIdxs,:]
	pForMATFREE.theta_phi_rad=copy(worker_theta_phi);
	put!(pFor,pForMATFREE)
	return pFor;
end

function myDump(mc::Vector,Dc,iter,pInv,pMis,method="",ndips=0,noiseAngle=0,noiseTranslation=0,useVisual=false,plotFun=modfunForPlotting)
		ntup = tuple((pInv.MInv.n)...);
		fullMc = plotFun(mc)[1];
		fullMc = reshape(IactPlot*fullMc[:],ntup);
		writedlm(string("twoHand1PCParams_iter_",iter,".dat"),convert(Array{Float64},mc));
		writedlm(string("twoHand1PCVolume_iter_",iter,"_",ntup,".dat"),convert(Array{Float16},fullMc));
end

pInv = getInverseParam(Mesh,modfun,regfun,alpha,mref,boundsLow,boundsHigh,
                     maxStep = 2e-1,pcgMaxIter=cgit,pcgTol=pcgTol,
					 minUpdate = 0.01, maxIter = maxit,HesPrec=HesPrec);

#### Projected Gauss Newton
mc = m0;
pInv.maxIter = 100;
Dc = nothing;
mc,Dc,flag,his = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);

pInv.maxIter = 10;
myDump(mc,0,1,pInv,0,"",0,noiseAnglesDeg,noiseTrans,invertVis);

nRBF=nRBF;
times=0;

if !(method == MATBased || method == MATFree)
	outerIter = 300;
	grad_u = nothing;
	new_RBF_location = [];
	for iterNum=2:outerIter
		println("Inversion using ", nRBF," basis functions")
		new_nRBF =10;
		global times = times + 1;

		# Dc - the data using our method.
		global Dc = Dc;
		u = modfunForPlotting(mc)[1];
		#I = find((u.>0.3) .& (u.<0.6));
		#Ii = P[randperm(length(P))[1:new_nRBF]];
		#Ii = findall(x -> x >0.3 && x<0.7,u);
		#Lets compute the gradient of the misfit function:
		locateRBFwithGrads = true;
		if(locateRBFwithGrads)
			#(m11,theta_phi1,b1) 	= splitRBFparamAndRotationsTranslations(mc,nRBF,ndips,numParamOfRBF)
			#updateAngles(pForDip_MATFree,theta_phi1);
			pMis_Free = getMisfitParam(pForDip_MATFree, Wd_Dip, dobsDirect, misfun, Iact_free,sback);
			computeMisfit(u,pMis_Free);
			gradMis = computeGradMisfit(u,Dc,pMis_Free)
			gradMis = gradMis.^2;
			println("size gradmis:",size(gradMis))
			amax = argmax(gradMis); vmax = maximum(gradMis);
			sp = sortperm(gradMis,rev = true);
			Ii = sp[1:new_nRBF];
			writedlm(string("gradMisfit",".txt"),convert(Array{Float64},gradMis[:]))
		end

		if length(Ii) > new_nRBF
		Ii = Ii[randperm(length(Ii))[1:new_nRBF]];
		end

   println("Ii size:",size(Ii))
    while(Ii[1] in new_RBF_location )#|| Ii.-1 in new_RBF_location || Ii.+1 in new_RBF_location)
      #Ii = findall(x -> x .>0.3 && x.<0.7,u);
	  #Ii = P[randperm(length(P))[1:new_nRBF]];
	  #Ii = Ii[randperm(length(Ii))[1:new_nRBF]];
       tmp = sp[1:outerIter*new_nRBF];
		Ii = tmp[randperm(outerIter*new_nRBF)[1:new_nRBF]];

    end
		global new_RBF_location = append!(new_RBF_location,Ii);
		indices = getindex.(ind2subv(n_tup,Ii),[1 2 3 ]);
		I1 = indices[:,1]; I2= indices[:,2]; I3 = indices[:,3];

		X1 = I1*Mesh.h[1] .- 0.5*Mesh.h[1]; .+ Mesh.domain[1];
		X2 = I2*Mesh.h[2] .- 0.5*Mesh.h[2]; .+ Mesh.domain[3];
		X3 = I3*Mesh.h[3] .- 0.5*Mesh.h[3]; .+ Mesh.domain[5];

		println("X1:",X1);println("X2:",X2);println("X3:",X3);

		tmp  = mc;
		n_mc_old = numParamOfRBF*nRBF;
		n_mc_new = numParamOfRBF*nRBF + length(Ii)*numParamOfRBF;
		mc_new = zeros(length(mc) + length(Ii)*numParamOfRBF);
		if isRBF10
			mref = pInv.mref[:,1];
		else
			mref = pInv.mref;
		end
		mref_new = zeros(Float32,length(mc) + length(Ii)*numParamOfRBF);

		mc_new[1:n_mc_old] = mc[1:n_mc_old];

		mc_new[(n_mc_new+1):end] = mc[(n_mc_old+1):end];
		mref_new[(n_mc_new+1):end] = mref[(n_mc_old+1):end];

		mc_new[(n_mc_old+1):numParamOfRBF:n_mc_new] .= 0.0; #ParamLevelSet.centerHeavySide - ParamLevelSet.deltaHeavySide;
		mc_new[(n_mc_old+numParamOfRBF-2):numParamOfRBF:n_mc_new] .= X1;
		mc_new[(n_mc_old+numParamOfRBF-1):numParamOfRBF:n_mc_new] .= X2;
		mc_new[(n_mc_old+numParamOfRBF):numParamOfRBF:n_mc_new] .= X3;

		pInv.boundsHigh = zeros(Float32,length(mc_new)) .+ b;
		boundsLow = zeros(Float32,length(mc_new)) .+ a;
		boundsLow[1:n_mc_old] = pInv.boundsLow[1:n_mc_old];
		if isRBF10
			mc_new[(n_mc_old+2):10:n_mc_new] .= 3.0;
			mc_new[(n_mc_old+5):10:n_mc_new] .= 3.0;
			mc_new[(n_mc_old+7):10:n_mc_new] .= 3.0;
			boundsLow[(n_mc_old+2):10:n_mc_new] .= 0.05;
			boundsLow[(n_mc_old+5):10:n_mc_new] .= 0.05;
			boundsLow[(n_mc_old+7):10:n_mc_new] .= 0.05;
		else
			mc_new[(n_mc_old+2):5:n_mc_new] .= 1.0;
			boundsLow[(n_mc_old+2):5:n_mc_new] .= 0.05;
		end
		pInv.boundsLow = boundsLow;
		mref_new[1:n_mc_new] = mc_new[1:n_mc_new];

		global mc = mc_new;
		#mref_new[end-9:end] .= 0.0;
		#II = speye(Float32,length(mc_new));
		len = Int64(length(mc_new));
		IIs = sparse(1.0I,len,len);
		pInv.regularizer = (m, mref, M)->TikhonovReg(m,mref,M,IIs);
		global nRBF += new_nRBF;
		if isRBF10
			spd_reg = (m,mref,M) -> RBF_SPD_regularization(m,mref,nRBF);
			pInv.regularizer  = [pInv.regularizer;spd_reg];
			pInv.mref 			  = [mref_new zeros(length(mref_new))];
			#pInv.alpha .*= 0.8;
			pInv.alpha[1] *= 0.8;
			pInv.alpha[2] *= 0.8;
		else
			pInv.mref = mref_new;
			pInv.alpha *= 0.8;
		end
		#mc[end-9:end] .= 0.0;
		mc,Dc,flag,his = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);
		myDump(mc,0,iterNum,pInv,0,methodName,0,noiseAnglesDeg,noiseTrans,invertVis);

		if iterNum==outerIter
			pInv.alpha *= 0.1;
			pInv.maxIter = 200;
			mc, = projGN(mc,pInv,pMisRFs,solveGN=projGNexplicit);
			myDump(mc,0,iterNum+1,pInv,0,methodName,0,noiseAnglesDeg,noiseTrans,invertVis);
		end

	end

		meshForPlotting = getRegularMesh(Mesh.domain,Mesh.n*1);
		ntup = tuple((meshForPlotting.n)...);
		if isSimple
			modfunForPlotting = (m)->ParamLevelSetModelFunc(meshForPlotting,m;computeJacobian = 0,
				sigma = getDefaultHeavySide() ,bf = 1, numParamOfRBF = numParamOfRBF);
		else
			modfunForPlotting = (m)->ParamLevelSetModelFunc(meshForPlotting,m[1:(numParamOfRBF*nRBF)];computeJacobian = 0,
				sigma = getDefaultHeavySide() ,bf = 1, numParamOfRBF = numParamOfRBF);
		end
		fullMc = modfunForPlotting(mc)[1];
		writedlm(string("oneHand160_noise0_Volume_final_",ntup,"_from_",tuple((Mesh.n)...),".dat"),convert(Array{Float16},fullMc));
		#m = convert(Array{Float16},m);
		#writedlm(string("Fandisk200Direct_Volume_orig.dat"),m);
		#writedlm(string("Fandisk200Direct_RBF_centers.dat"),mCenters);

end

max_stars_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_issues_repo_name: BGUCompSci/ShapeReconstructionPaLS

max_forks_repo_name: BGUCompSci/ShapeReconstructionPaLS


text:
export PlotRay

#rayfil = "Case0_Bellhop_f10.ray"

function PlotRay(rayfil;
    xlabs = "Range, m",
    ylabs = "Height, m",
    xlim = (0,10_000),
    ylim = (-1000,0))

    # parse functions
    function parse_numbers(s)
        pieces = split(s, ' ', keepempty=false)
        map(pieces) do piece
            parse(Float64, piece)
        end
    end

    function parse_int(s)
        pieces = split(s, ' ', keepempty=false)
        map(pieces) do piece
            parse(Int, piece)
        end
    end

    io = open(rayfil,"r")
    lines = readlines(io)
    close(io)

    TITLE   =   lines[1]
    FREQ    =   parse_numbers(lines[2])[1]
    Nsxyz   =   parse_numbers(lines[3])
    NBeamAngles = parse_int(lines[4])

    DEPTHT = parse_numbers(lines[5])
    DEPTHB = parse_numbers(lines[6])

    Type1        = lines[7]

    Nsx    = Nsxyz[1]
    Nsy    = Nsxyz[2]
    Nsz    = Nsxyz[3]

    Nalpha = NBeamAngles[1]
    Nbeta  = NBeamAngles[2]

    # Extract title
    m = match(r"BELLHOP- .*? ", TITLE)
    TITLE = m.match

    #cmp(Type1, " 'rz'") == 0

    ### Extract ray data
    df = DataFrame(Nalpha = Float64[],
                    TopBnc = Int[],
                    BotBnc = Int[],
                    r = Float64[],
                    z = Float64[])

    st = 8
    function ray_extract!(df,lines,st)
        for ibeam = 1:Nalpha
            alpha0 = parse_numbers(lines[st])
            nsteps = parse_int(lines[st+1])[1]

            NumTopBnc = parse_int(lines[st+1])[2]
            NumBotBnc = parse_int(lines[st+1])[3]

            for i = st+2 : st+2+nsteps-1
                ray = parse_numbers(lines[i])
                push!(df, (alpha0[1], NumTopBnc[1],NumBotBnc[1], ray[1],ray[2]))
            end

            st += 2 + nsteps
        end
    end

    ray_extract!(df,lines,st)

    #df.z = 1000 .- df.z
    # group data
    gdf = groupby(df, :Nalpha)

    # Plot ray
    #gr()
    p = plot()
    xlabel!(xlabs)
    ylabel!(ylabs)
    xlims!(xlim[1],xlim[2])
    ylims!(ylim[1],ylim[2])
    yflip!(true)

    for i=1:length(gdf)
        if (gdf[i].TopBnc[1] == 1) & (gdf[i].BotBnc[1] == 1) # hit both boundaries
            plot!(p,gdf[i].r,gdf[i].z,
                legend = false,
                color = "#e41a1c",
                α = 0.5)
        elseif (gdf[i].BotBnc[1] == 1) # hit bottom
            plot!(p,gdf[i].r,gdf[i].z,
                legend = false,
                color = "#377eb8",
                α = 0.5)

        elseif (gdf[i].TopBnc[1] == 1)  # hit top
            plot!(p,gdf[i].r,gdf[i].z,
                legend = false,
                color = "#4daf4a",
                α = 0.5)
        else
            plot!(p,gdf[i].r,gdf[i].z,
                legend = false,
                color = "#984ea3",
                α = 0.5)
        end

    end
    current()

    #return df

end

#savefig(p,"plot.pdf")


text:
module OpenStreetMapXSim

using DataFrames
using StatsBase
using SparseArrays
using OpenStreetMapX
using Serialization
using Distributed
using Dates
using Nanocsv
using Random
using Printf

export start_location
export destination_location!
export select_route
export node_statistics, stats_aggregator!
export run_simulation
export reduce_results

include("types.jl")
include("start_location.jl")
include("destination_location.jl")
include("routing_module.jl")
include("statistics_aggregator.jl")
include("run.jl")
include("reduce.jl")

end


text:
"""
This docstring is attached to the [`Showcase`](@ref) module itself.

The [`EXPORTS`](@ref) abbreviation creates a bulleted list of all the exported names:

$(EXPORTS)

Similarly, the [`IMPORTS`](@ref) abbreviation lists all the imported modules:

$(IMPORTS)

The [`README`](@ref) can be used to include the `README.md` file in a docstring. The content
between the horizontal lines is spliced by the abbreviation:

---

$(README)

---

The [`LICENSE`](@ref) abbreviation can be used in the same way for the `LICENSE.md` file.
"""
module Showcase
using DocStringExtensions

"""
This docstring is attached to an [empty function
definitions](https://docs.julialang.org/en/v1/manual/methods/#Empty-generic-functions-1).
The [`METHODLIST`](@ref) abbreviation allows you to list all the methods though:

$(METHODLIST)
"""
function foo end

"""
This docstring is attached to a method that uses default values for some positional
arguments: `foo(x::Int, y=3)`.

As this effectively means that there are two different methods taking different numbers of
arguments, the [`SIGNATURES`](@ref) abbreviation produces the following result:

$(SIGNATURES)

---

The [`TYPEDSIGNATURES`](@ref) abbreviation can be used to also get the types of the
variables in the function signature:

$(TYPEDSIGNATURES)

---

The [`FUNCTIONNAME`](@ref) abbreviation can be used to directly include the name of the
function in the docstring (e.g. here: $(FUNCTIONNAME)). This can be useful when writing your
own type signatures:

    $(FUNCTIONNAME)(x, ...)
"""
foo(x::Int, y=3) = nothing

"""
A different method for [`$(FUNCTIONNAME)`](@ref). [`SIGNATURES`](@ref) abbreviation:

$(SIGNATURES)

And the [`TYPEDSIGNATURES`](@ref) abbreviation:

$(TYPEDSIGNATURES)
"""
foo(x::AbstractString) = nothing

## Methods with type parameters

"""
A method for [`$(FUNCTIONNAME)`](@ref), with type parameters. Original declaration:

```julia
bar(x::AbstractArray{T}, y::T) where {T <: Integer} = nothing
```

And the result from [`TYPEDSIGNATURES`](@ref) abbreviation:

$(TYPEDSIGNATURES)

For comparison, [`SIGNATURES`](@ref) abbreviation:

$(SIGNATURES)
"""
bar(x::AbstractArray{T}, y::T) where {T <: Integer} = nothing

"""
A method for [`$(FUNCTIONNAME)`](@ref), with type parameters. Original declaration:

```julia
bar(x::AbstractArray{T}, ::String) where {T <: Integer} = x
```

And the result from [`TYPEDSIGNATURES`](@ref) abbreviation:

$(TYPEDSIGNATURES)

For comparison, [`SIGNATURES`](@ref) abbreviation:

$(SIGNATURES)
"""
bar(x::AbstractArray{T}, ::String) where {T <: Integer} = x

"""
A method for [`$(FUNCTIONNAME)`](@ref), with type parameters. Original declaration:

```julia
bar(x::AbstractArray{T}, y::U) where {T <: Integer, U <: AbstractString} = 0
```

And the result from [`TYPEDSIGNATURES`](@ref) abbreviation:

$(TYPEDSIGNATURES)

For comparison, [`SIGNATURES`](@ref) abbreviation:

$(SIGNATURES)
"""
bar(x::AbstractArray{T}, y::U) where {T <: Integer, U <: AbstractString} = 0

"""
The [`TYPEDEF`](@ref) abbreviation includes the type signature:

$(TYPEDEF)

---

The [`FIELDS`](@ref) abbreviation creates a list of all the fields of the type.
If the fields has a docstring attached, that will also get included.

$(FIELDS)

---

[`TYPEDFIELDS`](@ref) also adds in types for the fields:

$(TYPEDFIELDS)
"""
struct Foobar{T <: AbstractString}
    "Docstring for the `x` field."
    x :: Nothing
    y :: T # y is missing a docstring
    "Docstring for the `z` field."
    z :: Vector{T}
end

export foo, Foobar

end


text:
using DelimitedFiles
using LinearAlgebra
using Clustering
using StatsBase
using Distributed
using DataFrames
using CSV
import ClusterManagers
using LinearAlgebra
using Distances
using Random
using Base
using ArgParse

function ParseCommands()

    s = ArgParseSettings()

    @add_arg_table s begin
        "--mutation_file", "-f"
            help = "File with the matrix count"
            required = true
        "--min_signatures_to_extract", "-s"
            help = "Minimum number of signatures to extract"
            arg_type = Int
            required = true
        "--max_signatures_to_extract", "-m"
            help = "Maximum number of signatures to extract"
            arg_type = Int
            required = true
        "--outpath", "-o"
            help = "Folder where the results will be saved"
            required = true
        "--iterations", "-i"
            help = "Total NMF iterations performed"
            arg_type = Int
            default = 1024
        "--workers", "-n"
            help = "Number of slurm parallel process or local cores"
            arg_type = Int
            default = 1
        "--weak", "-w"
            help = "Whether to remove weak mutation types. If unsure, use default"
            arg_type = Int
            default = 0
        "--slurm"
            help = "Whether to run the pipeline using SLURM"
            action = :store_true

    end
    return parse_args(s)
end

parsed_args = ParseCommands()
numberWorkers = parsed_args["workers"]
println("Adding workers...")
if parsed_args["slurm"]
    ClusterManagers.addprocs_slurm(numberWorkers)
else
    addprocs(numberWorkers)
end

# get the path where the script is located
@eval @everywhere path_root=@__DIR__

# add all the necessary scripts
include(string(path_root, "/src/removeWeak.jl"))
include(string(path_root,"/src/addWeak.jl"))
include(string(path_root,"/src/filterOutIterations.jl"))
include(string(path_root,"/src/evaluateStability.jl"))
include(string(path_root,"/src/deconvoluteMutCounts.jl"))
include(string(path_root,"/src/extractSignaturesParallelNoArray.jl"))
include(string(path_root,"/src/removeAllSingleSignatures.jl"))
include(string(path_root,"/src/removeOneSig.jl"))
include(string(path_root,"/src/evalSingleSample.jl"))
include(string(path_root,"/src/stabilities.jl"))

# mutation file
mutation_file = parsed_args["mutation_file"]

# number of signatures to extract
min_numberProcessesToExtract = parsed_args["min_signatures_to_extract"]
max_numberProcessesToExtract = parsed_args["max_signatures_to_extract"]

# total number of iterations
totalIterations = parsed_args["iterations"]

# outpath where they will be saved
outpath = parsed_args["outpath"]

# define whether is a Pan whole-genome, so weak types won't be removed
isWeak = parsed_args["weak"]
remove_weak_muttypes=0.001
if isWeak==0
    remove_weak_muttypes=0.00
end

# add seed so that the extraction is reproducible
Random.seed!(123456)

totalReplicates = 100
removeLastPercentage = 0.07

# Load input
genomesOriginal, sample_names = readdlm(mutation_file, '\t', Int, '\n'; header=true)
sample_header = vec(Symbol.(sample_names))

# Deconvolute mutation counts
genomes, sort_idx, n_types_to_remove = removeWeak(genomesOriginal, remove_weak_muttypes)

# start processing
name_ttype = splitext(basename(mutation_file))[1]

@info "Processing $name_ttype"
@info "Outpath $outpath"

if parsed_args["slurm"]
    @info "Slurm mode"
end

# Bootstrap variables
@eval @everywhere parallelGenomes=$genomes
@eval @everywhere originalGenomes=$genomesOriginal
@eval @everywhere parallelTotalIterations=$totalIterations
@eval @everywhere OutPath=$outpath

@eval @everywhere mutfile = $name_ttype
@eval @everywhere initial_mutfile = $mutation_file
@eval @everywhere SavePath = joinpath(OutPath, mutfile)
mkpath(SavePath)

@everywhere parallelGenomesColsSum = sum(parallelGenomes, dims=1)
@everywhere parallelNormGenomes = parallelGenomes ./ repeat(parallelGenomesColsSum, size(parallelGenomes, 1), 1 )

for numberProcessesToExtract = min_numberProcessesToExtract : max_numberProcessesToExtract

    @eval @everywhere parallelTotalProcesses=$numberProcessesToExtract

    @info "Extracting $numberProcessesToExtract mutational signatures for $totalIterations iterations using $numberWorkers cores"

    #signature extraction
    @time processes, exposures, processStab, processStabAvg, exposures_fitting, avgReconstructionError, avgReconstructionErrorPercentage, avgStability, reconError, reconErrorPercentage, simError = deconvoluteMutCounts(numberProcessesToExtract, totalIterations, sort_idx, n_types_to_remove)

    # save outfiles for each signature extracted
    file = string(SavePath, "/exposures_", parallelTotalProcesses)
    CSV.write(file, DataFrame(exposures), header = sample_header, delim = "\t")

    file = string(SavePath, "/processes_",  parallelTotalProcesses)
    CSV.write(file, DataFrame(processes), delim = "\t")

    file = string(SavePath, "/processesStabAvg_", parallelTotalProcesses)
    processes_stab = convert(DataFrame, processStabAvg)
    CSV.write(file, processes_stab, delim = "\t")

    file = string(SavePath, "/exposures_fitting_", parallelTotalProcesses)
    CSV.write(file, DataFrame(exposures_fitting), header = sample_header, delim = "\t")

    file = string(SavePath, "/reconError_", parallelTotalProcesses)
    CSV.write(file, DataFrame(reconError'), header = sample_header, delim="\t")

    file = string(SavePath, "/reconErrorPercentage_", parallelTotalProcesses)
    CSV.write(file, DataFrame(reconErrorPercentage'), header = sample_header, delim="\t")

    file = string(SavePath, "/simError_", parallelTotalProcesses)
    CSV.write(file, DataFrame(simError'), header = sample_header, delim="\t")

    file = string(SavePath, "/avgReconstructionError_", parallelTotalProcesses)
    writedlm(file, avgReconstructionError,  "\t")

    file = string(SavePath, "/avgReconstructionErrorPercentage_", parallelTotalProcesses)
    writedlm(file, avgReconstructionErrorPercentage,  "\t")

    file = string(SavePath, "/avgStability_", parallelTotalProcesses)
    writedlm(file, avgStability,  "\t")
end

println("Removing workers...")
if numberWorkers>1
    for i in workers()
        rmprocs(i)
    end
end

max_stars_repo_name: oriolpich/sigprofilerjulia

max_issues_repo_name: oriolpich/sigprofilerjulia

max_forks_repo_name: oriolpich/sigprofilerjulia


text:
"""
    BlockMatrix{T} <: AbstractArray{T, 2}

Implicit block matrix of equally-sized matrices with elements of type `T`.

# Constructors
    BlockMatrix(blocks)
    BlockMatrix{T}(blocks)
    BlockMatrix{T}(rows, cols[, blocks...])
    BlockMatrix(rows, cols, first_block[, more_blocks...])

Create a new block matrix from the given matrix of matrices. Alternatively, create the
block matrix from a sequence of matrices, given the number of `rows` and `cols`.

# Examples

```jldoctest; setup = :(using ImplicitArrays)
julia> A = [1 2; 3 4]; B = zeros(Int, 2, 2);

julia> C = [x * y for x in 2:3, y in 3:4];

julia> BlockMatrix(1, 4, A, B, B, C)
2×8 BlockMatrix{Int64}:
 1  2  0  0  0  0  6   8
 3  4  0  0  0  0  9  12
```

!!! note
    This type is not suitable for high-performance purposes, as matrix blocks are
    stored in an abstract container to allow the composition of arbitrary matrix types.
    For more details, please refer to the
    [performance tips](https://docs.julialang.org/en/v1/manual/performance-tips/)
    in the official Julia documentation.
"""
struct BlockMatrix{T} <: AbstractArray{T, 2}
    "Dimensions of the (equally-sized) matrix blocks"
    bdims::NTuple{2, Int}
    "Matrix blocks"
    blocks::Array{AbstractArray{T, 2}, 2}

    function BlockMatrix{T}(
        blocks::Array{AbstractArray{T, 2}, 2}
    ) where T
        bdims = (1, 1)
        for (i, block) ∈ enumerate(blocks)
            i > 1 && @assert(
                size(block) == bdims,
                "invalid block dimensions (block $i expected $bdims, got $(size(block)))"
            )
            bdims = size(block)
        end
        new(bdims, blocks)
    end
end

@inline BlockMatrix(
    blocks::Array{AbstractArray{T, 2}, 2}
) where T = BlockMatrix{T}(blocks)

function BlockMatrix{T}(
    rows::Int,
    cols::Int,
    blocks::Vararg{AbstractArray{T, 2}}
) where T
    @assert(
        rows * cols == length(blocks),
        "invalid number of blocks given (expected $(rows * cols), got $(length(blocks)))"
    )
    B = Array{AbstractArray{T, 2}, 2}(undef, rows, cols)
    for (i, b) in enumerate(blocks)
        B[Int(ceil(i/cols)), (i - 1) % cols + 1] = b
    end
    BlockMatrix{T}(B)
end

@inline BlockMatrix(
    rows::Int,
    cols::Int,
    block::AbstractArray{T, 2},
    blocks::Vararg{AbstractArray{T, 2}}
) where T = BlockMatrix{T}(rows, cols, block, blocks...)

@inline Base.size(
    M::BlockMatrix{T}
) where T = M.bdims .* size(M.blocks)

@inline Base.getindex(
    M::BlockMatrix{T},
    I::Vararg{Int, 2}
) where T = Base.getindex(
    M.blocks[Int(ceil(I[1]/M.bdims[1])), Int(ceil(I[2]/M.bdims[2]))],
    (I[1] - 1) % M.bdims[1] + 1,
    (I[2] - 1) % M.bdims[2] + 1
)

@inline Base.setindex!(
    M::BlockMatrix{T},
    v::Any,
    I::Vararg{Int, 2}
) where T = Base.setindex!(
    M.blocks[Int(ceil(I[1]/M.bdims[1])), Int(ceil(I[2]/M.bdims[2]))],
    v,
    (I[1] - 1) % M.bdims[1] + 1,
    (I[2] - 1) % M.bdims[2] + 1
)


text:
module TBfromk2r
    using Dierckx
    using SparseArrays
    const eps = 1e-7

    const σ0 = [
    1 0
    0 1
    ]

    const σx = [
    0 1
    1 0
    ]

    const σy = [
    0 -im
    im 0
    ]

    const σz = [
    1 0
    0 -1
    ]

    export surfaceHamiltonian,σx,σy,σz,σ0

    function surfaceHamiltonian(Hk,norb;numhop=3,L=32,kpara="kx",BC="OBC")
        function calc_H(k)
            if kpara == "kx"
                fx(kperp) = Hk([k,kperp])[:,:]
                mat_H = make_Hamiltonian(fx,[numhop],L,norb,BC)
            elseif kpara == "ky"
                fy(kperp) = Hk([kperp,k])[:,:]
                mat_H = make_Hamiltonian(fy,[numhop],L,norb,BC)
            else
                println("Error! not supported.")
                exit()
            end

            return mat_H
        end
        k ->  calc_H(k)
    end

    function make_Hamiltonian(Hk,numhop,L,norb,BC)
        Vol = prod(L)
        dim = length(L)
        mat_H = spzeros(ComplexF64,Vol*norb,Vol*norb)
        #println(Vol*norb)

        mat_t = make_hop(Hk,numhop,norb)

        for i = 1:Vol
            for jorb = 1:norb
                for iorb = 1:norb
                    ii = (iorb-1)*Vol + i
                    #println(ii,"\t",mat_t[iorb,jorb])
                    isite = calc_i2site(i,L)

                    for hop in mat_t[iorb,jorb]
                        jsite = isite .+ hop[1]
                        if BC == "PBC"
                            for il=1:dim
                                jsite[il] += ifelse(jsite[il] < 1,L[il],0)
                                jsite[il] += ifelse(jsite[il] > L[il],-L[il],0)
                            end
                        end
                        flag = false

                        if dim ==1
                            if 1 <= jsite[1] <= L[1]
                                flag = true
                            end
                        elseif dim == 2
                            if 1 <= jsite[1] <= L[1] && 1 <= jsite[2] <= L[2]
                                flag = true
                            end
                        end

                        if flag
                            j = calc_r2j(jsite,L)

                            jj = (jorb-1)*Vol + j

                            mat_H[ii,jj] = hop[2]
                        end
                    end
                end
            end
        end
        return mat_H

    end

    function calc_i2site(i,L)
        dim = length(L)
        if dim == 1
            isite = [i]
        elseif dim ==2
            ix = (i-1) % L[1] + 1
            iy = div(i-ix,L[1])+ 1
            isite = [ix,iy]
        end
        return  isite
    end

    function calc_r2j(hopping,L)
        dim = length(L)
        if dim == 1
            j = hopping[1]
        elseif dim == 2
            j = (hopping[2]-1)*L[1]+hopping[1]
        end
        return j
    end

    function make_hop(Hk,n::Array{Int,1},norb)
        mat_t  = Array{Any}(undef, norb, norb)

        for i=1:norb
            for j=1:norb
                f(k) = Hk(k)[i,j]
                hop_t = calc_coeff(f,n)
                mat_t[i,j] = hop_t
                #println("hoppings matrix at ($i,$j) element")
                #for hop in hop_t
                #    println(hop[1],"\t",hop[2])
                #end
            end
        end

        #println(mat_t)
        return mat_t

    end

    function calc_coeff(fk,n::Array{Int,1})
        dim = length(n)

        #println(dim)
        if dim == 1
            ck = calc_coeff(fk,n[1],1)
        elseif dim == 2
            ck = calc_coeff(fk,n[1],n[2],1,1)
        else
            println("not supported yet: dimension is $dim")
            return
        end
        return ck
    end

    function calc_coeff(fk,nkx::Int,ax::Number;L=256)
        numhop = 2nkx + 1
        ck = zeros(ComplexF64,numhop)

        hopping = range(-nkx, stop=nkx)
        kmax_x = 2π/ax
        kx = range(0.0,stop=kmax_x ,length=L)

        for i=1:numhop
            f(k) = fk(k)*exp(im*hopping[i]*k)
            mat_f = make_f(f,L,kmax_x)./(2π)
            spline_re = Spline1D(kx, real.(mat_f))
            spline_im = Spline1D(kx, imag.(mat_f))
            integ_re = integrate(spline_re, 0, kmax_x)
            integ_im = integrate(spline_im, 0, kmax_x)

            if abs(integ_re) > eps
                ck[i] += integ_re
            end
            if abs(integ_im) > eps
                ck[i] += im*integ_im
            end

        end

        hop_t = []
        for ix=1:numhop

            if abs(ck[ix]) > eps
                push!(hop_t,
                    (((hopping[ix]),ck[ix]))
                )
            end

        end

        return hop_t #ck ./(2π)
    end

    function make_f(f,L,kmax_x,kmax_y)
        mat_f = zeros(ComplexF64,L,L)
        kx = range(0.0,stop=kmax_x,length=L)
        ky = range(0.0,stop=kmax_y,length=L)

        for ikx =1:L
            for iky=1:L
                mat_f[ikx,iky] = f([kx[ikx],ky[iky]])#f(kx[ikx],ky[iky])
            end
        end

        return mat_f
    end

    function make_f(f,L,kmax_x)
        mat_f = zeros(ComplexF64,L)
        kx = range(0.0,stop=kmax_x,length=L)

        for ikx =1:L
            mat_f[ikx] = f(kx[ikx])#f(kx[ikx],ky[iky])
        end

        return mat_f
    end

    function calc_coeff(fk,nkx::Int,nky::Int,ax::Number,ay::Number;L=256)
        numhop_x = 2nkx + 1
        numhop_y = 2nky + 1
        ck = zeros(ComplexF64,numhop_x,numhop_y)
        hopping_x = range(-nkx, stop=nkx)
        hopping_y = range(-nky, stop=nky)
        kmax_x = 2π/ax
        kmax_y = 2π/ay
        kx = range(0.0,stop=kmax_x ,length=L)
        ky = range(0.0,stop=kmax_y,length=L)
        hop_t = []

        for ix=1:numhop_x
            for iy=1:numhop_y
                f(k) = fk(k)*exp(im*hopping_x[ix]*k[1])*exp(im*hopping_y[iy]*k[2])
                mat_f = make_f(f,L,kmax_x,kmax_y)./(2π)^2
                spline_re = Spline2D(kx, ky, real.(mat_f))
                spline_im = Spline2D(kx, ky, imag.(mat_f))
                integ_re = integrate(spline_re, 0, kmax_x, 0, kmax_y)
                integ_im = integrate(spline_im, 0, kmax_x, 0, kmax_y)
                if abs(integ_re) > eps
                    ck[ix,iy] += integ_re
                end
                if abs(integ_im) > eps
                    ck[ix,iy] += im*integ_im
                end
                #ck[ix,iy] = integ_re + im*integ_im
                #println("$(hopping_x[ix]), $(hopping_y[iy]), $(ck[ix,iy])")
            end
        end

        hop_t = []
        for ix=1:numhop_x
            for iy=1:numhop_y
                if abs(ck[ix,iy]) > eps
                    push!(hop_t,
                        (((hopping_x[ix],hopping_y[iy]),ck[ix,iy]))
                    )
                end
            end
        end

        return hop_t

    end

end


text:
"""
  punch_holes_2D(centers, radius, xpoints, ypoints)

...
# Arguments

  - `centers::Union{Vector}`: the vector containing the punch centers
  - `radius::Vector`: the tuple containing the punch radii
  - `Nx::Int64`: the number of points in the x direction
  - `Ny::Int64`: the number of points in the y direction
...

...
# Outputs
  - `absolute_indices::Vector{Int64}`: vector containing the indices of coordinates
  inside the punch
...

"""
function punch_holes_2D(centers, radius::Union{T, Vector{T}},
                        Nx, Ny) where T<:Number
    clen = (typeof(centers) <: Vector) ? length(centers) : 1
    rad  = (typeof(radius) <: Vector) ? radius : radius * ones(clen)
    masking_data_points = []
    absolute_indices = Int64[]
    for a = 1:clen
        c = centers[a]
        count = 1
        for j = 1:Ny
            for h = 1:Nx
                if (((h-c[1]))^2 + ((j-c[2]))^2  <= radius[a]^2)
                    append!(masking_data_points,[(h,j)])
                    append!(absolute_indices, count)
                end
                count = count + 1
            end
        end
    end
    return absolute_indices
end

"""
  punch_holes_3D(centers, radius, xpoints, ypoints, zpoints)

...
# Arguments

  - `centers::Vector{T}`: the vector containing the centers of the punches
  - `radius::Float64`: the radius of the punch
  - `Nx::Int64`: the number of points in the x-direction, this code is hard-coded to start from one.
  - `Ny::Int64`: the number of points in the y-direction
  - `Nz::Int64`: the number of points in the z-direction
...

...
# Outputs
  - `absolute_indices::Vector{Int64}`: vector containing the indices of coordinates
  inside the punch
...

"""
function punch_holes_3D(centers, radius, Nx, Ny, Nz)
    clen = length(centers)
    masking_data_points = []
    absolute_indices = Int64[]
    for a = 1:clen
        c = centers[a]
        count = 1
        for i = 1:Nz
            for j = 1:Ny
                for h = 1:Nx
                    if((h-c[1])^2 + (j-c[2])^2 + (i - c[3])^2 <= radius^2)
                        append!(masking_data_points, [(h, j, i)])
                        append!(absolute_indices, count)
                    end
                    count = count +1
                end
            end
        end
    end
    return absolute_indices
end

"""
  punch_3D_cart(center, radius, xpoints, ypoints, zpoints; <kwargs>)

...
# Arguments

  - `center::Tuple{T}`: the tuple containing the center of a round punch
  - `radius::Union{Tuple{Float64},Float64}`: the radii/radius of the punch
  - `x::Vector{T} where T<:Real`: the vector containing the x coordinate
  - `y::Vector{T} where T<:Real`: the vector containing the y coordinate
  - `z::Vector{T} where T<:Real`: the vector containing the z coordinate

# Optional

  - `linear::Bool = false` can return a linear index
...

...
# Outputs
  - `inds::Vector{Int64}`: vector containing the indices of coordinates
  inside the punch
...

"""
function punch_3D_cart(center, radius, x, y, z; linear = false)
    radius_x, radius_y, radius_z = (typeof(radius) <: Tuple) ? radius :
                                                (radius, radius, radius)
    inds = filter(i -> (((x[i[1]]-center[1])/radius_x)^2
                        + ((y[i[2]]-center[2])/radius_y)^2
                        + ((z[i[3]] - center[3])/radius_z)^2 <= 1.0),
                  CartesianIndices((1:length(x), 1:length(y), 1:length(z))))
    (length(inds) == 0) && error("Empty punch.")
    if linear == false
      return inds
    else
      return LinearIndices(zeros(length(x), length(y), length(z)))[inds]
    end
end

function bounding_box(x::Array{CartesianIndex{3},1})
    #for i in 1:p
    xa = map(t -> t[1], x); xa_min, xa_max = (minimum(xa), maximum(xa))
    xb = map(t -> t[2], x); xb_min, xb_max = (minimum(xb), maximum(xb))
    xc = map(t -> t[3], x); xc_min, xc_max = (minimum(xc), maximum(xc))
    # This requires julia v 1.6
    # return CartesianIndices((xa_min:xa_max, xb_min:xb_max, xc_min:xc_max))
    list = [CartesianIndex(x,y,z) for x in xa_min:xa_max for y in xb_min:xb_max for z in xc_min:xc_max]
    return list
end

intersect_box(A, B) = maximum([minimum(A), minimum(B)]):minimum([maximum(A), maximum(B)])

my_floor(x) = (x>0) ? floor(x) : -floor(abs(x))

"""

Generate a list of centers

# Example

```<julia-repl>
xmin = 0
xmax = 6
ymin = 0
ymax = 8
zmin = 0
zmax = 8

centers = LaplaceInterpolation.center_list('A', xmin, xmax, ymin, ymax, zmin,
zmax)
```
"""
function center_list(symm, Qh_min, Qh_max, Qk_min, Qk_max, Ql_min, Ql_max)
    hs = my_floor(Qh_min):my_floor(Qh_max)
    ks = my_floor(Qk_min):my_floor(Qk_max)
    ls = my_floor(Ql_min):my_floor(Ql_max)
    hkl = [(h,k,l) for h in hs for k in ks for l in ls]
    if symm == 'P'
        centers = filter(i -> ~P(i...), hkl)
    elseif symm == 'A'
        centers = filter(i -> ~A(i...), hkl)
    elseif symm == 'B'
        centers = filter(i -> ~B(i...), hkl)
    elseif symm == 'C'
        centers = filter(i -> ~C(i...), hkl)
    elseif symm == 'I'
        centers = filter(i -> ~Ii(i...), hkl)
    elseif symm == 'F'
        centers = filter(i -> ~F(i...), hkl)
    elseif symm == 'R'
        centers = filter(i -> ~R(i...), hkl)
    else
        centers = hkl
    end
    return centers
end

function center_check(symm)
    if symm == 'P'
        tof = (h, k, l) -> ~P(h, k, l)
    elseif symm == 'A'
        tof = (h, k, l) -> ~A(h, k, l)
    elseif symm == 'B'
        tof = (h, k, l) -> ~B(h, k, l)
    elseif symm == 'C'
        tof = (h, k, l)-> ~C(h, k, l)
    elseif symm == 'I'
        tof = (h, k, l) -> ~Ii(h, k, l)
    elseif symm == 'F'
        tof = (h, k, l) -> ~F(h, k, l)
    elseif symm == 'R'
        tof = (h, k, l) -> ~R(h, k, l)
    else
        tof = (h, k, l) -> true
    end
    return tof
end

# Systematic absences

function P(h,k,l)
    return false
end

function A(h,k,l)
    return (k+l)%2 != 0
end

function B(h,k,l)
    return (h+l)%2 != 0
end

function C(h,k,l)
    return (h+k)%2 != 0
end

function Ii(h,k,l)
    return (h+k+l)%2 != 0
end

function F(h,k,l)
    return ((h+k)%2!=0)||((h+l)%2!=0)||((k+l)%2!=0)
end

function R(h,k,l)
    return (-h+k+l)%3 != 0
end

max_stars_repo_name: rayosborn/nxpeaks

max_issues_repo_name: rayosborn/nxpeaks

max_forks_repo_name: axmas-anl/nxrefine