text:
import numpy as np

from PyQt5 import QtWidgets, QtCore

from PyQt5.QtWidgets import QMainWindow, QWidget, QPushButton, QAction, \
    QLabel, QTabWidget, QVBoxLayout, QHBoxLayout, QFileDialog, QLineEdit, \
    QComboBox, QCheckBox, QGroupBox, QScrollArea

from matplotlib.backends.backend_qt5agg import NavigationToolbar2QT \
    as NavigationToolbar

from sculptor.specfit import SpecFit
from sculptor.specmodel import fitting_methods
from sculptor.masksmodels import mask_presets
from sculptor.specmodelwidget import SpecModelWidget
from sculptor.specfitcanvas import SpecFitCanvas
from sculptor.menu_dialogs import ResampleWindow, EmceeWindow, NormalizeWindow
from sculptor.menu_dialogs import ResampleWindow, EmceeWindow, NormalizeWindow

def update_float_variable_from_linedit(lineEdit, variable,
                                       expression='{:.2f}'):

    try:
        new_value = float(lineEdit.text())
        lineEdit.setText(expression.format(new_value))
        return new_value
    except:
        print('[INFO] Input value is not convertable to float.')
        return variable

class SpecFitGui(QMainWindow):
    """

    """

    def __init__(self, spectrum=None, redshift=0):

        QMainWindow.__init__(self)

        self.resize(1150, 824)
        self.setWindowTitle("Sculptor - Interactive Modelling of astronomic electromagnetic spectra")

        # Add the SpecFit class to the GUI
        self.specfit = SpecFit(spectrum=spectrum, redshift=redshift)

        # Initialize class variables
        self.x_pos_a = 0
        self.x_pos_b = 0
        self.y_pos_a = 0
        self.y_pos_b = 0

        # Random seed variable

        # Resampling variables
        self.nsamples = 100
        self.resample_seed = 1234
        self.save_result_plots = True
        self.resample_foldername = '.'

        # MCMC variables

        # Setup the main menu
        # Add exit action
        exitAction = QAction('&Exit', self)
        exitAction.setShortcut('Ctrl+Q')
        exitAction.setStatusTip('Exit application')
        exitAction.triggered.connect(self.close)
        # Add save action
        saveAction = QAction('&Save', self)
        saveAction.setShortcut('Ctrl+S')
        saveAction.setStatusTip('Save to SpecFit Folder')
        saveAction.triggered.connect(self.save)
        # Add load action
        loadAction = QAction('&Load', self)
        loadAction.setShortcut('Ctrl+L')
        loadAction.setStatusTip('Load from SpecFit Folder')
        loadAction.triggered.connect(self.load)

        # Spectrum actions
        # Add import spectrum action
        importIrafSpecAction = QAction('Import IRAF spectrum', self)
        # loadAction.setShortcut('Ctrl+L')
        importIrafSpecAction.setStatusTip('Import IRAF Spectrum. Overrides '
                                          'current spectrum')
        importIrafSpecAction.triggered.connect(lambda: self.import_spectrum(
            mode='IRAF'))

        importPypeitSpecAction = QAction('Import PypeIt spectrum', self)
        # loadAction.setShortcut('Ctrl+L')
        importPypeitSpecAction.setStatusTip('Import PypeIt Spectrum. Overrides '
                                          'current spectrum')
        importPypeitSpecAction.triggered.connect(lambda: self.import_spectrum(
            mode='PypeIt'))

        importSodSpecAction = QAction('Import SpecOneD spectrum', self)
        # loadAction.setShortcut('Ctrl+L')
        importSodSpecAction.setStatusTip('Import SpecOneD Spectrum. Overrides '
                                          'current spectrum')
        importSodSpecAction.triggered.connect(lambda: self.import_spectrum(
            mode='SpecOneD'))

        importSDSSSpecAction = QAction('Import SDSS spectrum', self)
        # loadAction.setShortcut('Ctrl+L')
        importSDSSSpecAction.setStatusTip('Import SDSS Spectrum. Overrides '
                                         'current spectrum')
        importSDSSSpecAction.triggered.connect(lambda: self.import_spectrum(
            mode='SDSS'))

        normalizeByErrorAction = QAction('Normalize by error', self)
        normalizeByErrorAction.setStatusTip('Normalize the flux density, '
                                            'flux density, and object model '
                                            'numerical values by the median '
                                            'value of the flux density error.'
                                            ' The physical flux properties will'
                                            ' be unchanged.')
        normalizeByErrorAction.triggered.connect(
            lambda: self.normalize_spectrum_by_error())

        normalizeByFactorAction = QAction('Open normalization menu', self)
        normalizeByFactorAction.setStatusTip('Normalize the flux density, '
                                             'flux density, and object model '
                                             'numerical values by/to a '
                                             'specified factor. The physical '
                                             'flux properties will be '
                                             'unchanged.')
        normalizeByFactorAction.triggered.connect(
            lambda: self.open_normalize_dialog())

        # SpecModel Actions
        addSpecModelAction = QAction('Add SpecModel', self)
        addSpecModelAction.setStatusTip('Add a SpecModel to the fit')
        addSpecModelAction.triggered.connect(self.add_specmodel)

        removeSpecModelAction = QAction('Remove current SpecModel', self)
        removeSpecModelAction.setStatusTip('Remove current SpecModel from fit')
        removeSpecModelAction.triggered.connect(self.remove_current_spec_model)

        resetSpecModelAction = QAction('Remove all SpecModels', self)
        resetSpecModelAction.setStatusTip('Remove all SpecModels from fit')
        resetSpecModelAction.triggered.connect(self.remove_all_models)

        # Fit Actions

        runResampleAction = QAction('Run resample and fit', self)
        runResampleAction.setStatusTip('Resample the spectrum, fit all models '
                                       'and save the posterior parameter '
                                       'distributions.')
        runResampleAction.triggered.connect(self.resample_and_fit)

        setResampleAction = QAction('Set resample and fit parameters', self)
        setResampleAction.setStatusTip('Set the resample and fit parameters')
        setResampleAction.triggered.connect(self.open_resample_dialog)

        setEmceeAction = QAction('Set MCMC parameters', self)
        setEmceeAction.setStatusTip('Set the MCMC fit parameters')
        setEmceeAction.triggered.connect(self.open_emcee_dialog)

        self.statusBar()

        mainMenu = self.menuBar()
        mainMenu.setNativeMenuBar(False)
        fileMenu = mainMenu.addMenu('&File')
        # fileMenu.addSeparator()
        fileMenu.addAction(loadAction)
        fileMenu.addAction(saveAction)
        fileMenu.addAction(exitAction)

        spectrumMenu = mainMenu.addMenu('&Spectrum')
        spectrumMenu.addAction(importIrafSpecAction)
        spectrumMenu.addAction(importPypeitSpecAction)
        spectrumMenu.addAction(importSodSpecAction)
        spectrumMenu.addAction(importSDSSSpecAction)
        spectrumMenu.addSeparator()
        spectrumMenu.addAction(normalizeByErrorAction)
        spectrumMenu.addAction(normalizeByFactorAction)

        specModelMenu = mainMenu.addMenu('&SpecModel')
        specModelMenu.addAction(addSpecModelAction)
        specModelMenu.addAction(removeSpecModelAction)
        specModelMenu.addAction(resetSpecModelAction)

        fitMenu = mainMenu.addMenu('&Fit')
        fitMenu.addAction(setEmceeAction)
        fitMenu.addAction(setResampleAction)
        fitMenu.addAction(runResampleAction)

        # Setup the central widget and and the main tab widget
        self.centralwidget = QWidget()
        self.mainVLayout = QVBoxLayout(self.centralwidget)
        self.mainHLayout = QHBoxLayout()

        self.setCentralWidget(self.centralwidget)
        self.tabWidget = QTabWidget(self.centralwidget)
        self.tabWidget.setTabShape(QTabWidget.Rounded)
        self.tabWidget.setUsesScrollButtons(False)
        self.tabWidget.setDocumentMode(False)
        self.tabWidget.currentChanged.connect(self.tabchanged)

        self.mainHLayout.addWidget(self.tabWidget)
        self.mainHLayout.setStretch(1, 2)

        self.mainVLayout.addLayout(self.mainHLayout)

        self.initialize_main_tab()

        # Initialize the Key Press Event for the SpecFitCanvas
        self.gcid = self.specFitCanvas.mpl_connect('key_press_event',
                                                   self.on_press)

        self.show()

    def initialize_main_tab(self):

        # Initialize the fit tab
        self.fittab = QWidget()
        self.hLayoutFit = QHBoxLayout(self.fittab)

        self.vLayoutFitProperties = QVBoxLayout()

        # Add positional input
        self.boxRegionSelect = QGroupBox('Region Select')
        self.labelDispersionPos = QLabel('Dispersion region (Shift + a/d)')
        self.leXposA = QLineEdit('{:.2f}'.format(self.specfit.xlim[0]))
        self.leXposB = QLineEdit('{:.2f}'.format(self.specfit.xlim[1]))
        self.hLayoutXpos = QHBoxLayout()
        self.hLayoutXpos.addWidget(self.leXposA)
        self.hLayoutXpos.addWidget(self.leXposB)
        self.labelFluxPos = QLabel('Flux region (Shift + w/s)')
        self.leYposA = QLineEdit('{:.2e}'.format(self.specfit.ylim[0]))
        self.leYposB = QLineEdit('{:.2e}'.format(self.specfit.ylim[1]))
        self.hLayoutYpos = QHBoxLayout()
        self.hLayoutYpos.addWidget(self.leYposA)
        self.hLayoutYpos.addWidget(self.leYposB)

        pos_le = [self.leXposA, self.leXposB, self.leYposA, self.leYposB]
        for le in pos_le:
            le.returnPressed.connect(self.update_region_from_ui)

        self.buttonSetX = QPushButton('Set dispersion range (x)')
        self.buttonSetX.clicked.connect(self.set_plot_dispersion)
        self.buttonSetY = QPushButton('Set flux range (y)')
        self.buttonSetY.clicked.connect(self.set_plot_fluxden)
        self.buttonResetPlot = QPushButton('Reset plot (r)')
        self.buttonResetPlot.clicked.connect(self.reset_plot_region)

        self.hLayoutRegionButtons = QHBoxLayout()
        self.hLayoutRegionButtons.addWidget(self.buttonSetX)
        self.hLayoutRegionButtons.addWidget(self.buttonSetY)
        self.hLayoutRegionButtons.addWidget(self.buttonResetPlot)

        self.vLayoutBoxRegionSelect = QVBoxLayout(self.boxRegionSelect)
        self.vLayoutBoxRegionSelect.addWidget(self.labelDispersionPos)
        self.vLayoutBoxRegionSelect.addLayout(self.hLayoutXpos)
        self.vLayoutBoxRegionSelect.addWidget(self.labelFluxPos)
        self.vLayoutBoxRegionSelect.addLayout(self.hLayoutYpos)
        self.vLayoutBoxRegionSelect.addLayout(self.hLayoutRegionButtons)

        # Masking
        self.boxMaskSelect = QGroupBox('Masking')
        self.vLayoutBoxMaskSelect = QVBoxLayout(self.boxMaskSelect)

        self.boxMaskPreset = QComboBox()
        for mask_preset_key in mask_presets.keys():
            self.boxMaskPreset.addItem(mask_preset_key)

        self.buttonLoadPreset = QPushButton('Load mask preset')
        self.buttonLoadPreset.clicked.connect(lambda:
                                              self.load_mask_preset())

        self.hLayoutMaskPreset = QHBoxLayout()
        self.hLayoutMaskPreset.addWidget(self.boxMaskPreset)
        self.hLayoutMaskPreset.addWidget(self.buttonLoadPreset)

        self.buttonMask = QPushButton('Mask (m)')
        self.buttonMask.clicked.connect(lambda: self.update_mask(mode='mask'))
        self.buttonUnmask = QPushButton('Unmask (u)')
        self.buttonUnmask.clicked.connect(lambda: self.update_mask(mode='unmask'))
        self.hLayoutMaskAction = QHBoxLayout()
        self.buttonResetMask = QPushButton('Reset mask (Shift + r)')
        self.buttonResetMask.clicked.connect(lambda: self.reset_mask())

        self.hLayoutMaskAction.addWidget(self.buttonMask)
        self.hLayoutMaskAction.addWidget(self.buttonUnmask)
        self.hLayoutMaskAction.addWidget(self.buttonResetMask)

        self.vLayoutBoxMaskSelect.addLayout(self.hLayoutMaskAction)
        self.vLayoutBoxMaskSelect.addLayout(self.hLayoutMaskPreset)

        # Super parameters
        # Add global parameter
        self.boxSuperParam = QGroupBox('Super parameters')
        self.vLayoutSuperParam = QVBoxLayout(self.boxSuperParam)
        self.hLayoutAddSuperParam = QHBoxLayout()
        self.hLayoutDelSuperParam = QHBoxLayout()
        self.leSuperParamName = QLineEdit('super_param')
        self.leSuperParamName.setMaxLength(20)
        self.buttonAddSuperParam = QPushButton('Add super parameter')
        self.buttonAddSuperParam.clicked.connect(self.add_super_param)

        self.boxSelectSuperParam = QComboBox()
        self.buttonDelSuperParam = QPushButton('Remove super parameter')
        self.buttonDelSuperParam.clicked.connect(self.remove_super_param)

        self.hLayoutAddSuperParam.addWidget(self.leSuperParamName)
        self.hLayoutAddSuperParam.addWidget(self.buttonAddSuperParam)
        self.hLayoutDelSuperParam.addWidget(self.boxSelectSuperParam)
        self.hLayoutDelSuperParam.addWidget(self.buttonDelSuperParam)

        self.vLayoutSuperParam.addLayout(self.hLayoutAddSuperParam)
        self.vLayoutSuperParam.addLayout(self.hLayoutDelSuperParam)

        # Redshift parameter
        self.labelRedshift = QLabel('Redshift')
        self.leRedshift = QLineEdit()
        if self.specfit.redshift is not None:
            self.leRedshift.setText('{}'.format(self.specfit.redshift))
        else:
            self.leRedshift.setText('None')
        self.leRedshift.setMaxLength(20)
        self.leRedshift.returnPressed.connect(self.update_specfit_from_ui)

        self.hLayoutRedshift = QHBoxLayout()
        self.hLayoutRedshift.addWidget(self.labelRedshift)
        self.hLayoutRedshift.addWidget(self.leRedshift)

        # Fitting
        self.buttonFit = QPushButton('Fit all')
        self.buttonFit.clicked.connect(self.fit)
        self.buttonFitSaveResults = QPushButton('Fit all + Save results')
        self.buttonFitSaveResults.clicked.connect(lambda: self.fit(
            save_results=True))
        self.boxFittingMethod = QComboBox()
        for key in fitting_methods.keys():
            self.boxFittingMethod.addItem(key)
        self.boxFittingMethod.currentTextChanged.connect(
            self.update_fitting_method)

        # Build vLayoutFitProperties
        self.vLayoutFitProperties.addWidget(self.boxRegionSelect)
        self.vLayoutFitProperties.addWidget(self.boxMaskSelect)
        self.vLayoutFitProperties.addWidget(self.boxSuperParam)
        self.vLayoutFitProperties.addLayout(self.hLayoutRedshift)
        self.vLayoutFitProperties.addWidget(self.boxFittingMethod)
        self.vLayoutFitProperties.addWidget(self.buttonFit)
        self.vLayoutFitProperties.addWidget(self.buttonFitSaveResults)

        # Add SpecFitCanvas
        self.specFitCanvas = SpecFitCanvas(self.specfit)
        self.toolbar = NavigationToolbar(self.specFitCanvas, self)
        self.vLayoutCanvas = QVBoxLayout()
        self.vLayoutCanvas.addWidget(self.specFitCanvas)
        self.vLayoutCanvas.addWidget(self.toolbar)

        # Set the ClickFocus active on the SpecFitCanvas
        self.specFitCanvas.setFocusPolicy(QtCore.Qt.ClickFocus)
        self.specFitCanvas.setFocus()

        # Super parameter scroll area widget
        self.build_super_param_widget()

        # Build hLayoutFit
        self.hLayoutFit.addLayout(self.vLayoutFitProperties)
        self.hLayoutFit.addLayout(self.vLayoutCanvas)
        self.hLayoutFit.setStretch(1, 2)

        # Add super parameter tab widget
        self.vLayoutCanvas.addWidget(self.ScArea)
        self.vLayoutCanvas.setStretch(0, 2)

        self.tabWidget.addTab(self.fittab, "Fit")

    def build_super_param_widget(self):

        superParamWidget = QWidget()
        hLayoutModel = QHBoxLayout(superParamWidget)

        superParams = self.specfit.super_params

        lineditlist = []
        varyboxlist = []
        widgetlist = []

        for jdx, param in enumerate(superParams):

            widgetlist = []

            groupBoxParam = QGroupBox(param)
            vLayoutGroupBoxParam = QVBoxLayout(groupBoxParam)

            label = QLabel(param)
            linedit = QLineEdit('{:.4E}'.format(superParams[param].value))
            linedit.setMaxLength(20)
            expr_linedit = QLineEdit('{}'.format(superParams[param].expr))
            expr_linedit.setMaxLength(20)
            min_label = QLabel("min")
            min_linedit = QLineEdit('{:.4E}'.format(superParams[param].min))
            min_linedit.setMaxLength(20)
            max_label = QLabel("max")
            max_linedit = QLineEdit('{:.4E}'.format(superParams[param].max))
            max_linedit.setMaxLength(20)
            vary_checkbox = QCheckBox("vary")
            vary_checkbox.setChecked(superParams[param].vary)

            widgetlist.extend(
                [label, linedit, expr_linedit, min_label, min_linedit,
                 max_label, max_linedit, vary_checkbox])
            lineditlist.extend(
                [linedit, expr_linedit, min_linedit, max_linedit])
            varyboxlist.append(vary_checkbox)

            for w in widgetlist:
                vLayoutGroupBoxParam.addWidget(w)

            hLayoutModel.addWidget(groupBoxParam)

            vary_checkbox.stateChanged.connect(self.update_super_params_from_ui)

        if widgetlist:
            self.super_params_widgetlist = widgetlist
            self.super_params_lineditlist = lineditlist
            self.super_params_varybox_list = varyboxlist

            # Activate input for lineEdit lists
            for l in lineditlist:
                l.returnPressed.connect(self.update_super_params_from_ui)

        # Add model parameters in scroll area
        if hasattr(self, 'ScArea'):
            self.vLayoutCanvas.removeWidget(self.ScArea)
            self.deleteItemsOfLayout(self.ScArea.layout())
        else:
            self.ScArea = QScrollArea()
            self.ScArea.setLayout(QHBoxLayout())
            self.ScArea.setWidgetResizable(True)
            self.ScArea.setHorizontalScrollBarPolicy(
                QtCore.Qt.ScrollBarAlwaysOn)

        self.ScArea.setWidget(superParamWidget)

        self.vLayoutCanvas.addWidget(self.ScArea)

    def on_press(self, event):

        if event.key == 'm':

            self.update_mask(mode='mask')

        elif event.key == 'u':

            self.update_mask(mode='unmask')

        elif event.key == 'f':

            self.fit()

        elif event.key == 'A':

            self.x_pos_a = event.xdata
            self.leXposA.setText('{:.2f}'.format(self.x_pos_a))

        elif event.key == 'D':

            self.x_pos_b = event.xdata
            self.leXposB.setText('{:.2f}'.format(self.x_pos_b))

        elif event.key == 'W':

            self.y_pos_b = event.ydata
            self.leYposB.setText('{:.2e}'.format(self.y_pos_b))

        elif event.key == 'S':

            self.y_pos_a = event.ydata
            self.leYposA.setText('{:.2e}'.format(self.y_pos_a))

        # Zoom into X-axis region
        elif event.key == 'x':

            self.set_plot_dispersion()

        elif event.key == 'y':

            self.set_plot_fluxden()

        # Reset the region values
        elif event.key == 'r':

            self.reset_plot_region()

        # Full reset (region + mask)
        elif event.key == 'R':
            self.reset_mask()

    # --------------------------------------------------------------------------
    # Region Actions
    # --------------------------------------------------------------------------

    def set_plot_dispersion(self):

        self.update_region_from_ui()
        self.specfit.xlim = [self.x_pos_a, self.x_pos_b]
        self.update_specfit_plot()

    def set_plot_fluxden(self):

        self.update_region_from_ui()
        self.specfit.ylim = [self.y_pos_a, self.y_pos_b]
        self.update_specfit_plot()

    def reset_plot_region(self):

        self.reset_region()
        self.update_specfit_plot()

    def reset_region(self):

        if self.specfit.spec is not None:
            self.x_pos_a = min(self.specfit.spec.dispersion)
            self.x_pos_b = max(self.specfit.spec.dispersion)
            self.y_pos_a = min(self.specfit.spec.fluxden)
            self.y_pos_b = max(self.specfit.spec.fluxden)

        else:
            self.x_pos_a = 0
            self.x_pos_b = 1
            self.y_pos_a = 0
            self.y_pos_b = 1

        # Excluded, to get back to previous region
        # self.leXposA.setText('{:.2f}'.format(self.x_pos_a))
        # self.leXposB.setText('{:.2f}'.format(self.x_pos_b))
        # self.leYposA.setText('{:.2E}'.format(self.y_pos_a))
        # self.leYposB.setText('{:.2E}'.format(self.y_pos_b))

        self.specfit.xlim = [self.x_pos_a, self.x_pos_b]
        self.specfit.ylim = [self.y_pos_a, self.y_pos_b]

    def update_region_from_ui(self):

        self.x_pos_a = update_float_variable_from_linedit(self.leXposA,
                                                          self.x_pos_a)
        self.x_pos_b = update_float_variable_from_linedit(self.leXposB,
                                                          self.x_pos_b)
        self.y_pos_a = update_float_variable_from_linedit(self.leYposA,
                                                          self.y_pos_a,
                                                          expression='{:.2E}')
        self.y_pos_b = update_float_variable_from_linedit(self.leYposB,
                                                          self.y_pos_b,
                                                          expression='{:.2E}')

        self.specFitCanvas.setFocus()

    # --------------------------------------------------------------------------
    # Mask Actions
    # --------------------------------------------------------------------------

    def update_mask(self, mode='mask'):

        self.x_pos_a = float(self.leXposA.text())
        self.x_pos_b = float(self.leXposB.text())

        self.mask_region(mode)

        self.update_specfit_plot()

    def mask_region(self, mode='mask'):

        if self.specfit.spec is not None:

            mask_between = np.sort(np.array([self.x_pos_a,
                                             self.x_pos_b]))

            spec = self.specfit.spec
            lo_index = np.argmin(np.abs(spec.dispersion - mask_between[0]))
            up_index = np.argmin(np.abs(spec.dispersion - mask_between[1]))
            if mode == 'mask':
                self.specfit.spec.mask[lo_index:up_index] = False
            elif mode == 'unmask':
                self.specfit.spec.mask[lo_index:up_index] = True

        self.specfit.update_specmodel_spectra()

    def reset_mask(self):

        self.reset_region()
        self.specfit.spec.reset_mask()
        self.update_specfit_plot()

        self.specfit.update_specmodel_spectra()

    def load_mask_preset(self):

        mask_preset_key = self.boxMaskPreset.currentText()

        mask_preset = mask_presets[mask_preset_key]

        if mask_preset['rest_frame']:
            one_p_z = 1 + self.specfit.redshift
        else:
            one_p_z = 1

        for mask_range in mask_preset['mask_ranges']:
            wave_a = mask_range[0] * one_p_z
            wave_b = mask_range[1] * one_p_z

            if self.specfit.spec is not None:
                self.specfit.spec.mask_between([wave_a, wave_b])

        self.update_specfit_plot()

    def update_specfit_from_ui(self):

        # update redshift
        self.specfit.redshift = float(self.leRedshift.text())
        self.specfit.update_specmodels()
        self.leRedshift.setText('{:.4f}'.format(self.specfit.redshift))

        self.specFitCanvas.plot(self.specfit)
        self.specFitCanvas.setFocus()

    # --------------------------------------------------------------------------
    # Super Param Actions
    # --------------------------------------------------------------------------

    def update_boxSelectSuperParam(self):
        """ Only used when loading a model to update the Combobox"""

        for param in self.specfit.super_params:
            if self.boxSelectSuperParam.findText(param) < 0:
                self.boxSelectSuperParam.addItem(param)

    def rebuild_super_params_widget(self):

        layout = self.ScArea.layout()
        self.deleteItemsOfLayout(layout)

        self.build_super_param_widget()

    def add_super_param(self):

        param_name = self.leSuperParamName.text()
        # Check if item already in QComboBox
        if self.boxSelectSuperParam.findText(param_name) < 0:
            # Add item to QComboBox
            self.boxSelectSuperParam.addItem(param_name)
            # Add param to SpecModel
            self.specfit.add_super_param(param_name)

        self.rebuild_super_params_widget()
        self.update_specmodelwidgets_for_global_params()

    def remove_super_param(self):

        param_name = self.boxSelectSuperParam.currentText()
        # Remove global param from SpecModel
        self.specfit.remove_super_param(param_name)
        # Remove item from QComboBox
        idx = self.boxSelectSuperParam.findText(param_name)
        self.boxSelectSuperParam.removeItem(idx)

        self.rebuild_super_params_widget()
        self.update_specmodelwidgets_for_global_params()

    def update_specmodelwidgets_for_global_params(self):

        for idx in range(self.tabWidget.count() - 1):
            self.tabWidget.setCurrentIndex(idx + 1)
            specmodel_widget = self.tabWidget.currentWidget()
            specmodel_widget.rebuild_global_params_tab()
            specmodel_widget.rebuild_model_tabs()

        self.tabWidget.setCurrentIndex(0)
        self.specFitCanvas.setFocus()

    # --------------------------------------------------------------------------
    # Spectrum Actions
    # --------------------------------------------------------------------------

    def normalize_spectrum_by_error(self):

        self.specfit.normalize_spectrum_by_error()
        current_widget = self.tabWidget.currentWidget()
        if self.tabWidget.currentIndex() > 0:
            current_widget.reset_plot_region()
        else:
            self.reset_plot_region()

    def normalize_spectrum_to_factor(self, factor):

        self.specfit.normalize_spectrum_to_factor(factor)
        current_widget = self.tabWidget.currentWidget()
        if self.tabWidget.currentIndex() > 0:
            current_widget.reset_plot_region()
        else:
            self.reset_plot_region()

    def normalize_spectrum_by_factor(self, factor):

        self.specfit.normalize_spectrum_by_factor(factor)
        current_widget = self.tabWidget.currentWidget()
        if self.tabWidget.currentIndex() > 0:
            current_widget.reset_plot_region()
        else:
            self.reset_plot_region()

    # --------------------------------------------------------------------------
    # Fit Actions
    # --------------------------------------------------------------------------

    def update_fitting_method(self, value):

        self.specfit.fitting_method = value

    def fit(self, save_results=False):

        if len(self.specfit.specmodels) > 0 and self.specfit.spec is not None:

            if save_results:
                # Select save folder
                foldername = str(
                    QFileDialog.getExistingDirectory(self, "Select Directory"))
                if foldername:
                    self.specfit.fit(save_results=save_results,
                                     foldername=foldername)
                    self.specFitCanvas.fig.savefig(foldername+'/SpecFit.png')
                else:
                    self.statusBar().showMessage('No SAVE folder selected.')
            else:
                self.specfit.fit()

            self.rebuild_super_params_widget()
            self.update_specfit_plot()

    def resample_and_fit(self):

        if len(self.specfit.specmodels) == 0:
            self.statusBar().showMessage('[ERROR] No SpecModel found. Cannot '
                                         'resample and fit spectrum.')

        else:

            self.specfit.resample(n_samples=self.nsamples,
                                  foldername=self.resample_foldername,
                                  save_result_plots=self.save_result_plots,
                                  seed=self.resample_seed)

    def open_resample_dialog(self):

        self.resampleDialog = ResampleWindow(self)
        self.resampleDialog.show()

    def open_emcee_dialog(self):

        self.emceeDialog = EmceeWindow(self)
        self.emceeDialog.show()

    def open_normalize_dialog(self):
        self.normalizeDialog = NormalizeWindow(self)
        self.normalizeDialog.show()

    # --------------------------------------------------------------------------
    # File Actions
    # --------------------------------------------------------------------------

    def closeEvent(self, event):
        """This function modifies the standard closing of the GUI

        Parameters
        ----------
        event : event
            The closing event.
        """

        result = QtWidgets.QMessageBox.question(self,
                      "Exit Dialog",
                      "Are you sure you want to exit ?",
                      QtWidgets.QMessageBox.Yes| QtWidgets.QMessageBox.No)
        event.ignore()

        if result == QtWidgets.QMessageBox.Yes:
            event.accept()

    def save(self):

        # Select save folder
        foldername = str(
            QFileDialog.getExistingDirectory(self, "Select Directory"))

        if foldername:
            self.specfit.save(foldername)
        else:
            self.statusBar().showMessage('[WARNING] No save directory '
                                         'selected.')

    def load(self):

        # Select load folder
        foldername = str(
            QFileDialog.getExistingDirectory(self, "Select Directory"))

        if foldername:
            # Remove all SpecModels
            self.remove_all_models()

            self.specfit.load(foldername)

            # Update dispersion and flux density limits from specfit
            self.leXposA.setText('{:.2f}'.format(self.specfit.xlim[0]))
            self.leXposB.setText('{:.2f}'.format(self.specfit.xlim[1]))
            self.leYposA.setText('{:.2e}'.format(self.specfit.ylim[0]))
            self.leYposB.setText('{:.2e}'.format(self.specfit.ylim[1]))

            # Rebuild TabWidget with loaded SpecModels
            self.add_specmodels_in_specfit()

            self.tabWidget.setCurrentIndex(0)
            # self.reset_region()
            self.update_specfit_plot()
            self.update_boxSelectSuperParam()
            self.rebuild_super_params_widget()

            # Update redshift value from specfit
            self.leRedshift.setText('{}'.format(self.specfit.redshift))

            # Update fitting method from specfit
            index = self.boxFittingMethod.findText(self.specfit.fitting_method,
                                                   QtCore.Qt.MatchFixedString)
            if index >= 0:
                self.boxFittingMethod.setCurrentIndex(index)

        else:
            self.statusBar().showMessage('[WARNING] No directory to '
                                         'load a model from selected.')

    # --------------------------------------------------------------------------
    # Spectrum Actions
    # --------------------------------------------------------------------------

    def import_spectrum(self, mode='IRAF'):
        """ Import spectrum """

        # Select spectrum
        fileName, fileFilter = QFileDialog.getOpenFileName(self, "Import "
                                                               "spectrum")

        if fileName:
            self.specfit.import_spectrum(fileName, filetype=mode)

            # Re-plot main canvas -> later update fittab
            self.reset_region()
            self.specFitCanvas.plot(self.specfit)

            # Re-plot all SpecModel tabs
            for idx in range(self.tabWidget.count()-1):

                self.tabWidget.setCurrentIndex(idx+1)
                specmodel_widget = self.tabWidget.currentWidget()
                specmodel_widget.reset_plot_region()
                specmodel_widget.update_specmodel_plot()

    # --------------------------------------------------------------------------
    # SpecModel Actions
    # --------------------------------------------------------------------------

    def add_specmodel(self):

        # Add SpecModel to SpecFit
        self.specfit.add_specmodel()

        specmodel_widget = SpecModelWidget(self, self.specfit.specmodels[-1])

        self.add_specmodel_tab(specmodel_widget)

    def add_specmodel_tab(self, specmodel_widget, specmodel_name='SpecModel'):

        index = self.tabWidget.addTab(specmodel_widget, specmodel_name)
        self.tabWidget.setCurrentIndex(index)

    def add_specmodels_in_specfit(self):

        for specmodel in self.specfit.specmodels:
            specmodel_widget = SpecModelWidget(self, specmodel)
            self.add_specmodel_tab(specmodel_widget,
                                   specmodel_name=specmodel.name)

            specmodel_widget.rebuild_model_tabs()
            specmodel_widget.update_boxSelectGlobalParam()

    def remove_all_models(self):

        for tabindex in reversed(range(self.tabWidget.count())[1:]):
            self.remove_spec_model(tabindex)

    def remove_current_spec_model(self):

        # Get current SpecModel index from TabWidget
        tabindex = self.tabWidget.currentIndex()
        self.remove_spec_model(tabindex)

    #     TODO UPDATE ALL TABS WITH REGARD TO THE PROPAGATED SPECTRUM

    def remove_spec_model(self, index):

        if index > 0:

            # Remove SpecModel from SpecFit
            self.specfit.delete_specmodel(index-1)
            # Update the TabWidget accordingly
            specmodel_widget = self.tabWidget.currentWidget()
            specmodel_widget.deleteLater()
            self.tabWidget.removeTab(index)

        else:
            self.statusBar().showMessage('Current tab does not contain a '
                                         'SpecModel')

    # --------------------------------------------------------------------------
    # Update Actions (SpecModelWidget, SpecModel, TabWidget)
    # --------------------------------------------------------------------------

    def update_specfit_plot(self):
        self.specFitCanvas.plot(self.specfit)

    def tabchanged(self):

        # Get current tab widget
        idx = self.tabWidget.currentIndex()

        # Update the tab
        if idx == 0:
            self.specFitCanvas.plot(self.specfit)
            self.update_specfit_plot()
            self.specFitCanvas.setFocus()

        elif idx > 0:
            self.specfit.update_specmodel_spectra()
            self.update_specmodel_tab(idx)

    def update_specmodel_tab(self, index):

        self.tabWidget.setCurrentIndex(index)

        # UPDATE SpecMODEL???

        specmodel_widget = self.tabWidget.currentWidget()
        # specmodel_widget.reset_plot_region()
        specmodel_widget.update_specmodel_plot()
        specmodel_widget.specModelCanvas.setFocus()

    def update_super_params_from_ui(self):

        for jdx, param in enumerate(self.specfit.super_params):

            new_value = float(
                self.super_params_lineditlist[jdx * 4 + 0].text())
            new_expr = self.super_params_lineditlist[jdx * 4 + 1].text()
            new_min = float(
                self.super_params_lineditlist[jdx * 4 + 2].text())
            new_max = float(
                self.super_params_lineditlist[jdx * 4 + 3].text())
            new_vary = self.super_params_varybox_list[jdx].isChecked()

            if new_expr == 'None':
                new_expr = None

            for specmodel in self.specfit.specmodels:

                # Set the new parameter values
                specmodel.global_params[param].set(value=new_value,
                                                       expr=new_expr,
                                                       min=new_min,
                                                       max=new_max,
                                                       vary=new_vary)

        self.update_specmodelwidgets_for_global_params()

    def deleteItemsOfLayout(self, layout):
         if layout is not None:
             while layout.count():
                 item = layout.takeAt(0)
                 widget = item.widget()
                 if widget is not None:
                     widget.setParent(None)
                 else:
                     self.deleteItemsOfLayout(item.layout())

max_stars_repo_name: jtschindler/sculptor

max_issues_repo_name: jtschindler/sculptor

max_forks_repo_name: jtschindler/sculptor


text:
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 19 13:28:23 2016

@author: bing
"""

import numpy as np
import matplotlib.pyplot as plt

import matplotlib as mpl

#plt.style.use('ggplot')

font = {'family' : 'Times New Roman',
        'weight' : 'normal',
        'size'   : '12'}

mpl.rc('font', **font)  # pass in the font dict as kwargs

#dat = np.genfromtxt('sobol/random.dat')
#
#fig, ((ax1,ax2), (ax3, ax4)) = plt.subplots(ncols=2,nrows=2)
#
#ax4.plot(dat[:,0],dat[:,1],'o',markersize=3)
#
#dat2 = np.genfromtxt('ran1/normal.dat')
#
#ax3.plot(dat2[:,0],dat2[:,1],'o',markersize=3)
#
#x,y = np.genfromtxt('random.dat',unpack=True)
#
#ax2.plot(x,y,'o',markersize=3)
#
#x,y = np.genfromtxt('ran2.dat',unpack=True)
#
#ax1.plot(x,y,'o',markersize=3)

x = np.array([10000, 20000, 100000 , 200000, 500000])
x = 1./np.sqrt(x)

y = np.array([1.0489922729470280, 1.0151243791486708, 1.0081307660222230,
     1.0036717974790175, 1.0014285885836964])
y = y-1.

z = np.array([1.0145262113813669,0.98850412318143310, 1.0001542763276223,0.99939065597754062,
                1.0000680719061754  ])
z = np.abs(z-1.)

fig, ax = plt.subplots()
ax.plot(x,y,'g-o',label='Normal')
ax.plot(x,z,'r-s',label='Uniform')
ax.set_xlabel('$1/\sqrt{N}$')
ax.set_ylabel('$\gamma -1 $')
ax.set_xlim(0,0.012)
ax.set_ylim(-0.002,0.05)

plt.legend(loc=0)

#ax.set_title('$\int e^{-x_1^2/2}dx$')
plt.show()
plt.grid(False)
plt.savefig('/home/bing/quasiMonteCarlo/integration.pdf')

ext: py

reason: import numpy


text:
import torch
import torch.nn as nn
import torch.nn.functional as F
from torch.nn.utils.rnn import pack_padded_sequence, pad_packed_sequence
from torch.nn import Conv1d, ReLU
from utils import common_util
import numpy as np

class ConnectionistTemporalClassification(nn.Module):
    def __init__(self, model_config):
        super(ConnectionistTemporalClassification, self).__init__()

        self.lstm = nn.LSTM(model_config.input_size,
                            model_config.lstm_dim,
                            num_layers=1, bidirectional=True, batch_first=True)

        self.model_config = model_config

        self.clf = nn.Linear(2 * model_config.lstm_dim, model_config.num_tags + 1)

        common_util.init_lstm_wt(self.lstm)
        common_util.init_linear_wt(self.hidden_layer)

    def forward(self, mfcc, length, **kwargs):
        lengths = length.view(-1).tolist()
        packed = pack_padded_sequence(mfcc, lengths, batch_first=True, enforce_sorted=False)
        output, hidden = self.lstm(packed)
        lstm_feats, _ = pad_packed_sequence(output, batch_first=True)  # h dim = B x t_k x n
        lstm_feats = lstm_feats.contiguous()
        logits = self.clf(lstm_feats)

        return logits

    def get_loss(self, logits, phone, length, **kwargs):
        b, t_k, d = list(logits.size())
        loss = 0
        for i in range(b):
            curr_len = length[i].cpu().data.numpy()
            curr_logit = logits[i].cpu().data.numpy()
            curr_phone = phone[i].cpu().data.numpy()
            loss += self.get_ctc_loss_single(curr_len, curr_logit, curr_phone)
        return loss

    '''
    For odd index return blank label
    otherwise map the index to ground truth phone
    index
    '''
    def get_phone_id(self, s, phone):
        if s % 2 == 1:
            return 0
        else:
            return phone[s // 2 - 1] + 1
    def has_same_label(self, s, phone):
        idx = s // 2
        return idx > 1 and phone[s // 2 - 1] == phone[s // 2 - 2]

    def get_ctc_loss_single(self, y, T, phone):
        num_phone = len(phone)
        S = 2*num_phone + 1
        alpha = np.zeros(T + 1, S + 1)
        t, s = 1, 1
        alpha[t, s] = y[t-1, self.get_phone_id(s, phone)]
        t, s = 1, 2
        alpha[t, s] = y[t - 1, self.get_phone_id(s, phone)]
        for t in range(2, T + 1):
            for s in range(1, S + 1):
                y_ = y[t - 1, self.get_phone_id(s, phone)]
                # blank or same labels
                if s % 2 == 1 or self.has_same_label(s, phone):
                    alpha[t, s] = (alpha[t-1, s] + alpha[t-1, s-1]) * y_
                else:
                    alpha[t, s] = (alpha[t-1, s] + alpha[t-1, s-1] + alpha[t-1, s-2]) * y_
        return alpha[T, S] + alpha[T, S - 1]

    def path_to_str(self, path):
        new_path = [path[0]] + [path[i] for i in range(1, len(path)) if path[i-1] != path[i]]
        return [p -1 for p in new_path if p > 0]

    def best_path_decode(self, logits, length):
        logits = logits.cpu().data.numpy()
        length = length.cpu().data.numpy()
        out_paths = []
        for y, curr_len in zip(logits, length):
            path = np.argmax(y[:curr_len, :], axis=1)
            out_paths.append(self.path_to_str(path))
        return out_paths

    def prefix_search_decode(self, logits, length):
        logprobs = logits.log_softmax(dim=2)
        logprobs = logprobs.cpu().data.numpy()
        length = length.cpu().data.numpy()

        out_paths = []
        for logprob, curr_len in zip(logprobs, length):
            beams = []


text:
# -*- coding: utf-8 -*-
"""chapter5.ipynb

Automatically generated by Colaboratory.

Original file is located at
    https://colab.research.google.com/drive/1FZ4JinLGstOWo72o-UtE1Q0A_V8k6unA
"""

# Chapter5 Biophysical ML
# we will explore in depth the problem of predicting how
# small drug-like molecules bind to a protein of interest in the human body.

# Our goal therefore is to design learning algorithms that can effectively predict when a
# given molecule is going to interact with a given protein. How can we do this?

#

# Commented out IPython magic to ensure Python compatibility.
##setup tensorflow v1
# %tensorflow_version 1.x

## this will install anaconda and deepchem, will add path, execution will take sometime

!wget -c https://repo.anaconda.com/archive/Anaconda3-2019.10-Linux-x86_64.sh
!chmod +x Anaconda3-2019.10-Linux-x86_64.sh
!bash ./Anaconda3-2019.10-Linux-x86_64.sh -b -f -p /usr/local
!conda install -y -c deepchem -c rdkit -c conda-forge -c omnia deepchem-gpu=2.3.0
import sys
sys.path.append('/usr/local/lib/python3.7/site-packages/')

############## check deepchem installation by
import deepchem as dc
dc.__version__ #should match with the installed

#need to reinstall deepchem bcz of an openmm error

!curl -Lo deepchem_installer.py https://raw.githubusercontent.com/deepchem/deepchem/master/scripts/colab_install.py

import deepchem_installer

# Commented out IPython magic to ensure Python compatibility.
# %time deepchem_installer.install(additional_packages=['mdtraj'], version='2.3.0')

import deepchem as dc

grid_featurizer = dc.feat.RdkitGridFeaturizer(voxel_width=2.0,feature_types=['hbond', 'salt_bridge', 'pi_stack',
'cation_pi', 'ecfp', 'splif'], sanitize=True, flatten=True)

# do not really need to use this, feature parameter
# in the next function will take care of it

#pdbbind dataset is large ~ 2GB, takes time

# bcz of an openmm error, I have to manually download
import deepchem as dc
from deepchem.utils import download_url

import os

data_dir = dc.utils.get_data_dir()
dataset_file = os.path.join(data_dir, "pdbbind_core_df.csv.gz")

if not os.path.exists(dataset_file):
    print('File does not exist. Downloading file...')
    download_url("https://s3-us-west-1.amazonaws.com/deepchem.io/datasets/pdbbind_core_df.csv.gz")
    print('File downloaded...')

raw_dataset = dc.utils.save.load_from_disk(dataset_file)

#raw_dataset

grid_featurizer = dc.feat.RdkitGridFeaturizer(
    voxel_width=16.0, feature_types=["ecfp", "splif", "hbond", "pi_stack", "cation_pi", "salt_bridge"],
    ecfp_power=5, splif_power=5, parallel=True, flatten=True, sanitize=True)

compound_featurizer = dc.feat.CircularFingerprint(size=128)

#load data
pdbbind_tasks, (train_dataset, valid_dataset, test_dataset), transformers = dc.molnet.load_pdbbind_grid(
    featurizer="ECFP", subset="refined")

train_dataset.get_data_shape(),valid_dataset.get_data_shape(),test_dataset.get_data_shape()

# now comes model
from sklearn.ensemble import RandomForestRegressor

#RSF model
sklearn_model = RandomForestRegressor(n_estimators=100)
model = dc.models.SklearnModel(sklearn_model)
model.fit(train_dataset)

#MLP model
n_features = train_dataset.X.shape[1]
model_mlp = dc.models.MultitaskRegressor(n_tasks=len(pdbbind_tasks),n_features=n_features,
layer_sizes=[2000, 1000],dropouts=0.5,learning_rate=0.0003)
model_mlp.fit(train_dataset, nb_epoch=250)

metric = dc.metrics.Metric(dc.metrics.pearson_r2_score)

#checking the RF model
print("Evaluating model")
train_scores = model.evaluate(train_dataset, [metric], transformers)
test_scores = model.evaluate(test_dataset, [metric], transformers)
print("Train scores ", train_scores)
print("Test scores ", test_scores)

#RF train r2: 0.89423
#RF test r2: 0.5099
#overfitting evident

#checking the MLP model
print("Evaluating model")
train_scores = model_mlp.evaluate(train_dataset, [metric], transformers)
test_scores = model_mlp.evaluate(test_dataset, [metric], transformers)
print("Train scores ", train_scores)
print("Test scores ", test_scores)

#MLP train r2: 0.918
#MLP test r2: 0.369
#overfitting more than RF

#This is different than what the book says, in book RF is more ovefit

test_dataset.y.shape

rf_predicted_test = model.predict(test_dataset)

rf_predicted_test.shape

# Commented out IPython magic to ensure Python compatibility.
#plots
# %matplotlib inline

import matplotlib
import numpy as np
import matplotlib.pyplot as plt

rf_predicted_test = model.predict(test_dataset)
mlp_predicted_test = model_mlp.predict(test_dataset)
rf_true_test = test_dataset.y

plt.figure(figsize=(10,5))
plt.subplot(1,2,1)
plt.scatter(rf_predicted_test, rf_true_test)
plt.xlabel('Predicted pIC50s')
plt.ylabel('True IC50')
plt.title(r'RF predicted IC50 vs. True IC50')
plt.xlim([-2.5, 2.5])
plt.ylim([-2.5, 2.5])
plt.plot([-2.5, 2.5], [-2.5, 2.5], color='k')
plt.text(-2.5, 2, 'RF R-squared = %0.2f' % 0.5099)
plt.subplot(1,2,2)
plt.scatter(mlp_predicted_test, rf_true_test)
plt.xlabel('Predicted pIC50s')
plt.ylabel('True IC50')
plt.title(r'MLP predicted IC50 vs. True IC50')
plt.xlim([-2.5, 2.5])
plt.ylim([-2.5, 2.5])
plt.plot([-2.5, 2.5], [-2.5, 2.5], color='k')
plt.text(-2.5, 2, 'MLP R-squared = %0.2f' % 0.369)
plt.show()

#chapter 5 completed

max_stars_repo_name: xinformatics/DeepLearningLifeSciences

max_issues_repo_name: xinformatics/DeepLearningLifeSciences

max_forks_repo_name: xinformatics/DeepLearningLifeSciences


text:
import numpy as np
from chordkit.hearing_models import cbw_volk as cbw, bark_zwicker as bark
from chordkit.pair_constants import SETHARES_CONSTANTS as sc, AUDITORY_CONSTANTS as ac, pair_volume, pair_distance
from chordkit.chord_utils import MergedSpectrum

# This file contains both the individual pairwise models used for assessing the
# roughness of partial pairs and the summing function that adds up all such
# pairwise contributions.

###################
# PAIRWISE MODELS #
###################

# Different pairwise roughness evaluation functions. Each has the same signature
# so that they can be called by the same line of code.

def helmholtz_roughness_pair(x_hz, *, x_p=1, ref_p=1, options={}):
    # Note that x_p and ref_p are now 1-indexed in the function call.

    # The following constants have not been given values by Helmholtz.
    # They are chosen here somewhat arbitrarily to make the graph resemble his.
    bPrime1 = 1
    bPrime2 = 1
    beta = 0.3 # Trial and error; Helmholtz doesn't specify

    # Take partial ref_p from the reference pitch
    ref_hz = options['ref'][ref_p - 1]

    delta = ((float(x_hz) / ref_hz) - 1) / 2
    theta = 15.0 / ref_hz

    s = 4 * bPrime1 * bPrime2 * (beta ** 2) / (beta ** 2 + (2 * np.pi * delta) ** 2)

    r = s * ((2 * theta * delta * x_p) ** 2) / ((theta ** 2 + (x_p * delta) ** 2) ** 2)

    return r

# Returns the value of the Sethares "sensory dissonance" (roughness)
# function at frequency x_hz with respect to ref_hz.
# Based on Sethares 1997, as implemented in Giordano 2015.
def sethares_roughness_pair(x_hz, ref_hz, v_x, v_ref, *, options = {
    'original': False,
    'amp_type': 'MIN',
    'cutoff': False,
}):
    s = sc['s_star'] / (sc['s1'] * min([x_hz, ref_hz]) + sc['s2'])

    if options['original'] == True:
        # This line not based on the Sethares 1993 paper, but on the implementation
        # on his website, https://sethares.engr.wisc.edu/comprog.html
        options['amp_type'] = 'MIN'

    v12 = pair_volume(v_x, v_ref, options['amp_type'])

    # Sethares' MATLAB implementation (but not the published paper)
    scaling = 5

    # The following scaling factor is introduced to ensure that the maximum
    # value obtained by the pairwise roughness function is approximately 1.
    # The numerator comes from Sethares' MATLAB implementation; the
    # denominator is my additional modification.
    if options['original'] == False:
        scaling = 5 / 0.8986

    # distance = pair_distance(x_hz, ref_hz)
    distance = np.abs(x_hz - ref_hz)

    # Cutoff: Following Hutchinson and Knopoff 1978, cuts off the Sethares
    # roughness function at 1.2 CBW, which prevents too-remote partials from
    # contributing to the final score. (Check whether this also prevents drift
    # for larger intervals.)
    # Sethares' original does not use this cutoff.
    if options['cutoff'] == True:
        cbw_limit = 1.2 * cbw(max[x_hz, ref_hz]) / 2
        if distance < ac['slow_beat_limit'] or distance >= cbw_limit:
            v12 = 0

    try:
        return v12 * scaling * (np.exp(-sc['a'] * s * distance) - np.exp(-sc['b'] * s * distance))
    except OverflowError:
        print(f"Overflow in computing roughness: a == {sc['a']}, b == {sc['b']}, s == {s}, distance == {distance}")

# Returns roughness contribution of two partials, based on an indicator
# function on the critical bandwidth, scaled to the amplitude of the partial.
def cbw_roughness_pair(x_hz, ref_hz, v_x, v_ref, options={ 'amp_type': 'MIN' }):
    cbw_limit = cbw(max([x_hz, ref_hz])) / 2
    distance = pair_distance(x_hz, ref_hz)

    if distance >= 15 and distance < cbw_limit:
        return pair_volume(v_x, v_ref, options['amp_type'])
    else:
        return 0

# From Bigand, Parncutt, and Lerdahl 1996. Note that the volumes are not used.
def parncutt_roughness_pair(x_hz, ref_hz, v_x, v_ref, options = {}):
    # Parameters asserted in BPL 1996 paper
    a = 0.25
    i_factor = 2

    x_bark = bark(x_hz)
    ref_bark = bark(ref_hz)
    distance = abs(x_bark - ref_bark)

    if distance < 1.2:
        return ((np.exp(1)/a) * distance * np.exp(-distance / a)) ** i_factor
    else:
        return 0

###################
# SUMMATION MODEL #
###################

# The following function is based on Sethares' 1993 model, which linearly sums
# all pairwise contributions to roughness. It can use any of the other
# pairwise roughness evaluation functions.
def roughness_complex(
    spectrum: MergedSpectrum,
    function_type: str = 'SETHARES',
    rough_limit: float = 0.1,
    *,
    options={
        'amp_type': 'MIN',
        'cutoff': False,
        'original': False,
        'show_partials': False
    }
):
    n = len(spectrum.partials['hz'])
    rough_partials = []
    rough_vals = np.zeros((n, n))

    if function_type.upper() == 'SETHARES':
        pair_assess = sethares_roughness_pair
    elif function_type.upper() == 'CBW':
        pair_assess = cbw_roughness_pair
    elif function_type.upper() == 'PARNCUTT':
        pair_assess = parncutt_roughness_pair
    elif function_type.upper() == 'HELMHOLTZ':
        pair_assess = helmholtz_roughness_pair
    else:
        raise ValueError(f'Invalid assessment function type: {function_type.upper()}')

    # Assess all pairs for roughness

    # Helmholtz's function works differently from the others.
    # In this case, `spectrum` is NOT a composite spectrum of the whole chord:
    # It contains exactly one pitch, which varies by transposition.
    # The second, reference pitch is stored under options['ref'] and is fixed.
    # This method will not work for cardinality > 2.
    if function_type.upper() == 'HELMHOLTZ':
        m = len(spectrum.partials['hz'])
        n = len(options['ref'])
        for i in range(m):
            for j in range(n):
                rough_vals[i][j] = helmholtz_roughness_pair(
                    spectrum.partials['hz'][i],
                    x_p=i + 1, # p is 1-indexed
                    ref_p=j + 1,
                    options=options
                )
                if options['show_partials'] == True and rough_vals[i][j] > rough_limit:
                    rough_partials.append((i, j))

    # Pairwise evaluations, as used by Hutchinson/Knopoff, Sethares, Parncutt)
    else:
        for i in range(n - 1):
            for j in range(i + 1, n):
                rough_vals[i][j] = pair_assess(
                    spectrum.partials['hz'][i],
                    spectrum.partials['hz'][j],
                    spectrum.partials['amp'][i],
                    spectrum.partials['amp'][j],
                    options=options
                )
                if options['show_partials'] == True and rough_vals[i][j] > rough_limit:
                    rough_partials.append((i, j))

    if options['show_partials']:
        return {
            'roughness': np.sum(rough_vals),
            'rough_partials': rough_partials
        }

    return np.sum(rough_vals)

max_stars_repo_name: bradley-gersh/chordkit

max_issues_repo_name: bradley-gersh/chordkit

max_forks_repo_name: bradley-gersh/chordkit


text:
# -*- coding: utf-8 -*-#
'''
# Name:         Softplus
# Description:
# Author:       super
# Date:         2020/5/24
'''

import numpy as np
import matplotlib.pyplot as plt

class CSoftplus(object):
    def forward(self, z):
        a = np.log(1 + np.exp(z))
        return a

    def backward(self, z, a, delta):
        p = np.exp(z)
        da = p / (1 + p)
        dz = np.multiply(delta, da)
        return da, dz


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

# Copyright (c) 2016-2019 by University of Kassel and Fraunhofer Institute for Energy Economics
# and Energy System Technology (IEE), Kassel. All rights reserved.

import numpy as np
import pandas as pd
from collections import UserDict

from pandapower.auxiliary import _select_is_elements_numba, _add_ppc_options, _add_auxiliary_elements
from pandapower.pd2ppc import _pd2ppc
from pandapower.pypower.idx_brch import branch_cols
from pandapower.pypower.idx_bus import bus_cols
import pandapower.pypower.idx_bus as idx_bus
from pandapower.pypower.makeYbus import makeYbus
from pandapower.pf.run_newton_raphson_pf import _run_dc_pf
from pandapower.run import rundcpp
from pandapower.build_branch import get_is_lines
from pandapower.create import create_buses, create_line_from_parameters

from pandapower.estimation.util import estimate_voltage_vector
from pandapower.estimation.idx_bus import *
from pandapower.estimation.idx_brch import *
from pandapower.estimation.results import _copy_power_flow_results

try:
    import pplog as logging
except ImportError:
    import logging
std_logger = logging.getLogger(__name__)
ZERO_INJECTION_STD_DEV = 0.001

def _initialize_voltage(net, init, calculate_voltage_angles):
    v_start, delta_start = None, None
    if init == 'results':
        v_start, delta_start = 'results', 'results'
    elif init == 'slack':
        res_bus = estimate_voltage_vector(net)
        v_start = res_bus.vm_pu.values
        if calculate_voltage_angles:
            delta_start = res_bus.va_degree.values
    elif init != 'flat':
        raise UserWarning("Unsupported init value. Using flat initialization.")
    return v_start, delta_start

def _init_ppc(net, v_start, delta_start, calculate_voltage_angles):
    # select elements in service and convert pandapower ppc to ppc
    net._options = {}
    _add_ppc_options(net, check_connectivity=False, init_vm_pu=v_start, init_va_degree=delta_start,
                     trafo_model="pi", mode="pf", enforce_q_lims=False,
                     calculate_voltage_angles=calculate_voltage_angles, switch_rx_ratio=2,
                     recycle=dict(_is_elements=False, ppc=False, Ybus=False))
    net["_is_elements"] = _select_is_elements_numba(net)
    _add_auxiliary_elements(net)
    ppc, ppci = _pd2ppc(net)

    # do dc power flow for phase shifting transformers
    if np.any(net.trafo.shift_degree):
        vm_backup = ppci["bus"][:, 7].copy()
        pq_backup = ppci["bus"][:, [2, 3]].copy()
        ppci["bus"][:, [2, 3]] = 0.
        ppci = _run_dc_pf(ppci)
        ppci["bus"][:, 7] = vm_backup
        ppci["bus"][:, [2, 3]] = pq_backup

    return ppc, ppci

def _add_measurements_to_ppci(net, ppci, zero_injection):
    """
    Add pandapower measurements to the ppci structure by adding new columns
    :param net: pandapower net
    :param ppci: generated ppci
    :return: ppc with added columns
    """
    meas = net.measurement.copy(deep=False)
    if meas.empty:
        raise Exception("No measurements are available in pandapower Network! Abort estimation!")

    # Convert side from string to bus id
    meas["side"] = meas.apply(lambda row:
                              net['line'].at[row["element"], row["side"]+"_bus"] if
                              row["side"] in ("from", "to") else
                              net[row["element_type"]].at[row["element"], row["side"]+'_bus'] if
                              row["side"] in ("hv", "mv", "lv") else row["side"], axis=1)

    # convert p, q, i measurement to p.u., u already in p.u.
    meas.loc[meas.measurement_type=="p", ["value", "std_dev"]] /= ppci["baseMVA"]
    meas.loc[meas.measurement_type=="q", ["value", "std_dev"]] /= ppci["baseMVA"]

    if not meas.query("measurement_type=='i'").empty:
        meas_i_mask = (meas.measurement_type=='i')
        base_i_ka = ppci["baseMVA"] / net.bus.loc[(meas.side.fillna(meas.element))[meas_i_mask].values,
                        "vn_kv"].values
        meas.loc[meas_i_mask, "value"] /= base_i_ka/np.sqrt(3)
        meas.loc[meas_i_mask, "std_dev"] /= base_i_ka/np.sqrt(3)

    if not meas.query("(measurement_type=='ia' )| (measurement_type=='va')").empty:
        meas_dg_mask = (meas.measurement_type=='ia')|(meas.measurement_type=='va')
        meas.loc[meas_dg_mask, "value"] = np.deg2rad(meas.loc[meas_dg_mask, "value"])
        meas.loc[meas_dg_mask, "std_dev"] = np.deg2rad(meas.loc[meas_dg_mask, "std_dev"])

    # Get elements mapping from pandapower to ppc
    map_bus = net["_pd2ppc_lookups"]["bus"]
    meas_bus = meas[(meas['element_type'] == 'bus')]
    if (map_bus[meas_bus['element'].values.astype(int)] >= ppci["bus"].shape[0]).any():
        std_logger.warning("Measurement defined in pp-grid does not exist in ppci! Will be deleted!")
        meas_bus = meas_bus[map_bus[meas_bus['element'].values.astype(int)] < ppci["bus"].shape[0]]

    # mapping to dict instead of np array ensures good performance for large indices
    # (e.g., 999999999 requires a large np array even if there are only 2 buses)
    # downside is loop comprehension to access the map
    map_line, map_trafo, map_trafo3w = None, None, None
    branch_mask = ppci['internal']['branch_is']
    if "line" in net["_pd2ppc_lookups"]["branch"]:
        map_line = {line_ix: br_ix for line_ix, br_ix in
                    zip(net.line.index, range(*net["_pd2ppc_lookups"]["branch"]["line"])) if branch_mask[br_ix]}

    if "trafo" in net["_pd2ppc_lookups"]["branch"]:
        trafo_ix_start, trafo_ix_end = net["_pd2ppc_lookups"]["branch"]["trafo"]
        trafo_ix_offset = np.sum(~branch_mask[:trafo_ix_start])
        trafo_ix_start, trafo_ix_end = trafo_ix_start - trafo_ix_offset, trafo_ix_end - trafo_ix_offset
        map_trafo = {trafo_ix: br_ix for trafo_ix, br_ix in
                     zip(net.trafo.index, range(trafo_ix_start, trafo_ix_end))
                     if branch_mask[br_ix+trafo_ix_offset]}

    if "trafo3w" in net["_pd2ppc_lookups"]["branch"]:
        trafo3w_ix_start, trafo3w_ix_end = net["_pd2ppc_lookups"]["branch"]["trafo3w"]
        trafo3w_ix_offset = np.sum(~branch_mask[:trafo3w_ix_start])
        num_trafo3w = net.trafo3w.shape[0]
        trafo3w_ix_start, trafo3w_ix_end = trafo3w_ix_start - trafo3w_ix_offset, trafo3w_ix_end - trafo3w_ix_offset
        map_trafo3w = {trafo3w_ix: {'hv': br_ix, 'mv': br_ix+num_trafo3w, 'lv': br_ix+2*num_trafo3w}
                        for trafo3w_ix, br_ix in
                       zip(net.trafo3w.index, range(trafo3w_ix_start, trafo3w_ix_start+num_trafo3w))
                       if branch_mask[br_ix+trafo3w_ix_offset]}

    # set measurements for ppc format
    # add 9 columns to ppc[bus] for Vm, Vm std dev, P, P std dev, Q, Q std dev,
    # pandapower measurement indices V, P, Q
    bus_append = np.full((ppci["bus"].shape[0], bus_cols_se), np.nan, dtype=ppci["bus"].dtype)

    v_measurements = meas_bus[(meas_bus.measurement_type == "v")]
    if len(v_measurements):
        bus_positions = map_bus[v_measurements.element.values.astype(int)]
        bus_append[bus_positions, VM] = v_measurements.value.values
        bus_append[bus_positions, VM_STD] = v_measurements.std_dev.values
        bus_append[bus_positions, VM_IDX] = v_measurements.index.values

    va_measurements = meas_bus[(meas_bus.measurement_type == "va")]
    if len(va_measurements):
        bus_positions = map_bus[va_measurements.element.values.astype(int)]
        bus_append[bus_positions, VA] = va_measurements.value.values
        bus_append[bus_positions, VA_STD] = va_measurements.std_dev.values
        bus_append[bus_positions, VA_IDX] = va_measurements.index.values

    p_measurements = meas_bus[(meas_bus.measurement_type == "p")]
    if len(p_measurements):
        bus_positions = map_bus[p_measurements.element.values.astype(int)]
        unique_bus_positions = np.unique(bus_positions)
        if len(unique_bus_positions) < len(bus_positions):
            std_logger.debug("P Measurement duplication will be automatically merged!")
            for bus in unique_bus_positions:
                p_meas_on_bus = p_measurements.iloc[np.argwhere(bus_positions==bus).ravel(), :]
                bus_append[bus, P] = p_meas_on_bus.value.sum()
                bus_append[bus, P_STD] = p_meas_on_bus.std_dev.max()
                bus_append[bus, P_IDX] = p_meas_on_bus.index[0]
        else:
            bus_append[bus_positions, P] = p_measurements.value.values
            bus_append[bus_positions, P_STD] = p_measurements.std_dev.values
            bus_append[bus_positions, P_IDX] = p_measurements.index.values

    q_measurements = meas_bus[(meas_bus.measurement_type == "q")]
    if len(q_measurements):
        bus_positions = map_bus[q_measurements.element.values.astype(int)]
        unique_bus_positions = np.unique(bus_positions)
        if len(unique_bus_positions) < len(bus_positions):
            std_logger.debug("Q Measurement duplication will be automatically merged!")
            for bus in unique_bus_positions:
                q_meas_on_bus = q_measurements.iloc[np.argwhere(bus_positions==bus).ravel(), :]
                bus_append[bus, Q] = q_meas_on_bus.value.sum()
                bus_append[bus, Q_STD] = q_meas_on_bus.std_dev.max()
                bus_append[bus, Q_IDX] = q_meas_on_bus.index[0]
        else:
            bus_positions = map_bus[q_measurements.element.values.astype(int)]
            bus_append[bus_positions, Q] = q_measurements.value.values
            bus_append[bus_positions, Q_STD] = q_measurements.std_dev.values
            bus_append[bus_positions, Q_IDX] = q_measurements.index.values

    #add zero injection measurement and labels defined in parameter zero_injection
    bus_append = _add_zero_injection(net, ppci, bus_append, zero_injection)
    # add virtual measurements for artificial buses, which were created because
    # of an open line switch. p/q are 0. and std dev is 1. (small value)
    new_in_line_buses = np.setdiff1d(np.arange(ppci["bus"].shape[0]), map_bus[map_bus >= 0])
    bus_append[new_in_line_buses, 2] = 0.
    bus_append[new_in_line_buses, 3] = 1.
    bus_append[new_in_line_buses, 4] = 0.
    bus_append[new_in_line_buses, 5] = 1.

    # add 15 columns to mpc[branch] for Im_from, Im_from std dev, Im_to, Im_to std dev,
    # P_from, P_from std dev, P_to, P_to std dev, Q_from, Q_from std dev,  Q_to, Q_to std dev,
    # pandapower measurement index I, P, Q
    branch_append = np.full((ppci["branch"].shape[0], branch_cols_se),
                            np.nan, dtype=ppci["branch"].dtype)

    if map_line is not None:
        i_measurements = meas[(meas.measurement_type == "i") & (meas.element_type == "line") &\
                              meas.element.isin(map_line)]
        if len(i_measurements):
            meas_from = i_measurements[(i_measurements.side.values.astype(int) ==
                                        net.line.from_bus[i_measurements.element]).values]
            meas_to = i_measurements[(i_measurements.side.values.astype(int) ==
                                      net.line.to_bus[i_measurements.element]).values]
            ix_from = [map_line[l] for l in meas_from.element.values.astype(int)]
            ix_to = [map_line[l] for l in meas_to.element.values.astype(int)]
            branch_append[ix_from, IM_FROM] = meas_from.value.values
            branch_append[ix_from, IM_FROM_STD] = meas_from.std_dev.values
            branch_append[ix_from, IM_FROM_IDX] = meas_from.index.values
            branch_append[ix_to, IM_TO] = meas_to.value.values
            branch_append[ix_to, IM_TO_STD] = meas_to.std_dev.values
            branch_append[ix_to, IM_TO_IDX] = meas_to.index.values

        ia_measurements = meas[(meas.measurement_type == "ia") & (meas.element_type == "line") & \
                                  meas.element.isin(map_line)]
        if len(ia_measurements):
            dg_meas_from = ia_measurements[(ia_measurements.side.values.astype(int) ==
                                            net.line.from_bus[ia_measurements.element]).values]
            dg_meas_to = ia_measurements[(ia_measurements.side.values.astype(int) ==
                                          net.line.to_bus[ia_measurements.element]).values]
            ix_from = [map_line[l] for l in dg_meas_from.element.values.astype(int)]
            ix_to = [map_line[l] for l in dg_meas_to.element.values.astype(int)]
            branch_append[ix_from, IA_FROM] = dg_meas_from.value.values
            branch_append[ix_from, IA_FROM_STD] = dg_meas_from.std_dev.values
            branch_append[ix_from, IA_FROM_IDX] = dg_meas_from.index.values
            branch_append[ix_to, IA_TO] = dg_meas_to.value.values
            branch_append[ix_to, IA_TO_STD] = dg_meas_to.std_dev.values
            branch_append[ix_to, IA_TO_IDX] = dg_meas_to.index.values

        p_measurements = meas[(meas.measurement_type == "p") & (meas.element_type == "line") &
                              meas.element.isin(map_line)]
        if len(p_measurements):
            meas_from = p_measurements[(p_measurements.side.values.astype(int) ==
                                        net.line.from_bus[p_measurements.element]).values]
            meas_to = p_measurements[(p_measurements.side.values.astype(int) ==
                                      net.line.to_bus[p_measurements.element]).values]
            ix_from = [map_line[l] for l in meas_from.element.values.astype(int)]
            ix_to = [map_line[l] for l in meas_to.element.values.astype(int)]
            branch_append[ix_from, P_FROM] = meas_from.value.values
            branch_append[ix_from, P_FROM_STD] = meas_from.std_dev.values
            branch_append[ix_from, P_FROM_IDX] = meas_from.index.values
            branch_append[ix_to, P_TO] = meas_to.value.values
            branch_append[ix_to, P_TO_STD] = meas_to.std_dev.values
            branch_append[ix_to, P_TO_IDX] = meas_to.index.values

        q_measurements = meas[(meas.measurement_type == "q") & (meas.element_type == "line") &
                              meas.element.isin(map_line)]
        if len(q_measurements):
            meas_from = q_measurements[(q_measurements.side.values.astype(int) ==
                                        net.line.from_bus[q_measurements.element]).values]
            meas_to = q_measurements[(q_measurements.side.values.astype(int) ==
                                      net.line.to_bus[q_measurements.element]).values]
            ix_from = [map_line[l] for l in meas_from.element.values.astype(int)]
            ix_to = [map_line[l] for l in meas_to.element.values.astype(int)]
            branch_append[ix_from, Q_FROM] = meas_from.value.values
            branch_append[ix_from, Q_FROM_STD] = meas_from.std_dev.values
            branch_append[ix_from, Q_FROM_IDX] = meas_from.index.values
            branch_append[ix_to, Q_TO] = meas_to.value.values
            branch_append[ix_to, Q_TO_STD] = meas_to.std_dev.values
            branch_append[ix_to, Q_TO_IDX] = meas_to.index.values

    # TODO review in 2019 -> is this a use case? create test with switches on lines
    # determine number of lines in ppci["branch"]
    # out of service lines and lines with open switches at both ends are not in the ppci
    # _is_elements = net["_is_elements"]
    # if "line" not in _is_elements:
    #     get_is_lines(net)
    # lines_is = _is_elements['line']
    # bus_is_idx = _is_elements['bus_is_idx']
    # slidx = (net["switch"]["closed"].values == 0) \
    #         & (net["switch"]["et"].values == "l") \
    #         & (np.in1d(net["switch"]["element"].values, lines_is.index)) \
    #         & (np.in1d(net["switch"]["bus"].values, bus_is_idx))
    # ppci_lines = len(lines_is) - np.count_nonzero(slidx)

    if map_trafo is not None:
        i_tr_measurements = meas[(meas.measurement_type == "i") & (meas.element_type == "trafo") &
                                 meas.element.isin(map_trafo)]
        if len(i_tr_measurements):
            meas_from = i_tr_measurements[(i_tr_measurements.side.values.astype(int) ==
                                           net.trafo.hv_bus[i_tr_measurements.element]).values]
            meas_to = i_tr_measurements[(i_tr_measurements.side.values.astype(int) ==
                                         net.trafo.lv_bus[i_tr_measurements.element]).values]
            ix_from = [map_trafo[t] for t in meas_from.element.values.astype(int)]
            ix_to = [map_trafo[t] for t in meas_to.element.values.astype(int)]
            branch_append[ix_from, IM_FROM] = meas_from.value.values
            branch_append[ix_from, IM_FROM_STD] = meas_from.std_dev.values
            branch_append[ix_from, IM_FROM_IDX] = meas_from.index.values
            branch_append[ix_to, IM_TO] = meas_to.value.values
            branch_append[ix_to, IM_TO_STD] = meas_to.std_dev.values
            branch_append[ix_to, IM_TO_IDX] = meas_to.index.values

        i_tr_dg_measurements = meas[(meas.measurement_type == "ia") & (meas.element_type == "trafo") &
                                     meas.element.isin(map_trafo)]
        if len(i_tr_dg_measurements):
            dg_meas_from = i_tr_dg_measurements[(i_tr_dg_measurements.side.values.astype(int) ==
                                           net.trafo.hv_bus[i_tr_dg_measurements.element]).values]
            dg_meas_to = i_tr_dg_measurements[(i_tr_dg_measurements.side.values.astype(int) ==
                                         net.trafo.lv_bus[i_tr_dg_measurements.element]).values]
            ix_from = [map_trafo[t] for t in dg_meas_from.element.values.astype(int)]
            ix_to = [map_trafo[t] for t in dg_meas_to.element.values.astype(int)]
            branch_append[ix_from, IA_FROM] = dg_meas_from.value.values
            branch_append[ix_from, IA_FROM_STD] = dg_meas_from.std_dev.values
            branch_append[ix_from, IA_FROM_IDX] = dg_meas_from.index.values
            branch_append[ix_to, IA_TO] = dg_meas_to.value.values
            branch_append[ix_to, IA_TO_STD] = dg_meas_to.std_dev.values
            branch_append[ix_to, IA_TO_IDX] = dg_meas_to.index.values

        p_tr_measurements = meas[(meas.measurement_type == "p") & (meas.element_type == "trafo") &
                                 meas.element.isin(map_trafo)]
        if len(p_tr_measurements):
            meas_from = p_tr_measurements[(p_tr_measurements.side.values.astype(int) ==
                                           net.trafo.hv_bus[p_tr_measurements.element]).values]
            meas_to = p_tr_measurements[(p_tr_measurements.side.values.astype(int) ==
                                         net.trafo.lv_bus[p_tr_measurements.element]).values]
            ix_from = [map_trafo[t] for t in meas_from.element.values.astype(int)]
            ix_to = [map_trafo[t] for t in meas_to.element.values.astype(int)]
            branch_append[ix_from, P_FROM] = meas_from.value.values
            branch_append[ix_from, P_FROM_STD] = meas_from.std_dev.values
            branch_append[ix_from, P_FROM_IDX] = meas_from.index.values
            branch_append[ix_to, P_TO] = meas_to.value.values
            branch_append[ix_to, P_TO_STD] = meas_to.std_dev.values
            branch_append[ix_to, P_TO_IDX] = meas_to.index.values

        q_tr_measurements = meas[(meas.measurement_type == "q") & (meas.element_type == "trafo") &
                                 meas.element.isin(map_trafo)]
        if len(q_tr_measurements):
            meas_from = q_tr_measurements[(q_tr_measurements.side.values.astype(int) ==
                                           net.trafo.hv_bus[q_tr_measurements.element]).values]
            meas_to = q_tr_measurements[(q_tr_measurements.side.values.astype(int) ==
                                         net.trafo.lv_bus[q_tr_measurements.element]).values]
            ix_from = [map_trafo[t] for t in meas_from.element.values.astype(int)]
            ix_to = [map_trafo[t] for t in meas_to.element.values.astype(int)]
            branch_append[ix_from, Q_FROM] = meas_from.value.values
            branch_append[ix_from, Q_FROM_STD] = meas_from.std_dev.values
            branch_append[ix_from, Q_FROM_IDX] = meas_from.index.values
            branch_append[ix_to, Q_TO] = meas_to.value.values
            branch_append[ix_to, Q_TO_STD] = meas_to.std_dev.values
            branch_append[ix_to, Q_TO_IDX] = meas_to.index.values

    # Add measurements for trafo3w
    if map_trafo3w is not None:
        i_tr3w_measurements = meas[(meas.measurement_type == "i") & (meas.element_type == "trafo3w") &
                                   meas.element.isin(map_trafo3w)]
        if len(i_tr3w_measurements):
            meas_hv = i_tr3w_measurements[(i_tr3w_measurements.side.values.astype(int) ==
                                           net.trafo3w.hv_bus[i_tr3w_measurements.element]).values]
            meas_mv = i_tr3w_measurements[(i_tr3w_measurements.side.values.astype(int) ==
                                           net.trafo3w.mv_bus[i_tr3w_measurements.element]).values]
            meas_lv = i_tr3w_measurements[(i_tr3w_measurements.side.values.astype(int) ==
                                           net.trafo3w.lv_bus[i_tr3w_measurements.element]).values]
            ix_hv = [map_trafo3w[t]['hv'] for t in meas_hv.element.values.astype(int)]
            ix_mv = [map_trafo3w[t]['mv'] for t in meas_mv.element.values.astype(int)]
            ix_lv = [map_trafo3w[t]['lv'] for t in meas_lv.element.values.astype(int)]
            branch_append[ix_hv, IM_FROM] = meas_hv.value.values
            branch_append[ix_hv, IM_FROM_STD] = meas_hv.std_dev.values
            branch_append[ix_hv, IM_FROM_IDX] = meas_hv.index.values
            branch_append[ix_mv, IM_TO] = meas_mv.value.values
            branch_append[ix_mv, IM_TO_STD] = meas_mv.std_dev.values
            branch_append[ix_mv, IM_TO_IDX] = meas_mv.index.values
            branch_append[ix_lv, IM_TO] = meas_lv.value.values
            branch_append[ix_lv, IM_TO_STD] = meas_lv.std_dev.values
            branch_append[ix_lv, IM_TO_IDX] = meas_lv.index.values

    # Add degree_measurements for trafo3w
    if map_trafo3w is not None:
        i_tr3w_dg_measurements = meas[(meas.measurement_type == "ia") & (meas.element_type == "trafo3w") &
                                   meas.element.isin(map_trafo3w)]
        if len(i_tr3w_dg_measurements):
            dg_meas_hv = i_tr3w_dg_measurements[(i_tr3w_dg_measurements.side.values.astype(int) ==
                                           net.trafo3w.hv_bus[i_tr3w_dg_measurements.element]).values]
            dg_meas_mv = i_tr3w_dg_measurements[(i_tr3w_dg_measurements.side.values.astype(int) ==
                                           net.trafo3w.mv_bus[i_tr3w_dg_measurements.element]).values]
            dg_meas_lv = i_tr3w_dg_measurements[(i_tr3w_dg_measurements.side.values.astype(int) ==
                                           net.trafo3w.lv_bus[i_tr3w_dg_measurements.element]).values]
            ix_hv = [map_trafo3w[t]['hv'] for t in dg_meas_hv.element.values.astype(int)]
            ix_mv = [map_trafo3w[t]['mv'] for t in dg_meas_mv.element.values.astype(int)]
            ix_lv = [map_trafo3w[t]['lv'] for t in dg_meas_lv.element.values.astype(int)]
            branch_append[ix_hv, IA_FROM] = dg_meas_hv.value.values
            branch_append[ix_hv, IA_FROM_STD] = dg_meas_hv.std_dev.values
            branch_append[ix_hv, IA_FROM_IDX] = dg_meas_hv.index.values
            branch_append[ix_mv, IA_TO] = dg_meas_mv.value.values
            branch_append[ix_mv, IA_TO_STD] = dg_meas_mv.std_dev.values
            branch_append[ix_mv, IA_TO_IDX] = dg_meas_mv.index.values
            branch_append[ix_lv, IA_TO] = dg_meas_lv.value.values
            branch_append[ix_lv, IA_TO_STD] = dg_meas_lv.std_dev.values
            branch_append[ix_lv, IA_TO_IDX] = dg_meas_lv.index.values

        p_tr3w_measurements = meas[(meas.measurement_type == "p") & (meas.element_type == "trafo3w") &
                                   meas.element.isin(map_trafo3w)]
        if len(p_tr3w_measurements):
            meas_hv = p_tr3w_measurements[(p_tr3w_measurements.side.values.astype(int) ==
                                           net.trafo3w.hv_bus[p_tr3w_measurements.element]).values]
            meas_mv = p_tr3w_measurements[(p_tr3w_measurements.side.values.astype(int) ==
                                           net.trafo3w.mv_bus[p_tr3w_measurements.element]).values]
            meas_lv = p_tr3w_measurements[(p_tr3w_measurements.side.values.astype(int) ==
                                           net.trafo3w.lv_bus[p_tr3w_measurements.element]).values]
            ix_hv = [map_trafo3w[t]['hv'] for t in meas_hv.element.values.astype(int)]
            ix_mv = [map_trafo3w[t]['mv'] for t in meas_mv.element.values.astype(int)]
            ix_lv = [map_trafo3w[t]['lv'] for t in meas_lv.element.values.astype(int)]
            branch_append[ix_hv, P_FROM] = meas_hv.value.values
            branch_append[ix_hv, P_FROM_STD] = meas_hv.std_dev.values
            branch_append[ix_hv, P_FROM_IDX] = meas_hv.index.values
            branch_append[ix_mv, P_TO] = meas_mv.value.values
            branch_append[ix_mv, P_TO_STD] = meas_mv.std_dev.values
            branch_append[ix_mv, P_TO_IDX] = meas_mv.index.values
            branch_append[ix_lv, P_TO] = meas_lv.value.values
            branch_append[ix_lv, P_TO_STD] = meas_lv.std_dev.values
            branch_append[ix_lv, P_TO_IDX] = meas_lv.index.values

        q_tr3w_measurements = meas[(meas.measurement_type == "q") & (meas.element_type == "trafo3w") &
                                   meas.element.isin(map_trafo3w)]
        if len(q_tr3w_measurements):
            meas_hv = q_tr3w_measurements[(q_tr3w_measurements.side.values.astype(int) ==
                                           net.trafo3w.hv_bus[q_tr3w_measurements.element]).values]
            meas_mv = q_tr3w_measurements[(q_tr3w_measurements.side.values.astype(int) ==
                                           net.trafo3w.mv_bus[q_tr3w_measurements.element]).values]
            meas_lv = q_tr3w_measurements[(q_tr3w_measurements.side.values.astype(int) ==
                                           net.trafo3w.lv_bus[q_tr3w_measurements.element]).values]
            ix_hv = [map_trafo3w[t]['hv'] for t in meas_hv.element.values.astype(int)]
            ix_mv = [map_trafo3w[t]['mv'] for t in meas_mv.element.values.astype(int)]
            ix_lv = [map_trafo3w[t]['lv'] for t in meas_lv.element.values.astype(int)]
            branch_append[ix_hv, Q_FROM] = meas_hv.value.values
            branch_append[ix_hv, Q_FROM_STD] = meas_hv.std_dev.values
            branch_append[ix_hv, Q_FROM_IDX] = meas_hv.index.values
            branch_append[ix_mv, Q_TO] = meas_mv.value.values
            branch_append[ix_mv, Q_TO_STD] = meas_mv.std_dev.values
            branch_append[ix_mv, Q_TO_IDX] = meas_mv.index.values
            branch_append[ix_lv, Q_TO] = meas_lv.value.values
            branch_append[ix_lv, Q_TO_STD] = meas_lv.std_dev.values
            branch_append[ix_lv, Q_TO_IDX] = meas_lv.index.values

    # Check append or update
    if ppci["bus"].shape[1] == bus_cols:
        ppci["bus"] = np.hstack((ppci["bus"], bus_append))
    else:
        ppci["bus"][:, bus_cols: bus_cols+bus_cols_se] = bus_append

    if ppci["branch"].shape[1] == branch_cols:
        ppci["branch"] = np.hstack((ppci["branch"], branch_append))
    else:
        ppci["branch"][:, branch_cols: branch_cols+branch_cols_se] = branch_append
    return ppci

def _add_zero_injection(net, ppci, bus_append, zero_injection):
    """
    Add zero injection labels to the ppci structure and add virtual measurements to those buses
    :param net: pandapower net
    :param ppci: generated ppci
    :param bus_append: added columns to the ppci bus with zero injection label
    :param zero_injection: parameter to control which bus to be identified as zero injection
    :return bus_append: added columns
    """
    bus_append[:, ZERO_INJ_FLAG] = False
    if zero_injection is not None:
        # identify aux bus to zero injection
        if net._pd2ppc_lookups['aux']:
            aux_bus_lookup = np.concatenate([v for k,v in net._pd2ppc_lookups['aux'].items() if k != 'xward'])
            aux_bus = net._pd2ppc_lookups['bus'][aux_bus_lookup]
            bus_append[aux_bus, ZERO_INJ_FLAG] = True

        if isinstance(zero_injection, str):
            if zero_injection == 'auto':
                # identify bus without elements and pq measurements as zero injection
                zero_inj_bus_mask = (ppci["bus"][:, 1] == 1) & (ppci["bus"][:, 2:6]==0).all(axis=1) &\
                    np.isnan(bus_append[:, P:(Q_STD+1)]).all(axis=1)
                bus_append[zero_inj_bus_mask, ZERO_INJ_FLAG] = True
            elif zero_injection != "aux_bus":
                raise UserWarning("zero injection parameter is not correctly initialized")
        elif hasattr(zero_injection, '__iter__'):
            zero_inj_bus = net._pd2ppc_lookups['bus'][zero_injection]
            bus_append[zero_inj_bus, ZERO_INJ_FLAG] = True

        zero_inj_bus = np.argwhere(bus_append[:, ZERO_INJ_FLAG]).ravel()
        bus_append[zero_inj_bus, P] = 0
        bus_append[zero_inj_bus, P_STD] = ZERO_INJECTION_STD_DEV
        bus_append[zero_inj_bus, Q] = 0
        bus_append[zero_inj_bus, Q_STD] = ZERO_INJECTION_STD_DEV
    return bus_append

def _build_measurement_vectors(ppci, update_meas_only=False):
    """
    Building measurement vector z, pandapower to ppci measurement mapping and covariance matrix R
    :param ppci: generated ppci which contains the measurement columns
    :param branch_cols: number of columns in original ppci["branch"] without measurements
    :param bus_cols: number of columns in original ppci["bus"] without measurements
    :return: both created vectors
    """
    p_bus_not_nan = ~np.isnan(ppci["bus"][:, bus_cols + P])
    p_line_f_not_nan = ~np.isnan(ppci["branch"][:, branch_cols + P_FROM])
    p_line_t_not_nan = ~np.isnan(ppci["branch"][:, branch_cols + P_TO])
    q_bus_not_nan = ~np.isnan(ppci["bus"][:, bus_cols + Q])
    q_line_f_not_nan = ~np.isnan(ppci["branch"][:, branch_cols + Q_FROM])
    q_line_t_not_nan = ~np.isnan(ppci["branch"][:, branch_cols + Q_TO])
    v_bus_not_nan = ~np.isnan(ppci["bus"][:, bus_cols + VM])
    v_degree_bus_not_nan = ~np.isnan(ppci["bus"][:, bus_cols + VA])
    i_line_f_not_nan = ~np.isnan(ppci["branch"][:, branch_cols + IM_FROM])
    i_line_t_not_nan = ~np.isnan(ppci["branch"][:, branch_cols + IM_TO])
    i_degree_line_f_not_nan = ~np.isnan(ppci["branch"][:, branch_cols + IA_FROM])
    i_degree_line_t_not_nan = ~np.isnan(ppci["branch"][:, branch_cols + IA_TO])
    # piece together our measurement vector z
    z = np.concatenate((ppci["bus"][p_bus_not_nan, bus_cols + P],
                        ppci["branch"][p_line_f_not_nan, branch_cols + P_FROM],
                        ppci["branch"][p_line_t_not_nan, branch_cols + P_TO],
                        ppci["bus"][q_bus_not_nan, bus_cols + Q],
                        ppci["branch"][q_line_f_not_nan, branch_cols + Q_FROM],
                        ppci["branch"][q_line_t_not_nan, branch_cols + Q_TO],
                        ppci["bus"][v_bus_not_nan, bus_cols + VM],
                        ppci["bus"][v_degree_bus_not_nan, bus_cols + VA],
                        ppci["branch"][i_line_f_not_nan, branch_cols + IM_FROM],
                        ppci["branch"][i_line_t_not_nan, branch_cols + IM_TO],
                        ppci["branch"][i_degree_line_f_not_nan, branch_cols + IA_FROM],
                        ppci["branch"][i_degree_line_t_not_nan, branch_cols + IA_TO]
                        )).real.astype(np.float64)
    if not update_meas_only:
        # conserve the pandapower indices of measurements in the ppci order
        pp_meas_indices = np.concatenate((ppci["bus"][p_bus_not_nan, bus_cols + P_IDX],
                                          ppci["branch"][p_line_f_not_nan, branch_cols + P_FROM_IDX],
                                          ppci["branch"][p_line_t_not_nan, branch_cols + P_TO_IDX],
                                          ppci["bus"][q_bus_not_nan, bus_cols + Q_IDX],
                                          ppci["branch"][q_line_f_not_nan, branch_cols + Q_FROM_IDX],
                                          ppci["branch"][q_line_t_not_nan, branch_cols + Q_TO_IDX],
                                          ppci["bus"][v_bus_not_nan, bus_cols + VM_IDX],
                                          ppci["bus"][v_degree_bus_not_nan, bus_cols + VA_IDX],
                                          ppci["branch"][i_line_f_not_nan, branch_cols + IM_FROM_IDX],
                                          ppci["branch"][i_line_t_not_nan, branch_cols + IM_TO_IDX],
                                          ppci["branch"][i_degree_line_f_not_nan, branch_cols + IA_FROM_IDX],
                                          ppci["branch"][i_degree_line_t_not_nan, branch_cols + IA_TO_IDX]
                                          )).real.astype(int)
        # Covariance matrix R
        r_cov = np.concatenate((ppci["bus"][p_bus_not_nan, bus_cols + P_STD],
                                ppci["branch"][p_line_f_not_nan, branch_cols + P_FROM_STD],
                                ppci["branch"][p_line_t_not_nan, branch_cols + P_TO_STD],
                                ppci["bus"][q_bus_not_nan, bus_cols + Q_STD],
                                ppci["branch"][q_line_f_not_nan, branch_cols + Q_FROM_STD],
                                ppci["branch"][q_line_t_not_nan, branch_cols + Q_TO_STD],
                                ppci["bus"][v_bus_not_nan, bus_cols + VM_STD],
                                ppci["bus"][v_degree_bus_not_nan, bus_cols + VA_STD],
                                ppci["branch"][i_line_f_not_nan, branch_cols + IM_FROM_STD],
                                ppci["branch"][i_line_t_not_nan, branch_cols + IM_TO_STD],
                                ppci["branch"][i_degree_line_f_not_nan, branch_cols + IA_FROM_STD],
                                ppci["branch"][i_degree_line_t_not_nan, branch_cols + IA_TO_STD]
                                )).real.astype(np.float64)
        meas_mask = np.concatenate([p_bus_not_nan,
                                    p_line_f_not_nan,
                                    p_line_t_not_nan,
                                    q_bus_not_nan,
                                    q_line_f_not_nan,
                                    q_line_t_not_nan,
                                    v_bus_not_nan,
                                    v_degree_bus_not_nan,
                                    i_line_f_not_nan,
                                    i_line_t_not_nan,
                                    i_degree_line_f_not_nan,
                                    i_degree_line_t_not_nan])
        any_i_meas = np.any(np.r_[i_line_f_not_nan, i_line_t_not_nan])
        any_degree_meas = np.any(np.r_[v_degree_bus_not_nan,
                                       i_degree_line_f_not_nan,
                                       i_degree_line_t_not_nan])
        return z, pp_meas_indices, r_cov, meas_mask, any_i_meas, any_degree_meas
    else:
        return z

def pp2eppci(net, v_start=None, delta_start=None, calculate_voltage_angles=True, zero_injection="aux_bus",
             ppc=None, eppci=None):
    # initialize result tables if not existent
    _copy_power_flow_results(net)
    if isinstance(eppci, ExtendedPPCI):
        eppci.data = _add_measurements_to_ppci(net, eppci.data, zero_injection)
        eppci.update_meas()
        return net, ppc, eppci
    else:
        # initialize ppc
        ppc, ppci = _init_ppc(net, v_start, delta_start, calculate_voltage_angles)

        # add measurements to ppci structure
        # Finished converting pandapower network to ppci
        ppci = _add_measurements_to_ppci(net, ppci, zero_injection)
        return net, ppc, ExtendedPPCI(ppci)

class ExtendedPPCI(UserDict):
    def __init__(self, ppci):
        self.data = ppci

        # Measurement relevant parameters
        self.z = None
        self.r_cov = None
        self.pp_meas_indices = None
        self.non_nan_meas_mask = None
        self.non_nan_meas_selector = None
        self.any_i_meas = False
        self.any_degree_meas = False
        self._initialize_meas()

        # check slack bus
        self.non_slack_buses = np.argwhere(ppci["bus"][:, idx_bus.BUS_TYPE] != 3).ravel()
        self.non_slack_bus_mask = (ppci['bus'][:, idx_bus.BUS_TYPE] != 3).ravel()
        self.num_non_slack_bus = np.sum(self.non_slack_bus_mask)
        self.delta_v_bus_mask = np.r_[self.non_slack_bus_mask,
                                      np.ones(self.non_slack_bus_mask.shape[0], dtype=bool)].ravel()
        self.delta_v_bus_selector = np.flatnonzero(self.delta_v_bus_mask)

        # Initialize state variable
        self.v_init = ppci["bus"][:, idx_bus.VM]
        self.delta_init = np.radians(ppci["bus"][:, idx_bus.VA])
        self.E_init = np.r_[self.delta_init[self.non_slack_bus_mask], self.v_init]
        self.v = self.v_init.copy()
        self.delta = self.delta_init.copy()
        self.E = self.E_init.copy()

    def _initialize_meas(self):
        # calculate relevant vectors from ppci measurements
        self.z, self.pp_meas_indices, self.r_cov, self.non_nan_meas_mask, self.any_i_meas, self.any_degree_meas =\
             _build_measurement_vectors(self, update_meas_only=False)
        self.non_nan_meas_selector = np.flatnonzero(self.non_nan_meas_mask)

    def update_meas(self):
        self.z = _build_measurement_vectors(self, update_meas_only=True)

    @property
    def V(self):
        return self.v * np.exp(1j * self.delta)

    def reset(self):
        self.v, self.delta, self.E = self.v_init.copy(), self.delta_init.copy(), self.E_init.copy()

    def update_E(self, E):
        self.E = E
        self.v = E[self.num_non_slack_bus:]
        self.delta[self.non_slack_buses] = E[:self.num_non_slack_bus]

    def E2V(self, E):
        self.update_E(E)
        return self.V

    def get_Y(self):
        # Using recycled version if available
        if "Ybus" in self["internal"] and self["internal"]["Ybus"].size:
            Ybus, Yf, Yt = self["internal"]['Ybus'], self["internal"]['Yf'], self["internal"]['Yt']
        else:
            ## build admittance matrices
            Ybus, Yf, Yt = makeYbus(self['baseMVA'], self['bus'], self['branch'])
            self["internal"]['Ybus'], self["internal"]['Yf'], self["internal"]['Yt'] = Ybus, Yf, Yt
        return Ybus, Yf, Yt

max_stars_repo_name: MDiesing/pandapower

max_issues_repo_name: MDiesing/pandapower

max_forks_repo_name: MDiesing/pandapower