Based on population and identified peak load data, the according grid
    topology for residential sector is determined and attached to the grid graph

    Parameters
    ----------
    lvgd : LVGridDistrictDing0
        Low-voltage grid district object
340961
def build_residential_branches(lvgd):
    

    # Choice of typified lv model grid depends on population within lv
    # grid district. If no population is given, lv grid is omitted and
    # load is represented by lv station's peak load
    if lvgd.population > 0 \
            and lvgd.peak_load_residential > 0:
        model_grid = select_grid_model_residential(lvgd)

        build_lv_graph_residential(lvgd, model_grid)

    # no residential load but population
    elif lvgd.population > 0 \
            and lvgd.peak_load_residential == 0:
        logger.warning(
            '{} has population but no residential load. '
            'No grid is created.'.format(
                repr(lvgd)))

    # residential load but no population
    elif lvgd.population == 0 \
            and lvgd.peak_load_residential > 0:
        logger.warning(
            '{} has no population but residential load. '
            'No grid is created and thus this load is '
            'missing in overall balance!'.format(
                repr(lvgd)))

    else:
        logger.info(
            '{} has got no residential load. '
            'No grid is created.'.format(
                repr(lvgd)))
18521
def get_discrete_grid(self):
        
        sets_grid = []
        for d in self.space:
            if d.type == 'discrete':
                sets_grid.extend([d.domain]*d.dimensionality)
        return np.array(list(itertools.product(*sets_grid)))
92056
def get_grid_spatial_dimensions(self, variable):
        """"""

        data = self.open_dataset(self.service).variables[variable.variable]
        dimensions = list(data.dimensions)
        return data.shape[dimensions.index(variable.x_dimension)], data.shape[dimensions.index(variable.y_dimension)]
1158968
def peak_generation(self):
        
        cum_peak_generation = 0

        for lv_grid_district in self._lv_grid_districts:
            cum_peak_generation += lv_grid_district.lv_grid.station().peak_generation

        return cum_peak_generation
136356
def validate_grid_districts(self):
        #TODO: check docstring
        

        msg_invalidity = []
        invalid_mv_grid_districts = []

        for grid_district in self.mv_grid_districts():

            # there's only one node (MV station) => grid is empty
            if len(grid_district.mv_grid._graph.nodes()) == 1:
                invalid_mv_grid_districts.append(grid_district)
                msg_invalidity.append('MV Grid District {} seems to be empty ' \
                                      'and ' \
                                      'was removed'.format(grid_district))
            # there're only aggregated load areas
            elif all([lvla.is_aggregated for lvla in
                      grid_district.lv_load_areas()]):
                invalid_mv_grid_districts.append(grid_district)
                msg_invalidity.append("MV Grid District {} contains only " \
                                 "aggregated Load Areas and was removed" \
                                 "".format(grid_district))

        for grid_district in invalid_mv_grid_districts:
            self._mv_grid_districts.remove(grid_district)

        logger.warning("\n".join(msg_invalidity))
        logger.info('=====> MV Grids validated')
        return msg_invalidity
1211991
def add_station(self, lv_station):
        """"""
        if not isinstance(lv_station, LVStationDing0):
            raise Exception('Given LV station is not a LVStationDing0 object.')
        if self._station is None:
            self._station = lv_station
            self.graph_add_node(lv_station)
            self.grid_district.lv_load_area.mv_grid_district.mv_grid.graph_add_node(lv_station)
96216
def _vmomentsurfacemassGrid(self,n,m,grid):
        """"""
        if len(grid.df.shape) == 3: tlist= True
        else: tlist= False
        if tlist: 
            nt= grid.df.shape[2]
            out= []
            for ii in range(nt):
                out.append(nu.dot(grid.vRgrid**n,nu.dot(grid.df[:,:,ii],grid.vTgrid**m))*\
                    (grid.vRgrid[1]-grid.vRgrid[0])*(grid.vTgrid[1]-grid.vTgrid[0]))
            return nu.array(out)
        else:
            return nu.dot(grid.vRgrid**n,nu.dot(grid.df,grid.vTgrid**m))*\
                (grid.vRgrid[1]-grid.vRgrid[0])*(grid.vTgrid[1]-grid.vTgrid[0])
117706
def spatialDomainNoGrid(self):
    
    self.w = np.zeros(self.xw.shape)

    if self.Debug:
      print("w = ")
      print(self.w.shape)
    
    for i in range(len(self.q)):
      # More efficient if we have created some 0-load points
      # (e.g., for where we want output)
      if self.q[i] != 0:
        dist = np.abs(self.xw - self.x[i])
        self.w -= self.q[i] * self.coeff * np.exp(-dist/self.alpha) * \
          ( np.cos(dist/self.alpha) + np.sin(dist/self.alpha) )
142028
def get_mv_impedance(grid):
    

    omega = 2 * math.pi * 50

    mv_grid = grid.grid_district.lv_load_area.mv_grid_district.mv_grid
    edges = mv_grid.find_path(grid._station, mv_grid._station, type='edges')
    r_mv_grid = sum([e[2]['branch'].type['R'] * e[2]['branch'].length / 1e3
                     for e in edges])
    x_mv_grid = sum([e[2]['branch'].type['L'] / 1e3 * omega * e[2][
        'branch'].length / 1e3 for e in edges])

    return [r_mv_grid, x_mv_grid]
1154865
def _compute_grid(self):
        
        widths = []
        for c in range(self.cols):
            c_widths = []
            for r in range(self.rows):
                subplot = self.subplots.get((r, c), None)
                nsubplots = 1 if subplot is None else len(subplot.layout)
                c_widths.append(2 if nsubplots > 1 else 1)
            widths.append(max(c_widths))

        heights = []
        for r in range(self.rows):
            r_heights = []
            for c in range(self.cols):
                subplot = self.subplots.get((r, c), None)
                nsubplots = 1 if subplot is None else len(subplot.layout)
                r_heights.append(2 if nsubplots > 2 else 1)
            heights.append(max(r_heights))

        # Generate empty grid
        rows = sum(heights)
        cols = sum(widths)
        grid = [[None]*cols for _ in range(rows)]

        return grid
1050549
void putIntoCanonical( SquareGrid grid ) {

		boolean rowOdd = grid.rows%2 == 1;
		boolean colOdd = grid.columns%2 == 1;

		if( colOdd == rowOdd ) {
			// if odd and square then 4 solutions.  Otherwise just two solution that are on
			// opposite sides on the grid
			if( rowOdd && grid.rows == grid.columns ) {
				int best = -1;
				double bestDistance = Double.MAX_VALUE;
				for (int i = 0; i < 4; i++) {
					SquareNode n = grid.getCornerByIndex(i);
					double d = n.center.normSq();
					if( d < bestDistance ) {
						best = i;
						bestDistance = d;
					}
				}

				for (int i = 0; i < best; i++) {
					tools.rotateCCW(grid);
				}
			} else {
				double first = grid.get(0,0).center.normSq();
				double last = grid.getCornerByIndex(2).center.normSq();

				if( last < first ) {
					tools.reverse(grid);
				}
			}
		}
		// if only one is odd then there is a unique solution.  Since uber is already in a legit
		// configuration nothing needs ot be done
	}
581649
def dSbr_dV(Yf, Yt, V, buses, branches):
    
    nl = len(branches)
    nb = len(V)

    f = matrix([l.from_bus._i for l in branches])
    t = matrix([l.to_bus._i for l in branches])

    # Compute currents.
    If = Yf * V
    It = Yt * V

    Vnorm = div(V, abs(V))

    diagVf = spdiag(V[f])
    diagIf = spdiag(If)
    diagVt = spdiag(V[t])
    diagIt = spdiag(It)
    diagV = spdiag(V)
    diagVnorm = spdiag(Vnorm)

    ibr = range(nl)
    size = (nl, nb)
    # Partial derivative of S w.r.t voltage phase angle.
    dSf_dVa = 1j * (conj(diagIf) *
        spmatrix(V[f], ibr, f, size) - diagVf * conj(Yf * diagV))

    dSt_dVa = 1j * (conj(diagIt) *
        spmatrix(V[t], ibr, t, size) - diagVt * conj(Yt * diagV))

    # Partial derivative of S w.r.t. voltage amplitude.
    dSf_dVm = diagVf * conj(Yf * diagVnorm) + conj(diagIf) * \
        spmatrix(Vnorm[f], ibr, f, size)

    dSt_dVm = diagVt * conj(Yt * diagVnorm) + conj(diagIt) * \
        spmatrix(Vnorm[t], ibr, t, size)

    # Compute power flow vectors.
    Sf = mul(V[f], conj(If))
    St = mul(V[t], conj(It))

    return dSf_dVa, dSf_dVm, dSt_dVa, dSt_dVm, Sf, St
1005006
def peak_load(self):
        
        if self._peak_load is None:
            self._peak_load = sum(
                [_.peak_load.sum()
                 for _ in self.graph.nodes_by_attribute('load')])

        return self._peak_load
764085
def connect_generators(self, debug=False):
        

        self._graph = mv_connect.mv_connect_generators(self.grid_district, self._graph, debug)
1059049
def governor(self, Xgov, Pgov, Vgov):
        
        governors = self.governors
        omegas = 2 * pi * self.freq

        F = zeros(Xgov.shape)

        typ1 = [g.generator._i for g in governors if g.model == CONST_POWER]
        typ2 = [g.generator._i for g in governors if g.model == GENERAL_IEEE]

        # Governor type 1: constant power
        F[typ1, 0] = 0

        # Governor type 2: IEEE general speed-governing system
        Pm = Xgov[typ2, 0]
        P = Xgov[typ2, 1]
        x = Xgov[typ2, 2]
        z = Xgov[typ2, 3]

        K = Pgov[typ2, 0]
        T1 = Pgov[typ2, 1]
        T2 = Pgov[typ2, 2]
        T3 = Pgov[typ2, 3]
        Pup = Pgov[typ2, 4]
        Pdown = Pgov[typ2, 5]
        Pmax = Pgov[typ2, 6]
        Pmin = Pgov[typ2, 7]
        P0 = Pgov[typ2, 8]

        omega = Vgov[typ2, 0]

        dx = K * (-1 / T1 * x + (1 - T2 / T1) * (omega - omegas))
        dP = 1 / T1 * x + T2 / T1 * (omega - omegas)

        y = 1 / T3 * (P0 - P - Pm)

        y2 = y

        if sum(flatnonzero(y > Pup)) >= 1:
            y2 = (1 - flatnonzero(y > Pup)) * y2 + flatnonzero(y > Pup) * Pup
        if sum(flatnonzero(y < Pdown)) >= 1:
            y2 = (1 - flatnonzero(y<Pdown)) * y2 + flatnonzero(y<Pdown) * Pdown

        dz = y2

        dPm = y2

        if sum(flatnonzero(z > Pmax)) >= 1:
            dPm = (1 - flatnonzero(z > Pmax)) * dPm + flatnonzero(z > Pmax) * 0
        if sum(flatnonzero(z < Pmin)) >= 1:
            dPm = (1 - flatnonzero(z < Pmin)) * dPm + flatnonzero(z < Pmin) * 0

        F[typ2, :] = c_[dPm, dP, dx, dz]

        # Governor type 3:

        # Governor type 4:

        return F
404764
def Grid(
    pos=(0, 0, 0),
    normal=(0, 0, 1),
    sx=1,
    sy=1,
    c="g",
    bc="darkgreen",
    lw=1,
    alpha=1,
    resx=10,
    resy=10,
):
    
    ps = vtk.vtkPlaneSource()
    ps.SetResolution(resx, resy)
    ps.Update()
    poly0 = ps.GetOutput()
    t0 = vtk.vtkTransform()
    t0.Scale(sx, sy, 1)
    tf0 = vtk.vtkTransformPolyDataFilter()
    tf0.SetInputData(poly0)
    tf0.SetTransform(t0)
    tf0.Update()
    poly = tf0.GetOutput()
    axis = np.array(normal) / np.linalg.norm(normal)
    theta = np.arccos(axis[2])
    phi = np.arctan2(axis[1], axis[0])
    t = vtk.vtkTransform()
    t.PostMultiply()
    t.RotateY(theta * 57.3)
    t.RotateZ(phi * 57.3)
    tf = vtk.vtkTransformPolyDataFilter()
    tf.SetInputData(poly)
    tf.SetTransform(t)
    tf.Update()
    pd = tf.GetOutput()
    actor = Actor(pd, c=c, bc=bc, alpha=alpha)
    actor.GetProperty().SetRepresentationToWireframe()
    actor.GetProperty().SetLineWidth(lw)
    actor.SetPosition(pos)
    settings.collectable_actors.append(actor)
    return actor
1030525
def plot_mv_voltages(self, **kwargs):
        
        if self.network.pypsa is not None:
            try:
                v_res = self.network.results.v_res()
            except:
                logging.warning("Voltages `pfa_v_mag_pu` from power flow "
                                "analysis must be available to plot them.")
                return
            plots.mv_grid_topology(
                self.network.pypsa, self.network.config,
                timestep=kwargs.get('timestep', None),
                node_color='voltage',
                filename=kwargs.get('filename', None),
                grid_district_geom=kwargs.get('grid_district_geom', True),
                background_map=kwargs.get('background_map', True),
                voltage=v_res,
                limits_cb_nodes=kwargs.get('limits_cb_nodes', None),
                xlim=kwargs.get('xlim', None), ylim=kwargs.get('ylim', None),
                title=kwargs.get('title', ''))
        else:
            logging.warning("pypsa representation of MV grid needed to "
                            "plot voltages.")