def simultaneous_backward_reduction(scenarios, number_of_reduced_scenarios, probability=None):
    """Simultaneous backward reduction algorithm

    Parameters
    ----------
    scenarios : numpy.array
        Contain the input scenarios.
        The columns representing the individual scenarios
        The rows are the vector of values in each scenario
    number_of_reduced_scenarios : int
        final number of scenarios that
        the reduced scenarios contain.
        If number of scenarios is equal to or greater than the input scenarios,
        then the original input scenario set is returned as the reduced set
    probability : numpy.array (default=None)
        probability is a numpy.array with length equal to number of scenarios.
        if probability is not defined, all scenarios get equal probabilities

    Returns
    -------
    reduced_scenarios : numpy.array
        reduced set of scenarios
    reduced_probability : numpy.array
        probability of reduced set of scenarios
    reduced_scenario_set : list
        scenario numbers of reduced set of scenarios

    Example
    -------

    Scenario reduction can be performed as shown below::

        >>> import numpy as np
        >>> import random
        >>> scenarios = np.array([[random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)]])
        >>> import psst.scenario
        >>> reduced_scenarios, reduced_probability, reduced_scenario_numbers = psst.scenario.simultaneous_backward_reduction(scenarios, probability, 2)
    """

    print("Running simultaneous backward reduction algorithm")

    number_of_scenarios = scenarios.shape[1]
    logger.debug("Input number of scenarios = %d", number_of_scenarios)
    # if probability is not defined assign equal probability to all scenarios
    if probability is None:
        probability = np.array([1/number_of_scenarios for i in range(0, number_of_scenarios)])

    # initialize z, c and J
    z = np.array([np.inf for i in range(0, number_of_scenarios)])
    c = np.zeros((number_of_scenarios, number_of_scenarios))
    J = []

    # no reduction necessary
    if number_of_reduced_scenarios >= number_of_scenarios:
        return(scenarios, probability, J)

    """compute the distance of scenario pairs"""

    for scenario_k in range(0, number_of_scenarios):
        for scenario_j in range(0, number_of_scenarios):
            c[scenario_k, scenario_j] = distance(scenarios[:, scenario_k], scenarios[:, scenario_j])

    for scenario_l in range(0, number_of_scenarios):
        lowest_value = np.inf
        for scenario_j in range(0, number_of_scenarios):
            if scenario_l == scenario_j:
                continue
            lowest_value = min(lowest_value, c[scenario_l, scenario_j])

        c[scenario_l, scenario_l] = lowest_value
        z[scenario_l] = probability[scenario_l]*c[scenario_l, scenario_l]

    J.append(np.argmin(z))

    for _ in range(0, number_of_scenarios - number_of_reduced_scenarios - 1):

        for scenario_l in range(0, number_of_scenarios):
            for scenario_k in range(0, number_of_scenarios):
                if scenario_k in J or scenario_k == scenario_l:
                    if scenario_l not in J:
                        lowest_value = np.inf
                        for scenario_j in range(0, number_of_scenarios):
                            if scenario_j not in J and scenario_j != scenario_l:
                                lowest_value = min(lowest_value, c[scenario_k, scenario_j])

                        c[scenario_k, scenario_l] = lowest_value

        for scenario_l in range(0, number_of_scenarios):
            if scenario_l not in J:
                summation = 0

                for scenario_k in range(0, number_of_scenarios):
                    if scenario_k in J or scenario_k == scenario_l:
                        summation = summation + probability[scenario_k]*c[scenario_k, scenario_l]

                z[scenario_l] = summation

        J.append(np.argmin([item if i not in J else np.inf for i, item in enumerate(z)]))

    reduced_scenario_set = []
    for scenario_number in range(0, number_of_scenarios):
        if scenario_number not in J:
            reduced_scenario_set.append(scenario_number)

    reduced_probability = []

    reduced_probability = copy.deepcopy(probability)
    for deleted_scenario_number in J:
        lowest_value = np.inf

        # find closest scenario_number
        for scenario_j in reduced_scenario_set:
            if c[deleted_scenario_number, scenario_j] < lowest_value:
                closest_scenario_number = scenario_j
                lowest_value = c[deleted_scenario_number, scenario_j]

        reduced_probability[closest_scenario_number] = reduced_probability[closest_scenario_number] + reduced_probability[deleted_scenario_number]

    reduced_scenarios = copy.deepcopy(scenarios[:, reduced_scenario_set])
    reduced_probability = reduced_probability[reduced_scenario_set]
    return(reduced_scenarios, reduced_probability, reduced_scenario_set)
Simultaneous backward reduction algorithm

    Parameters
    ----------
    scenarios : numpy.array
        Contain the input scenarios.
        The columns representing the individual scenarios
        The rows are the vector of values in each scenario
    number_of_reduced_scenarios : int
        final number of scenarios that
        the reduced scenarios contain.
        If number of scenarios is equal to or greater than the input scenarios,
        then the original input scenario set is returned as the reduced set
    probability : numpy.array (default=None)
        probability is a numpy.array with length equal to number of scenarios.
        if probability is not defined, all scenarios get equal probabilities

    Returns
    -------
    reduced_scenarios : numpy.array
        reduced set of scenarios
    reduced_probability : numpy.array
        probability of reduced set of scenarios
    reduced_scenario_set : list
        scenario numbers of reduced set of scenarios

    Example
    -------

    Scenario reduction can be performed as shown below::

        >>> import numpy as np
        >>> import random
        >>> scenarios = np.array([[random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)],
        >>>         [random.randint(500,1000) for i in range(0,24)]])
        >>> import psst.scenario
        >>> reduced_scenarios, reduced_probability, reduced_scenario_numbers = psst.scenario.simultaneous_backward_reduction(scenarios, probability, 2)