Compare revisions

from geopandas import GeoDataFrame

from ...misc.utils import discrete_sinusoid_matching_integral

from ...misc.utils import create_profile_using_time_weighted_state

from ...misc.utils import generate_manual_state_correlated_profile

from ...misc.utils import generate_state_correlated_profile, generate_profile

from .bbr import label_bbr_entrance_id, label_bbr_housing_area

# *****************************************************************************

    number_intervals: int,

    time_interval_durations: list,

    bdg_specific_demand: dict,

    bdg_ratio_min_max: dict,

    bdg_min_max_ratio: dict,

    bdg_demand_phase_shift: dict = None,

    key_osm_entr_id: str = label_osm_entrance_id,

    key_bbr_entr_id: str = label_bbr_entrance_id,

    avg_state: list = None,

    state_correlates_with_output: bool = False,

    deviation_gain: float = 1.0,

    solver: str = 'glpk',

    **kwargs

) -> dict:

    # initialise dict for each building entrance

        # for each building

        for building_index in building_indexes:

            # get relevant data

            # base_load_avg_ratio = 0.3

            # specific_demand = 107 # kWh/m2/year

            area = gdf_buildings.loc[building_index][label_bbr_housing_area]

            # estimate its demand

            # generate the profile            

            if type(avg_state) == type(None):

                # ignore states

                heat_demand_profiles.append(

                        discrete_sinusoid_matching_integral(

                            bdg_specific_demand[building_index] * area,

                            time_interval_durations=time_interval_durations,

                            min_to_max_ratio=bdg_ratio_min_max[building_index],

                            min_max_ratio=bdg_min_max_ratio[building_index],

                            phase_shift_radians=(

                                bdg_demand_phase_shift[building_index]

                            ),

                    )

                )

            elif type(deviation_gain) == type(None):

                # states matter but the gain must be determined

                # heat_demand_profiles.append(

                #     np.array(

                #         create_profile_using_time_weighted_state(

                #             integration_result=(

                #                 bdg_specific_demand[building_index] * area

                #             ),

                #             avg_state=avg_state,

                #             time_interval_durations=time_interval_durations,

                #             min_max_ratio=bdg_min_max_ratio[building_index],

                #             state_correlates_with_output=state_correlates_with_output,

                #         )

                #     )

                # )

                heat_demand_profiles.append(

                    np.array(

                        generate_state_correlated_profile(

                            integration_result=(

                                bdg_specific_demand[building_index] * area

                            ), 

                            states=avg_state, 

                            time_interval_durations=time_interval_durations,

                            states_correlate_profile=state_correlates_with_output,

                            min_max_ratio=bdg_min_max_ratio[building_index],

                            solver=solver

                            )

                        )

                    )

            else:

                # states matter

                # states matter and the gain is predefined

                heat_demand_profiles.append(

                    np.array(

                        create_profile_using_time_weighted_state(

                        generate_manual_state_correlated_profile(

                            integration_result=(

                                bdg_specific_demand[building_index] * area

                            ), 

                            avg_state=avg_state,

                            states=avg_state, 

                            time_interval_durations=time_interval_durations, 

                            min_to_max_ratio=bdg_ratio_min_max[building_index],

                            state_correlates_with_output=state_correlates_with_output,

                            deviation_gain=deviation_gain

                            )

                        )

                    )

    # return

    return demand_dict

# *****************************************************************************

# *****************************************************************************

# TODO: allow reusing the gain

def heat_demand_profiles(

    gdf_osm: GeoDataFrame,

    gdf_buildings: GeoDataFrame,

    time_interval_durations: list,

    assessments: list,

    annual_heat_demand: dict,

    air_temperature: dict = None, 

    reuse_deviation_gain: bool = True,

    **kwargs

) -> dict:    

    # calculate the total area

    total_area = total_heating_area(gdf_osm, gdf_buildings)

    # initialise data dict

    heat_demand_dict = {}

    # for each building entrance

    for osm_index in gdf_osm.index:

        # initialise dict for each building entrance

        bdg_entr_dict = {}

        # find the indexes for each building leading to the curr. cons. point

        building_indexes = gdf_buildings[

            gdf_buildings[label_bbr_entrance_id] == gdf_osm.loc[osm_index][label_osm_entrance_id]

        ].index

        for q in assessments:

            # define the specific heat demand

            specific_demand = annual_heat_demand[q]/total_area

            # 

            # initialise dict for each building consumption point

            bdg_profile_list = []

            # for each building

            for building_index in building_indexes:

                # get relevant data

                area = gdf_buildings.loc[building_index][label_bbr_housing_area]

                # handle states

                if type(air_temperature) != type(None):

                    kwargs['states'] = air_temperature[q]

                # append the profile for each building to the list

                _profile = generate_profile(

                    integration_result=specific_demand*area, 

                    time_interval_durations=time_interval_durations, 

                    **kwargs

                    # min_max_ratio, 

                    # states, 

                    # states_correlate_profile, 

                    # solver, 

                    # deviation_gain

                    )

                bdg_profile_list.append(np.array(_profile))

            # aggregate profiles for the same building entrance

            bdg_entr_profile = (

                [] 

                if len(bdg_profile_list) == 0 else 

                sum(profile for profile in bdg_profile_list)

                )

            # store the demand profile

            bdg_entr_dict[q] = bdg_entr_profile

        # add to the main dict

        heat_demand_dict[osm_index] = bdg_entr_dict

    return heat_demand_dict, total_area

# *****************************************************************************

# *****************************************************************************

def total_heating_area(

    gdf_osm: GeoDataFrame,

            area += gdf_buildings.loc[building_index][label_bbr_housing_area]

    return area

# *****************************************************************************

# *****************************************************************************

 No newline at end of file

def find_simplifiable_paths(

    network: nx.MultiDiGraph,

    excluded_nodes: list,

    protected_nodes: list,

    ignore_self_loops: bool = False,

    consider_reversed_edges: bool = False,

    include_both_directions: bool = False,

    ----------

    network : nx.MultiDiGraph

        The object describing the graph.

    excluded_nodes : list

    protected_nodes : list

        A list of keys for nodes that cannot be in any straight path.

    ignore_self_loops : bool, optional

        If True, paths including self-loops can still be straight. If False,

            node_key

            for node_key in network.nodes()

            # the node cannot be among those excluded

            if node_key not in excluded_nodes

            if node_key not in protected_nodes

            # the node has to be linked to two other nodes other than itself

            if len(set(neighbours(network, node_key, ignore_self_loops=True))) == 2

            # exclude nodes with self-loops if desired:

# *****************************************************************************

def remove_self_loops(network: nx.MultiDiGraph):

def remove_self_loops(network: nx.MultiDiGraph) -> list:

    """

    Removes self-loops from a directed graph defined in a MultiDiGraph object.

    Returns

    -------

    generator-expression

    list

        The keys to the nodes whose self-loops were removed.

    """

def remove_dead_ends(

    network: nx.MultiDiGraph, keepers: tuple = None, max_iterations: int = 1

    network: nx.MultiDiGraph, 

    protected_nodes: list = None, 

    max_iterations: int = 1

) -> list:

    """

    Removes dead ends (non-cyclical branches) from a directed graph.

    ----------

    network : nx.MultiDiGraph

        The object describing the directed graph.

    keepers : tuple, optional

    protected_nodes : list, optional

        A list of keys for the nodes that are not to be considered for removal.

        The default is None, which means all nodes are under consideration.

    max_iterations : int, optional

    """

    if type(keepers) == type(None):

        keepers = []

    if type(protected_nodes) == type(None):

        protected_nodes = []

    # while true

    nodes_removed = []

        target_nodes = [

            node_key

            for node_key in network.nodes()

            if node_key not in keepers

            if node_key not in protected_nodes

            # if it has at most one neighbour other than itself

            if len(set(gis_iden.neighbours(network, node_key, ignore_self_loops=True)))

            <= 1

def remove_longer_parallel_edges(

    network: nx.MultiDiGraph, ignore_edge_directions: bool = False

    network: nx.MultiDiGraph, 

    distance_key: str = osm.KEY_OSMNX_LENGTH, 

    ignore_edge_directions: bool = True,

    protected_edges: list = None

) -> list:

    """

    Removes longer parallel edges from the network.

    Parallel edges are those connecting the same nodes in the same direction. If

    there are parallel edges between any given pair of nodes, the longer ones

    will be removed. If desired, edge directions can be ignored. In that case,

    only the shortest edge between the same pair of nodes will be retained.

    Parallel edges refer to multiple edges connecting the same nodes. If there 

    are parallel edges between any given pair of nodes, only the shortest one 

    will be retained. If desired, edge directions can be ignored. By default,

    edge directions are taken into account when selecting parallel edges.

    Parameters

    ----------

    network : nx.MultiDiGraph

        The object describing the graph.

    distance_key : str, optional

        The key used to obtain distances. The default is osm.KEY_OSMNX_LENGTH.

    ignore_edge_directions : bool, optional

        If True, edge directions are ignored. The default is True.

    protected_edges : list, optional

        A list of edges that should be retained. The default is None, which

        means all edges are elligible.

    Returns

    -------

        A list of the edges removed.

    """

    # redundancy: having more than one edge between two nodes

    # solution: remove the longest one, leave the shortest one

    # for each node pair

    if type(protected_edges) == type(None):

        # default: empty list

        protected_edges = []

    edges_removed = []

    for node_one in network.nodes():

        for node_two in network.nodes():

                list_edges = gis_iden.get_edges_from_a_to_b(

                    network, node_start=node_one, node_end=node_two

                )

            # if none exist, skip

            if len(list_edges) == 0:

                continue

            # otherwise, find out which is the shortest one

            # otherwise, sort them by distance (shorter distances first)

            # note: protected edges are considered in the sorting too

            sorted_edges = sorted(

                (network.edges[edge_key][osm.KEY_OSMNX_LENGTH], edge_key)

                (network.edges[edge_key][distance_key], edge_key)

                for edge_key in list_edges

            )

            network.remove_edges_from(edge_tuple[1] for edge_tuple in sorted_edges[1:])

            edges_removed.extend(edge_tuple[1] for edge_tuple in sorted_edges[1:])

            # edges to be removed (drop protected edges here)

            edges_for_removal = tuple(

                edge_tuple[1] 

                for edge_tuple in sorted_edges[1:]

                if edge_tuple[1] not in protected_edges

                )

            # remove all but the shortest edge

            network.remove_edges_from(edges_for_removal)

            # update the list of edges removed

            edges_removed.extend(edges_for_removal)

    # return statement

    return edges_removed

# *****************************************************************************

# *****************************************************************************

def merge_points_into_linestring(

    line: LineString,

    points: tuple or list,

from pandas import MultiIndex, Series

from numpy import float64, int64

from geopandas import GeoDataFrame, read_file

from shapely.geometry import Point

from shapely.geometry import Point, LineString

import contextily as cx

from shapely import intersects

# local, internal

def simplify_network(

    network: MultiDiGraph,

    protected_nodes: list,

    protected_nodes: list = None,

    protected_edges: list = None,

    dead_end_probing_depth: int = 5,

    remove_opposite_parallel_edges: bool = False,

    ignore_edge_directions: bool = True,

    update_street_count_per_node: bool = True,

    transform_roundabouts: bool = False,

    max_number_iterations: int = 5,

    **roundabout_conditions

):

    """

    network : MultiDiGraph

        The object describing the network.

    protected_nodes : list

        A list with the keys for the nodes that must be preserved.

    dead_end_probing_depth: int

        The maximum number of nodes a dead end can have to be detectable.

    remove_opposite_parallel_edges : bool, optional

        If True, longer parallel edges in opposite directions are also removed.

        The default is False.

        The keys for the nodes that must be preserved.

    protected_edges : list

        The keys for the edges that must be preserved.

    dead_end_probing_depth : int, optional

        The maximum number of nodes a dead end can have to be detectable. The 

        default is 5.

    ignore_edge_directions : bool, optional

        If True, direction is ignored in the search for parallel edges and

        simplifiable paths. The default is True.

    update_street_count_per_node : bool, optional

        If True, updates the street count on each node. The default is True.

    transform_roundabouts : bool, optional

        If True, roundabouts are to be transformed. The default is False.

    max_number_iterations : int, optional

        The maximum number of iterations. The default is 5.

    **roundabout_conditions : keyword and value pairs

        The conditions used to define which roundabouts are simplified.

    """

    # 1) remove dead ends (tends to create straight paths)

    gis_mod.remove_dead_ends(

        network, protected_nodes, max_iterations=dead_end_probing_depth

    if type(protected_nodes) == type(None):

        protected_nodes = []

    if type(protected_edges) == type(None):

        protected_edges = []

    else:

        # if there are protected edges, then the nodes involved in those edges

        # must also be preserved, otherwise they can be removed too

        protected_nodes.extend(

            # set(aaa for aa in a for aaa in aa[0:-1])

            set(nn for nnn in protected_edges for nn in nnn[0:-1])

            )

    iteration_counter = 0

    while iteration_counter < max_number_iterations:

        # remove self loops (can create straight paths and dead ends)

        looping_node_keys = gis_mod.remove_self_loops(network)

        # remove longer parallel edges (can create dead ends and straight paths)

        edge_keys = gis_mod.remove_longer_parallel_edges(

            network, 

            ignore_edge_directions=ignore_edge_directions,

            protected_edges=protected_edges,

        )

    # 2) remove longer parallel edges (tends to create straight paths)

    gis_mod.remove_longer_parallel_edges(

        network, ignore_edge_directions=remove_opposite_parallel_edges

        # remove dead ends (can create straight paths)

        node_keys = gis_mod.remove_dead_ends(

            network, 

            protected_nodes=protected_nodes, 

            max_iterations=dead_end_probing_depth

        )

    # 3) remove self loops (tends to create straight paths and dead ends)

    gis_mod.remove_self_loops(network)

    # 4) join segments (can create self-loops)

    simplifiable_paths = gis_iden.find_simplifiable_paths(network, protected_nodes)

    for path in simplifiable_paths:

        # join segments (can create self-loops and parallel edges)

        paths = gis_iden.find_simplifiable_paths(

            network, 

            protected_nodes=protected_nodes,

            consider_reversed_edges=ignore_edge_directions)

        for path in paths:

            gis_mod.replace_path(network, path)

    # 4) remove self loops (tends to create straight paths and dead ends)

    gis_mod.remove_self_loops(network)

    # 5) transform roundabouts into crossroads (can create straight paths)

        # update iteration counter

        iteration_counter += 1

        # check if it makes sense to break out of the loop

        if (len(looping_node_keys) == 0 and 

            len(edge_keys) == 0 and 

            len(paths) == 0 and 

            len(node_keys) == 0):

            # no self-loops

            # no edges were removed

            # no paths were simplified

            # no nodes were removed

            break

    # transform roundabouts into crossroads (can create straight paths)

    if transform_roundabouts:

        list_roundabout_nodes = gis_iden.find_roundabouts(network, **roundabout_conditions)

        gis_mod.transform_roundabouts_into_crossroads(network, list_roundabout_nodes)

    # 6) update street count

    # update street count

    if update_street_count_per_node:

        gis_calc.update_street_count(network)

    # return statement

    return node_path

# *****************************************************************************

# *****************************************************************************

def create_edge_geometry(

        network: MultiDiGraph, 

        edge_key,

        x_key = osm.KEY_OSMNX_X,

        y_key = osm.KEY_OSMNX_Y) -> LineString:

    "Returns a newly-created geometry for a given edge."

    return LineString(

        [(network.nodes[edge_key[0]][x_key],

          network.nodes[edge_key[0]][y_key]),

         (network.nodes[edge_key[1]][x_key],

          network.nodes[edge_key[1]][y_key])]

        )

# *****************************************************************************

# *****************************************************************************

def find_overlapping_edges(

        network: MultiDiGraph,

        excluded_edges: list = None

        ) -> list:

    """

    Returns a list of key pairs for edges whose geometries overlap.

    Parameters

    ----------

    network : MultiDiGraph

        The object describing the network.

    excluded_edges : list, optional

        A list of edges that should not be considered. The default is None, in

        which case all edges in the network object will be considered.

    Returns

    -------

    list

        A list containing key pairs for overlapping edges.

    """

    # check if there are excluded edges

    if type(excluded_edges) == type(None):

        excluded_edges = list()

    # initialise the list of edges to check

    edges = list(

        edge_key 

        for edge_key in network.edges(keys=True)

        if edge_key not in excluded_edges

        )

    visited_edges = []

    out = []

    # for each edge

    for edge_key in edges:

        # remove the current edge so it is not considered again

        visited_edges.append(edge_key)

        # for each other edge

        for other_edge_key in edges:

            # for each other edge

            # skip edges having nodes in common

            # this will also skip identical edges

            if edge_key[0] in other_edge_key[0:2] or edge_key[1] in other_edge_key[0:2]:

                # has nodes in common, skip

                continue

            # skip edges that have already been considered in the first loop

            if other_edge_key in visited_edges:

                # this edge has already been tested against all other edges

                continue

            # first edge                

            if osm.KEY_OSMNX_GEOMETRY in network.edges[edge_key]:

                first_geo = network.edges[edge_key][osm.KEY_OSMNX_GEOMETRY]

            else:

                first_geo = create_edge_geometry(network, edge_key)

            # second edge

            if osm.KEY_OSMNX_GEOMETRY in network.edges[other_edge_key]:

                second_geo = network.edges[other_edge_key][osm.KEY_OSMNX_GEOMETRY] 

            else:

                second_geo = create_edge_geometry(network, other_edge_key)

            # check if they intersect

            if intersects(first_geo, second_geo):

                # they do, add tuple of the edges to the output

                out.append((edge_key, other_edge_key))

    # return tuples of overlapping edges

    return out

# *****************************************************************************

# *****************************************************************************

# *****************************************************************************

# standard

import uuid

import math

from statistics import mean

# local, external

import pyomo.environ as pyo

# *****************************************************************************

# *****************************************************************************

def discrete_sinusoid_matching_integral(

    integration_result: float,

    time_interval_durations: list,

    min_to_max_ratio: float,

    min_max_ratio: float,

    phase_shift_radians: float = None,

) -> list:

    """

    where:

    a = b*(1-min_to_max_ratio)/(1+min_to_max_ratio)

    a = b*(1-min_max_ratio)/(1+min_max_ratio)

    b = integration_result/integration_period

        The result of integrating the sinusoidal function for one period.

    time_interval_durations : list

        The time interval durations for each sample.

    min_to_max_ratio : float

    min_max_ratio : float

        The ratio between the minimum and maximum values of the function.

    phase_shift_radians : float, optional

        The phase shift for the sinusoidal function. The default is None, which

    b = integration_result / integration_period

    a = b * (1 - min_to_max_ratio) / (1 + min_to_max_ratio)

    a = b * (1 - min_max_ratio) / (1 + min_max_ratio)

    alpha = 2 * math.pi / integration_period

    By default, the profiles are synched: the highest sample in one is placed

    in the same position as the highest sample in the other; the second highest

    sample is placede in the same position as the second highest sample in the

    sample is placed in the same position as the second highest sample in the

    other profile; and so on. Alternatively, the profiles can be synched in

    reverse: the highest sample in one profile is placed in the same position

    as the lowest sample in the other; and so on and so forth.

def create_profile_using_time_weighted_state(

    integration_result: float,

    avg_state: list,

    states: list,

    time_interval_durations: list,

    min_to_max_ratio: float,

    state_correlates_with_output: bool = True,

    min_max_ratio: float,

    states_correlate_profile: bool = True,

) -> list:

    """

    Returns a profile that approximates a sinusoidal function in discrete time.

    where:

    a = b*(1-min_to_max_ratio)/(1+min_to_max_ratio)

    a = b*(1-min_max_ratio)/(1+min_max_ratio)

    b = integration_result/integration_period

    ----------

    integration_result : float

        The result of integrating the sinusoidal function for one period.

    avg_state : list

    states : list

        The average state during each time interval.

    time_interval_durations : list

        The time interval durations for each sample.

    min_to_max_ratio : float

    min_max_ratio : float

        The ratio between the minimum and maximum values of the function.

    state_correlates_with_output : bool, optional

    states_correlate_profile : bool, optional

        If True, the peak should happen when the state is at its highest point.

        If False, the peak should happen when the state is at its lowest point.

        The default is True.

    """

    if len(avg_state) != len(time_interval_durations):

    if len(states) != len(time_interval_durations):

        raise ValueError("The inputs are inconsistent.")

    period = sum(time_interval_durations)

    avg_time_interval_duration = mean(time_interval_durations)

    avg_state_weighted = [

    states_weighted = [

        (

            x_k * delta_k / avg_time_interval_duration

            if state_correlates_with_output

            if states_correlate_profile

            else -x_k * delta_k / avg_time_interval_duration

        )

        for delta_k, x_k in zip(time_interval_durations, avg_state)

        for delta_k, x_k in zip(time_interval_durations, states)

    ]

    # find the peak