text:
#!/usr/bin/env python
# coding: utf-8

# In[23]:

from path import Path
from fastai.vision import *
import numpy as np
import pandas as p
import matplotlib.cm as cmx
import matplotlib.colors as mcolors
from cycler import cycler
from typing import List, Union, Tuple, Callable, Optional

# In[20]:

class StdConv(nn.Module):
    def __init__(self, nin: int, nout: int, filter_size: int=3, stride: int=2, padding: int=1, dropout: float=0.1):
        super().__init__()
        self.conv = nn.Conv2d(nin, nout, filter_size, stride, padding)
        self.bn = nn.BatchNorm2d(nout)
        self.drop = nn.Dropout(dropout)
    def forward(self, x):
        return self.drop(self.bn(F.relu(self.conv(x))))

def flatten_conv(x: torch.Tensor, k: int):
    bs,nf,gx,gy = x.size()
    x = x.permute(0,2,3,1).contiguous()
    return x.view(bs, -1, nf//k)

def conv_params(sz_in: int, sz_out: int):
    """Solves for filter_size, padding and stride per the following equation,
       sz_out = (sz_in - filter_size + 2*padding) / stride + 1
       Attempts to find a solution by iterating over various filter_size, stride and padding
       in that order. If no solution is found, raises an error
    """
    filter_size, stride, padding = [3,2,4,5], [2,1,3], [1,2,3]
    for f in filter_size:
        for s in stride:
            for p in padding:
                if ((sz_in - f + 2*p) / s + 1) == sz_out:
                    return (f, s, p)
    raise Exception("Unable to find valid parameters for {0} to {1} convolution".format(sz_in, sz_out))

class OutConv(nn.Module):
    def __init__(self, k: int, nin: int, nclasses: int, bias: float):
        super().__init__()
        self.k = k
        self.oconv1 = nn.Conv2d(nin, nclasses * self.k, 3, padding=1)
        self.oconv2 = nn.Conv2d(nin, 4*self.k, 3, padding=1)
        self.oconv1.bias.data.zero_().add_(bias)

    def forward(self, x: torch.Tensor):
        return [flatten_conv(self.oconv1(x), self.k), flatten_conv(self.oconv2(x), self.k)]

class SSD_MultiHead(nn.Module):
    def __init__(self, k: int, grids: List[int], num_classes: int, bias: float, dropout: float=0.1, device: str='cuda:0'):
        super().__init__()
        self.k, self.grids, self.num_classes = k, grids, num_classes
        self.drop = nn.Dropout(dropout)

        self.sconv0 = StdConv(512, 256, stride=1, dropout=dropout).to(device)
        self.sconvs = []
        self.outconvs = []

        for i in range(len(self.grids)):
            sz_in = 7 if i == 0 else self.grids[i-1] #the 0th conv is taking a 7x7 grid from resnet34 output
            sz_out = self.grids[i]
            filter_size, stride, padding = conv_params(sz_in, sz_out)

            self.sconvs.append(StdConv(256, 256, filter_size=filter_size, stride=stride, padding=padding, dropout=dropout).to(device))
            self.outconvs.append(OutConv(self.k, 256, self.num_classes, bias).to(device))

    def forward(self, x: torch.Tensor):
        x = self.drop(F.relu(x))
        x = self.sconv0(x)
        oc, ol = [], []
        for i in range(len(self.grids)):
            x = self.sconvs[i](x)
            oc_i, ol_i = self.outconvs[i](x)
            oc.append(oc_i)
            ol.append(ol_i)
        return [torch.cat(oc, dim=1),
                torch.cat(ol, dim=1)]

# In[21]:

def one_hot_embedding(labels: torch.Tensor, num_classes: int):
    """Returns a one hot embedded matrix
    """
    return torch.eye(num_classes)[labels.data.cpu()]

class BCE_Loss(nn.Module):
    """Binary cross entropy loss
    """
    def __init__(self, num_classes: int):
        super().__init__()
        self.num_classes = num_classes

    def forward(self, pred: torch.Tensor, targ: torch.Tensor):
        t = one_hot_embedding(targ, self.num_classes)
        t = t[:,1:].contiguous()
        x = pred[:, 1:].cpu()
        w = (self.get_weight(x, t))
        return F.binary_cross_entropy_with_logits(x, t, w, size_average=False)/(self.num_classes - 1)

    def get_weight(self, x: torch.Tensor, t: torch.Tensor):
        return None

class FocalLoss(BCE_Loss):
    """Focal loss modifies BCE_Loss by de-emphasizing well classified categories
       and emphasizing poorly classified categories. Setting gamma to 0 is equivalent
       to the standard BCE_Loss
    """
    def get_weight(self, x: torch.Tensor, t: torch.Tensor, alpha: float=.5, gamma: float=2.):
        p = x.sigmoid()
        pt = p*t + (1-p)*(1-t)
        w = alpha*t + (1-alpha)*(1-t)
        w = w * (1-pt).pow(gamma)
        return w.detach()

# In[27]:

class SingleShotDetector():
    """Single shot detector model for multiple object detection
       Takes databunch object and creates learner for training and inference
    """
    def __init__(self, data: ObjectItemList, anc_grids: List[int]=[4, 2, 1], anc_zooms: List[float]=[.7, 1., 1.3],
                       anc_ratios: List[Tuple[float, float]]=[(1., 1.), (1., 0.5), (0.5, 1.)],
                       base_arch: Callable=models.resnet34, dropout: float=0.1, bias: float=-4.,
                       loss_func: nn.Module=FocalLoss, thresh: float=.4, metrics: List[Callable]=[],
                       num_color: int=12, verbose: bool=False):

        ###Init Data###
        self.data = data
        self.num_classes = self.data.c

        ###Init Device###
        self.device = 'cuda:0'

        ###Init Anchors###
        #default settings give a 4x4, 2x2 and 1x1 grid for a total of 16 + 4 + 1 = 21 anchor points
        #3 zooms and 3 ratios for a total of 9 boxes per anchor point
        #9 x 21 = 189 total anchor boxes
        self.anc_grids, self.anc_zooms, self.anc_ratios = anc_grids, anc_zooms, anc_ratios
        self.anchor_scales = [(anz*i, anz*j) for anz in anc_zooms for (i, j) in anc_ratios] #len(zooms) * len(ratios)
        self.k = len(self.anchor_scales)
        self.anc_offsets = [1/(o*2) for o in self.anc_grids] #scale grids to [0,1]
        self.anc_x = np.concatenate([np.repeat(np.linspace(ao, 1-ao, ag), ag) for ao,ag in zip(self.anc_offsets, self.anc_grids)])
        self.anc_y = np.concatenate([np.tile(np.linspace(ao, 1-ao, ag), ag) for ao,ag in zip(self.anc_offsets, self.anc_grids)])
        self.anc_ctrs = np.repeat(np.stack([self.anc_x, self.anc_y], axis=1), self.k, axis=0)
        self.anc_sizes  =   np.concatenate([np.array([[o/ag,p/ag] for i in range(ag*ag) for o,p in self.anchor_scales]) for ag in self.anc_grids])
        self.grid_sizes = torch.tensor(np.concatenate([np.array([1/ag for i in range(ag*ag) for o,p in self.anchor_scales])for ag in self.anc_grids])).unsqueeze(1).to(self.device)
        self.anchors = torch.tensor(np.concatenate([self.anc_ctrs, self.anc_sizes], axis=1)).to(self.device)
        self.anchor_cnr = self._hw2corners(self.anchors[:,:2], self.anchors[:,2:]).to(self.device)

        ###Init Model###
        self.base_arch, self.dropout, self.bias, self.thresh = base_arch, dropout, bias, thresh
        self.loss_func = loss_func(self.num_classes)
        self.head = SSD_MultiHead(self.k, self.anc_grids, self.num_classes, self.bias, self.dropout, device=self.device)
        self.learn = cnn_learner(self.data, self.base_arch, custom_head=self.head, metrics=metrics)
        self.learn.loss_func = self._ssd_loss
        self.learn.model.to(self.device)

        ###Init Drawing###
        self.num_color = num_color
        cmap = self._get_cmap(self.num_color)
        self.color_list = [cmap(float(x)) for x in range(self.num_color)]

        ###Init Misc###
        self.verbose = verbose

    def show_anchors(self):
        """Plots anchor points on blank plot
        """
        fig2 = plt.figure()
        ax2 = fig2.add_subplot(111, aspect='equal')
        plt.scatter(self.anc_x, self.anc_y)
        plt.xlim(0, 1)
        plt.ylim(0, 1)

    def show_results(self, prb_thresh: float=.2, nms_thresh: float=.4, ignore_bg: bool=True):
        """Shows predictions on a batch of images
        """
        self.learn.model.eval()
        x,y = next(iter(self.data.valid_dl))
        batch = self.learn.model(x)
        b_clas,b_bb = batch
        #let's visualize the the model after some training
        fig, axes = plt.subplots(3, 4, figsize=(16, 12))
        for i, ax in enumerate(axes.flat):
            ima = to_np(self.data.denorm(x)[i].permute(1,2,0)) #denormalize x at idx and assign to ima var
            bbox,clas = self._get_y(y[0][i], y[1][i]) #get ground truth bbox and class
            a_ic = self._actn_to_bb(b_bb[i], self.anchors) #convert bounding box predictions into plottable format
            to_keep = self._nms(a_ic, b_clas[i].max(1)[1], b_clas[i].max(1)[0].sigmoid(), thresh=nms_thresh)
            self._torch_gt(ax, ima, a_ic[to_keep], b_clas[i].max(1)[1][to_keep], b_clas[i].max(1)[0].sigmoid()[to_keep], thresh=prb_thresh, ignore_bg=ignore_bg)

        plt.subplots_adjust(wspace=0.15, hspace=0.15)

    def _ssd_1_loss(self, b_c: torch.Tensor, b_bb: torch.Tensor, bbox: torch.Tensor, clas: torch.Tensor,
                    thresh: Optional[float]=None, verbose: Optional[bool]=None):
        """Computes loss for one image and set of bbox predictions
           Takes b_c and b_bb which are predicted classes and bounding boxes (output from NN)
           and bbox and clas which are ground truth bounding box and classes for image
        """
        thresh = self.thresh if not thresh else thresh
        verbose = self.verbose if not verbose else verbose
        bbox, clas = self._get_y(bbox, clas) #reformat gt bbox and gt clas
        if len(bbox) == 0: return 0, 0 #if there's no bbox, we can't use this in our loss
        a_ic = self._actn_to_bb(b_bb, self.anchors) #map activation b_bb to anchor points (defined above)
        overlaps = self._jaccard(bbox.data, self.anchor_cnr.data.float()) #get overlaps between gt bbox and anchor boxes
        gt_overlap, gt_idx = self._map_to_ground_truth(overlaps, verbose)
        gt_clas = clas[gt_idx]
        pos = gt_overlap > thresh
        pos_idx = torch.nonzero(pos)[:, 0] #get indices of anchor boxes with sufficient overlap with GT
        gt_clas[1-pos] = 0 #assign background class
        gt_bbox = bbox[gt_idx]
        loc_loss = ((a_ic[pos_idx] - gt_bbox[pos_idx]).abs()).mean() #get average absolute loss between coordinates
        clas_loss = self.loss_func(b_c, gt_clas) #get cross entropy loss on predicted classes
        return loc_loss, clas_loss

    def _ssd_loss(self, pred: torch.Tensor, targ0: torch.Tensor, targ1: torch.Tensor,
                        verbose: Optional[bool]=None):
        """For a batch of predictions and a batch of (targ0, targ1) where targ0 is ground truth bbox and
           targ1 is ground truth classes compute the loss. We do this by iterating over the batch and computing
           loss on each individual image and summing the results. I believe the losses must be placed on the CPU
           due to the way fastai is handling the loss function
        """
        verbose = self.verbose if not verbose else verbose
        lcs,lls=0.,0.
        for b_c,b_bb,bbox,clas in zip(*pred, targ0, targ1):
            if len(bbox) == 0:
                continue
            loc_loss, clas_loss = self._ssd_1_loss(b_c, b_bb, bbox, clas, verbose=verbose)
            lls += loc_loss
            lcs += clas_loss
        if verbose:
            print("loc: {0}, clas: {1}".format(lls.data[0], lcs.data[0]))
        return lls.cpu() + lcs.cpu() #need to revisit why cpu is necessary

    def kaggle_1_loss(self, pred_bbox: torch.Tensor, pred_clas: torch.Tensor,
                            gt_bbox: torch.Tensor, gt_clas: torch.Tensor):
        """Custom kaggle error metric, checks if any predicted bboxes
           match the class of any ground truth bboxes AND if those bboxes
           have an overlap >= 50%, if so returns 0 else 1. This acts on
           one image. Note that this selects the top 5 highest confidence bboxes
           that are not background for comparison
        """
        gt_bbox, gt_clas = self._get_y(gt_bbox, gt_clas) #reformat gt bbox and gt clas
        pred_bbox = self._actn_to_bb(pred_bbox, self.anchors) #map pred_bbox to anchor points
        pred_clas_maxprb, pred_clas_maxidx = pred_clas[:, 1:].max(dim=1) #get the highest probability and class for each bbox
        pred_clas_maxprb = torch.exp(pred_clas_maxprb) * (pred_clas_maxidx != 0).float() #if background is most likely class, ignore
        pred_top = np.argsort(pred_clas_maxprb)[:5] #get the 5 bboxes with highest probability
        clas_top = pred_clas_maxidx[pred_top] #top 5 classes
        bbox_top = pred_bbox[pred_top] #top 5 bboxes
        for i_pred,clas_pred in enumerate(clas_top):
            for i_gt,clas_gt in enumerate(gt_clas):
                if clas_pred !=  clas_gt:
                    continue
                overlap = self._jaccard(bbox_top[i_pred].unsqueeze(0), gt_bbox[i_gt].unsqueeze(0)).item()
                if overlap >= .5:
                    return 0
        return 1

    def kaggle_loss(self, preds: List[torch.Tensor], targ_bbox: torch.Tensor, targ_clas: torch.Tensor):
        """Custom kaggle error metric, checks if any predicted bboxes
           match the class of any ground truth bboxes AND if those bboxes
           have an overlap >= 50%, if so returns 0 else 1. This acts on a batch of images.
        """
        loss = 0
        for pred_clas,pred_bbox,gt_bbox,gt_clas in zip(*preds, targ_bbox, targ_clas):
            if len(gt_bbox) == 0:
                continue
            loss += self.kaggle_1_loss(pred_bbox, pred_clas, gt_bbox, gt_clas)
        return torch.Tensor([loss / targ_clas.shape[0]]).to(self.device)

    def _show_ground_truth(self, ax, im, bbox: torch.Tensor, clas: Optional[torch.Tensor]=None,
                                 prs: Optional[torch.Tensor]=None, thresh: float=.2, ignore_bg: bool=False):
        """Given ax, image, list of bboxes, list of classes, list of probabilities and threshold
           Displays image with ground truth broxes
        """
        bb = [self._bb_hw(o) for o in bbox.reshape(-1,4)]
        if prs is None:  prs  = [None]*len(bb)
        if clas is None: clas = [None]*len(bb)
        ax = self._show_img(im, ax=ax)
        for i,(b,c,pr) in enumerate(zip(bb, clas, prs)):
            if c == 0 and ignore_bg:
                continue
            if((b[2]>0) and (pr is None or pr > thresh)):
                self._draw_rect(ax, b, color=self.color_list[i%self.num_color])
                txt = f'{i}: '
                if c is not None: txt += str(self.data.classes[c])
                if pr is not None: txt += f' {pr:.2f}'
                self._draw_text(ax, b[:2], txt, color=self.color_list[i%self.num_color])

    def _torch_gt(self, ax, ima, bbox: torch.Tensor, clas: torch.Tensor,
                        prs: Optional[torch.Tensor]=None, thresh: float=.2,
                        ignore_bg: bool=False):
        """Show ground truth given image and torch bbox
        """
        return self._show_ground_truth(ax, ima, to_np((bbox*224).long()),
             to_np(clas), to_np(prs) if prs is not None else None, thresh, ignore_bg=ignore_bg)

    def _actn_to_bb(self, actn: torch.Tensor, anchors: torch.Tensor):
        """Mapping activation to anchor boxes
        """
        actn_bbs = torch.tanh(actn) #map activations to [-1,1]
        actn_centers = ((actn_bbs[:, :2]/2).float() * self.grid_sizes.float()) + anchors[:,:2].float() #center relative to anchor box
        actn_hw = ((actn_bbs[:, 2:]/2).float() + 1) * anchors[:,2:].float() #get height and width relative to anchor box
        return self._hw2corners(actn_centers, actn_hw)

    def _nms(self, bbox: torch.Tensor, clas:torch.Tensor, clas_prb:torch.Tensor, thresh: float=.4):
        """Takes bbox, class and class probability vectors, and threshold for IOU
           Suppresses bboxes overlapping bboxes of the same class with sufficient IOU, keeping higher
           probability bbox. Returns an index of bboxes to keep
        """
        areas = self._box_sz(bbox) #bounding box sizes
        prb_idx = clas_prb.argsort().tolist() #probability sorted indices (ascending)
        to_keep = []
        while len(prb_idx):
            idx = prb_idx.pop() #get highest probability index
            to_keep.append(idx)
            if not len(prb_idx):
                #we're done
                break
            #calculate IOU between idx box and all others remaining in prb_idx
            y1 = torch.max(bbox[idx][0], bbox[prb_idx][:, 0]) #get the largest of the bottom ys
            x1 = torch.max(bbox[idx][1], bbox[prb_idx][:, 1]) #get the largest of the leftmost xs
            y2 = torch.max(bbox[idx][2], bbox[prb_idx][:, 2]) #get the smallest of the top ys
            x2 = torch.max(bbox[idx][3], bbox[prb_idx][:, 3]) #get the smallest of the rightmost xs
            inter = (np.maximum(y2 - y1, 0) * np.maximum(x2 - x1, 0)).to(self.device) #compute ious
            union = areas[idx] + areas[prb_idx] - inter
            ious = inter / union
            same_class = (clas[idx] == clas[prb_idx]) #flag if same class as idx
            not_bg = (clas[prb_idx] != 0) #flag if not bg
            suppress = ((ious > thresh) * same_class * not_bg).tolist() #suppress with iou > thresh, same class as idx and not background
            prb_idx = [prb_idx[i] for i,v in enumerate(suppress) if v == 0]
        return to_keep

    @staticmethod
    def _show_img(im, figsize=None, ax=None):
        """Standard func to display image
        """
        if not ax: fig,ax = plt.subplots(figsize=figsize)
        ax.imshow(im)
        ax.get_xaxis().set_visible(False)
        ax.get_yaxis().set_visible(False)
        return ax

    @staticmethod
    def _bb_hw(a: torch.Tensor):
        """Reshapes bounding box from bottom left (x,y), top right (x,y) to
           bottom left y, bottom left x, height, width
        """
        return np.array([a[1],a[0],a[3]-a[1],a[2]-a[0]])

    @staticmethod
    def _draw_outline(o, lw):
        """Adds black border around line
        """
        o.set_path_effects([patheffects.Stroke(
            linewidth=lw, foreground='black'), patheffects.Normal()])

    @staticmethod
    def _draw_rect(ax, b, color='white'):
        """Draws bounding rectangle
        """
        patch = ax.add_patch(patches.Rectangle(b[:2], *b[-2:], fill=False, edgecolor=color, lw=2))
        SingleShotDetector._draw_outline(patch, 4)

    @staticmethod
    def _draw_text(ax, xy, txt, sz=14, color='white'):
        """Draws text inside bounding rectangle
        """
        text = ax.text(*xy, txt,
            verticalalignment='top', color=color, fontsize=sz, weight='bold')
        SingleShotDetector._draw_outline(text, 1)

    @staticmethod
    def _get_cmap(N):
        """Gets colors
        """
        color_norm  = mcolors.Normalize(vmin=0, vmax=N-1)
        return cmx.ScalarMappable(norm=color_norm, cmap='Set3').to_rgba

    @staticmethod
    def _get_y(bbox: torch.tensor, clas: torch.tensor):
        """Returns tensor with 0 rows removed
        """
        if len(bbox.shape) == 1:
            bbox.unsqueeze_(0)
        bbox = (bbox + 1.) / 2. #normalize, remove from here?
        idx_keep = ((bbox[:,2]-bbox[:,0])>0).nonzero()[:,0]
        return bbox[idx_keep], clas[idx_keep]

    @staticmethod
    def _hw2corners(ctr: torch.Tensor, hw: torch.Tensor):
        """Converts array of x,y and height,width to top-left, bottom-right
        """
        return torch.cat([ctr-hw/2, ctr+hw/2], dim=1)

    @staticmethod
    def _intersection(box_a: torch.Tensor, box_b: torch.Tensor):
        """Computes intersection between two boxes
        """
        max_xy = torch.min(box_a[:, None, 2:], box_b[None, :, 2:])
        min_xy = torch.max(box_a[:, None, :2], box_b[None, :, :2])
        inter = torch.clamp((max_xy - min_xy), min=0)
        return inter[:, :, 0] * inter[:, :, 1]

    @staticmethod
    def _box_sz(b: torch.Tensor):
        """Computes box size of b
        """
        return (b[:, 2] - b[:, 0]) * (b[:, 3] - b[:, 1]) #height x width

    @staticmethod
    def _jaccard(box_a: torch.Tensor, box_b: torch.Tensor):
        """Computes jaccard distance between two boxes
        """
        box_a = box_a.float()
        box_b = box_b.float()
        inter = SingleShotDetector._intersection(box_a, box_b)
        union = (SingleShotDetector._box_sz(box_a).unsqueeze(1) +                  SingleShotDetector._box_sz(box_b).unsqueeze(0) - inter)
        return inter/union

    @staticmethod
    def _map_to_ground_truth(overlaps: torch.Tensor, verbose: bool=False):
        """Takes an array of IOU overlaps between anchor boxes and the ground truth bounding boxes
           Yields an array of overlaps by anchor box and indices of the class with which the anchor box overlaps
        """
        prior_overlap, prior_idx = overlaps.max(1)
        if verbose: print(prior_overlap)
        gt_overlap, gt_idx = overlaps.max(0)
        gt_overlap[prior_idx] = 1.99
        for i,o in enumerate(prior_idx): gt_idx[o] = i
        return gt_overlap, gt_idx

# In[ ]:

max_stars_repo_name: akshairajendran/dl-models

max_issues_repo_name: akshairajendran/dl-models

max_forks_repo_name: akshairajendran/dl-models


text:
# -*- coding: utf-8 -*-

# Copyright (c) 2016 by University of Kassel and Fraunhofer Institute for Wind Energy and Energy
# System Technology (IWES), Kassel. All rights reserved. Use of this source code is governed by a
# BSD-style license that can be found in the LICENSE file.

import numpy as np
import pandas as pd
import copy
import numbers
from collections import defaultdict
try:
    import pplog as logging
except:
    import logging

logger = logging.getLogger(__name__)

from pandapower.auxiliary import get_indices, PandapowerNet
from pandapower.create import create_empty_network

# --- Information

def lf_info(net, numv=1, numi=2):
    """
    Prints some basic information of the results in a net.
    (max/min voltage, max trafo load, max line load)
    """
    logger.info("Max voltage")
    for _, r in net.res_bus.sort_values("vm_pu", ascending=False).iloc[:numv].iterrows():
        logger.info("  %s at busidx %s (%s)", r.vm_pu, r.name, net.bus.name.at[r.name])
    logger.info("Min voltage")
    for _, r in net.res_bus.sort_values("vm_pu").iloc[:numv].iterrows():
        logger.info("  %s at busidx %s (%s)", r.vm_pu, r.name, net.bus.name.at[r.name])
    logger.info("Max loading trafo")
    if net.res_trafo is not None:
        for _, r in net.res_trafo.sort_values("loading_percent", ascending=False).iloc[
                :numi].iterrows():
            logger.info("  %s loading at trafo %s (%s)", r.loading_percent, r.name,
                        net.trafo.name.at[r.name])
    logger.info("Max loading line")
    for _, r in net.res_line.sort_values("loading_percent", ascending=False).iloc[:numi].iterrows():
        logger.info("  %s loading at line %s (%s)", r.loading_percent, r.name,
                    net.line.name.at[r.name])

def switch_info(net, sidx):
    """
    Prints what buses and elements are connected by a certain switch.
    """
    switch_type = net.switch.at[sidx, "et"]
    bidx = net.switch.at[sidx, "bus"]
    bus_name = net.bus.at[bidx, "name"]
    eidx = net.switch.at[sidx, "element"]
    if switch_type == "b":
        bus2_name = net.bus.at[eidx, "name"]
        logger.info("Switch %u connects bus %u (%s) with bus %u (%s)" % (sidx, bidx, bus_name, eidx,
                                                                         bus2_name))
    elif switch_type == "l":
        line_name = net.line.at[eidx, "name"]
        logger.info("Switch %u connects bus %u (%s) with line %u (%s)" % (sidx, bidx, bus_name, eidx,
                                                                          line_name))
    elif switch_type == "t":
        trafo_name = net.trafo.at[eidx, "name"]
        logger.info("Switch %u connects bus %u (%s) with trafo %u (%s)" % (sidx, bidx, bus_name, eidx,
                                                                           trafo_name))

def overloaded_lines(net, max_load=100):
    """
    Returns the results for all lines with loading_percent > max_load or None, if
    there are none.
    """
    if net.converged:
        return net["res_line"].index[net["res_line"]["loading_percent"] > max_load]
    else:
        raise UserWarning("The last loadflow terminated erratically, results are invalid!")

def violated_buses(net, u_min, u_max):
    """
    Returns the results for all buses with if vm_pu is not within u_max and u_min or None, if
    there are none.
    """
    if net.converged:
        return net["bus"].index[(net["res_bus"]["vm_pu"] < u_min) |
                                (net["res_bus"]["vm_pu"] > u_max)]
    else:
        raise UserWarning("The last loadflow terminated erratically, results are invalid!")

def equal_nets(x, y, check_only_results=False, tol=1.e-14):
    """
    Compares two networks. The networks are considered equal
    if they share the same keys and values, except of the
    'et' (elapsed time) entry which differs depending on
    runtime conditions and entries stating with '_'
    """
    eq = True
    eq_count = 0
    if isinstance(x, dict) or isinstance(x, PandapowerNet) and \
            isinstance(y, dict) or isinstance(y, PandapowerNet):

        # for two networks make sure both have the same keys ...
        if len(set(x.keys()) - set(y.keys())) + len(set(y.keys()) - set(x.keys())) > 0:
            logger.info("Networks entries mismatch:", list(x.keys()), " - VS. - ", list(y.keys()))
            return False

        # ... and then iter through the keys, checking for equality for each table
        for k in [k for k in x.keys() if(k != 'et' and not k.startswith("_"))]:
            if check_only_results and not k.startswith("res_"):
                continue  # skip anything thats not a result table

            eq = equal_nets(x[k], y[k], check_only_results, tol)

            # counts the mismatches and creates a list to print
            if not eq:
                if eq_count == 0:
                    logger_str = "Mismatch(es) at table(s): %s" % k
                else:
                    logger_str = logger_str + ", %s" % k
                eq_count += 1

        if eq_count:
            logger.info(logger_str)
            return False
    else:
        if isinstance(x, pd.DataFrame) and isinstance(y, pd.DataFrame):
            # for two DataFrames eval if all entries are equal
            if len(x) != len(y):
                logger.info("The number of elements of differ.")
                return False

            if len(x.dtypes) != len(y.dtypes):
                logger.info("The number of columns of differ.")
                return False

            # we use numpy.allclose to grant a tolerance
            numerical_equal = np.allclose(x.select_dtypes(include=[np.number]),
                                          y.select_dtypes(include=[np.number]),
                                          atol=tol, equal_nan=True)

            # we use pandas equals for the rest, which also evaluates NaNs to be equal
            rest_equal = x.select_dtypes(exclude=[np.number]).equals(
                y.select_dtypes(exclude=[np.number]))

            eq &= numerical_equal & rest_equal
        else:
            # else items are assumed to be single values, compare them
            if isinstance(x, numbers.Number):
                eq &= np.isclose(x, y, atol=tol)
            else:
                eq &= x == y

            # Note: if we deal with HDFStores we might encounter a problem here:
            # empty DataFrames are stored as None to HDFStore so we might be
            # comparing None to an empty DataFrame
    return eq

# --- Simulation setup and preparations

def convert_format(net):
    """
    Converts old nets to new format to ensure consistency. Converted net is returned
    """
    from pandapower.std_types import add_basic_std_types, create_std_type, parameter_from_std_type
    from pandapower.run import reset_results
    if "std_types" not in net:
        net.std_types = {"line": {}, "trafo": {}, "trafo3w": {}}
        add_basic_std_types(net)

        import os
        import json
        path, file = os.path.split(os.path.realpath(__file__))
        linedb = os.path.join(path, "linetypes.json")
        if os.path.isfile(linedb):
            with open(linedb, 'r') as f:
                lt = json.load(f)
        else:
            lt = {}
        for std_type in net.line.std_type.unique():
            if std_type in lt:
                if "shift_degree" not in lt[std_type]:
                    lt[std_type]["shift_degree"] = 0
                create_std_type(net, lt[std_type], std_type, element="line")
        trafodb = os.path.join(path, "trafotypes.json")
        if os.path.isfile(trafodb):
            with open(trafodb, 'r') as f:
                tt = json.load(f)
        else:
            tt = {}
        for std_type in net.trafo.std_type.unique():
            if std_type in tt:
                create_std_type(
                    net, tt[std_type], std_type, element="trafo")

    net.trafo.tp_side.replace(1, "hv", inplace=True)
    net.trafo.tp_side.replace(2, "lv", inplace=True)
    net.trafo.tp_side = net.trafo.tp_side.where(pd.notnull(net.trafo.tp_side), None)
    net.trafo3w.tp_side.replace(1, "hv", inplace=True)
    net.trafo3w.tp_side.replace(2, "mv", inplace=True)
    net.trafo3w.tp_side.replace(3, "lv", inplace=True)
    net.trafo3w.tp_side = net.trafo3w.tp_side.where(pd.notnull(net.trafo3w.tp_side), None)

    net["bus"] = net["bus"].rename(
        columns={'voltage_level': 'vn_kv', 'bus_type': 'type', "un_kv": "vn_kv"})
    net["bus"]["type"].replace("s", "b", inplace=True)
    net["bus"]["type"].replace("k", "n", inplace=True)
    net["line"] = net["line"].rename(columns={'vf': 'df', 'line_type': 'type'})
    net["ext_grid"] = net["ext_grid"].rename(columns={"angle_degree": "va_degree",
                                                      "ua_degree": "va_degree", "sk_max_mva": "s_sc_max_mva", "sk_min_mva": "s_sc_min_mva"})
    net["line"]["type"].replace("f", "ol", inplace=True)
    net["line"]["type"].replace("k", "cs", inplace=True)
    net["trafo"] = net["trafo"].rename(columns={'trafotype': 'std_type', "type": "std_type",
                                                "un1_kv": "vn_hv_kv", "un2_kv": "vn_lv_kv",
                                                'vfe_kw': 'pfe_kw', "unh_kv": "vn_hv_kv",
                                                "unl_kv": "vn_lv_kv", 'trafotype': 'std_type',
                                                "type": "std_type", 'vfe_kw': 'pfe_kw',
                                                "uk_percent": "vsc_percent",
                                                "ur_percent": "vscr_percent",
                                                "vnh_kv": "vn_hv_kv", "vnl_kv": "vn_lv_kv"})
    net["trafo3w"] = net["trafo3w"].rename(columns={"unh_kv": "vn_hv_kv", "unm_kv": "vn_mv_kv",
                                                    "unl_kv": "vn_lv_kv",
                                                    "ukh_percent": "vsc_hv_percent",
                                                    "ukm_percent": "vsc_mv_percent",
                                                    "ukl_percent": "vsc_lv_percent",
                                                    "urh_percent": "vscr_hv_percent",
                                                    "urm_percent": "vscr_mv_percent",
                                                    "url_percent": "vscr_lv_percent",
                                                    "vnh_kv": "vn_hv_kv", "vnm_kv": "vn_mv_kv",
                                                    "vnl_kv": "vn_lv_kv", "snh_kv": "sn_hv_kv",
                                                    "snm_kv": "sn_mv_kv", "snl_kv": "sn_lv_kv"})
    net["switch"]["type"].replace("LS", "CB", inplace=True)
    net["switch"]["type"].replace("LTS", "LBS", inplace=True)
    net["switch"]["type"].replace("TS", "DS", inplace=True)
    if "name" not in net.switch.columns:
        net.switch["name"] = None
    net["switch"] = net["switch"].rename(columns={'element_type': 'et'})
    net["ext_grid"] = net["ext_grid"].rename(columns={'voltage': 'vm_pu', "u_pu": "vm_pu",
                                                      "sk_max": "sk_max_mva", "ua_degree": "va_degree"})
    if "in_service" not in net["ext_grid"].columns:
        net["ext_grid"]["in_service"] = 1
    if "shift_mv_degree" not in net["trafo3w"].columns:
        net["trafo3w"]["shift_mv_degree"] = 0
    if "shift_lv_degree" not in net["trafo3w"].columns:
        net["trafo3w"]["shift_lv_degree"] = 0
    parameter_from_std_type(net, "shift_degree", element="trafo", fill=0)
    if "gen" not in net:
        net["gen"] = pd.DataFrame(np.zeros(0, dtype=[("name", np.dtype(object)),
                                                     ("bus", "u4"),
                                                     ("p_kw", "f8"),
                                                     ("vm_pu", "f8"),
                                                     ("sn_kva", "f8"),
                                                     ("scaling", "f8"),
                                                     ("in_service", "i8"),
                                                     ("min_q_kvar", "f8"),
                                                     ("max_q_kvar", "f8"),
                                                     ("type", np.dtype(object))]))

    if "impedance" not in net:
        net["impedance"] = pd.DataFrame(np.zeros(0, dtype=[("name", np.dtype(object)),
                                                           ("from_bus", "u4"),
                                                           ("to_bus", "u4"),
                                                           ("r_pu", "f8"),
                                                           ("x_pu", "f8"),
                                                           ("sn_kva", "f8"),
                                                           ("in_service", 'bool')]))
    if "ward" not in net:
        net["ward"] = pd.DataFrame(np.zeros(0, dtype=[("name", np.dtype(object)),
                                                      ("bus", "u4"),
                                                      ("ps_kw", "u4"),
                                                      ("qs_kvar", "f8"),
                                                      ("pz_kw", "f8"),
                                                      ("qz_kvar", "f8"),
                                                      ("in_service", "f8")]))
    if "xward" not in net:
        net["xward"] = pd.DataFrame(np.zeros(0, dtype=[("name", np.dtype(object)),
                                                       ("bus", "u4"),
                                                       ("ps_kw", "u4"),
                                                       ("qs_kvar", "f8"),
                                                       ("pz_kw", "f8"),
                                                       ("qz_kvar", "f8"),
                                                       ("r_ohm", "f8"),
                                                       ("x_ohm", "f8"),
                                                       ("vm_pu", "f8"),
                                                       ("in_service", "f8")]))
    if "shunt" not in net:
        net["shunt"] = pd.DataFrame(np.zeros(0, dtype=[("bus", "u4"),
                                                       ("name", np.dtype(object)),
                                                       ("p_kw", "f8"),
                                                       ("q_kvar", "f8"),
                                                       ("scaling", "f8"),
                                                       ("in_service", "i8")]))

    if "parallel" not in net.line:
        net.line["parallel"] = 1
    if "_empty_res_bus" not in net:
        net2 = create_empty_network()
        for key, item in net2.items():
            if key.startswith("_empty"):
                net[key] = copy.copy(item)
        reset_results(net)

    for attribute in ['tp_st_percent', 'tp_pos', 'tp_mid', 'tp_min', 'tp_max']:
        if net.trafo[attribute].dtype == 'O':
            net.trafo[attribute] = pd.to_numeric(net.trafo[attribute])
    net["gen"] = net["gen"].rename(columns={"u_pu": "vm_pu"})
    for element, old, new in [("trafo", "unh_kv", "vn_hv_kv"),
                              ("trafo", "unl_kv", "vn_lv_kv"),
                              ("trafo", "uk_percent", "vsc_percent"),
                              ("trafo", "ur_percent", "vscr_percent"),
                              ("trafo3w", "unh_kv", "vn_hv_kv"),
                              ("trafo3w", "unm_kv", "vn_mv_kv"),
                              ("trafo3w", "unl_kv", "vn_lv_kv")]:
        for std_type, parameters in net.std_types[element].items():
            if old in parameters:
                net.std_types[element][std_type][new] = net.std_types[element][std_type].pop(old)
    net.version = 1.0
    if "f_hz" not in net:
        net["f_hz"] = 50.

    if "type" not in net.load.columns:
        net.load["type"] = None
    if "zone" not in net.bus:
        net.bus["zone"] = None
    for element in ["line", "trafo", "bus", "load", "sgen", "ext_grid"]:
        net[element].in_service = net[element].in_service.astype(bool)
    net.switch.closed = net.switch.closed.astype(bool)
    return net

def add_zones_to_elements(net, elements=["line", "trafo", "ext_grid", "switch"]):
    """
    Adds zones to elements, inferring them from the zones of buses they are
    connected to.
    """
    for element in elements:
        if element == "sgen":
            net["sgen"]["zone"] = net["bus"]["zone"].loc[net["sgen"]["bus"]].values
        elif element == "load":
            net["load"]["zone"] = net["bus"]["zone"].loc[net["load"]["bus"]].values
        elif element == "ext_grid":
            net["ext_grid"]["zone"] = net["bus"]["zone"].loc[net["ext_grid"]["bus"]].values
        elif element == "switch":
            net["switch"]["zone"] = net["bus"]["zone"].loc[net["switch"]["bus"]].values
        elif element == "line":
            net["line"]["zone"] = net["bus"]["zone"].loc[net["line"]["from_bus"]].values
            crossing = sum(net["bus"]["zone"].loc[net["line"]["from_bus"]].values !=
                           net["bus"]["zone"].loc[net["line"]["to_bus"]].values)
            if crossing > 0:
                logger.warn("There have been %i lines with different zones at from- and to-bus"
                            % crossing)
        elif element == "trafo":
            net["trafo"]["zone"] = net["bus"]["zone"].loc[net["trafo"]["hv_bus"]].values
            crossing = sum(net["bus"]["zone"].loc[net["trafo"]["hv_bus"]].values !=
                           net["bus"]["zone"].loc[net["trafo"]["lv_bus"]].values)
            if crossing > 0:
                logger.warn("There have been %i trafos with different zones at lv_bus and hv_bus"
                            % crossing)
        elif element == "impedance":
            net["impedance"]["zone"] = net["bus"]["zone"].loc[net["impedance"]["from_bus"]].values
            crossing = sum(net["bus"]["zone"].loc[net["impedance"]["from_bus"]].values !=
                           net["bus"]["zone"].loc[net["impedance"]["to_bus"]].values)
            if crossing > 0:
                logger.warn("There have been %i impedances with different zones at from_bus and "
                            "to_bus" % crossing)
        elif element == "shunt":
            net["shunt"]["zone"] = net["bus"]["zone"].loc[net["shunt"]["bus"]].values
        elif element == "ward":
            net["ward"]["zone"] = net["bus"]["zone"].loc[net["ward"]["bus"]].values
        elif element == "xward":
            net["xward"]["zone"] = net["bus"]["zone"].loc[net["xward"]["bus"]].values
        else:
            raise UserWarning("Unkown element %s" % element)

def create_continuous_bus_index(net):
    """
    Creates a continuous bus index starting at zero and replaces all
    references of old indices by the new ones.
    """
    new_bus_idxs = np.arange(len(net.bus))
    bus_lookup = dict(zip(net["bus"].index.values, new_bus_idxs))
    net.bus.index = new_bus_idxs

    for element, value in [("line", "from_bus"), ("line", "to_bus"), ("trafo", "hv_bus"),
                           ("trafo", "lv_bus"), ("sgen", "bus"), ("load", "bus"),
                           ("switch", "bus"), ("ward", "bus"), ("xward", "bus"),
                           ("impedance", "from_bus"), ("impedance", "to_bus"),
                           ("shunt", "bus"), ("ext_grid", "bus")]:
        net[element][value] = get_indices(net[element][value], bus_lookup)
    return net

def set_scaling_by_type(net, scalings, scale_load=True, scale_sgen=True):
    """
    Sets scaling of loads and/or sgens according to a dictionary
    mapping type to a scaling factor. Note that the type-string is case
    sensitive.
    E.g. scaling = {"pv": 0.8, "bhkw": 0.6}

    :param net:
    :param scaling: A dictionary containing a mapping from element type to
    scaling factor.
    :return:
    """
    if not isinstance(scalings, dict):
        raise UserWarning("The parameter scaling has to be a dictionary, "
                          "see docstring")

    def scaleit(what):
        et = net[what]
        et["scaling"] = [scale[t] or s for t, s in zip(et.type.values, et.scaling.values)]

    scale = defaultdict(lambda: None, scalings)
    if scale_load:
        scaleit("load")
    if scale_sgen:
        scaleit("sgen")

# --- Modify topology

def close_switch_at_line_with_two_open_switches(net):
    """
    Finds lines that have opened switches at both ends and closes one of them.
    Function is usually used when optimizing section points to
    prevent the algorithm from ignoring isolated lines.
    """
    nl = net.switch[(net.switch.et == 'l') & (net.switch.closed == 0)]
    for i, switch in nl.groupby("element"):
        if len(switch.index) > 1:  # find all lines that have open switches at both ends
            # and close on of them
            net.switch.at[switch.index[0], "closed"] = 1

def drop_inactive_elements(net):
    """
    Drops any elements not in service AND any elements connected to inactive
    buses.
    """
    # removes inactive lines and its switches and geodata
    inactive_lines = net["line"][net["line"]["in_service"] == False].index
    drop_lines(net, inactive_lines)

    # removes inactive buses safely by deleting all connected elements and
    # its geodata as well
    disconnected_buses = (set(net["bus"].index)
                          - set(net["line"]["from_bus"])
                          - set(net["line"]["to_bus"])
                          - set(net["impedance"]["from_bus"])
                          - set(net["impedance"]["to_bus"])
                          - set(net["trafo"]["hv_bus"])
                          - set(net["trafo"]["lv_bus"])
                          - set(net["trafo3w"]["hv_bus"])
                          - set(net["trafo3w"]["mv_bus"])
                          - set(net["trafo3w"]["lv_bus"])
                          - set(net["switch"]["bus"])
                          - set(net["switch"][net["switch"]["et"] == "b"]["element"]))
    inactive_buses = (disconnected_buses | set(net["bus"].query("in_service==0").index))
    drop_buses(net, inactive_buses, also_drop_elements=True)

    # finally drops any remaining elements out of service
    for element in net.keys():
        # by iterating over all tables
        try:
            try:
                # and try to drop any elements OOS
                net[element] = net[element].query("in_service==1")
            except pd.computation.ops.UndefinedVariableError:
                # except for tables that have no "in_service"
                continue
        except AttributeError:
            # and except if net[element] is no DataFrame, and has no attribute "query"
            continue

def drop_buses(net, buses, also_drop_elements=True):
    """
    Drops buses and by default safely drops all elements connected to them as well.
    """
    # Note: lines and trafos are deleted seperately since switches may be
    # deleted along with them

    # drop any lines connected to the bus
    droplines = net["line"][(net["line"]["from_bus"].isin(buses)) | (
        net["line"]["to_bus"].isin(buses))].index
    if not also_drop_elements and len(droplines) > 0:
        raise UserWarning("Can not drop buses, "
                          "there are still lines connected")
    elif len(droplines) > 0:
        drop_lines(net, droplines)

    # drop any trafo connected to the bus
    droptrafos = net["trafo"][(net["trafo"]["hv_bus"].isin(buses)) | (
        net["trafo"]["lv_bus"].isin(buses))].index
    if not also_drop_elements and len(droptrafos) > 0:
        raise UserWarning("Can not drop buses, "
                          "there are still trafos connected")
    elif len(droptrafos) > 0:
        drop_trafos(net, droptrafos)

    # drop any other elements connected to the bus
    for element, value in [("sgen", "bus"), ("load", "bus"), ("switch", "bus"),
                           ("ext_grid", "bus"), ("shunt", "bus"),
                           ("impedance", "from_bus"), ("impedance", "to_bus"),
                           ("ward", "bus"), ("xward", "bus")]:
        i = net[element][net[element][value].isin(buses)].index
        if not also_drop_elements and i:
            raise UserWarning("Can not drop buses, "
                              "there are still elements connected")
        net[element].drop(i, inplace=True)

    # drop busbus switches
    i = net["switch"].index[
        (net["switch"]["element"].isin(buses)) & (net["switch"]["et"] == "b")]
    net["switch"].drop(i, inplace=True)

    # finally drop buses and their geodata
    net["bus"].drop(buses, inplace=True)
    net["bus_geodata"].drop(set(buses) & set(
        net["bus_geodata"].index), inplace=True)

def drop_trafos(net, trafos):
    """
    Deletes all trafos and in the given list of indices and removes
    any switches connected to it.
    """
    # drop any switches
    i = net["switch"].index[
        (net["switch"]["element"].isin(trafos)) & (net["switch"]["et"] == "t")]
    net["switch"].drop(i, inplace=True)

    # drop the lines+geodata
    net["trafo"].drop(trafos, inplace=True)

def drop_lines(net, lines):
    """
    Deletes all lines and their geodata in the given list of indices and removes
    any switches connected to it.
    """
    # drop any switches
    i = net["switch"].index[
        (net["switch"]["element"].isin(lines)) & (net["switch"]["et"] == "l")]
    net["switch"].drop(i, inplace=True)

    # drop the lines+geodata
    net["line"].drop(lines, inplace=True)
    net["line_geodata"].drop(set(lines) & set(
        net["line_geodata"].index), inplace=True)

def fuse_buses(net, b1, b2, drop=True):
    """
    Reroutes any connections to buses in b2 to the given bus b1. Additionally drops the buses b2,
    if drop=True (default).
    """
    try:
        b2.__iter__
    except:
        b2 = [b2]

    for element, value in [("line", "from_bus"), ("line", "to_bus"), ("impedance", "from_bus"),
                           ("impedance", "to_bus"), ("trafo", "hv_bus"), ("trafo", "lv_bus"),
                           ("sgen", "bus"), ("load", "bus"),
                           ("switch", "bus"), ("ext_grid", "bus"),
                           ("ward", "bus"), ("xward", "bus"),
                           ("shunt", "bus")]:
        i = net[element][net[element][value].isin(b2)].index
        net[element].loc[i, value] = b1

    i = net["switch"][(net["switch"]["et"] == 'b') & (
        net["switch"]["element"].isin(b2))].index
    net["switch"].loc[i, "element"] = b1
    net["switch"].drop(net["switch"][(net["switch"]["bus"] == net["switch"]["element"]) &
                                     (net["switch"]["et"] == "b")].index, inplace=True)
    if drop:
        net["bus"].drop(b2, inplace=True)
    return net

def set_element_status(net, buses, in_service):
    """
    Sets buses and all elements connected to them out of service.
    """
    net.bus.loc[buses, "in_service"] = in_service

    lines = net.line[(net.line.from_bus.isin(buses)) |
                     (net.line.to_bus.isin(buses))].index
    net.line.loc[lines, "in_service"] = in_service

    trafo = net.trafo[(net.trafo.hv_bus.isin(buses)) |
                      (net.trafo.lv_bus.isin(buses))].index
    net.trafo.loc[trafo, "in_service"] = in_service

    impedances = net.impedance[(net.impedance.from_bus.isin(buses)) |
                               (net.impedance.to_bus.isin(buses))].index
    net.impedance.loc[impedances, "in_service"] = in_service

    loads = net.load[net.load.bus.isin(buses)].index
    net.load.loc[loads, "in_service"] = in_service

    sgens = net.sgen[net.sgen.bus.isin(buses)].index
    net.sgen.loc[sgens, "in_service"] = in_service

    wards = net.ward[net.ward.bus.isin(buses)].index
    net.ward.loc[wards, "in_service"] = in_service

    xwards = net.xward[net.xward.bus.isin(buses)].index
    net.xward.loc[xwards, "in_service"] = in_service

    shunts = net.shunt[net.shunt.bus.isin(buses)].index
    net.shunt.loc[shunts, "in_service"] = in_service

def select_subnet(net, buses, include_switch_buses=False, include_results=False,
                  keep_everything_else=False):
    """
    Selects a subnet by a list of bus indices and returns a net with all elements
    connected to them.
    """
    buses = set(buses)
    if include_switch_buses:
        # we add both buses of a connected line, the one selected is not switch.bus

        # for all line switches
        for _, s in net["switch"].query("et=='l'").iterrows():
            # get from/to-bus of the connected line
            fb = net["line"]["from_bus"].at[s["element"]]
            tb = net["line"]["to_bus"].at[s["element"]]
            # if one bus of the line is selected and its not the switch-bus, add the other bus
            if fb in buses and s["bus"] != fb:
                buses.add(tb)
            if tb in buses and s["bus"] != tb:
                buses.add(fb)

    p2 = create_empty_network()

    p2.bus = net.bus.loc[buses]
    p2.ext_grid = net.ext_grid[net.ext_grid.bus.isin(buses)]
    p2.load = net.load[net.load.bus.isin(buses)]
    p2.sgen = net.sgen[net.sgen.bus.isin(buses)]
    p2.gen = net.gen[net.gen.bus.isin(buses)]
    p2.shunt = net.shunt[net.shunt.bus.isin(buses)]
    p2.ward = net.ward[net.ward.bus.isin(buses)]
    p2.xward = net.xward[net.xward.bus.isin(buses)]

    p2.line = net.line[(net.line.from_bus.isin(buses)) & (net.line.to_bus.isin(buses))]
    p2.trafo = net.trafo[(net.trafo.hv_bus.isin(buses)) & (net.trafo.lv_bus.isin(buses))]
    p2.impedance = net.impedance[(net.impedance.from_bus.isin(buses)) &
                                 (net.impedance.to_bus.isin(buses))]

    if include_results:
        for table in net.keys():
            if net[table] is None:
                continue
            elif table == "res_bus":
                p2[table] = net[table].loc[buses]
            elif table.startswith("res_"):
                p2[table] = net[table].loc[p2[table.split("res_")[1]].index]
    if "bus_geodata" in net:
        p2["bus_geodata"] = net["bus_geodata"].loc[net["bus_geodata"].index.isin(buses)]
    if "line_geodata" in net:
        lines = p2.line.index
        p2["line_geodata"] = net["line_geodata"].loc[net["line_geodata"].index.isin(lines)]

    # switches
    si = [i for i, s in net["switch"].iterrows()
          if s["bus"] in buses and
          ((s["et"] == "b" and s["element"] in p2["bus"].index) or
           (s["et"] == "l" and s["element"] in p2["line"].index) or
           (s["et"] == "t" and s["element"] in p2["trafo"].index))]
    p2["switch"] = net["switch"].loc[si]
    # return a PandapowerNet
    if keep_everything_else:
        newnet = copy.deepcopy(net)
        newnet.update(p2)
        return PandapowerNet(newnet)
    return PandapowerNet(p2)

# --- item/element selections

def get_element_index(net, element, name, exact_match=True, regex=False):
    """
    Returns the element identified by the element string and a name.

    INPUT:
      **net** - pandapower network

      **element** - line indices of lines that are considered. If None, all lines in the network are \
                  considered.

      **name** - if True, switches are also considered as connected elements.

    OUTPUT:

      **connections** - pandapower Series with number of connected elements for each bus. If a \
                      selection of lines is given, only buses connected to these lines are in \
                      the returned series.
    """
    if exact_match:
        idx = net[element][net[element]["name"] == name].index
        if len(idx) == 0:
            raise UserWarning("There is no %s with name %s" % (element, name))
        if len(idx) > 1:
            raise UserWarning("Duplicate %s names for %s" % (element, name))
    else:
        return net[element][net[element]["name"].str.contains(name, regex=regex)].index
    return idx[0]

def next_bus(net, bus, element_id, et='line', **kwargs):
    """
    Returns the index of the second bus an element is connected to, given a
    first one. E.g. the from_bus given the to_bus of a line.
    """
    # for legacy compliance
    if "element_type" in kwargs:
        et = kwargs["element_type"]

    if et == 'line':
        bc = ["from_bus", "to_bus"]
    elif et == 'trafo':
        bc = ["hv_bus", "lv_bus"]
    elif et == "switch":
        bc = ["bus", "element"]
    else:
        raise Exception("unknown element type")
    nb = list(net[et].loc[element_id, bc].values)
    nb.remove(bus)
    return nb[0]

def get_connected_elements(net, element, buses, respect_switches=True, respect_in_service=False):
    """
     Returns elements connected to a given bus.

     INPUT:

        **net** (PandapowerNet)

        **element** (string, name of the element table)

        **buses** (single integer or iterable of ints)

     OPTIONAL:

        **respect_switches** (boolean, True)    - True: open switches will be respected
                                                  False: open switches will be ignored
        **respect_in_service** (boolean, False) - True: in_service status of connected lines will be
                                                        respected
                                                  False: in_service status will be ignored
     RETURN:

        **cl** (set) - Returns connected lines.

    """

    if not hasattr(buses, "__iter__"):
        buses = [buses]

    if element in ["line", "l"]:
        element = "l"
        element_table = net.line
        connected_elements = set(net.line.index[net.line.from_bus.isin(buses)
                                                | net.line.to_bus.isin(buses)])

    elif element in ["trafo"]:
        element = "t"
        element_table = net.trafo
        connected_elements = set(net["trafo"].index[(net.trafo.hv_bus.isin(buses))
                                                    | (net.trafo.lv_bus.isin(buses))])
    elif element in ["trafo3w", "t3w"]:
        element = "t3w"
        element_table = net.trafo3w
        connected_elements = set(net["trafo3w"].index[(net.trafo3w.hv_bus.isin(buses))
                                                      | (net.trafo3w.mv_bus.isin(buses))
                                                      | (net.trafo3w.lv_bus.isin(buses))])
    elif element == "impedance":
        element_table = net.impedance
        connected_elements = set(net["impedance"].index[(net.impedance.from_bus.isin(buses))
                                                        | (net.impedance.to_bus.isin(buses))])
    elif element == "load":
        element_table = net.load
        connected_elements = set(element_table.index[(element_table.bus.isin(buses))])
    elif element == "sgen":
        element_table = net.sgen
        connected_elements = set(element_table.index[(element_table.bus.isin(buses))])
    elif element == "ward":
        element_table = net.ward
        connected_elements = set(element_table.index[(element_table.bus.isin(buses))])
    elif element == "shunt":
        element_table = net.shunt
        connected_elements = set(element_table.index[(element_table.bus.isin(buses))])
    elif element == "xward":
        element_table = net.xward
        connected_elements = set(element_table.index[(element_table.bus.isin(buses))])
    elif element == "ext_grid":
        element_table = net.ext_grid
        connected_elements = set(element_table.index[(element_table.bus.isin(buses))])
    elif element == "gen":
        element_table = net.gen
        connected_elements = set(element_table.index[(element_table.bus.isin(buses))])
    else:
        raise UserWarning("Unknown element! ", element)

    if respect_switches and element in ["l", "t", "t3w"]:
        open_switches = get_connected_switches(net, buses, consider=element, status="open")
        if open_switches:
            open_and_connected = net.switch.loc[net.switch.index.isin(open_switches)
                                                & net.switch.element.isin(connected_elements)].index
            connected_elements -= set(net.switch.element[open_and_connected])

    if respect_in_service:
        connected_elements -= set(element_table[element_table.in_service == False].index)

    return connected_elements

def get_connected_buses(net, buses, consider=("l", "s", "t"), respect_switches=True, respect_in_service=False):
    """
     Returns buses connected to given buses. The source buses will NOT be returned.

     INPUT:

        **net** (PandapowerNet)

        **buses** (single integer or iterable of ints)

     OPTIONAL:

        **respect_switches** (boolean, True)        - True: open switches will be respected
                                                      False: open switches will be ignored
        **respect_in_service** (boolean, False)     - True: in_service status of connected buses
                                                            will be respected
                                                            False: in_service status will be
                                                            ignored
        **consider** (iterable, ("l", "s", "t"))    - Determines, which types of connections will
                                                      be will be considered.
                                                      l: lines
                                                      s: switches
                                                      t: trafos
     RETURN:

        **cl** (set) - Returns connected buses.

    """
    if not hasattr(buses, "__iter__"):
        buses = [buses]

    cb = set()
    if "l" in consider:
        cl = get_connected_elements(net, "line", buses, respect_switches, respect_in_service)
        cb |= set(net.line[net.line.index.isin(cl)].from_bus)
        cb |= set(net.line[net.line.index.isin(cl)].to_bus)

    if "s" in consider:
        cs = get_connected_switches(net, buses, consider='b',
                                    status="closed" if respect_switches else "all")
        cb |= set(net.switch[net.switch.index.isin(cs)].element)
        cb |= set(net.switch[net.switch.index.isin(cs)].bus)

    if "t" in consider:
        ct = get_connected_elements(net, "trafo", buses, respect_switches, respect_in_service)
        cb |= set(net.trafo[net.trafo.index.isin(ct)].lv_bus)
        cb |= set(net.trafo[net.trafo.index.isin(cl)].hv_bus)

    if respect_in_service:
        cb -= set(net.bus[~net.bus.in_service].index)

    return cb - set(buses)

def get_connected_buses_at_element(net, element, et, respect_in_service=False):
    """
     Returns buses connected to a given line, switch or trafo. In case of a bus switch, two buses
     will be returned, else one.

     INPUT:

        **net** (PandapowerNet)

        **element** (integer)

        **et** (string)                             - Type of the source element:
                                                      l: line
                                                      s: switch
                                                      t: trafo

     OPTIONAL:

        **respect_in_service** (boolean, False)     - True: in_service status of connected buses
                                                            will be respected
                                                      False: in_service status will be ignored
     RETURN:

        **cl** (set) - Returns connected switches.

    """

    cb = set()
    if et == 'l':
        cb.add(net.line.from_bus.at[element])
        cb.add(net.line.to_bus.at[element])

    elif et == 's':
        cb.add(net.switch.bus.at[element])
        if net.switch.et.at[element] == 'b':
            cb.add(net.switch.element.at[element])
    elif et == 't':
        cb.add(net.trafo.hv_bus.at[element])
        cb.add(net.trafo.lv_bus.at[element])

    if respect_in_service:
        cb -= set(net.bus[~net.bus.in_service].index)

    return cb

def get_connected_switches(net, buses, consider=('b', 'l', 't'), status="all"):
    """
    Returns switches connected to given buses.

    INPUT:

        **net** (PandapowerNet)

        **buses** (single integer or iterable of ints)

    OPTIONAL:

        **respect_switches** (boolean, True)        - True: open switches will be respected
                                                     False: open switches will be ignored

        **respect_in_service** (boolean, False)     - True: in_service status of connected
                                                            buses will be respected

                                                      False: in_service status will be ignored
        **consider** (iterable, ("l", "s", "t"))    - Determines, which types of connections
                                                      will be will be considered.
                                                      l: lines
                                                      s: switches
                                                      t: trafos

        **status** (string, ("all", "closed", "open"))    - Determines, which switches will
                                                            be considered
    RETURN:

       **cl** (set) - Returns connected buses.

    """

    if not hasattr(buses, "__iter__"):
        buses = [buses]

    if status == "closed":
        switch_selection = net.switch.closed == True
    elif status == "open":
        switch_selection = net.switch.closed == False
    elif status == "all":
        switch_selection = np.full(len(net.switch), True, dtype=bool)
    else:
        logger.warn("Unknown switch status \"%s\" selected! "
                    "Selecting all switches by default." % status)

    cs = set()
    if 'b' in consider:
        cs |= set(net['switch'].index[(net['switch']['bus'].isin(buses)
                                       | net['switch']['element'].isin(buses))
                                      & (net['switch']['et'] == 'b')
                                      & switch_selection])
    if 'l' in consider:
        cs |= set(net['switch'].index[(net['switch']['bus'].isin(buses))
                                      & (net['switch']['et'] == 'l')
                                      & switch_selection])

    if 't' in consider:
        cs |= set(net['switch'].index[net['switch']['bus'].isin(buses)
                                      & (net['switch']['et'] == 't')
                                      & switch_selection])

    return cs

ext: py

max_stars_repo_name: lucassm/cigre-montecarlo

max_issues_repo_name: lucassm/cigre-montecarlo

max_forks_repo_name: lucassm/cigre-montecarlo

reason: import numpy


text:
import numpy as np
from skimage.transform import resize
from skimage.io import imread

def resize_image(img):
    im = resize(img, (120, 160, 3))
    im_arr = im.reshape((120, 160, 3))
    return im_arr

print('Preparing data')

X = []
Y = []

image_files = np.loadtxt('data/inputs.csv', delimiter=';', dtype=str, usecols=(0,))
joystick_values = np.loadtxt('data/inputs.csv', delimiter=';', usecols=(1,2))

Y.append(joystick_values)
Y = np.concatenate(Y)

for image_file in image_files:
    image = imread(image_file)
    vec = resize_image(image)
    X.append(vec)

X = np.asarray(X)

print('Saving data')

np.save("data/X", X)
np.save("data/Y", Y)

print('Done!')

max_stars_repo_name: mudziok/TensorKart

max_issues_repo_name: mudziok/TensorKart

max_forks_repo_name: mudziok/TensorKart


text:
"""
Defines a PDE class whose right hand side is given as a string

.. codeauthor:: David Zwicker <david.zwicker@ds.mpg.de>
"""

import re
from collections import OrderedDict
from typing import Any, Callable, Dict

import numpy as np

from ..fields import FieldCollection, VectorField
from ..fields.base import DataFieldBase, FieldBase, OptionalArrayLike
from ..grids.boundaries.axes import BoundariesData
from ..pdes.base import PDEBase
from ..tools.docstrings import fill_in_docstring
from ..tools.numba import jit, nb

class PDE(PDEBase):
    """PDE defined by mathematical expressions

    Attributes:
        variables (tuple):
            The name of the variables (i.e., fields) in the order they are
            expected to appear in the `state`.
        diagnostics (dict):
            Additional diagnostic information that might help with analyzing
            problems, e.g., when :mod:`sympy` cannot parse or :mod`numba` cannot
            compile a function.
    """

    @fill_in_docstring
    def __init__(
        self,
        rhs: "OrderedDict[str, str]",
        noise: OptionalArrayLike = 0,
        bc: BoundariesData = "natural",
        bc_ops: "OrderedDict[str, BoundariesData]" = None,
        user_funcs: Dict[str, Any] = None,
    ):
        """
        Warning:
            {WARNING_EXEC}

        Args:
            rhs (OrderedDict):
                The expressions defining the evolution rate. The dictionary keys define
                the name of the fields whose evolution is considered, while the values
                specify their evolution rate as a string that can be parsed by
                :mod:`sympy`. These expression may contain variables (i.e., the fields
                themselves, spatial coordinates of the grid, and `t` for the time),
                standard local mathematical operators defined by sympy, and the
                operators defined in the :mod:`pde` package. Note that operators need to
                be specified with their full name, i.e., `laplace` for a scalar
                Laplacian and `vector_laplace` for a Laplacian operating on a vector
                field. Moreover, the dot product between two vector fields can be
                denoted by using `dot(field1, field2)` in the expression, while an outer
                product is calculated using `outer(field1, field2)`.
            noise (float or :class:`numpy.ndarray`):
                Magnitude of additive Gaussian white noise. The default value of zero
                implies deterministic partial differential equations will be solved.
                Different noise magnitudes can be supplied for each field in coupled
                PDEs.
            bc:
                Boundary conditions for the operators used in the expression. The
                conditions here are applied to all operators that do not have a
                specialized condition given in `bc_ops`.
                {ARG_BOUNDARIES}
            bc_ops (dict):
                Special boundary conditions for some operators. The keys in this
                dictionary specify where the boundary condition will be applied.
                The keys follow the format "VARIABLE:OPERATOR", where VARIABLE specifies
                the expression in `rhs` where the boundary condition is applied to the
                operator specified by OPERATOR. For both identifiers, the wildcard
                symbol "*" denotes that all fields and operators are affected,
                respectively.
            user_funcs (dict, optional):
                A dictionary with user defined functions that can be used in the
                expressions in `rhs`.

        Note:
            The order in which the fields are given in `rhs` defines the order in which
            they need to appear in the `state` variable when the evolution rate is
            calculated. Note that `dict` keep the insertion order since Python version
            3.7, so a normal dictionary can be used to define the equations.
        """
        from sympy.core.function import AppliedUndef

        from ..tools.expressions import ScalarExpression

        super().__init__(noise=noise)

        # validate input
        if not isinstance(rhs, OrderedDict):
            rhs = OrderedDict(rhs)
        if "t" in rhs:
            raise RuntimeError("`t` is not allowed as a variable since it denotes time")

        # turn the expression strings into sympy expressions
        self._rhs_expr, self._operators = {}, {}
        explicit_time_dependence = False
        complex_valued = False
        for var, rhs_str in rhs.items():
            rhs_expr = ScalarExpression(rhs_str, user_funcs=user_funcs)
            self._rhs_expr[var] = rhs_expr

            if rhs_expr.depends_on("t"):
                explicit_time_dependence = True
            if rhs_expr.complex:
                complex_valued = True

            # determine undefined functions in the expression
            self._operators[var] = {
                func.__class__.__name__
                for func in rhs_expr._sympy_expr.atoms(AppliedUndef)
                if func.__class__.__name__ not in rhs_expr.user_funcs
            }

        # set public instance attributes
        self.rhs = rhs
        self.variables = tuple(rhs.keys())
        self.explicit_time_dependence = explicit_time_dependence
        self.complex_valued = complex_valued
        operators = frozenset().union(*self._operators.values())  # type: ignore

        # setup boundary conditions
        if bc_ops is None:
            bcs = {"*:*": bc}
        else:
            bcs = OrderedDict(bc_ops)
            if "*:*" in bcs and bc != "natural":
                self._logger.warning("Two default boundary conditions.")
            bcs["*:*"] = bc  # append default boundary conditions

        self.bcs: Dict[str, Any] = {}
        for key_str, value in bcs.items():
            # split on . and :
            parts = re.split(r"\.|:", key_str)
            if len(parts) == 1:
                if len(self.variables):
                    key = f"{self.variables[0]}:{key_str}"
                else:
                    raise ValueError(
                        f'Boundary condition "{key_str}" is ambiguous. Use format '
                        '"VARIABLE:OPERATOR" instead.'
                    )
            elif len(parts) == 2:
                key = ":".join(parts)
            else:
                raise ValueError(f'Cannot parse boundary condition "{key_str}"')
            if key in self.bcs:
                self._logger.warning(f"Two boundary conditions for key {key}")
            self.bcs[key] = value

        # save information for easy inspection
        self.diagnostics = {
            "variables": self.variables,
            "explicit_time_dependence": explicit_time_dependence,
            "complex_valued_rhs": complex_valued,
            "operators": operators,
        }
        self._cache: Dict[str, Dict[str, Any]] = {}

    @property
    def expressions(self) -> Dict[str, str]:
        """ show the expressions of the PDE """
        return {k: v.expression for k, v in self._rhs_expr.items()}

    def _prepare_cache(
        self, state: FieldBase, backend: str = "numpy"
    ) -> Dict[str, Any]:
        """prepare the expression by setting internal variables in the cache

        Note that the expensive calculations in this method are only carried
        out if the state attributes change.

        Args:
            state (:class:`~pde.fields.FieldBase`):
                The field describing the state of the PDE
        """
        # check the cache
        cache = self._cache.get(backend, {})
        if state.attributes == cache.get("state.attributes", None):
            return cache  # this cache was already prepared
        cache = self._cache[backend] = {}  # clear cache, if there was any

        # check whether the state is compatible with the PDE
        num_fields = len(self.variables)
        self.diagnostics["num_fields"] = num_fields
        if isinstance(state, FieldCollection):
            if num_fields != len(state):
                raise ValueError(
                    f"Expected {num_fields} fields in state, but got {len(state)} ones"
                )
        elif isinstance(state, DataFieldBase):
            if num_fields != 1:
                raise ValueError(
                    f"Expected {num_fields} fields in state, but got only one"
                )
        else:
            raise ValueError(f"Unknown state class {state.__class__.__name__}")

        # obtain functions used in the expression
        ops_general = {}

        # create special operators if necessary
        if "dot" in self.diagnostics["operators"]:  # type: ignore
            # add dot product between two vector fields. This can for instance
            # appear when two gradients of scalar fields need to be multiplied
            ops_general["dot"] = VectorField(state.grid).make_dot_operator()

        if "inner" in self.diagnostics["operators"]:  # type: ignore
            # synonym for dot product operator
            ops_general["inner"] = VectorField(state.grid).make_dot_operator()

        if "outer" in self.diagnostics["operators"]:  # type: ignore
            # synonym for dot product operator
            ops_general["outer"] = VectorField(state.grid).make_outer_prod_operator()

        # obtain the python functions for the rhs
        rhs_funcs = []
        for var in self.variables:
            ops = ops_general.copy()

            # obtain the (differential) operators
            for func in self._operators[var]:
                if func in ops:
                    continue
                # determine boundary conditions
                for bc_key, bc in self.bcs.items():
                    bc_var, bc_func = bc_key.split(":")
                    var_match = bc_var == var or bc_var == "*"
                    func_match = bc_func == func or bc_func == "*"
                    if var_match and func_match:
                        break  # found a matching boundary condition
                else:
                    raise RuntimeError(
                        "Could not find suitable boundary condition for function "
                        f"`{func}` applied in equation for `{var}`"
                    )

                try:
                    ops[func] = state.grid.get_operator(func, bc=bc)
                except ValueError:
                    self._logger.info(
                        "Operator %s was not defined, so we assume that sympy knows it",
                        func,
                    )

            # obtain the function to calculate the right hand side
            expr = self._rhs_expr[var]
            signature = self.variables + ("t",)

            def _get_expr_func(signature):
                """ helper function obtaining expression and checking the signature """
                extra_vars = set(expr.vars) - set(signature)
                if extra_vars:
                    extra_vars_str = ", ".join(sorted(extra_vars))
                    raise RuntimeError(
                        f"Undefined variable in expression: {extra_vars_str}"
                    )
                expr.vars = signature
                if backend == "numpy":
                    result = expr._get_function(single_arg=False, user_funcs=ops)
                elif backend == "numba":
                    result = expr._get_function(
                        single_arg=False, user_funcs=ops, prepare_compilation=True
                    )
                else:
                    raise ValueError(f"Unsupported backend {backend}")
                return result

            if any(expr.depends_on(c) for c in state.grid.axes):
                # expression has a spatial dependence, too

                # determine and check the signature
                inner_func = jit(_get_expr_func(signature + tuple(state.grid.axes)))

                # inject the spatial coordinates into the expression for the rhs
                coords_tuple = tuple(  # @UnusedVariable
                    state.grid.cell_coords[..., i] for i in range(state.grid.num_axes)
                )

                def rhs_func(*args):
                    """ wrapper that inserts the spatial variables """
                    return inner_func(*args, *coords_tuple)

            else:
                # expression only depends on the actual variables
                rhs_func = _get_expr_func(signature)

            rhs_funcs.append(rhs_func)

        cache["rhs_funcs"] = rhs_funcs

        # add extra information for field collection
        if isinstance(state, FieldCollection):
            # isscalar be False even if start == stop (e.g. vector fields)
            isscalar = tuple(field.rank == 0 for field in state)
            starts = tuple(slc.start for slc in state._slices)
            stops = tuple(slc.stop for slc in state._slices)

            def get_data_tuple(state_data):
                """ helper for turning state_data into a tuple of field data """
                return tuple(
                    state_data[starts[i]]
                    if isscalar[i]
                    else state_data[starts[i] : stops[i]]
                    for i in range(num_fields)
                )

            cache["get_data_tuple"] = get_data_tuple

        # store the attributes in the cache, which allows to later circumvent
        # calculating the quantities above again. Note that this has to be the
        # last expression of the method, so the cache is only valid when the
        # prepare function worked successfully
        cache["state_attributes"] = state.attributes
        return cache

    def evolution_rate(self, state: FieldBase, t: float = 0) -> FieldBase:
        """evaluate the right hand side of the PDE

        Args:
            state (:class:`~pde.fields.FieldBase`):
                The field describing the state of the PDE
            t (float):
                The current time point

        Returns:
            :class:`~pde.fields.FieldBase`:
            Field describing the evolution rate of the PDE
        """
        cache = self._prepare_cache(state, backend="numpy")

        if isinstance(state, DataFieldBase):
            rhs = cache["rhs_funcs"][0]
            return state.copy(data=rhs(state.data, t))

        elif isinstance(state, FieldCollection):
            result = state.copy()
            for i in range(len(state)):
                data_tpl = cache["get_data_tuple"](state.data)
                result[i].data = cache["rhs_funcs"][i](*data_tpl, t)
            return result

        else:
            raise TypeError(f"Unsupported field {state.__class__.__name__}")

    def _make_pde_rhs_numba_coll(
        self, state: FieldCollection, cache: Dict[str, Any]
    ) -> Callable:
        """create the compiled rhs if `state` is a field collection

        Args:
            state (:class:`~pde.fields.FieldCollection`):
                An example for the state defining the grid and data types
            cache (dict):
                Cached information that will be used in the function. The cache is
                populated by :meth:`PDE._prepare_cache`.

        Returns:
            A function with signature `(state_data, t)`, which can be called
            with an instance of :class:`numpy.ndarray` of the state data and
            the time to obtained an instance of :class:`numpy.ndarray` giving
            the evolution rate.
        """
        num_fields = len(state)
        data_shape = state.data.shape
        rhs_list = tuple(jit(cache["rhs_funcs"][i]) for i in range(num_fields))

        starts = tuple(slc.start for slc in state._slices)
        stops = tuple(slc.stop for slc in state._slices)

        # In the future, the following should be possible:
        #         @jit
        #         def evolution_rate(state_data, t, out):
        #             """ evolve all agents explicitly """
        #             for i in nb.literal_unroll(range(num_fields)):
        #                 out[i] = rhs_list[i](*state_data, t)
        #         return evolver

        get_data_tuple = cache["get_data_tuple"]

        def chain(i=0, inner=None):
            """ recursive helper function for applying all rhs """
            # run through all functions
            rhs = rhs_list[i]

            if inner is None:
                # the innermost function does not need to call a child
                @jit
                def wrap(data_tpl, t, out):
                    out[starts[i] : stops[i]] = rhs(*data_tpl, t)

            else:
                # all other functions need to call one deeper in the chain
                @jit
                def wrap(data_tpl, t, out):
                    inner(data_tpl, t, out)
                    out[starts[i] : stops[i]] = rhs(*data_tpl, t)

            if i < num_fields - 1:
                # there are more items in the chain
                return chain(i + 1, inner=wrap)
            else:
                # this is the outermost function
                @jit
                def evolution_rate(state_data: np.ndarray, t: float = 0):
                    out = np.empty(data_shape)
                    with nb.objmode():
                        data_tpl = get_data_tuple(state_data)
                        wrap(data_tpl, t, out)
                    return out

                return evolution_rate

        # compile the recursive chain
        return chain()  # type: ignore

    def _make_pde_rhs_numba(self, state: FieldBase) -> Callable:
        """create a compiled function evaluating the right hand side of the PDE

        Args:
            state (:class:`~pde.fields.FieldBase`):
                An example for the state defining the grid and data types

        Returns:
            A function with signature `(state_data, t)`, which can be called
            with an instance of :class:`numpy.ndarray` of the state data and
            the time to obtained an instance of :class:`numpy.ndarray` giving
            the evolution rate.
        """
        cache = self._prepare_cache(state, backend="numba")

        if isinstance(state, DataFieldBase):
            return jit(cache["rhs_funcs"][0])  # type: ignore
        elif isinstance(state, FieldCollection):
            return self._make_pde_rhs_numba_coll(state, cache)
        else:
            raise TypeError(f"Unsupported field {state.__class__.__name__}")

max_stars_repo_name: volpatto/py-pde

max_issues_repo_name: volpatto/py-pde

max_forks_repo_name: volpatto/py-pde

reason: import numpy,from sympy


text:
from __future__ import print_function
import argparse
import torch
import torch.nn as nn
import torch.nn.functional as F
import torch.optim as optim
from torchvision import datasets, transforms
from torch.optim.lr_scheduler import StepLR
import csv
from tensorboardX import SummaryWriter
from datetime import datetime
import math
import numpy as np

now = datetime.now()
current_time = now.strftime("%H:%M:%S")
writer = SummaryWriter(logdir=f'ae_log/ae_log_{current_time}')

activation = nn.LeakyReLU
# authors use this initializer, but it doesn't seem essential
def Initializer(layers, slope=0.2):
    for layer in layers:
        if hasattr(layer, 'weight'):
            w = layer.weight.data
            std = 1/np.sqrt((1 + slope**2) * np.prod(w.shape[:-1]))
            w.normal_(std=std)
        if hasattr(layer, 'bias'):
            layer.bias.data.zero_()

def Encoder(scales, depth, latent, colors):
    layers = []
    layers.append(nn.Conv2d(colors, depth, 1, padding=1))
    kp = depth
    for scale in range(scales):
        k = depth << scale
        layers.extend([nn.Conv2d(kp, k, 3, padding=1), activation()])
        layers.extend([nn.Conv2d(k, k, 3, padding=1), activation()])
        layers.append(nn.AvgPool2d(2))
        kp = k
    k = depth << scales
    layers.extend([nn.Conv2d(kp, k, 3, padding=1), activation()])
    layers.append(nn.Conv2d(k, latent, 3, padding=1))
    Initializer(layers)
    return nn.Sequential(*layers)

def Decoder(scales, depth, latent, colors):
    layers = []
    kp = latent
    for scale in range(scales - 1, -1, -1):
        k = depth << scale
        layers.extend([nn.Conv2d(kp, k, 3, padding=1), activation()])
        layers.extend([nn.Conv2d(k, k, 3, padding=1), activation()])
        layers.append(nn.Upsample(scale_factor=2))
        kp = k
    layers.extend([nn.Conv2d(kp, depth, 3, padding=1), activation()])
    layers.append(nn.Conv2d(depth, colors, 3, padding=1))
    Initializer(layers)
    return nn.Sequential(*layers)

def train(args, encoder, decoder, device, train_loader, optimizer, epoch):
    encoder.train()
    decoder.train()
    train_ae_loss = 0.0
    for batch_idx, (data, target) in enumerate(train_loader):
        data, target = data.to(device), target.to(device)
        optimizer.zero_grad()
        h = encoder(data)
        output = decoder(h)
        ae_loss = F.mse_loss(output, data, reduction='mean')
        loss = ae_loss
        loss.backward()
        optimizer.step()
        train_ae_loss += loss.item()
        if batch_idx % args.log_interval == 0:
            #print('Train Epoch: {} [{}/{} ({:.0f}%)]\tLoss: {:.6f}'.format(
            #    epoch, batch_idx * len(data), len(train_loader.dataset),
            #	100. * batch_idx / len(train_loader), loss.item()))
            if args.dry_run:
                break

    writer.add_scalar('data/train_ae_loss', train_ae_loss, epoch)
    print('Train set: Average ae loss: {:.8f})'.format(train_ae_loss))

def test(encoder, decoder, device, test_loader, epoch):
    encoder.eval()
    decoder.eval()
    test_ae_loss = 0
    correct = 0
    #ce_loss = nn.CrossEntropyLoss(reduction='mean')
    with torch.no_grad():
        for data, target in test_loader:
            data, target = data.to(device), target.to(device)
            h = encoder(data)
            output = decoder(h)
            test_ae_loss += F.mse_loss(output, data, reduction='mean').item()  # sum up batch loss
            #pred = output.argmax(dim=1, keepdim=True)  # get the index of the max log-probability
            #correct += pred.eq(target.view_as(pred)).sum().item()

    #test_loss /= len(test_loader.dataset)
    writer.add_scalar('data/test_ae_loss', test_ae_loss, epoch)
    print('Test set: Average ae loss: {:.8f})'.format(test_ae_loss))
    #print('\nTest set: Average loss: {:.4f}, Accuracy: {}/{} ({:.0f}%)\n'.format(
    #	test_loss, correct, len(test_loader.dataset),
    #	100. * correct / len(test_loader.dataset)))

def output_rep(encoder, decoder, device, train_loader, test_loader):
    encoder.eval()
    decoder.eval()
    fout = open('rep_ae_test.csv', 'w')
    writer = csv.writer(fout, delimiter=',')
    with torch.no_grad():
        for data, target in test_loader:
            data, target = data.to(device), target.to(device)
            output = encoder(data)
            output = output.view((output.size(0),-1))
            output = list(output.data.cpu().numpy())
            labels = list(target.data.cpu().numpy())
            for i in range(len(labels)):
                line = list(output[i]) + [labels[i]]
                writer.writerow(line)
    fout.close()

    fout = open('rep_ae_train.csv', 'w')
    writer = csv.writer(fout, delimiter=',')
    with torch.no_grad():
        for data, target in train_loader:
            data, target = data.to(device), target.to(device)
            output = encoder(data)
            output = output.view((output.size(0),-1))
            output = list(output.data.cpu().numpy())
            labels = list(target.data.cpu().numpy())
            for i in range(len(labels)):
                line = list(output[i]) + [labels[i]]
                writer.writerow(line)
    fout.close()

def main():
    # Training settings
    parser = argparse.ArgumentParser(description='PyTorch MNIST Example')
    parser.add_argument('--batch-size', type=int, default=64, metavar='N',
            help='input batch size for training (default: 64)')
    parser.add_argument('--test-batch-size', type=int, default=1000, metavar='N',
            help='input batch size for testing (default: 1000)')
    parser.add_argument('--epochs', type=int, default=100, metavar='N',
            help='number of epochs to train (default: 100)')
    parser.add_argument('--lr', type=float, default=0.0001, metavar='LR',
            help='learning rate (default: 0.0001)')
    parser.add_argument('--gamma', type=float, default=0.7, metavar='M',
            help='Learning rate step gamma (default: 0.7)')
    parser.add_argument('--no-cuda', action='store_true', default=False,
            help='disables CUDA training')
    parser.add_argument('--dry-run', action='store_true', default=False,
            help='quickly check a single pass')
    parser.add_argument('--seed', type=int, default=1, metavar='S',
            help='random seed (default: 1)')
    parser.add_argument('--log-interval', type=int, default=10, metavar='N',
            help='how many batches to wait before logging training status')
    parser.add_argument('--save-model', action='store_true', default=False,
            help='For Saving the current Model')
    parser.add_argument('--width', type=int, default=32,
            help='width (default: 32)')
    parser.add_argument('--latent_width', type=int, default=4,
            help='latent_width (default: 4)')
    parser.add_argument('--depth', type=int, default=16,
            help='depth (default: 16)')
    parser.add_argument('--advdepth', type=int, default=16,
            help='advdepth (default: 16)')
    parser.add_argument('--advweight', type=float, default=0.5,
            help='advweight (default: 0.5)')
    parser.add_argument('--reg', type=float, default=0.2,
            help='reg (default: 0.2)')
    parser.add_argument('--latent', type=int, default=2,
            help='latent (default: 2)')
    parser.add_argument('--colors', type=int, default=1,
            help='colors (default: 1)')

    args = parser.parse_args()
    use_cuda = not args.no_cuda and torch.cuda.is_available()

    torch.manual_seed(args.seed)

    device = torch.device("cuda" if use_cuda else "cpu")

    train_kwargs = {'batch_size': args.batch_size}
    test_kwargs = {'batch_size': args.test_batch_size}
    if use_cuda:
        cuda_kwargs = {'num_workers': 1,'pin_memory': True,'shuffle': True}
        train_kwargs.update(cuda_kwargs)
        test_kwargs.update(cuda_kwargs)

    transform=transforms.Compose([
            transforms.Resize((24,24)),
            transforms.ToTensor(),
            #transforms.Normalize((0.1307,), (0.3081,))
            ])
    dataset1 = datasets.MNIST('./data', train=True, download=True,
            transform=transform)
    dataset2 = datasets.MNIST('./data', train=False,
            transform=transform)
    train_loader = torch.utils.data.DataLoader(dataset1,**train_kwargs)
    test_loader = torch.utils.data.DataLoader(dataset2, **test_kwargs)

    scales = int(round(math.log(args.width // args.latent_width, 2)))
    encoder = Encoder(scales, args.depth, args.latent, args.colors).to(device)
    decoder = Decoder(scales, args.depth, args.latent, args.colors).to(device)

    optimizer = optim.Adam(list(encoder.parameters()) + list(decoder.parameters()), lr=args.lr, weight_decay=1e-5)

    for epoch in range(1, args.epochs + 1):
        print(f'epoch: {epoch}')
        train(args, encoder, decoder, device, train_loader, optimizer, epoch)
        #test(encoder, decoder, device, test_loader, epoch)

        if args.save_model:
            torch.save(encoder.state_dict(), "mnist_cnn.pt")

    output_rep(encoder, decoder, device, train_loader, test_loader)

if __name__ == '__main__':
        main()


text:
#!/usr/bin/python3

import io
import sys
import os
from pathlib import Path
import struct
import scipy.io
import math
import numpy as np
from numpy.linalg import inv

def write_calibration(output_directory):
  f = open("{}/calib.txt".format(output_directory), "w")

  f.write("Tr: 1.0 0.0 0.0 0.0 0.0 1.0 0.0 0.0 0.0 0.0 1.0 0.0")

  f.close()

def rotxyz(r, p, h):
  """ Compute rotation matrix about the XYZ-axes.

      R = rotxyz(rph) returns a 3x3 rotation matrix R where (r,p,h) is a 3-vector
      of Euler angles (roll, pitch, heading) measured in radians.
  """

  cr = math.cos(r); sr = math.sin(r)
  cp = math.cos(p); sp = math.sin(p)
  ch = math.cos(h); sh = math.sin(h)

  R = np.array([[ch*cp, (-sh*cr + ch*sp*sr), ( sh*sr + ch*sp*cr)], \
                  [sh*cp, ( ch*cr + sh*sp*sr), (-ch*sr + sh*sp*cr)], \
                  [-sp,         cp*sr,                cp*cr       ]])
  return R

def rotZ(alpha, N = 3):
  """ Get rotation matrix (of dim N x N) about z-axis with angle alpha in randians. """
  R = np.identity(N)
  R[0,0] = math.cos(alpha)
  R[0,1] = -math.sin(alpha)
  R[1,0] = math.sin(alpha)
  R[1,1] = math.cos(alpha)

  return R

def initialize_firstpose(filename):
  mat = scipy.io.loadmat(filename)
  xyzrph = mat["SCAN"]["X_wv"][0,0]
  # print(xyzrph)
  T = np.identity(4)
  T[0:3, 0:3] = rotxyz(*xyzrph[3:6])
  T[0:3, 3:4] = xyzrph[0:3]

  first_pose = np.linalg.inv(T)

  return first_pose

if __name__ == "__main__":
  if len(sys.argv) < 2:
    print("Convert Michigan Velodyne files to KITTI compatible binary format and folder structure.", file = sys.stderr)
    print("Usage: convert <source> <sequence-id>", file = sys.stderr)
    print("Where <source> corresponds to the root folder of the dataset.", file = sys.stderr)

    sys.exit(1)

  input_dir = Path(sys.argv[1]) / "SCANS"
  sequence_dir = Path("dataset/sequences/{}/".format(sys.argv[2]))
  poses_file = Path("dataset/poses/{}.txt".format(sys.argv[2]))

  skip_existing = False

  if sequence_dir.exists() or poses_file.exists():
    print("Warning: Squence directory already exists.")
    answer = input("Do you want to overwrite all files? [y/s/N] ")
    if answer != "y" and answer != "s":
      print("Stopping conversion.")
      sys.exit(0)
    if answer == "y": print("Overwritting files.")
    if answer == "s":
      print("Skipping existing files.")
      skip_existing = True

  else:
    print("-- Creating directories.")
    sequence_dir.mkdir(parents = True)
    (sequence_dir / "velodyne").mkdir(parents = True)
    poses_dir = Path("dataset/poses/")
    if not poses_dir.exists(): poses_dir.mkdir(parents = True)

  write_calibration(sequence_dir.as_posix())

  scan_files = sorted([f for f in input_dir.iterdir() if f.is_file()])
  print("Converting files: ", end="")
  sys.stdout.flush()
  progress = 0

  out_poses = None
  if skip_existing:
    out_poses = poses_file.open("a")
  else:
    out_poses = poses_file.open("w")

  first_pose = initialize_firstpose(scan_files[0].as_posix())

  # Convert in KITTI Velodyne frame. (x pointing forward)
  C = rotZ(math.radians(-90.0))

  for i, f in enumerate(scan_files):
      if 100 * i / len(scan_files) > progress:
        progress += 10
        print(".", end="")
        sys.stdout.flush()

      out_filename = "{}/velodyne/{:06d}.bin".format(sequence_dir.as_posix(), i)
      if Path(out_filename).exists() and skip_existing: continue

      try:

        mat = scipy.io.loadmat(f.as_posix())

        xyzrph = mat["SCAN"]["X_wv"][0,0]
        # print(xyzrph)
        T = np.identity(4)
        T[0:3, 0:3] = rotxyz(*xyzrph[3:6])
        T[0:3, 3:4] = xyzrph[0:3]

        T = np.dot(first_pose, T)

        pose_string = "{} {} {} {}".format(*T[0, :])
        pose_string += " {} {} {} {}".format(*T[1, :])
        pose_string += " {} {} {} {}\n".format(*T[2, :])
        # print(pose_string)
        out_poses.write(pose_string)

        out_file = open(out_filename, "wb")

        M, N = mat["SCAN"]["XYZ"][0,0].shape
        for k in range(N):
          values = np.dot(C, mat["SCAN"]["XYZ"][0,0][:,k])

          byte_values = struct.pack("ffff", values[0], values[1], values[2], 0.0)
          out_file.write(byte_values)

        out_file.close()

      except Exception as e:
        print("\n\nException while processing: {}:".format(f.as_posix()), file = sys.stderr)
        print("Skipping...")

      else:
        pass

  while progress < 100:
    progress += 10
    print(".", end="")
  print(" finished.")

  #print(mat)

reason: import numpy,from numpy,import scipy


text:
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import re

import numpy as np
from tensorflow.contrib import learn  # pylint: disable=g-bad-import-order

TOKENIZER_RE = re.compile(r"[A-Z]{2,}(?![a-z])|[A-Z][a-z]+(?=[A-Z])|[\'\w\-]+",
                          re.UNICODE)

def tokenizer_char(iterator):
    for value in iterator:
        yield list(value)

def tokenizer_word(iterator):
    for value in iterator:
        yield TOKENIZER_RE.findall(value)

class MyVocabularyProcessor(learn.preprocessing.VocabularyProcessor):
    def __init__(self,
                 max_document_length,
                 min_frequency=0,
                 vocabulary=None):
        tokenizer_fn = tokenizer_char
        sup = super(MyVocabularyProcessor, self)
        sup.__init__(max_document_length, min_frequency, vocabulary, tokenizer_fn)

    def transform(self, raw_documents):
        """Transform documents to word-id matrix.
        Convert words to ids with vocabulary fitted with fit or the one
        provided in the constructor.
        Args:
          raw_documents: An iterable which yield either str or unicode.
        Yields:
          x: iterable, [n_samples, max_document_length]. Word-id matrix.
        """
        for tokens in self._tokenizer(raw_documents):
            word_ids = np.zeros(self.max_document_length, np.int64)
            for idx, token in enumerate(tokens):
                if idx >= self.max_document_length:
                    break
                word_ids[idx] = self.vocabulary_.get(token)
            yield word_ids